Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
GSLNLSMinimizer.cxx
Go to the documentation of this file.
1// @(#)root/mathmore:$Id$
2// Author: L. Moneta Wed Dec 20 17:16:32 2006
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
7 * *
8 * *
9 **********************************************************************/
10
11// Implementation file for class GSLNLSMinimizer
12
14
17
18#include "Math/Error.h"
19#include "GSLMultiFit.h"
20#include "gsl/gsl_errno.h"
21
22
24//#include "Math/Derivator.h"
25
26#include <iostream>
27#include <iomanip>
28#include <cassert>
29
30namespace ROOT {
31
32 namespace Math {
33
34
35/// Internal class used by GSLNLSMinimizer to implement the transformation of the chi2
36/// function used by GSL Non-linear Least-square fitting
37/// The class is template on the FitMethodFunction type to support both gradient and non
38/// gradient functions
39template<class FMFunc>
41
42public:
43
44 FitTransformFunction(const FMFunc & f, std::unique_ptr<MinimTransformFunction> transFunc ) :
45 FMFunc( f.NDim(), f.NPoints() ),
46 fFunc(f),
47 fTransform(std::move(transFunc)),
48 fGrad( std::vector<double>(f.NDim() ) )
49 {
50 // constructor from a given FitMethodFunction and Transformation object.
51 // Ownership of the transformation object is passed to this class
52 }
53
55 }
56
57 // re-implement data element
58 virtual double DataElement(const double * x, unsigned i, double * g = nullptr, double * = nullptr, bool = false) const {
59 // transform from x internal to x external
60 const double * xExt = fTransform->Transformation(x);
61 if ( g == nullptr) return fFunc.DataElement( xExt, i );
62 // use gradient
63 double val = fFunc.DataElement( xExt, i, &fGrad[0]);
64 // transform gradient
65 fTransform->GradientTransformation( x, &fGrad.front(), g);
66 return val;
67 }
68
69
70 virtual IMultiGenFunction * Clone() const {
71 // not supported
72 return nullptr;
73 }
74
75 // dimension (this is number of free dimensions)
76 virtual unsigned int NDim() const {
77 return fTransform->NDim();
78 }
79
80 unsigned int NTot() const {
81 return fTransform->NTot();
82 }
83
84 // forward of transformation functions
85 const double * Transformation( const double * x) const { return fTransform->Transformation(x); }
86
87
88 /// inverse transformation (external -> internal)
89 void InvTransformation(const double * xext, double * xint) const { fTransform->InvTransformation(xext,xint); }
90
91 /// inverse transformation for steps (external -> internal) at external point x
92 void InvStepTransformation(const double * x, const double * sext, double * sint) const { fTransform->InvStepTransformation(x,sext,sint); }
93
94 ///transform gradient vector (external -> internal) at internal point x
95 void GradientTransformation(const double * x, const double *gext, double * gint) const { fTransform->GradientTransformation(x,gext,gint); }
96
97 void MatrixTransformation(const double * x, const double *cint, double * cext) const { fTransform->MatrixTransformation(x,cint,cext); }
98
99private:
100
101 // objects of this class are not meant for copying or assignment
104
105 virtual double DoEval(const double * x) const {
106 return fFunc( fTransform->Transformation(x) );
107 }
108
109 virtual double DoDerivative(const double * /* x */, unsigned int /*icoord*/) const {
110 // not used
111 throw std::runtime_error("FitTransformFunction::DoDerivative");
112 return 0;
113 }
114
116 const FMFunc & fFunc; // pointer to original fit method function
117 std::unique_ptr<MinimTransformFunction> fTransform; // pointer to transformation function
118 mutable std::vector<double> fGrad; // cached vector of gradient values
119
120};
121
122//________________________________________________________________________________
123/**
124 LSResidualFunc class description.
125 Internal class used for accessing the residuals of the Least Square function
126 and their derivatives which are estimated numerically using GSL numerical derivation.
127 The class contains a pointer to the fit method function and an index specifying
128 the i-th residual and wraps it in a multi-dim gradient function interface
129 ROOT::Math::IGradientFunctionMultiDim.
130 The class is used by ROOT::Math::GSLNLSMinimizer (GSL non linear least square fitter)
131
132 @ingroup MultiMin
133*/
134template<class Func>
136public:
137
138 //default ctor (required by CINT)
140 {}
141
142
143 LSResidualFunc(const Func & func, unsigned int i) :
144 fIndex(i),
145 fChi2(&func)
146 {}
147
148
149 // copy ctor
153 {
154 operator=(rhs);
155 }
156
157 // assignment
159 {
160 fIndex = rhs.fIndex;
161 fChi2 = rhs.fChi2;
162 return *this;
163 }
164
165 IMultiGenFunction * Clone() const override {
166 return new LSResidualFunc<Func>(*fChi2,fIndex);
167 }
168
169 unsigned int NDim() const override { return fChi2->NDim(); }
170
171 void Gradient( const double * x, double * g) const override {
172 double f0 = 0;
173 FdF(x,f0,g);
174 }
175
176 void FdF (const double * x, double & f, double * g) const override {
177 f = fChi2->DataElement(x,fIndex,g);
178 // for likelihood fits need to rescale g ??
179 // if (fChi2->Type() == Func::kPoissonLikelihood) {
180 // f *= -1;
181 // for (unsigned int i = 0; i < NDim(); i++)
182 // g[i] *= -1.; // /= f;
183 // }
184 }
185
186
187private:
188
189 double DoEval (const double * x) const override {
190 return fChi2->DataElement(x, fIndex, nullptr);
191 }
192
193 double DoDerivative(const double * /* x */, unsigned int /* icoord */) const override {
194 //this function should not be called by GSL
195 throw std::runtime_error("LSRESidualFunc::DoDerivative");
196 return 0;
197 }
198
199 unsigned int fIndex;
200 const Func * fChi2;
201};
202
203
204// GSLNLSMinimizer implementation
205
207{
208 // Constructor implementation : create GSLMultiFit wrapper object
209 const gsl_multifit_fdfsolver_type * gsl_type = nullptr; // use default type defined in GSLMultiFit
210 if (type == 1) gsl_type = gsl_multifit_fdfsolver_lmsder; // scaled lmder version
211 if (type == 2) gsl_type = gsl_multifit_fdfsolver_lmder; // unscaled version
212
213 fGSLMultiFit = new GSLMultiFit( gsl_type );
214
215 fEdm = -1;
216
217 // default tolerance and max iterations
219 if (niter <= 0) niter = 100;
220 SetMaxIterations(niter);
221
223 if (fLSTolerance <=0) fLSTolerance = 0.0001; // default internal value
224
226}
227
229 assert(fGSLMultiFit != nullptr);
230 delete fGSLMultiFit;
231}
232
233
234
236 // set the function to minimizer
237 // need to create vector of functions to be passed to GSL multifit
238 // support now only CHi2 implementation
239
240 // call base class method. It will clone the function and set number of dimensions
242 fNFree = NDim();
244}
245
246
248
249 if (ObjFunction() == nullptr) {
250 MATH_ERROR_MSG("GSLNLSMinimizer::Minimize","Function has not been set");
251 return false;
252 }
253 // check type of function (if it provides gradient)
254 auto fitFunc = (!fUseGradFunction) ? dynamic_cast<const ROOT::Math::FitMethodFunction *>(ObjFunction()) : nullptr;
255 auto fitGradFunc = (fUseGradFunction) ? dynamic_cast<const ROOT::Math::FitMethodGradFunction *>(ObjFunction()) : nullptr;
256 if (fitFunc == nullptr && fitGradFunc == nullptr) {
257 if (PrintLevel() > 0) std::cout << "GSLNLSMinimizer: Invalid function set - only FitMethodFunction types are supported" << std::endl;
258 return false;
259 }
260
261 if (fitGradFunc)
262 return DoMinimize<ROOT::Math::FitMethodGradFunction>(*fitGradFunc);
263 else
264 return DoMinimize<ROOT::Math::FitMethodFunction>(*fitFunc);
265}
266
267template<class Func>
268bool GSLNLSMinimizer::DoMinimize(const Func & fitFunc) {
269
270 unsigned int size = fitFunc.NPoints();
271 fNCalls = 0; // reset number of function calls
272
273 std::vector<LSResidualFunc<Func>> residualFuncs;
274 residualFuncs.reserve(size);
275
276 // set initial parameters of the minimizer
277 int debugLevel = PrintLevel();
278
279 assert (fGSLMultiFit != nullptr);
280
281 unsigned int npar = NPar();
282 unsigned int ndim = NDim();
283 if (npar == 0 || npar < ndim) {
284 MATH_ERROR_MSGVAL("GSLNLSMinimizer::Minimize","Wrong number of parameters",npar);
285 return false;
286 }
287
288 // set residual functions and check if a transformation is needed
289 std::vector<double> startValues;
290
291 // transformation need a grad function.
292 std::unique_ptr<MultiNumGradFunction> gradFunction;
293 std::unique_ptr<MinimTransformFunction> trFuncRaw;
294 if (!fUseGradFunction) {
295 gradFunction = std::make_unique<MultiNumGradFunction>(fitFunc);
296 trFuncRaw.reset(CreateTransformation(startValues, gradFunction.get()));
297 }
298 else {
299 // use pointer stored in BasicMinimizer
300 trFuncRaw.reset(CreateTransformation(startValues));
301 }
302 // need to transform in a FitTransformFunction which is set in the residual functions
303 std::unique_ptr<FitTransformFunction<Func>> trFunc;
304 if (trFuncRaw) {
305 //pass ownership of trFuncRaw to FitTransformFunction
306 trFunc = std::make_unique<FitTransformFunction<Func>>(fitFunc, std::move(trFuncRaw));
307 assert(npar == trFunc->NTot() );
308 for (unsigned int ires = 0; ires < size; ++ires) {
309 residualFuncs.emplace_back(LSResidualFunc<Func>(*trFunc, ires));
310 }
311 } else {
312 for (unsigned int ires = 0; ires < size; ++ires) {
313 residualFuncs.emplace_back( LSResidualFunc<Func>(fitFunc, ires) );
314 }
315 }
316
317 if (debugLevel >=1 ) std::cout <<"Minimize using GSLNLSMinimizer " << std::endl;
318
319
320 int iret = fGSLMultiFit->Set( residualFuncs, &startValues.front() );
321
322 if (iret) {
323 MATH_ERROR_MSGVAL("GSLNLSMinimizer::Minimize","Error setting the residual functions ",iret);
324 return false;
325 }
326
327 if (debugLevel >=1 ) std::cout <<"GSLNLSMinimizer: " << fGSLMultiFit->Name() << " - start iterating......... " << std::endl;
328
329 // start iteration
330 unsigned int iter = 0;
331 int status;
332 bool minFound = false;
333 do {
334 status = fGSLMultiFit->Iterate();
335
336 if (debugLevel >=1) {
337 std::cout << "----------> Iteration " << iter << " / " << MaxIterations() << " status " << gsl_strerror(status) << std::endl;
338 const double * x = fGSLMultiFit->X();
339 if (trFunc) x = trFunc->Transformation(x);
340 int pr = std::cout.precision(18);
341 std::cout << " FVAL = " << (fitFunc)(x) << std::endl;
342 std::cout.precision(pr);
343 std::cout << " X Values : ";
344 for (unsigned int i = 0; i < NDim(); ++i)
345 std::cout << " " << VariableName(i) << " = " << X()[i];
346 std::cout << std::endl;
347 }
348
349 if (status) break;
350
351 // check also the delta in X()
352 status = fGSLMultiFit->TestDelta( Tolerance(), Tolerance() );
353 if (status == GSL_SUCCESS) {
354 minFound = true;
355 }
356
357 // double-check with the gradient
358 int status2 = fGSLMultiFit->TestGradient( Tolerance() );
359 if ( minFound && status2 != GSL_SUCCESS) {
360 // check now edm
361 fEdm = fGSLMultiFit->Edm();
362 if (fEdm > Tolerance() ) {
363 // continue the iteration
364 status = status2;
365 minFound = false;
366 }
367 }
368
369 if (debugLevel >=1) {
370 std::cout << " after Gradient and Delta tests: " << gsl_strerror(status);
371 if (fEdm > 0) std::cout << ", edm is: " << fEdm;
372 std::cout << std::endl;
373 }
374
375 iter++;
376
377 }
378 while (status == GSL_CONTINUE && iter < MaxIterations() );
379
380 // check edm
381 fEdm = fGSLMultiFit->Edm();
382 if ( fEdm < Tolerance() ) {
383 minFound = true;
384 }
385
386 // save state with values and function value
387 const double * x = fGSLMultiFit->X();
388 if (x == nullptr) return false;
389 // apply transformation outside SetFinalValues(..)
390 // because trFunc is not a MinimTransformFunction but a FitTransFormFunction
391 if (trFunc) x = trFunc->Transformation(x);
393
394 SetMinValue( (fitFunc)(x) );
395 fStatus = status;
396 fNCalls = fitFunc.NCalls();
397 fErrors.resize(NDim());
398
399 // get errors from cov matrix
400 const double * cov = fGSLMultiFit->CovarMatrix();
401 if (cov) {
402
403 fCovMatrix.resize(ndim*ndim);
404
405 if (trFunc) {
406 trFunc->MatrixTransformation(x, fGSLMultiFit->CovarMatrix(), fCovMatrix.data() );
407 }
408 else {
409 std::copy(cov, cov + fCovMatrix.size(), fCovMatrix.begin() );
410 }
411
412 for (unsigned int i = 0; i < ndim; ++i)
413 fErrors[i] = std::sqrt(fCovMatrix[i*ndim + i]);
414 }
415
416 if (minFound) {
417
418 if (debugLevel >=1 ) {
419 std::cout << "GSLNLSMinimizer: Minimum Found" << std::endl;
420 int pr = std::cout.precision(18);
421 std::cout << "FVAL = " << MinValue() << std::endl;
422 std::cout << "Edm = " << fEdm << std::endl;
423 std::cout.precision(pr);
424 std::cout << "NIterations = " << iter << std::endl;
425 std::cout << "NFuncCalls = " << fitFunc.NCalls() << std::endl;
426 for (unsigned int i = 0; i < NDim(); ++i)
427 std::cout << std::setw(12) << VariableName(i) << " = " << std::setw(12) << X()[i] << " +/- " << std::setw(12) << fErrors[i] << std::endl;
428 }
429
430 return true;
431 }
432 else {
433 if (debugLevel >=0 ) {
434 std::cout << "GSLNLSMinimizer: Minimization did not converge: " << std::endl;
435 if (status == GSL_ENOPROG) // case status 27
436 std::cout << "\t iteration is not making progress towards solution" << std::endl;
437 else
438 std::cout << "\t failed with status " << status << std::endl;
439 }
440 if (debugLevel >=1 ) {
441 std::cout << "FVAL = " << MinValue() << std::endl;
442 std::cout << "Edm = " << fGSLMultiFit->Edm() << std::endl;
443 std::cout << "Niterations = " << iter << std::endl;
444 }
445 return false;
446 }
447 return false;
448}
449
450const double * GSLNLSMinimizer::MinGradient() const {
451 // return gradient (internal values)
452 return fGSLMultiFit->Gradient();
453}
454
455
456double GSLNLSMinimizer::CovMatrix(unsigned int i , unsigned int j ) const {
457 // return covariance matrix element
458 unsigned int ndim = NDim();
459 if ( fCovMatrix.empty()) return 0;
460 if (i > ndim || j > ndim) return 0;
461 return fCovMatrix[i*ndim + j];
462}
463
465 // return covariance matrix status = 0 not computed,
466 // 1 computed but is approximate because minimum is not valid, 3 is fine
467 if ( fCovMatrix.empty()) return 0;
468 // case minimization did not finished correctly
469 if (fStatus != GSL_SUCCESS) return 1;
470 return 3;
471}
472
473
474 } // end namespace Math
475
476} // end namespace ROOT
477
#define MATH_ERROR_MSGVAL(loc, txt, x)
Definition Error.h:109
#define MATH_ERROR_MSG(loc, str)
Definition Error.h:83
#define f(i)
Definition RSha256.hxx:104
#define g(i)
Definition RSha256.hxx:105
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
FitMethodFunction class Interface for objective functions (like chi2 and likelihood used in the fit) ...
virtual unsigned int NPar() const
total number of parameter defined
unsigned int NDim() const override
number of dimensions
void SetFinalValues(const double *x, const MinimTransformFunction *func=nullptr)
double MinValue() const override
return minimum function value
MinimTransformFunction * CreateTransformation(std::vector< double > &startValues, const ROOT::Math::IMultiGradFunction *func=nullptr)
void SetFunction(const ROOT::Math::IMultiGenFunction &func) override
set the function to minimize
const ROOT::Math::IMultiGenFunction * ObjFunction() const
return pointer to used objective function
const double * X() const override
return pointer to X values at the minimum
std::string VariableName(unsigned int ivar) const override
get name of variables (override if minimizer support storing of variable names)
Internal class used by GSLNLSMinimizer to implement the transformation of the chi2 function used by G...
FitTransformFunction & operator=(const FitTransformFunction &rhs)=delete
void InvTransformation(const double *xext, double *xint) const
inverse transformation (external -> internal)
std::unique_ptr< MinimTransformFunction > fTransform
void InvStepTransformation(const double *x, const double *sext, double *sint) const
inverse transformation for steps (external -> internal) at external point x
virtual IMultiGenFunction * Clone() const
void MatrixTransformation(const double *x, const double *cint, double *cext) const
void GradientTransformation(const double *x, const double *gext, double *gint) const
transform gradient vector (external -> internal) at internal point x
FitTransformFunction(const FitTransformFunction &rhs)=delete
virtual double DoEval(const double *x) const
FitTransformFunction(const FMFunc &f, std::unique_ptr< MinimTransformFunction > transFunc)
virtual double DataElement(const double *x, unsigned i, double *g=nullptr, double *=nullptr, bool=false) const
virtual unsigned int NDim() const
virtual double DoDerivative(const double *, unsigned int) const
const double * Transformation(const double *x) const
GSLMultiFit, internal class for implementing GSL non linear least square GSL fitting.
Definition GSLMultiFit.h:53
int TestGradient(double absTol) const
test gradient (ask from solver gradient vector)
int TestDelta(double absTol, double relTol) const
test using abs and relative tolerance |dx| < absTol + relTol*|x| for every component
const double * Gradient() const
gradient value at the minimum
std::string Name() const
const double * CovarMatrix() const
return covariance matrix of the parameters
const double * X() const
parameter values at the minimum
int Set(const std::vector< Func > &funcVec, const double *x)
set the solver parameters
double CovMatrix(unsigned int, unsigned int) const override
return covariance matrices elements if the variable is fixed the matrix is zero The ordering of the v...
int CovMatrixStatus() const override
return covariance matrix status
void SetFunction(const ROOT::Math::IMultiGenFunction &func) override
set the function to minimize
std::vector< double > fErrors
std::vector< double > fCovMatrix
~GSLNLSMinimizer() override
Destructor (no operations)
bool DoMinimize(const Func &f)
Internal method to perform minimization template on the type of method function.
bool Minimize() override
method to perform the minimization
GSLNLSMinimizer(int type=0)
Default constructor.
const double * MinGradient() const override
return pointer to gradient values at the minimum
ROOT::Math::GSLMultiFit * fGSLMultiFit
Documentation for the abstract class IBaseFunctionMultiDim.
Definition IFunction.h:61
Interface (abstract class) for multi-dimensional functions providing a gradient calculation.
Definition IFunction.h:168
LSResidualFunc class description.
double DoEval(const double *x) const override
LSResidualFunc< Func > & operator=(const LSResidualFunc< Func > &rhs)
void FdF(const double *x, double &f, double *g) const override
LSResidualFunc(const LSResidualFunc< Func > &rhs)
void Gradient(const double *x, double *g) const override
double DoDerivative(const double *, unsigned int) const override
LSResidualFunc(const Func &func, unsigned int i)
IMultiGenFunction * Clone() const override
Clone a function.
unsigned int NDim() const override
Retrieve the dimension of the function.
double Tolerance() const
absolute tolerance
Definition Minimizer.h:315
void SetMaxIterations(unsigned int maxiter)
set maximum iterations (one iteration can have many function calls)
Definition Minimizer.h:349
int fStatus
status of minimizer
Definition Minimizer.h:391
unsigned int MaxIterations() const
max iterations
Definition Minimizer.h:312
void SetPrintLevel(int level)
set print level
Definition Minimizer.h:343
int PrintLevel() const
minimizer configuration parameters
Definition Minimizer.h:306
Double_t x[n]
Definition legend1.C:17
Namespace for new Math classes and functions.
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.