mirror of
https://github.com/aportelli/LatAnalyze.git
synced 2025-04-11 03:20:46 +01:00
interface to GSL minimisers
This commit is contained in:
parent
8683721c5a
commit
9bf538dfca
@ -2,8 +2,7 @@
|
|||||||
#include <cmath>
|
#include <cmath>
|
||||||
#include <LatAnalyze/CompiledModel.hpp>
|
#include <LatAnalyze/CompiledModel.hpp>
|
||||||
#include <LatAnalyze/Io.hpp>
|
#include <LatAnalyze/Io.hpp>
|
||||||
#include <LatAnalyze/MinuitMinimizer.hpp>
|
#include <LatAnalyze/GslMinimizer.hpp>
|
||||||
#include <LatAnalyze/NloptMinimizer.hpp>
|
|
||||||
#include <LatAnalyze/Plot.hpp>
|
#include <LatAnalyze/Plot.hpp>
|
||||||
#include <LatAnalyze/XYStatData.hpp>
|
#include <LatAnalyze/XYStatData.hpp>
|
||||||
|
|
||||||
@ -48,26 +47,12 @@ int main(void)
|
|||||||
data.setYError(0, DVec::Constant(data.getYSize(), yErr));
|
data.setYError(0, DVec::Constant(data.getYSize(), yErr));
|
||||||
|
|
||||||
// set minimizers
|
// set minimizers
|
||||||
DVec init = DVec::Constant(2, 0.1);
|
DVec init = DVec::Constant(2, 0.1);
|
||||||
FitResult p;
|
FitResult p;
|
||||||
NloptMinimizer globalMin(NloptMinimizer::Algorithm::GN_CRS2_LM);
|
GslMinimizer min(GslMinimizer::Algorithm::bfgs2);
|
||||||
MinuitMinimizer localMin;
|
|
||||||
vector<Minimizer *> min{&globalMin, &localMin};
|
|
||||||
|
|
||||||
globalMin.setVerbosity(Minimizer::Verbosity::Normal);
|
|
||||||
globalMin.setPrecision(0.1);
|
|
||||||
globalMin.setMaxIteration(10000);
|
|
||||||
globalMin.useLowLimit(0);
|
|
||||||
globalMin.setLowLimit(0, 0.);
|
|
||||||
globalMin.useHighLimit(0);
|
|
||||||
globalMin.setHighLimit(0, 20.);
|
|
||||||
globalMin.useLowLimit(1);
|
|
||||||
globalMin.setLowLimit(1, 0.);
|
|
||||||
globalMin.useHighLimit(1);
|
|
||||||
globalMin.setHighLimit(1, 20.);
|
|
||||||
localMin.setVerbosity(Minimizer::Verbosity::Normal);
|
|
||||||
|
|
||||||
// fit
|
// fit
|
||||||
|
min.setVerbosity(Minimizer::Verbosity::Debug);
|
||||||
cout << "-- fit..." << endl;
|
cout << "-- fit..." << endl;
|
||||||
f.parName().setName(0, "m");
|
f.parName().setName(0, "m");
|
||||||
f.parName().setName(1, "A");
|
f.parName().setName(1, "A");
|
||||||
|
@ -1,6 +1,5 @@
|
|||||||
#include <LatAnalyze/CompiledModel.hpp>
|
#include <LatAnalyze/CompiledModel.hpp>
|
||||||
#include <LatAnalyze/MinuitMinimizer.hpp>
|
#include <LatAnalyze/GslMinimizer.hpp>
|
||||||
#include <LatAnalyze/NloptMinimizer.hpp>
|
|
||||||
#include <LatAnalyze/Plot.hpp>
|
#include <LatAnalyze/Plot.hpp>
|
||||||
#include <LatAnalyze/XYSampleData.hpp>
|
#include <LatAnalyze/XYSampleData.hpp>
|
||||||
|
|
||||||
@ -47,22 +46,9 @@ int main(void)
|
|||||||
data.assumeXExact(true, 1);
|
data.assumeXExact(true, 1);
|
||||||
|
|
||||||
// set minimizers
|
// set minimizers
|
||||||
DVec init = DVec::Constant(2, 0.1);
|
DVec init = DVec::Constant(2, 0.1);
|
||||||
SampleFitResult p;
|
SampleFitResult p;
|
||||||
NloptMinimizer globalMin(NloptMinimizer::Algorithm::GN_CRS2_LM);
|
GslMinimizer min(GslMinimizer::Algorithm::bfgs2);
|
||||||
MinuitMinimizer localMin;
|
|
||||||
vector<Minimizer *> min{&globalMin, &localMin};
|
|
||||||
|
|
||||||
globalMin.setPrecision(0.1);
|
|
||||||
globalMin.setMaxIteration(10000);
|
|
||||||
globalMin.useLowLimit(0);
|
|
||||||
globalMin.setLowLimit(0, 0.);
|
|
||||||
globalMin.useHighLimit(0);
|
|
||||||
globalMin.setHighLimit(0, 20.);
|
|
||||||
globalMin.useLowLimit(1);
|
|
||||||
globalMin.setLowLimit(1, 0.);
|
|
||||||
globalMin.useHighLimit(1);
|
|
||||||
globalMin.setHighLimit(1, 20.);
|
|
||||||
|
|
||||||
// fit
|
// fit
|
||||||
cout << "-- fit..." << endl;
|
cout << "-- fit..." << endl;
|
||||||
|
@ -41,6 +41,7 @@ CONST_EXC(Range, Logic("range error: " + msg, loc))
|
|||||||
CONST_EXC(Size, Logic("size error: " + msg, loc))
|
CONST_EXC(Size, Logic("size error: " + msg, loc))
|
||||||
// runtime errors
|
// runtime errors
|
||||||
CONST_EXC(Runtime, runtime_error(Env::msgPrefix + msg + ERR_SUFF))
|
CONST_EXC(Runtime, runtime_error(Env::msgPrefix + msg + ERR_SUFF))
|
||||||
|
CONST_EXC(Argument, Runtime("argument error: " + msg, loc))
|
||||||
CONST_EXC(Compilation, Runtime("compilation error: " + msg, loc))
|
CONST_EXC(Compilation, Runtime("compilation error: " + msg, loc))
|
||||||
CONST_EXC(Io, Runtime("IO error: " + msg, loc))
|
CONST_EXC(Io, Runtime("IO error: " + msg, loc))
|
||||||
CONST_EXC(Memory, Runtime("memory error: " + msg, loc))
|
CONST_EXC(Memory, Runtime("memory error: " + msg, loc))
|
||||||
|
@ -51,6 +51,7 @@ namespace Exceptions
|
|||||||
DECL_EXC(Size, Logic);
|
DECL_EXC(Size, Logic);
|
||||||
// runtime errors
|
// runtime errors
|
||||||
DECL_EXC(Runtime, std::runtime_error);
|
DECL_EXC(Runtime, std::runtime_error);
|
||||||
|
DECL_EXC(Argument, Runtime);
|
||||||
DECL_EXC(Compilation, Runtime);
|
DECL_EXC(Compilation, Runtime);
|
||||||
DECL_EXC(Io, Runtime);
|
DECL_EXC(Io, Runtime);
|
||||||
DECL_EXC(Memory, Runtime);
|
DECL_EXC(Memory, Runtime);
|
||||||
|
359
lib/GslMinimizer.cpp
Normal file
359
lib/GslMinimizer.cpp
Normal file
@ -0,0 +1,359 @@
|
|||||||
|
/*
|
||||||
|
* GslMinimizer.cpp, part of LatAnalyze
|
||||||
|
*
|
||||||
|
* Copyright (C) 2013 - 2016 Antonin Portelli
|
||||||
|
*
|
||||||
|
* LatAnalyze 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.
|
||||||
|
*
|
||||||
|
* LatAnalyze 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 LatAnalyze. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
*/
|
||||||
|
|
||||||
|
#include <LatAnalyze/GslMinimizer.hpp>
|
||||||
|
#include <LatAnalyze/includes.hpp>
|
||||||
|
#include <LatAnalyze/Math.hpp>
|
||||||
|
#include <gsl/gsl_multimin.h>
|
||||||
|
#include <gsl/gsl_blas.h>
|
||||||
|
|
||||||
|
using namespace std;
|
||||||
|
using namespace Latan;
|
||||||
|
|
||||||
|
/******************************************************************************
|
||||||
|
* GslMinimizer implementation *
|
||||||
|
******************************************************************************/
|
||||||
|
// constructor /////////////////////////////////////////////////////////////////
|
||||||
|
GslMinimizer::GslMinimizer(const Algorithm algorithm)
|
||||||
|
{
|
||||||
|
setAlgorithm(algorithm);
|
||||||
|
der_.setOrder(1, 1);
|
||||||
|
}
|
||||||
|
|
||||||
|
// access //////////////////////////////////////////////////////////////////////
|
||||||
|
GslMinimizer::Algorithm GslMinimizer::getAlgorithm(void) const
|
||||||
|
{
|
||||||
|
return algorithm_;
|
||||||
|
}
|
||||||
|
|
||||||
|
void GslMinimizer::setAlgorithm(const Algorithm algorithm)
|
||||||
|
{
|
||||||
|
algorithm_ = algorithm;
|
||||||
|
}
|
||||||
|
|
||||||
|
bool GslMinimizer::supportLimits(void) const
|
||||||
|
{
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
|
||||||
|
// test ////////////////////////////////////////////////////////////////////////
|
||||||
|
bool GslMinimizer::isDerAlgorithm(const Algorithm algorithm)
|
||||||
|
{
|
||||||
|
return (algorithm <= Algorithm::lastDerAlg);
|
||||||
|
}
|
||||||
|
|
||||||
|
// minimization ////////////////////////////////////////////////////////////////
|
||||||
|
const DVec & GslMinimizer::operator()(const DoubleFunction &f)
|
||||||
|
{
|
||||||
|
DVec &x = getState();
|
||||||
|
|
||||||
|
// resize minimizer state to match function number of arguments
|
||||||
|
if (f.getNArg() != x.size())
|
||||||
|
{
|
||||||
|
resize(f.getNArg());
|
||||||
|
}
|
||||||
|
|
||||||
|
// set function data
|
||||||
|
GslFuncData data;
|
||||||
|
|
||||||
|
der_.setFunction(f);
|
||||||
|
data.f = &f;
|
||||||
|
data.d = &der_;
|
||||||
|
|
||||||
|
// set initial position
|
||||||
|
gsl_vector *gslX = gsl_vector_alloc(getDim());
|
||||||
|
|
||||||
|
for (Index i = 0; i < getDim(); ++i)
|
||||||
|
{
|
||||||
|
gsl_vector_set(gslX, i, x(i));
|
||||||
|
}
|
||||||
|
|
||||||
|
// minimization
|
||||||
|
int status;
|
||||||
|
|
||||||
|
if (isDerAlgorithm(getAlgorithm()))
|
||||||
|
{
|
||||||
|
// set function
|
||||||
|
gsl_multimin_function_fdf gslFunc;
|
||||||
|
|
||||||
|
gslFunc.n = getDim();
|
||||||
|
gslFunc.f = &fWrapper;
|
||||||
|
gslFunc.df = &dfWrapper;
|
||||||
|
gslFunc.fdf = &fdfWrapper;
|
||||||
|
gslFunc.params = &data;
|
||||||
|
|
||||||
|
// create and set minimizer
|
||||||
|
const gsl_multimin_fdfminimizer_type *gslAlg;
|
||||||
|
gsl_multimin_fdfminimizer *gslMin;
|
||||||
|
|
||||||
|
switch (getAlgorithm())
|
||||||
|
{
|
||||||
|
case Algorithm::cgFR:
|
||||||
|
gslAlg = gsl_multimin_fdfminimizer_conjugate_fr;
|
||||||
|
break;
|
||||||
|
case Algorithm::cgPR:
|
||||||
|
gslAlg = gsl_multimin_fdfminimizer_conjugate_pr;
|
||||||
|
break;
|
||||||
|
case Algorithm::bfgs:
|
||||||
|
gslAlg = gsl_multimin_fdfminimizer_vector_bfgs;
|
||||||
|
break;
|
||||||
|
case Algorithm::bfgs2:
|
||||||
|
gslAlg = gsl_multimin_fdfminimizer_vector_bfgs2;
|
||||||
|
break;
|
||||||
|
case Algorithm::steepDesc:
|
||||||
|
gslAlg = gsl_multimin_fdfminimizer_vector_bfgs2;
|
||||||
|
break;
|
||||||
|
default:
|
||||||
|
LATAN_ERROR(Argument, "unknow GSL minization algorithm "
|
||||||
|
+ strFrom(static_cast<int>(getAlgorithm())));
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
gslMin = gsl_multimin_fdfminimizer_alloc(gslAlg, getDim());
|
||||||
|
|
||||||
|
// minimize
|
||||||
|
unsigned int pass = 0, it;
|
||||||
|
double dxRel;
|
||||||
|
|
||||||
|
do
|
||||||
|
{
|
||||||
|
pass++;
|
||||||
|
gsl_multimin_fdfminimizer_set(gslMin, &gslFunc, gslX, 0.01, 0.001);
|
||||||
|
if (getVerbosity() >= Verbosity::Normal)
|
||||||
|
{
|
||||||
|
cout << "========== GSL minimization, pass #" << pass;
|
||||||
|
cout << " ==========" << endl;
|
||||||
|
cout << "Algorithm: " << getAlgorithmName(getAlgorithm());
|
||||||
|
cout << endl;
|
||||||
|
cout << "Max eval.= " << getMaxIteration();
|
||||||
|
cout << " -- Precision= " << getPrecision() << endl;
|
||||||
|
printf("Starting f(x)= %.10e\n", f(x));
|
||||||
|
}
|
||||||
|
it = 0;
|
||||||
|
do
|
||||||
|
{
|
||||||
|
it++;
|
||||||
|
gsl_multimin_fdfminimizer_iterate(gslMin);
|
||||||
|
dxRel = gsl_blas_dnrm2(gslMin->dx)/gsl_blas_dnrm2(gslMin->x);
|
||||||
|
status = (dxRel < getPrecision()) ? GSL_SUCCESS : GSL_CONTINUE;
|
||||||
|
if (getVerbosity() >= Verbosity::Debug)
|
||||||
|
{
|
||||||
|
printf("iteration %4d: f= %.10e dxrel= %.10e eval= %d\n",
|
||||||
|
it, gslMin->f, dxRel, data.evalCount);
|
||||||
|
}
|
||||||
|
|
||||||
|
} while (status == GSL_CONTINUE and
|
||||||
|
(data.evalCount < getMaxIteration()));
|
||||||
|
if (getVerbosity() >= Verbosity::Normal)
|
||||||
|
{
|
||||||
|
printf("Found minimum %.10e at:\n", gslMin->f);
|
||||||
|
for (Index i = 0; i < x.size(); ++i)
|
||||||
|
{
|
||||||
|
printf("%8s= %.10e\n", f.varName().getName(i).c_str(),
|
||||||
|
gsl_vector_get(gslMin->x, i));
|
||||||
|
}
|
||||||
|
cout << "after " << data.evalCount << " evaluations" << endl;
|
||||||
|
cout << "Minimization ended with code " << status;
|
||||||
|
cout << endl;
|
||||||
|
}
|
||||||
|
data.evalCount = 0;
|
||||||
|
for (Index i = 0; i < getDim(); ++i)
|
||||||
|
{
|
||||||
|
gsl_vector_set(gslX, i, gsl_vector_get(gslMin->x, i));
|
||||||
|
}
|
||||||
|
} while (status != GSL_SUCCESS and (pass < getMaxPass()));
|
||||||
|
|
||||||
|
// deallocate GSL minimizer
|
||||||
|
gsl_multimin_fdfminimizer_free(gslMin);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
// set function
|
||||||
|
gsl_multimin_function gslFunc;
|
||||||
|
|
||||||
|
gslFunc.n = getDim();
|
||||||
|
gslFunc.f = &fWrapper;
|
||||||
|
gslFunc.params = &data;
|
||||||
|
|
||||||
|
// create and set minimizer
|
||||||
|
const gsl_multimin_fminimizer_type *gslAlg;
|
||||||
|
gsl_multimin_fminimizer *gslMin;
|
||||||
|
|
||||||
|
switch (getAlgorithm())
|
||||||
|
{
|
||||||
|
case Algorithm::simplex:
|
||||||
|
gslAlg = gsl_multimin_fminimizer_nmsimplex;
|
||||||
|
break;
|
||||||
|
case Algorithm::simplex2:
|
||||||
|
gslAlg = gsl_multimin_fminimizer_nmsimplex2;
|
||||||
|
break;
|
||||||
|
case Algorithm::simplex2R:
|
||||||
|
gslAlg = gsl_multimin_fminimizer_nmsimplex2rand;
|
||||||
|
break;
|
||||||
|
default:
|
||||||
|
LATAN_ERROR(Argument, "unknow GSL minization algorithm "
|
||||||
|
+ strFrom(static_cast<int>(getAlgorithm())));
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
gslMin = gsl_multimin_fminimizer_alloc(gslAlg, getDim());
|
||||||
|
|
||||||
|
// minimize
|
||||||
|
unsigned int pass = 0, it;
|
||||||
|
gsl_vector *step = gsl_vector_alloc(getDim());
|
||||||
|
double relSize;
|
||||||
|
|
||||||
|
gsl_vector_set_all(step, 0.01);
|
||||||
|
do
|
||||||
|
{
|
||||||
|
pass++;
|
||||||
|
gsl_multimin_fminimizer_set(gslMin, &gslFunc, gslX, step);
|
||||||
|
if (getVerbosity() >= Verbosity::Normal)
|
||||||
|
{
|
||||||
|
cout << "========== GSL minimization, pass #" << pass;
|
||||||
|
cout << " ==========" << endl;
|
||||||
|
cout << "Algorithm: " << getAlgorithmName(getAlgorithm());
|
||||||
|
cout << endl;
|
||||||
|
cout << "Max eval.= " << getMaxIteration();
|
||||||
|
cout << " -- Precision= " << getPrecision() << endl;
|
||||||
|
printf("Starting f(x)= %.10e\n", f(x));
|
||||||
|
}
|
||||||
|
it = 0;
|
||||||
|
do
|
||||||
|
{
|
||||||
|
it++;
|
||||||
|
gsl_multimin_fminimizer_iterate(gslMin);
|
||||||
|
relSize = Math::pow<2>(gslMin->size)/gsl_blas_dnrm2(gslMin->x);
|
||||||
|
status = (relSize < getPrecision()) ? GSL_SUCCESS
|
||||||
|
: GSL_CONTINUE;
|
||||||
|
if (getVerbosity() >= Verbosity::Debug)
|
||||||
|
{
|
||||||
|
printf("iteration %4d: f= %.10e relSize= %.10e eval= %d\n",
|
||||||
|
it, gslMin->fval, relSize, data.evalCount);
|
||||||
|
}
|
||||||
|
|
||||||
|
} while (status == GSL_CONTINUE and
|
||||||
|
(data.evalCount < getMaxIteration()));
|
||||||
|
if (getVerbosity() >= Verbosity::Normal)
|
||||||
|
{
|
||||||
|
printf("Found minimum %.10e at:\n", gslMin->fval);
|
||||||
|
for (Index i = 0; i < x.size(); ++i)
|
||||||
|
{
|
||||||
|
printf("%8s= %.10e\n", f.varName().getName(i).c_str(),
|
||||||
|
gsl_vector_get(gslMin->x, i));
|
||||||
|
}
|
||||||
|
cout << "after " << data.evalCount << " evaluations" << endl;
|
||||||
|
cout << "Minimization ended with code " << status;
|
||||||
|
cout << endl;
|
||||||
|
}
|
||||||
|
data.evalCount = 0;
|
||||||
|
for (Index i = 0; i < getDim(); ++i)
|
||||||
|
{
|
||||||
|
gsl_vector_set(gslX, i, gsl_vector_get(gslMin->x, i));
|
||||||
|
}
|
||||||
|
} while (status != GSL_SUCCESS and (pass < getMaxPass()));
|
||||||
|
|
||||||
|
// deallocate GSL minimizer
|
||||||
|
gsl_multimin_fminimizer_free(gslMin);
|
||||||
|
gsl_vector_free(step);
|
||||||
|
}
|
||||||
|
if (status != GSL_SUCCESS)
|
||||||
|
{
|
||||||
|
LATAN_WARNING("invalid minimum: maximum number of call reached");
|
||||||
|
}
|
||||||
|
|
||||||
|
// save final result
|
||||||
|
for (Index i = 0; i < getDim(); ++i)
|
||||||
|
{
|
||||||
|
x(i) = gsl_vector_get(gslX, i);
|
||||||
|
}
|
||||||
|
|
||||||
|
// deallocate GSL state and return
|
||||||
|
gsl_vector_free(gslX);
|
||||||
|
|
||||||
|
return x;
|
||||||
|
}
|
||||||
|
|
||||||
|
// function wrappers ///////////////////////////////////////////////////////////
|
||||||
|
double GslMinimizer::fWrapper(const gsl_vector *x, void *vdata)
|
||||||
|
{
|
||||||
|
GslFuncData &data = *static_cast<GslFuncData *>(vdata);
|
||||||
|
|
||||||
|
data.evalCount++;
|
||||||
|
|
||||||
|
return (*data.f)(x->data);
|
||||||
|
}
|
||||||
|
|
||||||
|
void GslMinimizer::dfWrapper(const gsl_vector *x, void *vdata, gsl_vector * df)
|
||||||
|
{
|
||||||
|
GslFuncData &data = *static_cast<GslFuncData *>(vdata);
|
||||||
|
const unsigned int n = data.f->getNArg();
|
||||||
|
|
||||||
|
for (unsigned int i = 0; i < n; ++i)
|
||||||
|
{
|
||||||
|
data.d->setDir(i);
|
||||||
|
gsl_vector_set(df, i, (*(data.d))(x->data));
|
||||||
|
}
|
||||||
|
data.evalCount += data.d->getNPoint()*n;
|
||||||
|
}
|
||||||
|
|
||||||
|
void GslMinimizer::fdfWrapper(const gsl_vector *x, void *vdata, double *f,
|
||||||
|
gsl_vector * df)
|
||||||
|
{
|
||||||
|
GslFuncData &data = *static_cast<GslFuncData *>(vdata);
|
||||||
|
const unsigned int n = data.f->getNArg();
|
||||||
|
|
||||||
|
for (unsigned int i = 0; i < n; ++i)
|
||||||
|
{
|
||||||
|
data.d->setDir(i);
|
||||||
|
gsl_vector_set(df, i, (*(data.d))(x->data));
|
||||||
|
}
|
||||||
|
*f = (*data.f)(x->data);
|
||||||
|
data.evalCount += data.d->getNPoint()*n + 1;
|
||||||
|
}
|
||||||
|
|
||||||
|
// algorithm names /////////////////////////////////////////////////////////////
|
||||||
|
string GslMinimizer::getAlgorithmName(const Algorithm algorithm)
|
||||||
|
{
|
||||||
|
switch (algorithm)
|
||||||
|
{
|
||||||
|
case Algorithm::cgFR:
|
||||||
|
return "Fletcher-Reeves conjugate gradient";
|
||||||
|
break;
|
||||||
|
case Algorithm::cgPR:
|
||||||
|
return "Polak-Ribiere conjugate gradient";
|
||||||
|
break;
|
||||||
|
case Algorithm::bfgs:
|
||||||
|
return "Broyden-Fletcher-Goldfarb-Shanno";
|
||||||
|
break;
|
||||||
|
case Algorithm::bfgs2:
|
||||||
|
return "improved Broyden-Fletcher-Goldfarb-Shanno";
|
||||||
|
break;
|
||||||
|
case Algorithm::steepDesc:
|
||||||
|
return "steepest descent";
|
||||||
|
break;
|
||||||
|
case Algorithm::simplex:
|
||||||
|
return "Nelder-Mead simplex";
|
||||||
|
break;
|
||||||
|
case Algorithm::simplex2:
|
||||||
|
return "improved Nelder-Mead simplex";
|
||||||
|
break;
|
||||||
|
case Algorithm::simplex2R:
|
||||||
|
return "improved Nelder-Mead simplex with random start";
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
86
lib/GslMinimizer.hpp
Normal file
86
lib/GslMinimizer.hpp
Normal file
@ -0,0 +1,86 @@
|
|||||||
|
/*
|
||||||
|
* GslMinimizer.hpp, part of LatAnalyze
|
||||||
|
*
|
||||||
|
* Copyright (C) 2013 - 2016 Antonin Portelli
|
||||||
|
*
|
||||||
|
* LatAnalyze 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.
|
||||||
|
*
|
||||||
|
* LatAnalyze 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 LatAnalyze. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
*/
|
||||||
|
|
||||||
|
#ifndef Latan_GslMinimizer_hpp_
|
||||||
|
#define Latan_GslMinimizer_hpp_
|
||||||
|
|
||||||
|
#include <LatAnalyze/Global.hpp>
|
||||||
|
#include <LatAnalyze/Derivative.hpp>
|
||||||
|
#include <LatAnalyze/Function.hpp>
|
||||||
|
#include <LatAnalyze/Minimizer.hpp>
|
||||||
|
#include <gsl/gsl_vector.h>
|
||||||
|
|
||||||
|
BEGIN_LATAN_NAMESPACE
|
||||||
|
|
||||||
|
/******************************************************************************
|
||||||
|
* interface to the GSL minimizers *
|
||||||
|
******************************************************************************/
|
||||||
|
class GslMinimizer: public Minimizer
|
||||||
|
{
|
||||||
|
public:
|
||||||
|
enum class Algorithm
|
||||||
|
{
|
||||||
|
cgFR = 1,
|
||||||
|
cgPR = 2,
|
||||||
|
bfgs = 3,
|
||||||
|
bfgs2 = 4,
|
||||||
|
steepDesc = 5,
|
||||||
|
lastDerAlg = 5,
|
||||||
|
simplex = 6,
|
||||||
|
simplex2 = 7,
|
||||||
|
simplex2R = 8
|
||||||
|
};
|
||||||
|
private:
|
||||||
|
struct GslFuncData
|
||||||
|
{
|
||||||
|
const DoubleFunction *f{nullptr};
|
||||||
|
Derivative *d{nullptr};
|
||||||
|
unsigned int evalCount{0};
|
||||||
|
};
|
||||||
|
public:
|
||||||
|
// constructor
|
||||||
|
explicit GslMinimizer(const Algorithm algorithm = defaultAlg_);
|
||||||
|
// destructor
|
||||||
|
virtual ~GslMinimizer(void) = default;
|
||||||
|
// access
|
||||||
|
Algorithm getAlgorithm(void) const;
|
||||||
|
void setAlgorithm(const Algorithm algorithm);
|
||||||
|
virtual bool supportLimits(void) const;
|
||||||
|
// minimization
|
||||||
|
virtual const DVec & operator()(const DoubleFunction &f);
|
||||||
|
private:
|
||||||
|
// test
|
||||||
|
static bool isDerAlgorithm(const Algorithm algorithm);
|
||||||
|
// function wrappers
|
||||||
|
static double fWrapper(const gsl_vector *x, void * params);
|
||||||
|
static void dfWrapper(const gsl_vector *x, void * params,
|
||||||
|
gsl_vector * df);
|
||||||
|
static void fdfWrapper(const gsl_vector *x, void *params, double *f,
|
||||||
|
gsl_vector * df);
|
||||||
|
// algorithm names
|
||||||
|
std::string getAlgorithmName(const Algorithm algorithm);
|
||||||
|
private:
|
||||||
|
Algorithm algorithm_;
|
||||||
|
static constexpr Algorithm defaultAlg_ = Algorithm::simplex2;
|
||||||
|
CentralDerivative der_;
|
||||||
|
};
|
||||||
|
|
||||||
|
END_LATAN_NAMESPACE
|
||||||
|
|
||||||
|
#endif // Latan_GslMinimizer_hpp_
|
@ -32,6 +32,7 @@ libLatAnalyze_la_SOURCES = \
|
|||||||
Function.cpp \
|
Function.cpp \
|
||||||
Global.cpp \
|
Global.cpp \
|
||||||
GslHybridRootFinder.cpp\
|
GslHybridRootFinder.cpp\
|
||||||
|
GslMinimizer.cpp \
|
||||||
GslQagsIntegrator.cpp \
|
GslQagsIntegrator.cpp \
|
||||||
Hdf5File.cpp \
|
Hdf5File.cpp \
|
||||||
Histogram.cpp \
|
Histogram.cpp \
|
||||||
@ -64,6 +65,7 @@ libLatAnalyze_la_HEADERS = \
|
|||||||
FitInterface.hpp \
|
FitInterface.hpp \
|
||||||
Global.hpp \
|
Global.hpp \
|
||||||
GslHybridRootFinder.hpp\
|
GslHybridRootFinder.hpp\
|
||||||
|
GslMinimizer.hpp \
|
||||||
GslQagsIntegrator.hpp \
|
GslQagsIntegrator.hpp \
|
||||||
Hdf5File.hpp \
|
Hdf5File.hpp \
|
||||||
Histogram.hpp \
|
Histogram.hpp \
|
||||||
|
@ -43,8 +43,15 @@ void Minimizer::resize(const Index dim)
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
#define checkSupport \
|
||||||
|
if (!supportLimits())\
|
||||||
|
{\
|
||||||
|
LATAN_ERROR(Implementation, "minimizer does not support limits");\
|
||||||
|
}
|
||||||
|
|
||||||
double Minimizer::getHighLimit(const Index i) const
|
double Minimizer::getHighLimit(const Index i) const
|
||||||
{
|
{
|
||||||
|
checkSupport;
|
||||||
if (i >= getDim())
|
if (i >= getDim())
|
||||||
{
|
{
|
||||||
LATAN_ERROR(Size, "invalid variable index");
|
LATAN_ERROR(Size, "invalid variable index");
|
||||||
@ -55,11 +62,14 @@ double Minimizer::getHighLimit(const Index i) const
|
|||||||
|
|
||||||
const DVec & Minimizer::getHighLimit(const PlaceHolder ph __dumb) const
|
const DVec & Minimizer::getHighLimit(const PlaceHolder ph __dumb) const
|
||||||
{
|
{
|
||||||
|
checkSupport;
|
||||||
|
|
||||||
return highLimit_;
|
return highLimit_;
|
||||||
}
|
}
|
||||||
|
|
||||||
double Minimizer::getLowLimit(const Index i) const
|
double Minimizer::getLowLimit(const Index i) const
|
||||||
{
|
{
|
||||||
|
checkSupport;
|
||||||
if (i >= getDim())
|
if (i >= getDim())
|
||||||
{
|
{
|
||||||
LATAN_ERROR(Size, "invalid variable index");
|
LATAN_ERROR(Size, "invalid variable index");
|
||||||
@ -70,11 +80,14 @@ double Minimizer::getLowLimit(const Index i) const
|
|||||||
|
|
||||||
const DVec & Minimizer::getLowLimit(const PlaceHolder ph __dumb) const
|
const DVec & Minimizer::getLowLimit(const PlaceHolder ph __dumb) const
|
||||||
{
|
{
|
||||||
|
checkSupport;
|
||||||
|
|
||||||
return lowLimit_;
|
return lowLimit_;
|
||||||
}
|
}
|
||||||
|
|
||||||
bool Minimizer::hasHighLimit(const Index i) const
|
bool Minimizer::hasHighLimit(const Index i) const
|
||||||
{
|
{
|
||||||
|
checkSupport;
|
||||||
if (i >= getDim())
|
if (i >= getDim())
|
||||||
{
|
{
|
||||||
LATAN_ERROR(Size, "invalid variable index");
|
LATAN_ERROR(Size, "invalid variable index");
|
||||||
@ -85,6 +98,7 @@ bool Minimizer::hasHighLimit(const Index i) const
|
|||||||
|
|
||||||
bool Minimizer::hasLowLimit(const Index i) const
|
bool Minimizer::hasLowLimit(const Index i) const
|
||||||
{
|
{
|
||||||
|
checkSupport;
|
||||||
if (i >= getDim())
|
if (i >= getDim())
|
||||||
{
|
{
|
||||||
LATAN_ERROR(Size, "invalid variable index");
|
LATAN_ERROR(Size, "invalid variable index");
|
||||||
@ -95,6 +109,7 @@ bool Minimizer::hasLowLimit(const Index i) const
|
|||||||
|
|
||||||
void Minimizer::setHighLimit(const Index i, const double l)
|
void Minimizer::setHighLimit(const Index i, const double l)
|
||||||
{
|
{
|
||||||
|
checkSupport;
|
||||||
if (i >= getDim())
|
if (i >= getDim())
|
||||||
{
|
{
|
||||||
resize(i + 1);
|
resize(i + 1);
|
||||||
@ -105,6 +120,7 @@ void Minimizer::setHighLimit(const Index i, const double l)
|
|||||||
|
|
||||||
void Minimizer::setHighLimit(const PlaceHolder ph __dumb, const DVec &l)
|
void Minimizer::setHighLimit(const PlaceHolder ph __dumb, const DVec &l)
|
||||||
{
|
{
|
||||||
|
checkSupport;
|
||||||
if (l.size() != getDim())
|
if (l.size() != getDim())
|
||||||
{
|
{
|
||||||
resize(l.size());
|
resize(l.size());
|
||||||
@ -115,6 +131,7 @@ void Minimizer::setHighLimit(const PlaceHolder ph __dumb, const DVec &l)
|
|||||||
|
|
||||||
void Minimizer::setLowLimit(const Index i, const double l)
|
void Minimizer::setLowLimit(const Index i, const double l)
|
||||||
{
|
{
|
||||||
|
checkSupport;
|
||||||
if (i >= getDim())
|
if (i >= getDim())
|
||||||
{
|
{
|
||||||
resize(i + 1);
|
resize(i + 1);
|
||||||
@ -125,6 +142,7 @@ void Minimizer::setLowLimit(const Index i, const double l)
|
|||||||
|
|
||||||
void Minimizer::setLowLimit(const PlaceHolder ph __dumb, const DVec &l)
|
void Minimizer::setLowLimit(const PlaceHolder ph __dumb, const DVec &l)
|
||||||
{
|
{
|
||||||
|
checkSupport;
|
||||||
if (l.size() != getDim())
|
if (l.size() != getDim())
|
||||||
{
|
{
|
||||||
resize(l.size());
|
resize(l.size());
|
||||||
@ -135,6 +153,7 @@ void Minimizer::setLowLimit(const PlaceHolder ph __dumb, const DVec &l)
|
|||||||
|
|
||||||
void Minimizer::useHighLimit(const Index i, const bool use)
|
void Minimizer::useHighLimit(const Index i, const bool use)
|
||||||
{
|
{
|
||||||
|
checkSupport;
|
||||||
if (i >= getDim())
|
if (i >= getDim())
|
||||||
{
|
{
|
||||||
resize(i + 1);
|
resize(i + 1);
|
||||||
@ -144,11 +163,13 @@ void Minimizer::useHighLimit(const Index i, const bool use)
|
|||||||
|
|
||||||
void Minimizer::useHighLimit(const PlaceHolder ph __dumb, const bool use)
|
void Minimizer::useHighLimit(const PlaceHolder ph __dumb, const bool use)
|
||||||
{
|
{
|
||||||
|
checkSupport;
|
||||||
hasHighLimit_.fill(use);
|
hasHighLimit_.fill(use);
|
||||||
}
|
}
|
||||||
|
|
||||||
void Minimizer::useLowLimit(const Index i, const bool use)
|
void Minimizer::useLowLimit(const Index i, const bool use)
|
||||||
{
|
{
|
||||||
|
checkSupport;
|
||||||
if (i >= getDim())
|
if (i >= getDim())
|
||||||
{
|
{
|
||||||
resize(i + 1);
|
resize(i + 1);
|
||||||
@ -158,6 +179,7 @@ void Minimizer::useLowLimit(const Index i, const bool use)
|
|||||||
|
|
||||||
void Minimizer::useLowLimit(const PlaceHolder ph __dumb, const bool use)
|
void Minimizer::useLowLimit(const PlaceHolder ph __dumb, const bool use)
|
||||||
{
|
{
|
||||||
|
checkSupport;
|
||||||
hasLowLimit_.fill(use);
|
hasLowLimit_.fill(use);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -56,6 +56,7 @@ public:
|
|||||||
virtual void useLowLimit(const Index i, const bool use = true);
|
virtual void useLowLimit(const Index i, const bool use = true);
|
||||||
virtual void useLowLimit(const PlaceHolder ph = _,
|
virtual void useLowLimit(const PlaceHolder ph = _,
|
||||||
const bool use = true);
|
const bool use = true);
|
||||||
|
virtual bool supportLimits(void) const = 0;
|
||||||
virtual unsigned int getMaxPass(void) const;
|
virtual unsigned int getMaxPass(void) const;
|
||||||
virtual void setMaxPass(const unsigned int maxPass);
|
virtual void setMaxPass(const unsigned int maxPass);
|
||||||
// minimization
|
// minimization
|
||||||
|
@ -47,6 +47,11 @@ void MinuitMinimizer::setAlgorithm(const Algorithm algorithm)
|
|||||||
algorithm_ = algorithm;
|
algorithm_ = algorithm;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
bool MinuitMinimizer::supportLimits(void) const
|
||||||
|
{
|
||||||
|
return true;
|
||||||
|
}
|
||||||
|
|
||||||
// minimization ////////////////////////////////////////////////////////////////
|
// minimization ////////////////////////////////////////////////////////////////
|
||||||
const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
|
const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
|
||||||
{
|
{
|
||||||
@ -72,13 +77,13 @@ const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
|
|||||||
}
|
}
|
||||||
switch (getAlgorithm())
|
switch (getAlgorithm())
|
||||||
{
|
{
|
||||||
case Algorithm::Migrad:
|
case Algorithm::migrad:
|
||||||
minuitAlg = kMigrad;
|
minuitAlg = kMigrad;
|
||||||
break;
|
break;
|
||||||
case Algorithm::Simplex:
|
case Algorithm::simplex:
|
||||||
minuitAlg = kSimplex;
|
minuitAlg = kSimplex;
|
||||||
break;
|
break;
|
||||||
case Algorithm::Combined:
|
case Algorithm::combined:
|
||||||
minuitAlg = kCombined;
|
minuitAlg = kCombined;
|
||||||
break;
|
break;
|
||||||
}
|
}
|
||||||
|
@ -35,9 +35,9 @@ class MinuitMinimizer: public Minimizer
|
|||||||
public:
|
public:
|
||||||
enum class Algorithm
|
enum class Algorithm
|
||||||
{
|
{
|
||||||
Migrad = 1,
|
migrad = 1,
|
||||||
Simplex = 2,
|
simplex = 2,
|
||||||
Combined = 3
|
combined = 3
|
||||||
};
|
};
|
||||||
public:
|
public:
|
||||||
// constructor
|
// constructor
|
||||||
@ -45,13 +45,14 @@ public:
|
|||||||
// destructor
|
// destructor
|
||||||
virtual ~MinuitMinimizer(void) = default;
|
virtual ~MinuitMinimizer(void) = default;
|
||||||
// access
|
// access
|
||||||
Algorithm getAlgorithm(void) const;
|
Algorithm getAlgorithm(void) const;
|
||||||
void setAlgorithm(const Algorithm algorithm);
|
void setAlgorithm(const Algorithm algorithm);
|
||||||
|
virtual bool supportLimits(void) const;
|
||||||
// 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;
|
static constexpr Algorithm defaultAlg_ = Algorithm::combined;
|
||||||
};
|
};
|
||||||
|
|
||||||
END_LATAN_NAMESPACE
|
END_LATAN_NAMESPACE
|
||||||
|
@ -44,6 +44,11 @@ void NloptMinimizer::setAlgorithm(const Algorithm algorithm)
|
|||||||
algorithm_ = algorithm;
|
algorithm_ = algorithm;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
bool NloptMinimizer::supportLimits(void) const
|
||||||
|
{
|
||||||
|
return true;
|
||||||
|
}
|
||||||
|
|
||||||
// minimization ////////////////////////////////////////////////////////////////
|
// minimization ////////////////////////////////////////////////////////////////
|
||||||
const DVec & NloptMinimizer::operator()(const DoubleFunction &f)
|
const DVec & NloptMinimizer::operator()(const DoubleFunction &f)
|
||||||
{
|
{
|
||||||
@ -75,7 +80,7 @@ const DVec & NloptMinimizer::operator()(const DoubleFunction &f)
|
|||||||
min.set_lower_bounds(lb);
|
min.set_lower_bounds(lb);
|
||||||
min.set_upper_bounds(hb);
|
min.set_upper_bounds(hb);
|
||||||
|
|
||||||
// minimise
|
// minimize
|
||||||
double res;
|
double res;
|
||||||
vector<double> vx(x.size());
|
vector<double> vx(x.size());
|
||||||
nlopt::result status;
|
nlopt::result status;
|
||||||
|
@ -52,8 +52,9 @@ public:
|
|||||||
// destructor
|
// destructor
|
||||||
virtual ~NloptMinimizer(void) = default;
|
virtual ~NloptMinimizer(void) = default;
|
||||||
// access
|
// access
|
||||||
Algorithm getAlgorithm(void) const;
|
Algorithm getAlgorithm(void) const;
|
||||||
void setAlgorithm(const Algorithm algorithm);
|
void setAlgorithm(const Algorithm algorithm);
|
||||||
|
virtual bool supportLimits(void) const;
|
||||||
// minimization
|
// minimization
|
||||||
virtual const DVec & operator()(const DoubleFunction &f);
|
virtual const DVec & operator()(const DoubleFunction &f);
|
||||||
private:
|
private:
|
||||||
|
@ -313,16 +313,19 @@ FitResult XYStatData::fit(vector<Minimizer *> &minimizer, const DVec &init,
|
|||||||
for (auto &m: minimizer)
|
for (auto &m: minimizer)
|
||||||
{
|
{
|
||||||
m->setInit(totalInit);
|
m->setInit(totalInit);
|
||||||
//// do not allow more than maxXsiDev std. deviations on the x-axis
|
if (m->supportLimits())
|
||||||
for (Index p = nPar; p < totalNPar; ++p)
|
|
||||||
{
|
{
|
||||||
double err;
|
//// do not allow more than maxXsiDev std. deviations on the x-axis
|
||||||
|
for (Index p = nPar; p < totalNPar; ++p)
|
||||||
err = sqrt(fitVar_.diagonal()(layout.totalYSize + p - nPar));
|
{
|
||||||
m->useLowLimit(p);
|
double err;
|
||||||
m->useHighLimit(p);
|
|
||||||
m->setLowLimit(p, totalInit(p) - maxXsiDev*err);
|
err = sqrt(fitVar_.diagonal()(layout.totalYSize + p - nPar));
|
||||||
m->setHighLimit(p, totalInit(p) + maxXsiDev*err);
|
m->useLowLimit(p);
|
||||||
|
m->useHighLimit(p);
|
||||||
|
m->setLowLimit(p, totalInit(p) - maxXsiDev*err);
|
||||||
|
m->setHighLimit(p, totalInit(p) + maxXsiDev*err);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
//// minimize and store results
|
//// minimize and store results
|
||||||
result = (*m)(chi2);
|
result = (*m)(chi2);
|
||||||
|
Loading…
x
Reference in New Issue
Block a user