convolution norm fixing

Author: AAAlvesJr

examples/convolution/convolute_functions.inl

@@ -47,7 +47,8 @@
 #include <hydra/functions/Ipatia.h>
 #include <hydra/device/System.h>
 #include <hydra/functions/ConvolutionFunctor.h>
-#include <hydra/Random.h>
+#include <hydra/GaussKronrodQuadrature.h>
+
 //hydra
 #if HYDRA_DEVICE_SYSTEM == CUDA
 #include <hydra/CuFFT.h>
@@ -112,9 +113,9 @@ int main(int argv, char** argc)
 
 	// gaussian
 	auto mean   = hydra::Parameter::Create( "mean").Value(0.0).Error(0.0001);
-	auto sigma  = hydra::Parameter::Create("sigma").Value(3.5).Error(0.0001);
+	auto sigma  = hydra::Parameter::Create("sigma").Value(1.0).Error(0.0001);
 
-	hydra::Gaussian<double> gaussian_kernel(mean,  sigma);
+	hydra::Gaussian<double> kernel(mean,  sigma);
 
 	//===========================
 	// signals
@@ -162,7 +163,7 @@ int main(int argv, char** argc)
 	auto start_d = std::chrono::high_resolution_clock::now();
 
 	hydra::convolute(hydra::device::sys, fft_backend,
-			signal, gaussian_kernel, min, max,  conv_result, true);
+			signal, kernel, min, max,  conv_result, true);
 
 	auto end_d = std::chrono::high_resolution_clock::now();
 
@@ -176,7 +177,22 @@ int main(int argv, char** argc)
 	/*
 	 * using the hydra::ConvolutionFunctor
 	 */
-    auto convoluton = hydra::make_convolution<double>(  hydra::device::sys,  fft_backend, signal, gaussian_kernel, min, max,  conv_result.size() );
+    auto convoluton = hydra::make_convolution<double>(  hydra::device::sys,  fft_backend,
+    		signal, kernel, min, max,  conv_result.size() );
+
+    hydra::GaussKronrodQuadrature<61,100, hydra::device::sys_t> GKQ61_Kernel(-25.0, 25.0);
+    hydra::GaussKronrodQuadrature<61,100, hydra::device::sys_t> GKQ61_Signal(min,  max);
+
+    auto kernel_int = GKQ61_Kernel.Integrate(kernel);
+    auto signal_int = GKQ61_Signal.Integrate(signal);
+    auto convol_int = GKQ61_Signal.Integrate(convoluton);
+
+    std::cout << "===========================================" << std::endl;
+    std::cout << "Kernel Integral: " <<kernel_int.first <<"+/-"<< kernel_int.second << std::endl;
+    std::cout << "Signal Integral: " <<signal_int.first <<"+/-"<< signal_int.second << std::endl;
+    std::cout << "Convol Integral: " <<convol_int.first <<"+/-"<< convol_int.second << std::endl;
+    std::cout << "===========================================" << std::endl;
+
 
 	//------------------------
 	//------------------------
@@ -193,7 +209,7 @@ int main(int argv, char** argc)
 		hist_convol_functor->SetBinContent(i, convoluton(hist_convol_functor->GetBinCenter(i)) );
 		hist_convol->SetBinContent(i, conv_result[i-1] );
 		hist_signal->SetBinContent(i, signal(hist_signal->GetBinCenter(i) ) );
-		hist_kernel->SetBinContent(i, gaussian_kernel( hist_kernel->GetBinCenter(i)));
+		hist_kernel->SetBinContent(i, kernel( hist_kernel->GetBinCenter(i)));
 	}
 #endif //_ROOT_AVAILABLE_
 

hydra/Convolution.h

@@ -79,7 +79,6 @@ convolute(detail::BackendPolicy<BACKEND> policy, detail::FFTPolicy<T, FFTBackend
 
 	int nsamples = std::forward<Iterable>(output).size();
 
-
 	T delta = (max - min)/(nsamples);
 
 	auto complex_buffer  = hydra_thrust::get_temporary_buffer<complex_type>(policy, nsamples+1);
@@ -88,20 +87,23 @@ convolute(detail::BackendPolicy<BACKEND> policy, detail::FFTPolicy<T, FFTBackend
 
 	//
 	auto counting_samples = range(0, 2*nsamples);
+
 	// sample kernel
 	auto kernel_sampler = hydra::detail::convolution::KernelSampler<Kernel>(kernel, nsamples, delta);
 
 	hydra_thrust::transform(policy, counting_samples.begin(), counting_samples.end(),
 			kernel_samples.first , kernel_sampler);
 
-	auto norm_factor = hydra_thrust::reduce(policy, kernel_samples.first, kernel_samples.first + kernel_samples.second );
+	//auto norm_factor = hydra_thrust::reduce(policy, kernel_samples.first, kernel_samples.first + kernel_samples.second );
 
 	// sample function
 	auto functor_sampler = hydra::detail::convolution::FunctorSampler<Functor>(functor, nsamples,  min, delta);
 
 	hydra_thrust::transform( policy, counting_samples.begin(), counting_samples.end(),
 			functor_samples.first, functor_sampler);
 
+	//norm_factor *= hydra_thrust::reduce(policy, functor_samples.first, functor_samples.first + functor_samples.second );
+
 	//transform kernel
 	auto fft_kernel = _RealToComplexFFT( kernel_samples.second );
 
@@ -131,19 +133,19 @@ convolute(detail::BackendPolicy<BACKEND> policy, detail::FFTPolicy<T, FFTBackend
 
 	//transform product back to real
 
-	auto fft_product = _ComplexToRealFFT( 2*nsamples );
+	auto fft_product = _ComplexToRealFFT( 2.0*nsamples );
 
 	fft_product.LoadInputData(complex_buffer.second, complex_buffer.first);
 	fft_product.Execute();
 
 	auto fft_product_output =  fft_product.GetOutputData();
 
-	T n = 2*nsamples*norm_factor;
+	T n = 2*nsamples;
 
 	auto normalize_fft =  detail::convolution::NormalizeFFT<T>(n);
 
     auto first = hydra_thrust::make_transform_iterator( fft_product_output.first,normalize_fft);
-    auto last  = hydra_thrust::make_transform_iterator(fft_product_output.first + nsamples+1,normalize_fft);
+    auto last  = hydra_thrust::make_transform_iterator(fft_product_output.first + nsamples,normalize_fft);
 
 	auto fft_product_range = make_range(first, last);
 

hydra/detail/Convolution.inl

@@ -70,12 +70,11 @@ struct KernelSampler
 
 	KernelSampler(Kernel const& kernel, int nsamples , double delta):
 		fDelta(delta),
-
 		fKernel(kernel)
 	{
 		fNZero = nsamples;
-		fNMin  = nsamples- nsamples/16;
-	    fNMax  = nsamples+ nsamples/16;
+		fNMin  = nsamples- nsamples/2;
+	    fNMax  = nsamples+ nsamples/2;
 	}
 
 	__hydra_host__ __hydra_device__
@@ -215,7 +214,7 @@ struct FunctorSampler
 
         double t =  index*fDelta + fMin;
 
-        if(  (t >= fMin) && ( t <=  (fNSamples)*fDelta + fMin)){
+        if(  (t >= fMin) && ( t < (fNSamples)*fDelta + fMin)){
 
 		    value = fFunctor(t);
 		}