diff --git a/math/mathcore/inc/Fit/FitUtil.h b/math/mathcore/inc/Fit/FitUtil.h index 6556376e01616d69d06adf18907d3069312d80cf..6be79c05af23126dfa6d1bd8b83f4fcbfb18ccfd 100644 --- a/math/mathcore/inc/Fit/FitUtil.h +++ b/math/mathcore/inc/Fit/FitUtil.h @@ -413,8 +413,8 @@ namespace FitUtil { vecCore::Load<T>(invErrorVec, invErrorptr); const T *x; + std::vector<T> xc; if(data.NDim() > 1) { - std::vector<T> xc; xc.resize(data.NDim()); xc[0] = x1; for (unsigned int j = 1; j < data.NDim(); ++j) @@ -565,8 +565,9 @@ namespace FitUtil { vecCore::Load<T>(x1, data.GetCoordComponent(i * vecSize, 0)); const T *x = nullptr; unsigned int ndim = data.NDim(); + std::vector<T> xc; if (ndim > 1) { - std::vector<T> xc(ndim); + xc.resize(ndim); xc[0] = x1; for (unsigned int j = 1; j < ndim; ++j) vecCore::Load<T>(xc[j], data.GetCoordComponent(i * vecSize, j)); @@ -964,8 +965,11 @@ namespace FitUtil { const T *x = nullptr; unsigned int ndim = data.NDim(); + // need to declare vector outside if statement + // otherwise pointer will be invalid + std::vector<T> xc; if (ndim > 1) { - std::vector<T> xc(ndim); + xc.resize(ndim); xc[0] = x1; for (unsigned int j = 1; j < ndim; ++j) vecCore::Load<T>(xc[j], data.GetCoordComponent(i * vecSize, j)); @@ -1150,8 +1154,9 @@ namespace FitUtil { const T *x = nullptr; unsigned ndim = data.NDim(); + std::vector<T> xc; if (ndim > 1) { - std::vector<T> xc(ndim); + xc.resize(ndim); xc[0] = x1; for (unsigned int j = 1; j < ndim; ++j) vecCore::Load<T>(xc[j], data.GetCoordComponent(i * vecSize, j)); @@ -1293,8 +1298,9 @@ namespace FitUtil { const T *x = nullptr; unsigned int ndim = data.NDim(); + std::vector<T> xc(ndim); if (ndim > 1) { - std::vector<T> xc(ndim); + xc.resize(ndim); xc[0] = x1; for (unsigned int j = 1; j < ndim; ++j) vecCore::Load<T>(xc[j], data.GetCoordComponent(i * vecSize, j)); diff --git a/math/mathcore/src/FitUtil.cxx b/math/mathcore/src/FitUtil.cxx index 5c260ee37b8edc9c47dd005cea7ec283816ad93a..aaff90833e49d059d8a1b178628f2b2ce699bbda 100644 --- a/math/mathcore/src/FitUtil.cxx +++ b/math/mathcore/src/FitUtil.cxx @@ -33,6 +33,8 @@ #include <numeric> //#include <memory> +#include "TROOT.h" + //#define DEBUG #ifdef DEBUG #define NSAMPLE 10 @@ -966,7 +968,7 @@ double FitUtil::EvaluateLogL(const IModelFunctionTempl<double> &func, const UnBi #ifdef USE_PARAMCACHE fval = func(x.data()); #else - fval = func(x.data(), p); + fval = func(x.data(), p); #endif // one -dim case @@ -975,7 +977,7 @@ double FitUtil::EvaluateLogL(const IModelFunctionTempl<double> &func, const UnBi #ifdef USE_PARAMCACHE fval = func(x); #else - fval = func(x, p); + fval = func(x, p); #endif } @@ -1151,9 +1153,10 @@ void FitUtil::EvaluateLogLGradient(const IModelFunction &f, const UnBinData &dat std::vector<double> pointContribution(npar); - const double * x = nullptr; + const double * x = nullptr; + std::vector<double> xc; if (data.NDim() > 1) { - std::vector<double> xc(data.NDim()); + xc.resize(data.NDim() ); for (unsigned int j = 0; j < data.NDim(); ++j) xc[j] = *data.GetCoordComponent(i, j); x = xc.data(); @@ -1163,6 +1166,17 @@ void FitUtil::EvaluateLogLGradient(const IModelFunction &f, const UnBinData &dat double fval = func(x, p); func.ParameterGradient(x, p, &gradFunc[0]); + +#ifdef DEBUG + { + R__LOCKGUARD(gROOTMutex); + if (i < 5 || (i > data.Size()-5) ) { + if (data.NDim() > 1) std::cout << i << " x " << x[0] << " y " << x[1] << " func " << fval + << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl; + else std::cout << i << " x " << x[0] << " gradient " << gradFunc[0] << " " << gradFunc[1] << " " << gradFunc[3] << std::endl; + } + } +#endif for (unsigned int kpar = 0; kpar < npar; ++kpar) { if (fval > 0) diff --git a/math/mathcore/test/testGradient.cxx b/math/mathcore/test/testGradient.cxx index 169455c59d619edddc6eda7518232b3fcf2cd593..d080cbd7b3613cc6788dc8a25961550e6f381b79 100644 --- a/math/mathcore/test/testGradient.cxx +++ b/math/mathcore/test/testGradient.cxx @@ -14,6 +14,8 @@ #include "HFitInterface.h" #include "TF1.h" #include "TH1.h" +#include "TF2.h" +#include "TH2.h" #include "TRandom3.h" @@ -32,10 +34,11 @@ // "Scalar", "Vectorial") // PolicyInfoStr points to a human-readable string describing // ExecutionPolicyType (e.g., "Serial", "Multithread") -template <typename U, ROOT::Fit::ExecutionPolicy V, const char *dataInfoStr, const char *policyInfoStr> +template <typename U, ROOT::Fit::ExecutionPolicy V, int W, const char *dataInfoStr, const char *policyInfoStr> struct GradientTestTraits { using DataType = U; static constexpr ROOT::Fit::ExecutionPolicy ExecutionPolicyType() { return V; }; + static constexpr int Dimensions() { return W; }; static void PrintTypeInfo(const std::string &fittingInfo) { @@ -43,6 +46,7 @@ struct GradientTestTraits { std::cout << "- Fitting type: " << fittingInfo << std::endl; std::cout << "- Data type: " << dataInfoStr << std::endl; std::cout << "- Execution policy: " << policyInfoStr << std::endl; + std::cout << "- Dimensions: " << Dimensions() << std::endl; std::cout << "-------------------------------------------" << std::endl; } }; @@ -57,51 +61,55 @@ char mthreadStr[] = "Multithread"; // Typedefs of GradientTestTraits for scalar (serial and multithreaded) // scenarios -using ScalarSerial = GradientTestTraits<Double_t, ROOT::Fit::ExecutionPolicy::kSerial, scalarStr, serialStr>; -using ScalarMultithread = GradientTestTraits<Double_t, ROOT::Fit::ExecutionPolicy::kMultithread, scalarStr, mthreadStr>; +using ScalarSerial1D = GradientTestTraits<Double_t, ROOT::Fit::ExecutionPolicy::kSerial, 1, scalarStr, serialStr>; +using ScalarMultithread1D = + GradientTestTraits<Double_t, ROOT::Fit::ExecutionPolicy::kMultithread, 1, scalarStr, mthreadStr>; +using ScalarSerial2D = GradientTestTraits<Double_t, ROOT::Fit::ExecutionPolicy::kSerial, 2, scalarStr, serialStr>; +using ScalarMultithread2D = + GradientTestTraits<Double_t, ROOT::Fit::ExecutionPolicy::kMultithread, 2, scalarStr, mthreadStr>; #ifdef R__HAS_VECCORE // Typedefs of GradientTestTraits for vectorial (serial and multithreaded) // scenarios -using VectorialSerial = GradientTestTraits<ROOT::Double_v, ROOT::Fit::ExecutionPolicy::kSerial, vectorStr, serialStr>; -using VectorialMultithread = - GradientTestTraits<ROOT::Double_v, ROOT::Fit::ExecutionPolicy::kMultithread, vectorStr, mthreadStr>; +using VectorialSerial1D = + GradientTestTraits<ROOT::Double_v, ROOT::Fit::ExecutionPolicy::kSerial, 1, vectorStr, serialStr>; +using VectorialMultithread1D = + GradientTestTraits<ROOT::Double_v, ROOT::Fit::ExecutionPolicy::kMultithread, 1, vectorStr, mthreadStr>; +using VectorialSerial2D = + GradientTestTraits<ROOT::Double_v, ROOT::Fit::ExecutionPolicy::kSerial, 2, vectorStr, serialStr>; +using VectorialMultithread2D = + GradientTestTraits<ROOT::Double_v, ROOT::Fit::ExecutionPolicy::kMultithread, 2, vectorStr, mthreadStr>; #endif -// Model function to test the gradient evaluation +// Interface abstracting the model function and its related data (number of parameters, parameters, the model function +// itself...) template <class T> -static T modelFunction(const T *data, const double *params) -{ - return params[0] * exp(-(*data + (-130.)) * (*data + (-130.)) / 2) + - params[1] * exp(-(params[2] * (*data * (0.01)) - params[3] * ((*data) * (0.01)) * ((*data) * (0.01)))); -} +struct Model { + virtual void FillModelData(ROOT::Fit::BinData *&data) = 0; + virtual void FillModelData(ROOT::Fit::UnBinData *&data) = 0; -// Helper class used to encapsulate the calls to the gradient interfaces, -// templated with a GradientTestTraits type -template <class T, class U> -struct GradientTestEvaluation { - // Types to instantiate Chi2FCN: the gradient function interface handles the - // parameters, so its base typedef - // has to be always a double; the parametric function interface is templated - // to test both serial and vectorial - // cases. - using GradFunctionType = ROOT::Math::IGradientFunctionMultiDimTempl<double>; - using BaseFunctionType = ROOT::Math::IParamMultiFunctionTempl<typename T::DataType>; + TF1 *fModelFunction; + unsigned fNumParams; + const Double_t *fParams; +}; - GradientTestEvaluation() +// Class storing a 1D model function, its parameters and a histogram with data from the function +template <class T> +struct Model1D : public Model<T> { + Model1D() { + Model<T>::fNumParams = fNumParams; + Model<T>::fParams = fParams; + // Create TF1 from model function and initialize the fit function std::stringstream streamTF1; streamTF1 << "f" << this; std::string nameTF1 = streamTF1.str(); - fModelFunction = new TF1(nameTF1.c_str(), modelFunction<typename T::DataType>, 100, 200, 4); - fModelFunction->SetParameters(fParams); - - fFitFunction = - new ROOT::Math::WrappedMultiTF1Templ<typename T::DataType>(*fModelFunction, fModelFunction->GetNdim()); + Model<T>::fModelFunction = new TF1(nameTF1.c_str(), Function, 100, 200, 4); + Model<T>::fModelFunction->SetParameters(Model<T>::fParams); // Assure the to-be-created histogram does not replace an old one std::stringstream streamTH1; @@ -114,33 +122,162 @@ struct GradientTestEvaluation { // Create TH1 and fill it with values from model function fNumPoints = 12801; + //fNumPoints = 11; fHistogram = new TH1D(nameTH1.c_str(), "Test random numbers", fNumPoints, 100, 200); gRandom->SetSeed(1); fHistogram->FillRandom(nameTF1.c_str(), 1000000); + } - // Fill (binned or unbinned) data - FillData(); + static T Function(const T *data, const double *params) + { + return params[0] * exp(-(*data + (-130.)) * (*data + (-130.)) / 2) + + params[1] * exp(-(params[2] * (*data * (0.01)) - params[3] * ((*data) * (0.01)) * ((*data) * (0.01)))); } // Fill binned data - template <class FitType = U> - typename std::enable_if<std::is_same<FitType, U>::value && std::is_same<FitType, ROOT::Fit::BinData>::value>::type - FillData() + void FillModelData(ROOT::Fit::BinData *&data) override { - fData = new ROOT::Fit::BinData(); - ROOT::Fit::FillData(*fData, fHistogram, fModelFunction); + data = new ROOT::Fit::BinData(fNumPoints, 1); + ROOT::Fit::FillData(*data, fHistogram, Model<T>::fModelFunction); } // Fill unbinned data - template <class FitType = U> - typename std::enable_if<std::is_same<FitType, U>::value && std::is_same<FitType, ROOT::Fit::UnBinData>::value>::type - FillData() + void FillModelData(ROOT::Fit::UnBinData *&data) override { - fData = new ROOT::Fit::UnBinData(fNumPoints); + data = new ROOT::Fit::UnBinData(fNumPoints, 1); TAxis *coords = fHistogram->GetXaxis(); for (unsigned i = 0; i < fNumPoints; i++) - fData->Add(coords->GetBinCenter(i)); + data->Add(coords->GetBinCenter(i)); + } + + TH1D *fHistogram; + static const unsigned fNumParams = 4; + const Double_t fParams[fNumParams] = {1, 1000, 7.5, 1.5}; + unsigned fNumPoints; +}; + +// Class storing a DD model function, its parameters and a histogram with data from the function +template <class T> +struct Model2D : public Model<T> { + Model2D() + { + Model<T>::fNumParams = fNumParams; + Model<T>::fParams = fParams; + + // Create TF1 from model function and initialize the fit function + std::stringstream streamTF2; + streamTF2 << "f" << this; + std::string nameTF2 = streamTF2.str(); + + Model<T>::fModelFunction = new TF2(nameTF2.c_str(), Function, -5, 5, -5, 5, fNumParams); + Model<T>::fModelFunction->SetParameters(Model<T>::fParams); + + // Assure the to-be-created histogram does not replace an old one + std::stringstream streamTH2; + streamTH2 << "h" << this; + std::string nameTH2 = streamTH2.str(); + + auto oldTH2 = gROOT->FindObject(nameTH2.c_str()); + if (oldTH2) + delete oldTH2; + + // Create TH1 and fill it with values from model function + fNumPoints = 801; + fHistogram = new TH2D(nameTH2.c_str(), "Test random numbers", fNumPoints, -5, 5, fNumPoints, -5, 5); + gRandom->SetSeed(1); + fHistogram->FillRandom(nameTF2.c_str(), 1000); + } + +#ifdef R__HAS_VECCORE + static T TemplatedGaus(T x, Double_t mean, Double_t sigma, Bool_t norm = false) + { + if (sigma == 0) + return 1.e30; + + T arg = (x - mean) / sigma; + + // for |arg| > 39 result is zero in double precision + vecCore::Mask_v<T> mask = !(arg < -39.0 || arg > 39.0); + + // Initialize the result to 0.0 + T res(0.0); + + // Compute the function only when the arg meets the criteria, using the mask computed before + vecCore::MaskedAssign<T>(res, mask, vecCore::math::Exp(-0.5 * arg * arg)); + + if (!norm) + return res; + + return res / (2.50662827463100024 * sigma); // sqrt(2*Pi)=2.50662827463100024 + } + + static T Function(const T *data, const double *params) + { + return params[0] * TemplatedGaus(data[0], params[1], params[2]) * TemplatedGaus(data[1], params[3], params[4]); + } + +#else + + static T Function(const T *data, const double *params) + { + return params[0] * TMath::Gaus(data[0], params[1], params[2]) * TMath::Gaus(data[1], params[3], params[4]); + } +#endif + + // Fill binned data + void FillModelData(ROOT::Fit::BinData *&data) override + { + data = new ROOT::Fit::BinData(fNumPoints, 2); + ROOT::Fit::FillData(*data, fHistogram, Model<T>::fModelFunction); + } + + // Fill unbinned data + void FillModelData(ROOT::Fit::UnBinData *&data) override + { + data = new ROOT::Fit::UnBinData(fNumPoints, 2); + + TAxis *xCoords = fHistogram->GetXaxis(); + TAxis *yCoords = fHistogram->GetYaxis(); + for (unsigned i = 0; i < fNumPoints; i++) + data->Add(xCoords->GetBinCenter(i), yCoords->GetBinCenter(i)); + } + + TH2D *fHistogram; + static const unsigned fNumParams = 5; + const Double_t fParams[fNumParams] = {300, 0., 2., 0., 3.}; + unsigned fNumPoints; +}; + +// Helper class used to encapsulate the calls to the gradient interfaces, +// templated with a GradientTestTraits type +template <class T, class U> +struct GradientTestEvaluation { + // Types to instantiate Chi2FCN: the gradient function interface handles the + // parameters, so its base typedef + // has to be always a double; the parametric function interface is templated + // to test both serial and vectorial + // cases. + using GradFunctionType = ROOT::Math::IGradientFunctionMultiDimTempl<double>; + using BaseFunctionType = ROOT::Math::IParamMultiFunctionTempl<typename T::DataType>; + + // Basic type to compare against + using ScalarSerial = + GradientTestTraits<double, ROOT::Fit::ExecutionPolicy::kSerial, T::Dimensions(), scalarStr, serialStr>; + + GradientTestEvaluation() + { + if (T::Dimensions() == 1) + fModel = new Model1D<typename T::DataType>(); + else if (T::Dimensions() == 2) + fModel = new Model2D<typename T::DataType>(); + + fNumParams = fModel->fNumParams; + fFitFunction = new ROOT::Math::WrappedMultiTF1Templ<typename T::DataType>(*(fModel->fModelFunction), + fModel->fModelFunction->GetNdim()); + + // Fill (binned or unbinned) data + fModel->FillModelData(fData); } virtual void SetFitter() = 0; @@ -151,22 +288,25 @@ struct GradientTestEvaluation { start = std::chrono::system_clock::now(); for (int i = 0; i < fNumRepetitions; i++) - fFitter->Gradient(fParams, solution); + fFitter->Gradient(fModel->fParams, solution); end = std::chrono::system_clock::now(); + // std::cout << "Gradient is : " << fFitter->NDim() << " "; + // for (unsigned int i = 0; i < fNumParams ; ++i) + // std::cout << " " << solution[i]; + // std::cout << std::endl; + + std::chrono::duration<Double_t> timeElapsed = end - start; return timeElapsed.count() / fNumRepetitions; + } - static const int fNumParams = 4; - static const int fNumRepetitions = 100; - - const Double_t fParams[fNumParams] = {1, 1000, 7.5, 1.5}; - unsigned int fNumPoints; + static const int fNumRepetitions = 2; - TF1 *fModelFunction; - TH1D *fHistogram; + unsigned fNumParams; + Model<typename T::DataType> *fModel; ROOT::Math::WrappedMultiTF1Templ<typename T::DataType> *fFitFunction; U *fData; ROOT::Fit::BasicFCN<GradFunctionType, BaseFunctionType, U> *fFitter; @@ -176,6 +316,7 @@ template <class T> struct Chi2GradientTestEvaluation : public GradientTestEvaluation<T, ROOT::Fit::BinData> { using typename GradientTestEvaluation<T, ROOT::Fit::BinData>::GradFunctionType; using typename GradientTestEvaluation<T, ROOT::Fit::BinData>::BaseFunctionType; + using typename GradientTestEvaluation<T, ROOT::Fit::BinData>::ScalarSerial; Chi2GradientTestEvaluation() { SetFitter(); } @@ -190,6 +331,7 @@ template <class T> struct PoissonLikelihoodGradientTestEvaluation : public GradientTestEvaluation<T, ROOT::Fit::BinData> { using typename GradientTestEvaluation<T, ROOT::Fit::BinData>::GradFunctionType; using typename GradientTestEvaluation<T, ROOT::Fit::BinData>::BaseFunctionType; + using typename GradientTestEvaluation<T, ROOT::Fit::BinData>::ScalarSerial; PoissonLikelihoodGradientTestEvaluation() { SetFitter(); } @@ -204,6 +346,7 @@ template <class T> struct LogLikelihoodGradientTestEvaluation : public GradientTestEvaluation<T, ROOT::Fit::UnBinData> { using typename GradientTestEvaluation<T, ROOT::Fit::UnBinData>::GradFunctionType; using typename GradientTestEvaluation<T, ROOT::Fit::UnBinData>::BaseFunctionType; + using typename GradientTestEvaluation<T, ROOT::Fit::UnBinData>::ScalarSerial; LogLikelihoodGradientTestEvaluation() { SetFitter(); } @@ -221,6 +364,8 @@ struct LogLikelihoodGradientTestEvaluation : public GradientTestEvaluation<T, RO // type. template <class T> class Chi2GradientTest : public ::testing::Test, public Chi2GradientTestEvaluation<T> { + using typename Chi2GradientTestEvaluation<T>::ScalarSerial; + protected: virtual void SetUp() { @@ -244,6 +389,8 @@ protected: // type. template <class T> class PoissonLikelihoodGradientTest : public ::testing::Test, public PoissonLikelihoodGradientTestEvaluation<T> { + using typename PoissonLikelihoodGradientTestEvaluation<T>::ScalarSerial; + protected: virtual void SetUp() { @@ -267,6 +414,8 @@ protected: // type. template <class T> class LogLikelihoodGradientTest : public ::testing::Test, public LogLikelihoodGradientTestEvaluation<T> { + using typename LogLikelihoodGradientTestEvaluation<T>::ScalarSerial; + protected: virtual void SetUp() { @@ -286,22 +435,25 @@ protected: // Types used by Google Test to instantiate the tests. #ifdef R__HAS_VECCORE # ifdef R__USE_IMT -typedef ::testing::Types<ScalarMultithread, VectorialSerial, VectorialMultithread> TestTypes; +typedef ::testing::Types<ScalarMultithread1D, VectorialSerial1D, VectorialMultithread1D, ScalarMultithread2D, + VectorialSerial2D, VectorialMultithread2D> + TestTypes; # else -typedef ::testing::Types<ScalarSerial, VectorialSerial, VectorialMultithread> TestTypes; +typedef ::testing::Types<VectorialSerial1D, VectorialSerial2D> TestTypes; # endif #else # ifdef R__USE_IMT -typedef ::testing::Types<ScalarSerial> TestTypes; +typedef ::testing::Types<ScalarMultithread1D, ScalarMultithread2D> TestTypes; # else -typedef ::testing::Types<ScalarSerial> TestTypes; +typedef ::testing::Types<> TestTypes; # endif #endif -TYPED_TEST_CASE(Chi2GradientTest, TestTypes); + +TYPED_TEST_CASE(LogLikelihoodGradientTest, TestTypes); // Test EvalChi2Gradient and outputs its speedup against the scalar serial case. -TYPED_TEST(Chi2GradientTest, Chi2Gradient) +TYPED_TEST(LogLikelihoodGradientTest, LogLikelihoodGradient) { Double_t solution[TestFixture::fNumParams]; @@ -316,10 +468,11 @@ TYPED_TEST(Chi2GradientTest, Chi2Gradient) } } -TYPED_TEST_CASE(PoissonLikelihoodGradientTest, TestTypes); + +TYPED_TEST_CASE(Chi2GradientTest, TestTypes); // Test EvalChi2Gradient and outputs its speedup against the scalar serial case. -TYPED_TEST(PoissonLikelihoodGradientTest, PoissonLikelihoodGradient) +TYPED_TEST(Chi2GradientTest, Chi2Gradient) { Double_t solution[TestFixture::fNumParams]; @@ -334,10 +487,10 @@ TYPED_TEST(PoissonLikelihoodGradientTest, PoissonLikelihoodGradient) } } -TYPED_TEST_CASE(LogLikelihoodGradientTest, TestTypes); +TYPED_TEST_CASE(PoissonLikelihoodGradientTest, TestTypes); // Test EvalChi2Gradient and outputs its speedup against the scalar serial case. -TYPED_TEST(LogLikelihoodGradientTest, LogLikelihoodGradient) +TYPED_TEST(PoissonLikelihoodGradientTest, PoissonLikelihoodGradient) { Double_t solution[TestFixture::fNumParams]; @@ -351,3 +504,4 @@ TYPED_TEST(LogLikelihoodGradientTest, LogLikelihoodGradient) EXPECT_NEAR(solution[i], TestFixture::fReferenceSolution[i], 1e-6); } } +