...

Author: AAAlvesJr

CHANGELOG.md

@@ -165,7 +165,40 @@ The filling functions can be called also with a specific backend policy and with
 
 #### Data fitting
 
-1. Adding support to simultaneous fit.
+1. Added support to multi-layered simultaneous fit of different models, over different datasets, deploying different parallelization strategies for each model. No categorization of the dataset is needed, but just to set up preliminarly the different component FCNs, that can be optimized in isolation or in the context of the simultaneous FCN. Simultaneous FCNs can be created via direct instantiation or using the convenience function `hydra::make_simultaneous_fcn(...)`, as shown in the following snippet.
+
+    ```cpp
+    ...
+    #include <hydra/LogLikelihoodFCN.h>
+    ...
+    
+    int main(int argv, char** argc)
+    {
+     ...
+     //=====================================================================================
+     //                                                           +----< fcn(model-x) 
+     //                           +----< simultaneous fcn 1 ----- |
+     //                           |                               +----< fcn(model-y)  
+     //   simultaneous fcn   <----+
+     //                           |                               +----< fcn(model-w)
+     //                           +----< simultaneous fcn 2 ------|
+     //                           |                               +----< fcn(model-z) 
+     //                           +----< fcn(model-v)
+     //=====================================================================================        
+     auto fcnX    = hydra::make_loglikehood_fcn(modelX, dataX);
+     auto fcnY    = hydra::make_loglikehood_fcn(modelY, dataY);
+     auto fcnW    = hydra::make_loglikehood_fcn(modelY, dataY);
+     auto fcnZ    = hydra::make_loglikehood_fcn(modelZ, dataZ);    
+     auto fcnV    = hydra::make_loglikehood_fcn(modelv, datav);
+     
+     auto sim_fcn1 = hydra::make_simultaneous_fcn(fcnx, fcny);
+     auto sim_fcn2 = hydra::make_simultaneous_fcn(fcnw, fcnz);
+     auto sim_fcn  = hydra::make_simultaneous_fcn(sim_fcn1, sim_fcn2, fcnV);
+     ...
+    }
+    ```
+    Moreover, the generic interface allows to build up a simultaneous FCN object by composing usual FCNs and simultaneous FCNs. An example of such new method can be found in `examples/fit/simultaneous_fit.inl`.
+
 2. Fitting of convoluted PDFs.
 
 #### General

CMakeLists.txt

@@ -18,7 +18,7 @@ SET(CMAKE_VERBOSE_MAKEFILE  ON)
 #check if compiler is C++14 compliant
 include(CheckCXXCompilerFlag)
 CHECK_CXX_COMPILER_FLAG("--std=c++14" COMPILER_SUPPORTS_CXX14)
-if(NOT COMPILER_SUPPORTS_CXX11)
+if(NOT COMPILER_SUPPORTS_CXX14)
  message(FATAL "The compiler ${CMAKE_CXX_COMPILER} has no C++14 support. Please use a different C++ compiler.")
 endif()
 

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>
@@ -114,7 +115,7 @@ int main(int argv, char** argc)
 	auto mean   = hydra::Parameter::Create( "mean").Value(0.0).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,28 +79,34 @@ 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);
 	auto kernel_samples  = hydra_thrust::get_temporary_buffer<T>(policy, 2*nsamples);
 	auto functor_samples = hydra_thrust::get_temporary_buffer<T>(policy, 2*nsamples);
 
 	//
-	auto counting_samples = range(0, 2*nsamples-1);
+
+	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 );
+
+
 	// 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 );
 
@@ -134,14 +140,17 @@ convolute(detail::BackendPolicy<BACKEND> policy, detail::FFTPolicy<T, FFTBackend
 
 	//transform product back to real
 
+
 	auto fft_product = _ComplexToRealFFT( complex_buffer.second );
 
+
 	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);
 

hydra/Random.h

@@ -199,6 +199,26 @@ struct is_callable: std::conditional<
 }  // namespace detail
 
 
