portelli/LatAnalyze
1
0
Fork 0
mirror of https://github.com/aportelli/LatAnalyze.git synced 2025-04-11 03:20:46 +01:00
Code Releases Activity

simplification of the Minuit minimiser interface, (successful) test of the new fit interface

Browse Source
This commit is contained in:
Antonin Portelli 2016-03-23 17:08:25 +00:00
parent ee8ed05b81
commit 9b9c86cf72
5 changed files with 168 additions and 185 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
examples/exFit.cpp
View File

@ -1,6 +1,7 @@
#include <iostream> #include <iostream>
#include <cmath> #include <cmath>
#include <LatAnalyze/CompiledModel.hpp> #include <LatAnalyze/CompiledModel.hpp>
#include <LatAnalyze/Io.hpp>
#include <LatAnalyze/MinuitMinimizer.hpp> #include <LatAnalyze/MinuitMinimizer.hpp>
#include <LatAnalyze/RandGen.hpp> #include <LatAnalyze/RandGen.hpp>
#include <LatAnalyze/XYStatData.hpp> #include <LatAnalyze/XYStatData.hpp>
@ -8,39 +9,50 @@
using namespace std; using namespace std;
using namespace Latan; using namespace Latan;
const Index nPoint = 30; const Index nPoint1 = 5, nPoint2 = 5;
const double xErr = .2, yErr = .1; const double xErr = .1, yErr = .1;
const double exactPar[2] = {0.5,5.0}, dx = 10.0/static_cast<double>(nPoint); const double exactPar[2] = {0.5,5.};
const double dx1 = 10.0/static_cast<double>(nPoint1);
const double dx2 = 5.0/static_cast<double>(nPoint2);
int main(void) int main(void)
{ {
// generate fake data // generate fake data
XYStatData data; XYStatData data;
RandGen rg; RandGen rg;
double x1_k, x2_k; double xBuf[2];
DoubleModel f([](const double *x, const double *p) DoubleModel f([](const double *x, const double *p)
{return p[1]*exp(-x[0]*p[0]);}, 1, 2); {return p[1]*exp(-x[0]*p[0])+x[1];}, 2, 2);
data.addXDim("x", nPoint); data.addXDim("x", nPoint1);
data.addXDim("off", nPoint2);
data.addYDim("y"); data.addYDim("y");
for (Index k = 0; k < nPoint; ++k) for (Index i1 = 0; i1 < nPoint1; ++i1)
{ {
x1_k = rg.gaussian(k*dx, xErr); xBuf[0] = i1*dx1;
x2_k = rg.gaussian(k*dx, xErr); data.x(i1, 0) = rg.gaussian(xBuf[0], xErr);
data.x(k) = x1_k; for (Index i2 = 0; i2 < nPoint2; ++i2)
data.y(k) = rg.gaussian(f(&x2_k, exactPar), yErr); {
printf("% 8e % 8e % 8e % 8e\n", data.x(k), data.y(k), xErr, yErr); xBuf[1] = i2*dx2;
data.x(i2, 1) = xBuf[1];
data.y(data.dataIndex(i1, i2)) = rg.gaussian(f(xBuf, exactPar),
yErr);
printf("% 8e % 8e % 8e % 8e % 8e\n", data.x(i1, 0), xErr,
data.x(i2, 1), data.y(i1), yErr);
}
} }
cout << endl; cout << endl;
data.setXError(0, DVec::Constant(nPoint, xErr)); data.setXError(0, DVec::Constant(data.getXSize(0), xErr));
data.setYError(0, DVec::Constant(nPoint, yErr)); data.assumeXExact(true, 1);
data.setYError(0, DVec::Constant(data.getYSize(), yErr));
cout << data << endl; cout << data << endl;
// fit // fit
DVec init = DVec::Constant(2, 0.5); DVec init = DVec::Constant(2, 0.1);
FitResult p; FitResult p;
MinuitMinimizer minimizer; MinuitMinimizer minimizer;
minimizer.setVerbosity(Minimizer::Verbosity::Normal);
p = data.fit(minimizer, init, f); p = data.fit(minimizer, init, f);
cout << "a= " << p(0) << " b= " << p(1) << endl; cout << "a= " << p(0) << " b= " << p(1) << endl;
cout << "chi^2/ndof= " << p.getChi2PerDof() << endl; cout << "chi^2/ndof= " << p.getChi2PerDof() << endl;

42
examples/exFitSample.cpp
View File

@ -8,41 +8,49 @@
using namespace std; using namespace std;
using namespace Latan; using namespace Latan;
const Index nPoint = 30; const Index nPoint1 = 10, nPoint2 = 10;
const Index nSample = 1000; const Index nSample = 1000;
const double xErr = .2, yErr = .1; const double xErr = .1, yErr = .1;
const double exactPar[2] = {0.5,5.0}, dx = 10.0/static_cast<double>(nPoint); const double exactPar[2] = {0.5,5.};
const double dx1 = 10.0/static_cast<double>(nPoint1);
const double dx2 = 5.0/static_cast<double>(nPoint2);
int main(void) int main(void)
{ {
// generate fake data // generate fake data
DMatSample x(nSample, nPoint, 1), y(nSample, nPoint, 1);
XYSampleData data(nSample); XYSampleData data(nSample);
RandGen rg; RandGen rg;
double x1_k, x2_k; double xBuf[2];
DoubleModel f([](const double *t, const double *p) DoubleModel f([](const double *x, const double *p)
{return p[1]*exp(-t[0]*p[0]);}, 1, 2); {return p[1]*exp(-x[0]*p[0])+x[1];}, 2, 2);
data.addXDim("x", nPoint); data.addXDim("x", nPoint1);
data.addXDim("off", nPoint2);
data.addYDim("y"); data.addYDim("y");
FOR_STAT_ARRAY(x, s) for (Index s = central; s < nSample; ++s)
{ {
for (Index k = 0; k < nPoint; ++k) for (Index i1 = 0; i1 < nPoint1; ++i1)
{ {
x1_k = rg.gaussian(k*dx, xErr); xBuf[0] = i1*dx1;
x2_k = rg.gaussian(k*dx, xErr); data.x(i1, 0)[s] = rg.gaussian(xBuf[0], xErr);
data.x(k)[s] = x1_k; for (Index i2 = 0; i2 < nPoint2; ++i2)
data.y(k)[s] = rg.gaussian(f(&x2_k, exactPar), yErr); {
xBuf[1] = i2*dx2;
data.x(i2, 1)[s] = xBuf[1];
data.y(data.dataIndex(i1, i2))[s] =
rg.gaussian(f(xBuf, exactPar), yErr);
}
} }
} }
data.assumeXExact(true, 1);
cout << data << endl; cout << data << endl;
// fit // fit
DVec init = DVec::Constant(2, 0.5); DVec init = DVec::Constant(2, 0.1);
DMat err; DMat err;
SampleFitResult p; SampleFitResult p;
MinuitMinimizer minimizer; MinuitMinimizer minimizer;
p = data.fit(minimizer, init, f); p = data.fit(minimizer, init, f);
err = p.variance().cwiseSqrt(); err = p.variance().cwiseSqrt();
cout << "a= " << p[central](0) << " +/- " << err(0) << endl; cout << "a= " << p[central](0) << " +/- " << err(0) << endl;

227
lib/MinuitMinimizer.cpp
View File

@ -1,7 +1,7 @@
/* /*
* MinuitMinimizer.cpp, part of LatAnalyze 3 * MinuitMinimizer.cpp, part of LatAnalyze 3
* *
* Copyright (C) 2013 - 2015 Antonin Portelli * Copyright (C) 2013 - 2016 Antonin Portelli
* *
* LatAnalyze 3 is free software: you can redistribute it and/or modify * LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by * it under the terms of the GNU General Public License as published by
@ -19,48 +19,18 @@
#include <LatAnalyze/MinuitMinimizer.hpp> #include <LatAnalyze/MinuitMinimizer.hpp>
#include <LatAnalyze/includes.hpp> #include <LatAnalyze/includes.hpp>
#include <Minuit2/Minuit2Minimizer.h>
#pragma GCC diagnostic ignored "-Wconversion" #include <Math/Functor.h>
#include <Minuit2/FCNBase.h>
#include <Minuit2/FunctionMinimum.h>
#include <Minuit2/MnMigrad.h>
#include <Minuit2/MnPrint.h>
#include <Minuit2/MnPlot.h>
#include <Minuit2/MnScan.h>
#include <Minuit2/MnSimplex.h>
#include <Minuit2/ScanMinimizer.h>
#include <Minuit2/SimplexMinimizer.h>
#include <Minuit2/VariableMetricMinimizer.h>
using namespace std; using namespace std;
using namespace ROOT;
using namespace Minuit2;
using namespace Latan; using namespace Latan;
static constexpr double initErr = 0.5; static constexpr double initErr = 0.1;
static constexpr unsigned int maxTry = 10u; static constexpr unsigned int maxTry = 10u;
/****************************************************************************** /******************************************************************************
* MinuitMinimizer implementation * * MinuitMinimizer implementation *
******************************************************************************/ ******************************************************************************/
// MinuitFunction constructor //////////////////////////////////////////////////
MinuitMinimizer::MinuitFunction::MinuitFunction(const DoubleFunction &f)
: f_(&f)
{}
// MinuitFunction minuit members ///////////////////////////////////////////////
double MinuitMinimizer::MinuitFunction::operator()
(const vector<double> &x) const
{
return (*f_)(x);
}
double MinuitMinimizer::MinuitFunction::Up(void) const
{
return 1.;
}
// constructors //////////////////////////////////////////////////////////////// // constructors ////////////////////////////////////////////////////////////////
MinuitMinimizer::MinuitMinimizer(const Algorithm algorithm) MinuitMinimizer::MinuitMinimizer(const Algorithm algorithm)
{ {
@ -79,126 +49,135 @@ MinuitMinimizer::Algorithm MinuitMinimizer::getAlgorithm(void) const
return algorithm_; return algorithm_;
} }
double MinuitMinimizer::getPrecision(void) const
{
LATAN_ERROR(Implementation,
"Minuit minimizer precision cannot be accessed");
return 0.;
}
void MinuitMinimizer::setAlgorithm(const Algorithm algorithm) void MinuitMinimizer::setAlgorithm(const Algorithm algorithm)
{ {
algorithm_ = algorithm; algorithm_ = algorithm;
} }
void MinuitMinimizer::setPrecision(const double precision __dumb)
{
LATAN_ERROR(Implementation,
"Minuit minimizer precision cannot be accessed");
}
// minimization //////////////////////////////////////////////////////////////// // minimization ////////////////////////////////////////////////////////////////
const DVec & MinuitMinimizer::operator()(const DoubleFunction &f) const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
{ {
DVec &x = getState(); using namespace ROOT;
Verbosity verbosity = getVerbosity(); using namespace Minuit2;
DVec &x = getState();
int printLevel;
EMinimizerType minuitAlg;
unique_ptr<Math::Minimizer> min;
// convert Latan parameters to Minuit parameters
switch (getVerbosity())
{
case Verbosity::Silent:
printLevel = 0;
break;
case Verbosity::Normal:
printLevel = 2;
break;
case Verbosity::Debug:
printLevel = 3;
break;
}
switch (getAlgorithm())
{
case Algorithm::Migrad:
minuitAlg = kMigrad;
break;
case Algorithm::Simplex:
minuitAlg = kSimplex;
break;
case Algorithm::Combined:
minuitAlg = kCombined;
break;
}
// resize minimizer state to match function number of arguments
if (f.getNArg() != x.size()) if (f.getNArg() != x.size())
{ {
resize(f.getNArg()); resize(f.getNArg());
} }
// set parameters // create and set minimizer
MinuitFunction minuitF(f); min.reset(new Minuit2Minimizer(minuitAlg));
MnUserParameters parameters; min->SetMaxFunctionCalls(getMaxIteration());
min->SetTolerance(getPrecision());
min->SetPrintLevel(printLevel);
// set function and variables
Math::Functor minuitF(f, x.size());
string name;
double val, step;
min->SetFunction(minuitF);
for (Index i = 0; i < x.size(); ++i) for (Index i = 0; i < x.size(); ++i)
{ {
parameters.Add("x_" + strFrom(i), x(i), initErr*fabs(x(i))); name = "x_" + strFrom(i);
if (hasLowLimit(i)&&hasHighLimit(i)) val = x(i);
step = (fabs(x(i)) != 0.) ? initErr*fabs(x(i)) : 1.;
if (hasHighLimit(i) and !hasLowLimit(i))
{ {
parameters.SetLimits(static_cast<unsigned int>(i), min->SetUpperLimitedVariable(i, name, val, step, getHighLimit(i));
getLowLimit(i), getHighLimit(i));
} }
else if (hasLowLimit(i)) else if (!hasHighLimit(i) and hasLowLimit(i))
{ {
parameters.SetLowerLimit(static_cast<unsigned int>(i), min->SetLowerLimitedVariable(i, name, val, step, getLowLimit(i));
getLowLimit(i));
} }
else if (hasHighLimit(i)) else if (hasHighLimit(i) and hasLowLimit(i))
{ {
parameters.SetUpperLimit(static_cast<unsigned int>(i), min->SetLimitedVariable(i, name, val, step, getLowLimit(i),
getHighLimit(i)); getHighLimit(i));
}
else
{
min->SetVariable(i, name, val, step);
} }
} }
// pre-minimization // minimize
MnMigrad preMinimizer(minuitF, parameters, 1); int status;
FunctionMinimum preMin = preMinimizer(); unsigned int n = 0;
if (verbosity >= Verbosity::Debug) if (getVerbosity() >= Verbosity::Normal)
{ {
cout << "-- MINUIT pre-minimization -----------------------------"; cout << "========== Minuit pre-minimization " << endl;
cout << preMin;
cout << "--------------------------------------------------------";
cout << endl;
} }
for (unsigned int i = 0; i < x.size(); ++i) min->SetStrategy(0);
min->Minimize();
do
{ {
parameters.SetValue(i, preMin.UserParameters().Value(i)); n++;
parameters.SetError(i, 2.*preMin.UserParameters().Error(i)); if (getVerbosity() >= Verbosity::Normal)
}
// minimization and output
unique_ptr<MnApplication> minimizer(nullptr);
if (algorithm_ == Algorithm::Migrad)
{
minimizer.reset(new MnMigrad(minuitF, parameters, 2));
}
else if (algorithm_ == Algorithm::Simplex)
{
minimizer.reset(new MnSimplex(minuitF, parameters, 2));
}
unsigned int iTry = 0;
FunctionMinimum min = (*minimizer)();
while ((!min.IsValid())&&(iTry < maxTry))
{
min = (*minimizer)();
iTry++;
}
if (!min.IsValid())
{
LATAN_WARNING("MINUIT library reported that minimization result is not valid");
}
for (unsigned int i = 0; i < x.size(); ++i)
{
x(i) = min.UserParameters().Value(i);
}
if (verbosity >= Verbosity::Normal)
{
cout << "-- MINUIT minimization ---------------------------------";
cout << min;
cout << "--------------------------------------------------------";
cout << endl;
}
if (verbosity >= Verbosity::Debug)
{
vector<pair<double, double>> scanRes;
MnPlot plot;
MnScan scanner(minuitF, preMin.UserParameters(), 2);
cout << "-- MINUIT scan -----------------------------------------";
cout << endl;
for (unsigned int i = 0; i < x.size(); i++)
{ {
cout << "Parameter x_" << i << endl; cout << "========== Minuit minimization, try #" << n << endl;
scanRes = scanner.Scan(i);
plot(scanRes);
} }
cout << "--------------------------------------------------------"; min->SetStrategy(2);
cout << endl; min->Minimize();
status = min->Status();
} while (status and (n < maxTry));
if (getVerbosity() >= Verbosity::Normal)
{
cout << "==============================" << endl;
}
switch (status)
{
case 1:
LATAN_WARNING("invalid minimum: covariance matrix was made positive");
break;
case 2:
LATAN_WARNING("invalid minimum: Hesse analysis is not valid");
break;
case 3:
LATAN_WARNING("invalid minimum: requested precision not reached");
break;
case 4:
LATAN_WARNING("invalid minimum: iteration limit reached");
break;
}
// save and return result
for (Index i = 0; i < x.size(); ++i)
{
x(i) = min->X()[i];
} }
return x; return x;

40
lib/MinuitMinimizer.hpp
View File

@ -1,7 +1,7 @@
/* /*
* MinuitMinimizer.hpp, part of LatAnalyze 3 * MinuitMinimizer.hpp, part of LatAnalyze 3
* *
* Copyright (C) 2013 - 2015 Antonin Portelli * Copyright (C) 2013 - 2016 Antonin Portelli
* *
* LatAnalyze 3 is free software: you can redistribute it and/or modify * LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by * it under the terms of the GNU General Public License as published by
@ -23,52 +23,36 @@
#include <LatAnalyze/Global.hpp> #include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Function.hpp> #include <LatAnalyze/Function.hpp>
#include <LatAnalyze/Minimizer.hpp> #include <LatAnalyze/Minimizer.hpp>
#include <Minuit2/FCNBase.h>
BEGIN_LATAN_NAMESPACE BEGIN_LATAN_NAMESPACE
/****************************************************************************** /******************************************************************************
* Minuit minimizer * * MinuitMinimizer *
******************************************************************************/ ******************************************************************************/
class MinuitMinimizer: public Minimizer class MinuitMinimizer: public Minimizer
{ {
public: public:
enum class Algorithm enum class Algorithm
{ {
Migrad = 1, Migrad = 1,
Simplex = 2 Simplex = 2,
}; Combined = 3
private:
class MinuitFunction: public ROOT::Minuit2::FCNBase
{
public:
// constructor
explicit MinuitFunction(const DoubleFunction &f);
// destructor
virtual ~MinuitFunction(void) = default;
// minuit members
virtual double operator()(const std::vector<double> &x) const;
virtual double Up(void) const;
private:
const DoubleFunction *f_;
}; };
public: public:
// constructors // constructor
MinuitMinimizer(const Algorithm algorithm = Algorithm::Migrad); MinuitMinimizer(const Algorithm algorithm = defaultAlg_);
explicit MinuitMinimizer(const Index dim, explicit MinuitMinimizer(const Index dim,
const Algorithm algorithm = Algorithm::Migrad); const Algorithm algorithm = defaultAlg_);
// destructor // destructor
virtual ~MinuitMinimizer(void) = default; virtual ~MinuitMinimizer(void) = default;
// access // access
virtual double getPrecision(void) const; Algorithm getAlgorithm(void) const;
Algorithm getAlgorithm(void) const; void setAlgorithm(const Algorithm algorithm);
virtual void setPrecision(const double precision);
void setAlgorithm(const Algorithm algorithm);
// minimization // minimization
virtual const DVec & operator()(const DoubleFunction &f); virtual const DVec & operator()(const DoubleFunction &f);
private: private:
Algorithm algorithm_; Algorithm algorithm_;
static constexpr Algorithm defaultAlg_ = Algorithm::Combined;
}; };
END_LATAN_NAMESPACE END_LATAN_NAMESPACE

2
lib/Solver.hpp
View File

@ -31,7 +31,7 @@ BEGIN_LATAN_NAMESPACE
class Solver class Solver
{ {
public: public:
static const unsigned int defaultMaxIteration = 1000u; static const unsigned int defaultMaxIteration = 10000u;
static constexpr double defaultPrec = 1.0e-7; static constexpr double defaultPrec = 1.0e-7;
public: public:
enum class Verbosity enum class Verbosity