From 4d489a740b5c0c49f15afe9711c16d0a51f2706a Mon Sep 17 00:00:00 2001
From: Lorenzo Moneta <Lorenzo.Moneta@cern.ch>
Date: Fri, 1 Nov 2013 14:32:32 +0100
Subject: [PATCH] Fix math core fit tests
---
math/mathcore/test/fit/WrapperRooPdf.h | 5 +-
math/mathcore/test/fit/testFitPerf.cxx | 9 ++--
math/mathcore/test/fit/testMinim.cxx | 62 +++++++++++++++++--------
math/mathcore/test/fit/testRooFit.cxx | 64 +++++++++++++++++---------
4 files changed, 96 insertions(+), 44 deletions(-)
diff --git a/math/mathcore/test/fit/WrapperRooPdf.h b/math/mathcore/test/fit/WrapperRooPdf.h
index 978b7c9f72d..2b8e07c5bf0 100644
--- a/math/mathcore/test/fit/WrapperRooPdf.h
+++ b/math/mathcore/test/fit/WrapperRooPdf.h
@@ -36,8 +36,10 @@ public:
fParams = fPdf->getParameters(obsList);
assert(fX!=0);
assert(fParams!=0);
+#ifdef DEBUG
fX->Print("v");
fParams->Print("v");
+#endif
}
@@ -58,9 +60,10 @@ public:
fParams = fPdf->getParameters(obsList);
assert(fX!=0);
assert(fParams!=0);
+#ifdef DEBUG
fX->Print("v");
fParams->Print("v");
-
+#endif
// // iterate on fX
// TIterator* itr = fX->createIterator() ;
// RooAbsArg* arg = 0;
diff --git a/math/mathcore/test/fit/testFitPerf.cxx b/math/mathcore/test/fit/testFitPerf.cxx
index 76a35e3b319..2ee6b0d3804 100644
--- a/math/mathcore/test/fit/testFitPerf.cxx
+++ b/math/mathcore/test/fit/testFitPerf.cxx
@@ -655,6 +655,9 @@ int FitUsingRooFit2(TTree * tree) {
else if (j > 1) {
pdf[j] = new RooProdPdf(pname.c_str(),pname.c_str(),RooArgSet(*g[j],*pdf[j-1]) );
}
+// else
+// pdf[0] = g[0];
+
}
@@ -667,7 +670,7 @@ int FitUsingRooFit2(TTree * tree) {
(pdf[N-1]->getVariables())->Print("v"); // print the parameters
std::cout << "\n\nDo the fit now \n\n";
#else
- int level = -1;
+ int level = -1;
bool save = false;
#endif
@@ -685,7 +688,7 @@ int FitUsingRooFit2(TTree * tree) {
#endif
- iret |= (result == 0);
+ iret |= (result != 0);
// free
for (int j = 0; j < N; ++j) {
@@ -693,7 +696,7 @@ int FitUsingRooFit2(TTree * tree) {
delete m[j];
delete s[j];
delete g[j];
- delete pdf[j];
+ if (j> 0) delete pdf[j];
}
if (iret != 0) return iret;
diff --git a/math/mathcore/test/fit/testMinim.cxx b/math/mathcore/test/fit/testMinim.cxx
index 2bf88a9c193..19117cd66e7 100644
--- a/math/mathcore/test/fit/testMinim.cxx
+++ b/math/mathcore/test/fit/testMinim.cxx
@@ -19,17 +19,19 @@
#include "TRandom3.h"
#include "TMath.h"
+#include "RVersion.h"
+
//#define DEBUG
int gNCall = 0;
int gNCall2 = 0;
-int gNmin = 1000;
-int gVerbose = 0;
+int gNmin = 1;
+int gVerbose = 1;
bool minos = false;
-double gAbsTolerance = 0.001;
+double gAbsTolerance = 5.E-6; // otherwise gsl_conjugate_PR fails
// Rosenbrok function to be minimize
@@ -62,11 +64,15 @@ public :
return new RosenBrockFunction();
}
- const double * TrueMinimum() const {
+ const double * TrueXMinimum() const {
fTrueMin[0] = 1;
fTrueMin[1] = 1;
return fTrueMin;
}
+
+ double TrueMinimum() const {
+ return 0;
+ }
private:
@@ -153,10 +159,11 @@ public :
}
- const double * TrueMinimum() const {
+ const double * TrueXMinimum() const {
return x0.GetMatrixArray();
}
+ double TrueMinimum() const { return 0; }
void Gradient (const double * x, double * g) const {
gNCall2++;
@@ -270,10 +277,12 @@ public :
return new ChebyQuadFunction(*this);
}
- const double * TrueMinimum() const {
+ const double * TrueXMinimum() const {
return &fTrueMin[0];
}
+ double TrueMinimum() const { return 0; }
+
// use equally spaced points
void StartPoints(double * x, double * s) {
for (unsigned int i = 0; i < fDim; ++i) {
@@ -367,14 +376,24 @@ public :
-const double * TrueMinimum(const ROOT::Math::IMultiGenFunction & func) {
+const double * TrueXMinimum(const ROOT::Math::IMultiGenFunction & func) {
+
+ const RosenBrockFunction * fRB = dynamic_cast< const RosenBrockFunction *> (&func);
+ if (fRB != 0) return fRB->TrueXMinimum();
+ const TrigoFletcherFunction * fTF = dynamic_cast< const TrigoFletcherFunction *> (&func);
+ if (fTF != 0) return fTF->TrueXMinimum();
+// const ChebyQuadFunction * fCQ = dynamic_cast< const ChebyQuadFunction *> (&func);
+// if (fCQ != 0) return fCQ->TrueXMinimum();
+ return 0;
+}
+double TrueMinimum(const ROOT::Math::IMultiGenFunction & func) {
const RosenBrockFunction * fRB = dynamic_cast< const RosenBrockFunction *> (&func);
if (fRB != 0) return fRB->TrueMinimum();
const TrigoFletcherFunction * fTF = dynamic_cast< const TrigoFletcherFunction *> (&func);
if (fTF != 0) return fTF->TrueMinimum();
// const ChebyQuadFunction * fCQ = dynamic_cast< const ChebyQuadFunction *> (&func);
-// if (fCQ != 0) return fCQ->TrueMinimum();
+// if (fCQ != 0) return fCQ->TrueXMinimum();
return 0;
}
@@ -399,10 +418,11 @@ int DoNewMinimization( const ROOT::Math::IMultiGenFunction & func, const double
if (func.NDim() >= 10) {
min->SetMaxFunctionCalls(1000000);
min->SetMaxIterations(100000);
- min->SetTolerance(0.001);
- if (func.NDim() >= 10) min->SetTolerance(0.01);
-
+ min->SetTolerance(0.01);
}
+ else
+ min->SetTolerance(0.001);
+
min->SetPrintLevel(gVerbose);
// check if func provides gradient
@@ -428,10 +448,11 @@ int DoNewMinimization( const ROOT::Math::IMultiGenFunction & func, const double
const double * xmin = min->X();
bool ok = true;
- const double * trueMin = TrueMinimum(func);
+ const double * trueMin = TrueXMinimum(func);
if (trueMin != 0) {
+ ok &= (std::fabs(minval - TrueMinimum(func) ) < gAbsTolerance );
for (unsigned int i = 0; i < func.NDim(); ++i)
- ok &= (std::fabs(xmin[i]-trueMin[i] ) < gAbsTolerance);
+ ok &= (std::fabs(xmin[i]-trueMin[i] ) < sqrt(gAbsTolerance));
}
if (!ok) iret = -2;
@@ -487,13 +508,15 @@ int DoOldMinimization( FCN func, TVirtualFitter * min, double &minval, double &
min->GetParameter(0,parName,parx,we,al,bl);
min->GetParameter(1,parName,pary,we,al,bl);
- bool ok = ( TMath::Abs(parx-1.) < gAbsTolerance &&
- TMath::Abs(pary-1.) < gAbsTolerance );
-
+ bool ok = ( TMath::Abs(parx-1.) < sqrt(gAbsTolerance) &&
+ TMath::Abs(pary-1.) < sqrt(gAbsTolerance) );
+
+
double errdef = 0;
int nvpar, nparx;
min->GetStats(minval,edm,errdef,nvpar,nparx);
+
if (!ok) iret = -2;
min->Clear(); // for further use
@@ -521,7 +544,6 @@ int testNewMinimizer( const ROOT::Math::IMultiGenFunction & func, const double *
return -1;
}
-
for (int i = 0; i < gNmin; ++i) {
gNCall = 0; gNCall2 = 0;
iret |= DoNewMinimization(func, x0, s0, min,minval,edm,&xmin[0]);
@@ -599,10 +621,14 @@ int testRosenBrock() {
// minimize using Rosenbrock Function
#ifndef DEBUG
- gNmin = 10000;
+ gNmin = 1;
#endif
+
+
+#if ROOT_VERSION_CODE < ROOT_VERSION(5,99,00)
iret |= testOldMinimizer(RosenBrock,"Minuit",2);
iret |= testOldMinimizer(RosenBrock,"Minuit2",2);
+#endif
RosenBrockFunction fRB;
double xRB[2] = { -1.,1.2};
diff --git a/math/mathcore/test/fit/testRooFit.cxx b/math/mathcore/test/fit/testRooFit.cxx
index 38a1cae57f7..4a7d0f993ef 100644
--- a/math/mathcore/test/fit/testRooFit.cxx
+++ b/math/mathcore/test/fit/testRooFit.cxx
@@ -37,11 +37,12 @@
#include "Math/WrappedParamFunction.h"
#include <cmath>
-const int N = 1; // n must be greater than 1
+const int N = 3; // n must be greater than 1
const int nfit = 1;
const int nEvents = 10000;
double iniPar[2*N];
+//#define DEBUG
typedef ROOT::Math::IParamMultiFunction Func;
@@ -93,10 +94,13 @@ void FillUnBinData(ROOT::Fit::UnBinData &d, TTree * tree ) {
for (int j = 0; j < N; ++j)
m[j] += vx[j];
}
+
+#ifdef DEBUG
std::cout << "average values of means :\n";
for (int j = 0; j < N; ++j)
std::cout << m[j]/n << " ";
std::cout << "\n";
+#endif
return;
@@ -134,14 +138,14 @@ class MultiGaussRooPdf {
g[j] = new RooGaussian(gname.c_str(),"gauss(x,mean,sigma)",*x[j],*m[j],*s[j]);
std::string pname = "prod_" + ROOT::Math::Util::ToString(j);
- if (j == 1) {
+ if (j == 0)
+ pdf[0] = g[0];
+ else if (j == 1) {
pdf[1] = new RooProdPdf(pname.c_str(),pname.c_str(),RooArgSet(*g[1],*g[0]) );
}
else if (j > 1) {
pdf[j] = new RooProdPdf(pname.c_str(),pname.c_str(),RooArgSet(*g[j],*pdf[j-1]) );
}
- else if (j == 0)
- pdf[0] = g[0];
}
@@ -164,7 +168,7 @@ class MultiGaussRooPdf {
delete m[j];
delete s[j];
delete g[j];
- delete pdf[j];
+ if (j> 0) delete pdf[j]; // no pdf allocated for j = 0
}
}
@@ -199,13 +203,13 @@ int FitUsingRooFit(TTree & tree, RooAbsPdf & pdf, RooArgSet & xvars) {
#ifdef DEBUG
- int level = 3;
+ int level = 2;
std::cout << "num entries = " << data.numEntries() << std::endl;
bool save = true;
- (pdf[N-1]->getVariables())->Print("v"); // print the parameters
+ pdf.getVariables()->Print("v"); // print the parameters
std::cout << "\n\nDo the fit now \n\n";
#else
- int level = 1;
+ int level = -1;
bool save = false;
#endif
@@ -226,6 +230,7 @@ int FitUsingRooFit(TTree & tree, RooAbsPdf & pdf, RooArgSet & xvars) {
w.Stop();
+ std::cout << "RooFit result " << std::endl;
RooArgSet * params = pdf.getParameters(xvars);
params->Print("v");
delete params;
@@ -248,7 +253,7 @@ int DoFit(TTree * tree, Func & func, bool debug = false, bool = false ) {
//std::cout << "Fit parameter 2 " << f.Parameters()[2] << std::endl;
ROOT::Fit::Fitter fitter;
- fitter.Config().MinimizerOptions().SetPrintLevel(2);
+ fitter.Config().MinimizerOptions().SetPrintLevel(0);
fitter.Config().SetMinimizer(MinType::name().c_str(),MinType::name2().c_str());
fitter.Config().MinimizerOptions().SetTolerance(1.); // to be consistent with RooFit
@@ -267,8 +272,7 @@ int DoFit(TTree * tree, Func & func, bool debug = false, bool = false ) {
std::cout << " Fit Failed " << std::endl;
return -1;
}
- if (debug)
- fitter.Result().Print(std::cout);
+ fitter.Result().Print(std::cout);
return 0;
}
@@ -301,13 +305,29 @@ int FitUsingNewFitter(FitObj * fitobj, Func & func, bool useGrad=false) {
return iret;
}
-double gausnorm(const double *x, const double *p) {
- //return p[0]*TMath::Gaus(x[0],p[1],p[2]);
- double invsig = 1./p[1];
- double tmp = (x[0]-p[0]) * invsig;
- const double sqrt_2pi = 1./std::sqrt(2.* 3.14159 );
- return std::exp(-0.5 * tmp*tmp ) * sqrt_2pi * invsig;
-}
+template <int N>
+struct GausNorm {
+ static double F(const double *x, const double *p) {
+ return ROOT::Math::normal_pdf(x[N-1], p[2*N-1], p[2*N-2] ) * GausNorm<N-1>::F(x,p);
+ }
+};
+template <>
+struct GausNorm<1> {
+
+ static double F(const double *x, const double *p) {
+ return ROOT::Math::normal_pdf(x[0], p[1], p[0] );
+ }
+};
+
+
+// double gausnorm(
+// return ROOT::Math::normal_pdf(x[0], p[1], p[0] );
+// // //return p[0]*TMath::Gaus(x[0],p[1],p[2]);
+// // double invsig = 1./p[1];
+// // double tmp = (x[0]-p[0]) * invsig;
+// // const double sqrt_2pi = 1./std::sqrt(2.* 3.14159 );
+// // return std::exp(-0.5 * tmp*tmp ) * sqrt_2pi * invsig;
+// }
@@ -346,10 +366,10 @@ int main() {
// in case of N = 1 do also a simple gauss fit
// using TF1 gausN
- if (N == 1) {
- ROOT::Math::WrappedParamFunction<> gausn(&gausnorm,2,iniPar,iniPar+1);
- FitUsingNewFitter<TMINUIT>(&tree,gausn);
- }
+// if (N == 1) {
+ ROOT::Math::WrappedParamFunction<> gausn(&GausNorm<N>::F,2*N,iniPar,iniPar+2*N);
+ FitUsingNewFitter<TMINUIT>(&tree,gausn);
+ //}
--
GitLab