+/**
+ * \ingroup random
+ *
+ * This functions reorder a dataset to produce a unweighted sample according to the weights
+ * [wbegin, wend]. The length of the range [wbegin, wend] should be equal or greater than
+ * the dataset size.
+ *
+ * @param policy parallel backend to perform the unweighting
+ * @param data_begin iterator pointing to the begin of the range of weights
+ * @param data_end iterator pointing to the begin of the range of weights
+ * @param weights_begin iterator pointing to the begin of the range of data
+ * @return
+ */
+template<typename RNG=default_random_engine, typename DerivedPolicy, typename IteratorData, typename IteratorWeight>
+typename std::enable_if<
+detail::random::is_iterator<IteratorData>::value && detail::random::is_iterator<IteratorWeight>::value,
+Range<IteratorData> >::type
+unweight( hydra_thrust::detail::execution_policy_base<DerivedPolicy>  const& policy, IteratorData data_begin, IteratorData data_end, IteratorWeight weights_begin);
+
+
 
 
 /**
@@ -276,6 +296,25 @@ typename std::enable_if<
 unweight( IterableData data, IterableWeight weights);
 
 
+/**
+ * \ingroup random
+ *
+ * This functions reorder a dataset to produce an unweighted sample according to @param functor .
+ *
+ * @param policy
+ * @param begin
+ * @param end
+ * @param functor
+ * @return the index of the last entry of the unweighted event.
+ */
+template<typename RNG=default_random_engine, typename Functor, typename Iterator, typename DerivedPolicy>
+typename std::enable_if<
+	detail::random::is_callable<Functor>::value && detail::random::is_iterator<Iterator>::value,
+	Range<Iterator>
+>::type
+unweight( hydra_thrust::detail::execution_policy_base<DerivedPolicy> const& policy, Iterator begin, Iterator end, Functor const& functor);
+
+
 /**
  * \ingroup random
  *

hydra/SeedRNG.h

@@ -0,0 +1,91 @@
+/*----------------------------------------------------------------------------
+ *
+ *   Copyright (C) 2016 - 2020 Antonio Augusto Alves Junior
+ *
+ *   This file is part of Hydra Data Analysis Framework.
+ *
+ *   Hydra is free software: you can redistribute it and/or modify
+ *   it under the terms of the GNU General Public License as published by
+ *   the Free Software Foundation, either version 3 of the License, or
+ *   (at your option) any later version.
+ *
+ *   Hydra is distributed in the hope that it will be useful,
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *   GNU General Public License for more details.
+ *
+ *   You should have received a copy of the GNU General Public License
+ *   along with Hydra.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ *---------------------------------------------------------------------------*/
+
+/*
+ * SeedRNG.h
+ *
+ *  Created on: 30/06/2020
+ *      Author: Antonio Augusto Alves Junior
+ */
+
+#ifndef SEEDRNG_H_
+#define SEEDRNG_H_
+
+
+#include <hydra/detail/Config.h>
+#include <stdint.h>
+
+
+namespace hydra {
+
+class SeedRNG
+{
+public:
+
+	 __hydra_host__ __hydra_device__
+	 SeedRNG(size_t state=1337 ):
+	fState(state)
+	{}
+
+	 __hydra_host__ __hydra_device__
+	 SeedRNG(SeedRNG const& other):
+	fState(other.GetState())
+	{}
+
+	 __hydra_host__ __hydra_device__
+	 inline SeedRNG& operator=(SeedRNG const& other)
+	{
+		if(this==&other) return *this;
+
+		fState=other.GetState();
+
+		return *this;
+	}
+
+	 __hydra_host__ __hydra_device__
+	size_t GetState() const {
+		return fState;
+	}
+
+	__hydra_host__ __hydra_device__
+	void SetState(size_t state) {
+		fState = state;
+	}
+
+	__hydra_host__ __hydra_device__
+	size_t operator()()
+	{
+		size_t z = (fState += 0x9e3779b97f4a7c15);
+		z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+		z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+		return z ^ (z >> 31);
+	}
+
+
+private:
+
+	size_t fState;
+
+};
+
+}  // namespace hydra
+
+#endif /* SEEDRNG_H_ */

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);
 		}

hydra/detail/Decays.inl

@@ -85,7 +85,7 @@ public:
 
 		hydra_thrust::uniform_real_distribution<double> uniDist(0.0,1.0);
 
-		return (  weight/fMaxWeight >= uniDist(randEng)) ;
+		return (  weight/fMaxWeight > uniDist(randEng)) ;
 	}
 
 	__hydra_host__ __hydra_device__
@@ -474,7 +474,7 @@ typedef hydra_thrust::transform_iterator<reweight_functor,iterator> reweight_ite
 typedef detail::FlagDaugthers< reweight_functor> tagger_type;
 
 	//number of events to trial
-	size_t ntrials = this->size();
+	size_t ntrials = fDecays.size();
 
 	//create iterators
 	hydra_thrust::counting_iterator < size_t > first(0);