diff --git a/examples/exFit.cpp b/examples/exFit.cpp index 850a142..2c36333 100644 --- a/examples/exFit.cpp +++ b/examples/exFit.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include #include @@ -8,39 +9,50 @@ using namespace std; using namespace Latan; -const Index nPoint = 30; -const double xErr = .2, yErr = .1; -const double exactPar[2] = {0.5,5.0}, dx = 10.0/static_cast(nPoint); +const Index nPoint1 = 5, nPoint2 = 5; +const double xErr = .1, yErr = .1; +const double exactPar[2] = {0.5,5.}; +const double dx1 = 10.0/static_cast(nPoint1); +const double dx2 = 5.0/static_cast(nPoint2); int main(void) { // generate fake data XYStatData data; RandGen rg; - double x1_k, x2_k; + double xBuf[2]; 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"); - for (Index k = 0; k < nPoint; ++k) + for (Index i1 = 0; i1 < nPoint1; ++i1) { - x1_k = rg.gaussian(k*dx, xErr); - x2_k = rg.gaussian(k*dx, xErr); - data.x(k) = x1_k; - 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[0] = i1*dx1; + data.x(i1, 0) = rg.gaussian(xBuf[0], xErr); + for (Index i2 = 0; i2 < nPoint2; ++i2) + { + 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; - data.setXError(0, DVec::Constant(nPoint, xErr)); - data.setYError(0, DVec::Constant(nPoint, yErr)); + data.setXError(0, DVec::Constant(data.getXSize(0), xErr)); + data.assumeXExact(true, 1); + data.setYError(0, DVec::Constant(data.getYSize(), yErr)); cout << data << endl; // fit - DVec init = DVec::Constant(2, 0.5); + DVec init = DVec::Constant(2, 0.1); FitResult p; MinuitMinimizer minimizer; + minimizer.setVerbosity(Minimizer::Verbosity::Normal); p = data.fit(minimizer, init, f); cout << "a= " << p(0) << " b= " << p(1) << endl; cout << "chi^2/ndof= " << p.getChi2PerDof() << endl; diff --git a/examples/exFitSample.cpp b/examples/exFitSample.cpp index 7c6ddd3..6fe5012 100644 --- a/examples/exFitSample.cpp +++ b/examples/exFitSample.cpp @@ -8,41 +8,49 @@ using namespace std; using namespace Latan; -const Index nPoint = 30; +const Index nPoint1 = 10, nPoint2 = 10; const Index nSample = 1000; -const double xErr = .2, yErr = .1; -const double exactPar[2] = {0.5,5.0}, dx = 10.0/static_cast(nPoint); +const double xErr = .1, yErr = .1; +const double exactPar[2] = {0.5,5.}; +const double dx1 = 10.0/static_cast(nPoint1); +const double dx2 = 5.0/static_cast(nPoint2); int main(void) { // generate fake data - DMatSample x(nSample, nPoint, 1), y(nSample, nPoint, 1); XYSampleData data(nSample); - RandGen rg; - double x1_k, x2_k; - DoubleModel f([](const double *t, const double *p) - {return p[1]*exp(-t[0]*p[0]);}, 1, 2); + RandGen rg; + double xBuf[2]; + DoubleModel f([](const double *x, const double *p) + {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"); - 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); - x2_k = rg.gaussian(k*dx, xErr); - data.x(k)[s] = x1_k; - data.y(k)[s] = rg.gaussian(f(&x2_k, exactPar), yErr); + xBuf[0] = i1*dx1; + data.x(i1, 0)[s] = rg.gaussian(xBuf[0], xErr); + for (Index i2 = 0; i2 < nPoint2; ++i2) + { + 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; // fit - DVec init = DVec::Constant(2, 0.5); + DVec init = DVec::Constant(2, 0.1); DMat err; SampleFitResult p; MinuitMinimizer minimizer; - + p = data.fit(minimizer, init, f); err = p.variance().cwiseSqrt(); cout << "a= " << p[central](0) << " +/- " << err(0) << endl; diff --git a/lib/MinuitMinimizer.cpp b/lib/MinuitMinimizer.cpp index 4d8bab6..7f5de0f 100644 --- a/lib/MinuitMinimizer.cpp +++ b/lib/MinuitMinimizer.cpp @@ -1,7 +1,7 @@ /* * 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 * it under the terms of the GNU General Public License as published by @@ -19,48 +19,18 @@ #include #include - -#pragma GCC diagnostic ignored "-Wconversion" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include +#include using namespace std; -using namespace ROOT; -using namespace Minuit2; using namespace Latan; -static constexpr double initErr = 0.5; +static constexpr double initErr = 0.1; 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 &x) const -{ - return (*f_)(x); -} - -double MinuitMinimizer::MinuitFunction::Up(void) const -{ - return 1.; -} - // constructors //////////////////////////////////////////////////////////////// MinuitMinimizer::MinuitMinimizer(const Algorithm algorithm) { @@ -79,126 +49,135 @@ MinuitMinimizer::Algorithm MinuitMinimizer::getAlgorithm(void) const return algorithm_; } -double MinuitMinimizer::getPrecision(void) const -{ - LATAN_ERROR(Implementation, - "Minuit minimizer precision cannot be accessed"); - - return 0.; -} - void MinuitMinimizer::setAlgorithm(const Algorithm algorithm) { algorithm_ = algorithm; } -void MinuitMinimizer::setPrecision(const double precision __dumb) -{ - LATAN_ERROR(Implementation, - "Minuit minimizer precision cannot be accessed"); -} - // minimization //////////////////////////////////////////////////////////////// const DVec & MinuitMinimizer::operator()(const DoubleFunction &f) { - DVec &x = getState(); - Verbosity verbosity = getVerbosity(); + using namespace ROOT; + using namespace Minuit2; + DVec &x = getState(); + int printLevel; + EMinimizerType minuitAlg; + unique_ptr 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()) { resize(f.getNArg()); } - // set parameters - MinuitFunction minuitF(f); - MnUserParameters parameters; + // create and set minimizer + min.reset(new Minuit2Minimizer(minuitAlg)); + 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) { - parameters.Add("x_" + strFrom(i), x(i), initErr*fabs(x(i))); - if (hasLowLimit(i)&&hasHighLimit(i)) + name = "x_" + strFrom(i); + val = x(i); + step = (fabs(x(i)) != 0.) ? initErr*fabs(x(i)) : 1.; + if (hasHighLimit(i) and !hasLowLimit(i)) { - parameters.SetLimits(static_cast(i), - getLowLimit(i), getHighLimit(i)); + min->SetUpperLimitedVariable(i, name, val, step, getHighLimit(i)); } - else if (hasLowLimit(i)) + else if (!hasHighLimit(i) and hasLowLimit(i)) { - parameters.SetLowerLimit(static_cast(i), - getLowLimit(i)); + min->SetLowerLimitedVariable(i, name, val, step, getLowLimit(i)); } - else if (hasHighLimit(i)) + else if (hasHighLimit(i) and hasLowLimit(i)) { - parameters.SetUpperLimit(static_cast(i), - getHighLimit(i)); + min->SetLimitedVariable(i, name, val, step, getLowLimit(i), + getHighLimit(i)); + } + else + { + min->SetVariable(i, name, val, step); } } - // pre-minimization - MnMigrad preMinimizer(minuitF, parameters, 1); - FunctionMinimum preMin = preMinimizer(); + // minimize + int status; + unsigned int n = 0; - if (verbosity >= Verbosity::Debug) + if (getVerbosity() >= Verbosity::Normal) { - cout << "-- MINUIT pre-minimization -----------------------------"; - cout << preMin; - cout << "--------------------------------------------------------"; - cout << endl; + cout << "========== Minuit pre-minimization " << endl; } - for (unsigned int i = 0; i < x.size(); ++i) + min->SetStrategy(0); + min->Minimize(); + do { - parameters.SetValue(i, preMin.UserParameters().Value(i)); - parameters.SetError(i, 2.*preMin.UserParameters().Error(i)); - } - - // minimization and output - unique_ptr 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> scanRes; - MnPlot plot; - MnScan scanner(minuitF, preMin.UserParameters(), 2); - - cout << "-- MINUIT scan -----------------------------------------"; - cout << endl; - for (unsigned int i = 0; i < x.size(); i++) + n++; + if (getVerbosity() >= Verbosity::Normal) { - cout << "Parameter x_" << i << endl; - scanRes = scanner.Scan(i); - plot(scanRes); + cout << "========== Minuit minimization, try #" << n << endl; } - cout << "--------------------------------------------------------"; - cout << endl; + min->SetStrategy(2); + 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; diff --git a/lib/MinuitMinimizer.hpp b/lib/MinuitMinimizer.hpp index f79873c..91b1bc8 100644 --- a/lib/MinuitMinimizer.hpp +++ b/lib/MinuitMinimizer.hpp @@ -1,7 +1,7 @@ /* * 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 * it under the terms of the GNU General Public License as published by @@ -23,52 +23,36 @@ #include #include #include -#include BEGIN_LATAN_NAMESPACE /****************************************************************************** - * Minuit minimizer * + * MinuitMinimizer * ******************************************************************************/ - class MinuitMinimizer: public Minimizer { public: enum class Algorithm { - Migrad = 1, - Simplex = 2 - }; -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 &x) const; - virtual double Up(void) const; - private: - const DoubleFunction *f_; + Migrad = 1, + Simplex = 2, + Combined = 3 }; public: - // constructors - MinuitMinimizer(const Algorithm algorithm = Algorithm::Migrad); + // constructor + MinuitMinimizer(const Algorithm algorithm = defaultAlg_); explicit MinuitMinimizer(const Index dim, - const Algorithm algorithm = Algorithm::Migrad); + const Algorithm algorithm = defaultAlg_); // destructor virtual ~MinuitMinimizer(void) = default; // access - virtual double getPrecision(void) const; - Algorithm getAlgorithm(void) const; - virtual void setPrecision(const double precision); - void setAlgorithm(const Algorithm algorithm); + Algorithm getAlgorithm(void) const; + void setAlgorithm(const Algorithm algorithm); // minimization virtual const DVec & operator()(const DoubleFunction &f); private: - Algorithm algorithm_; + Algorithm algorithm_; + static constexpr Algorithm defaultAlg_ = Algorithm::Combined; }; END_LATAN_NAMESPACE diff --git a/lib/Solver.hpp b/lib/Solver.hpp index 23601c5..1d82e27 100644 --- a/lib/Solver.hpp +++ b/lib/Solver.hpp @@ -31,7 +31,7 @@ BEGIN_LATAN_NAMESPACE class Solver { public: - static const unsigned int defaultMaxIteration = 1000u; + static const unsigned int defaultMaxIteration = 10000u; static constexpr double defaultPrec = 1.0e-7; public: enum class Verbosity