ROOT   6.10/09 Reference Guide
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
13 #include "Math/GSLNLSMinimizer.h"
14
17
18 #include "Math/Error.h"
19 #include "GSLMultiFit.h"
20 #include "gsl/gsl_errno.h"
21
22
23 #include "Math/FitMethodFunction.h"
24 //#include "Math/Derivator.h"
25
26 #include <iostream>
27 #include <iomanip>
28 #include <cassert>
29
30 namespace ROOT {
31
32  namespace Math {
33
34
35 // class to implement transformation of chi2 function
36 // in general could make template on the fit method function type
37
38 class FitTransformFunction : public FitMethodFunction {
39
40 public:
41
42  FitTransformFunction(const FitMethodFunction & f, const std::vector<EMinimVariableType> & types, const std::vector<double> & values,
43  const std::map<unsigned int, std::pair<double, double> > & bounds) :
44  FitMethodFunction( f.NDim(), f.NPoints() ),
45  fOwnTransformation(true),
46  fFunc(f),
47  fTransform(new MinimTransformFunction( new MultiNumGradFunction(f), types, values, bounds) ),
49  {
50  // constructor
51  // need to pass to MinimTransformFunction a new pointer which will be managed by the class itself
52  // pass a gradient pointer although it will not be used byb the class
53  }
54
55  FitTransformFunction(const FitMethodFunction & f, MinimTransformFunction *transFunc ) :
56  FitMethodFunction( f.NDim(), f.NPoints() ),
57  fOwnTransformation(false),
58  fFunc(f),
59  fTransform(transFunc),
61  {
62  // constructor from al already existing Transformation object. Ownership of the transformation onbect is passed to caller
63  }
64
65  ~FitTransformFunction() {
66  if (fOwnTransformation) {
67  assert(fTransform);
68  delete fTransform;
69  }
70  }
71
72  // re-implement data element
73  virtual double DataElement(const double * x, unsigned i, double * g = 0) const {
74  // transform from x internal to x external
75  const double * xExt = fTransform->Transformation(x);
76  if ( g == 0) return fFunc.DataElement( xExt, i );
78  double val = fFunc.DataElement( xExt, i, &fGrad[0]);
81  return val;
82  }
83
84
85  IMultiGenFunction * Clone() const {
86  // not supported
87  return 0;
88  }
89
90  // dimension (this is number of free dimensions)
91  unsigned int NDim() const {
92  return fTransform->NDim();
93  }
94
95  unsigned int NTot() const {
96  return fTransform->NTot();
97  }
98
99  // forward of transformation functions
100  const double * Transformation( const double * x) const { return fTransform->Transformation(x); }
101
102
103  /// inverse transformation (external -> internal)
104  void InvTransformation(const double * xext, double * xint) const { fTransform->InvTransformation(xext,xint); }
105
106  /// inverse transformation for steps (external -> internal) at external point x
107  void InvStepTransformation(const double * x, const double * sext, double * sint) const { fTransform->InvStepTransformation(x,sext,sint); }
108
109  ///transform gradient vector (external -> internal) at internal point x
110  void GradientTransformation(const double * x, const double *gext, double * gint) const { fTransform->GradientTransformation(x,gext,gint); }
111
112  void MatrixTransformation(const double * x, const double *cint, double * cext) const { fTransform->MatrixTransformation(x,cint,cext); }
113
114 private:
115
116  // objects of this class are not meant for copying or assignment
117  FitTransformFunction(const FitTransformFunction& rhs);
118  FitTransformFunction& operator=(const FitTransformFunction& rhs);
119
120  double DoEval(const double * x) const {
121  return fFunc( fTransform->Transformation(x) );
122  }
123
124  bool fOwnTransformation;
125  const FitMethodFunction & fFunc; // pointer to original fit method function
126  MinimTransformFunction * fTransform; // pointer to transformation function
128
129 };
130
131
132
133
134 // GSLNLSMinimizer implementation
135
137  //fNFree(0),
138  fSize(0),
139  fChi2Func(0)
140 {
141  // Constructor implementation : create GSLMultiFit wrapper object
142  const gsl_multifit_fdfsolver_type * gsl_type = 0; // use default type defined in GSLMultiFit
143  if (type == 1) gsl_type = gsl_multifit_fdfsolver_lmsder; // scaled lmder version
144  if (type == 2) gsl_type = gsl_multifit_fdfsolver_lmder; // unscaled version
145
146  fGSLMultiFit = new GSLMultiFit( gsl_type );
147
148  fEdm = -1;
149
150  // default tolerance and max iterations
152  if (niter <= 0) niter = 100;
153  SetMaxIterations(niter);
154
156  if (fLSTolerance <=0) fLSTolerance = 0.0001; // default internal value
157
159 }
160
162  assert(fGSLMultiFit != 0);
163  delete fGSLMultiFit;
164 }
165
166
167
169  // set the function to minimizer
170  // need to create vector of functions to be passed to GSL multifit
171  // support now only CHi2 implementation
172
173  // call base class method. It will clone the function and set ndimension
175  //need to check if function can be used
176  const ROOT::Math::FitMethodFunction * chi2Func = dynamic_cast<const ROOT::Math::FitMethodFunction *>(ObjFunction());
177  if (chi2Func == 0) {
178  if (PrintLevel() > 0) std::cout << "GSLNLSMinimizer: Invalid function set - only Chi2Func supported" << std::endl;
179  return;
180  }
181  fSize = chi2Func->NPoints();
182  fNFree = NDim();
183
184  // use vector by value
185  fResiduals.reserve(fSize);
186  for (unsigned int i = 0; i < fSize; ++i) {
187  fResiduals.push_back( LSResidualFunc(*chi2Func, i) );
188  }
189  // keep pointers to the chi2 function
190  fChi2Func = chi2Func;
191  }
192
194  // set the function to minimizer using gradient interface
195  // not supported yet, implemented using the other SetFunction
196  return SetFunction(static_cast<const ROOT::Math::IMultiGenFunction &>(func) );
197 }
198
199
201  // set initial parameters of the minimizer
202  int debugLevel = PrintLevel();
203
204
205  assert (fGSLMultiFit != 0);
206  if (fResiduals.size() != fSize || fChi2Func == 0) {
207  MATH_ERROR_MSG("GSLNLSMinimizer::Minimize","Function has not been set");
208  return false;
209  }
210
211  unsigned int npar = NPar();
212  unsigned int ndim = NDim();
213  if (npar == 0 || npar < ndim) {
214  MATH_ERROR_MSGVAL("GSLNLSMinimizer::Minimize","Wrong number of parameters",npar);
215  return false;
216  }
217
218  // set residual functions and check if a transformation is needed
219  std::vector<double> startValues;
220
221  // transformation need a grad function. Delegate fChi2Func to given object
223  MinimTransformFunction * trFuncRaw = CreateTransformation(startValues, gradFunction);
224  // need to transform in a FitTransformFunction which is set in the residual functions
225  std::unique_ptr<FitTransformFunction> trFunc;
226  if (trFuncRaw) {
227  trFunc.reset(new FitTransformFunction(*fChi2Func, trFuncRaw) );
228  //FitTransformationFunction *trFunc = new FitTransformFunction(*fChi2Func, trFuncRaw);
229  for (unsigned int ires = 0; ires < fResiduals.size(); ++ires) {
230  fResiduals[ires] = LSResidualFunc(*trFunc, ires);
231  }
232
233  assert(npar == trFunc->NTot() );
234  }
235
236  if (debugLevel >=1 ) std::cout <<"Minimize using GSLNLSMinimizer " << std::endl;
237
238 // // use a global step size = min (step vectors)
239 // double stepSize = 1;
240 // for (unsigned int i = 0; i < fSteps.size(); ++i)
241 // //stepSize += fSteps[i];
242 // if (fSteps[i] < stepSize) stepSize = fSteps[i];
243
244  int iret = fGSLMultiFit->Set( fResiduals, &startValues.front() );
245  if (iret) {
246  MATH_ERROR_MSGVAL("GSLNLSMinimizer::Minimize","Error setting the residual functions ",iret);
247  return false;
248  }
249
250  if (debugLevel >=1 ) std::cout <<"GSLNLSMinimizer: " << fGSLMultiFit->Name() << " - start iterating......... " << std::endl;
251
252  // start iteration
253  unsigned int iter = 0;
254  int status;
255  bool minFound = false;
256  do {
257  status = fGSLMultiFit->Iterate();
258
259  if (debugLevel >=1) {
260  std::cout << "----------> Iteration " << iter << " / " << MaxIterations() << " status " << gsl_strerror(status) << std::endl;
261  const double * x = fGSLMultiFit->X();
262  if (trFunc.get()) x = trFunc->Transformation(x);
263  int pr = std::cout.precision(18);
264  std::cout << " FVAL = " << (*fChi2Func)(x) << std::endl;
265  std::cout.precision(pr);
266  std::cout << " X Values : ";
267  for (unsigned int i = 0; i < NDim(); ++i)
268  std::cout << " " << VariableName(i) << " = " << X()[i];
269  std::cout << std::endl;
270  }
271
272  if (status) break;
273
274  // check also the delta in X()
275  status = fGSLMultiFit->TestDelta( Tolerance(), Tolerance() );
276  if (status == GSL_SUCCESS) {
277  minFound = true;
278  }
279
280  // double-check with the gradient
281  int status2 = fGSLMultiFit->TestGradient( Tolerance() );
282  if ( minFound && status2 != GSL_SUCCESS) {
283  // check now edm
284  fEdm = fGSLMultiFit->Edm();
285  if (fEdm > Tolerance() ) {
286  // continue the iteration
287  status = status2;
288  minFound = false;
289  }
290  }
291
292  if (debugLevel >=1) {
293  std::cout << " after Gradient and Delta tests: " << gsl_strerror(status);
294  if (fEdm > 0) std::cout << ", edm is: " << fEdm;
295  std::cout << std::endl;
296  }
297
298  iter++;
299
300  }
301  while (status == GSL_CONTINUE && iter < MaxIterations() );
302
303  // check edm
304  fEdm = fGSLMultiFit->Edm();
305  if ( fEdm < Tolerance() ) {
306  minFound = true;
307  }
308
309  // save state with values and function value
310  const double * x = fGSLMultiFit->X();
311  if (x == 0) return false;
312
313  SetFinalValues(x);
314
315  SetMinValue( (*fChi2Func)(x) );
316  fStatus = status;
317
318  fErrors.resize(NDim());
319
320  // get errors from cov matrix
321  const double * cov = fGSLMultiFit->CovarMatrix();
322  if (cov) {
323
324  fCovMatrix.resize(ndim*ndim);
325
326  if (trFunc.get() ) {
327  trFunc->MatrixTransformation(x, fGSLMultiFit->CovarMatrix(), &fCovMatrix[0] );
328  }
329  else {
330  std::copy(cov, cov + ndim*ndim, fCovMatrix.begin() );
331  }
332
333  for (unsigned int i = 0; i < ndim; ++i)
334  fErrors[i] = std::sqrt(fCovMatrix[i*ndim + i]);
335  }
336
337  if (minFound) {
338
339  if (debugLevel >=1 ) {
340  std::cout << "GSLNLSMinimizer: Minimum Found" << std::endl;
341  int pr = std::cout.precision(18);
342  std::cout << "FVAL = " << MinValue() << std::endl;
343  std::cout << "Edm = " << fEdm << std::endl;
344  std::cout.precision(pr);
345  std::cout << "NIterations = " << iter << std::endl;
346  std::cout << "NFuncCalls = " << fChi2Func->NCalls() << std::endl;
347  for (unsigned int i = 0; i < NDim(); ++i)
348  std::cout << std::setw(12) << VariableName(i) << " = " << std::setw(12) << X()[i] << " +/- " << std::setw(12) << fErrors[i] << std::endl;
349  }
350
351  return true;
352  }
353  else {
354  if (debugLevel >=1 ) {
355  std::cout << "GSLNLSMinimizer: Minimization did not converge" << std::endl;
356  std::cout << "FVAL = " << MinValue() << std::endl;
357  std::cout << "Edm = " << fGSLMultiFit->Edm() << std::endl;
358  std::cout << "Niterations = " << iter << std::endl;
359  }
360  return false;
361  }
362  return false;
363 }
364
365 const double * GSLNLSMinimizer::MinGradient() const {
366  // return gradient (internal values)
368 }
369
370
371 double GSLNLSMinimizer::CovMatrix(unsigned int i , unsigned int j ) const {
372  // return covariance matrix element
373  unsigned int ndim = NDim();
374  if ( fCovMatrix.size() == 0) return 0;
375  if (i > ndim || j > ndim) return 0;
376  return fCovMatrix[i*ndim + j];
377 }
378
380  // return covariance matrix status = 0 not computed,
381  // 1 computed but is approximate because minimum is not valid, 3 is fine
382  if ( fCovMatrix.size() == 0) return 0;
383  // case minimization did not finished correctly
384  if (fStatus != GSL_SUCCESS) return 1;
385  return 3;
386 }
387
388
389  } // end namespace Math
390
391 } // end namespace ROOT
392
Interface (abstract class) for multi-dimensional functions providing a gradient calculation.
Definition: IFunction.h:330
virtual const double * X() const
return pointer to X values at the minimum
void SetMaxIterations(unsigned int maxiter)
set maximum iterations (one iteration can have many function calls)
Definition: Minimizer.h:451
~GSLNLSMinimizer()
Destructor (no operations)
std::vector< double > fErrors
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
const ROOT::Math::FitMethodFunction * fChi2Func
virtual double CovMatrix(unsigned int, unsigned int) const
return covariance matrices elements if the variable is fixed the matrix is zero The ordering of the v...
GSLMultiFit, internal class for implementing GSL non linear least square GSL fitting.
Definition: GSLMultiFit.h:52
virtual int CovMatrixStatus() const
return covariance matrix status
int TestDelta(double absTol, double relTol) const
test using abs and relative tolerance |dx| < absTol + relTol*|x| for every component ...
Definition: GSLMultiFit.h:203
MinimTransformFunction class to perform a transformations on the variables to deal with fixed or limi...
virtual unsigned int NPoints() const
return the number of data points used in evaluating the function
MinimTransformFunction * CreateTransformation(std::vector< double > &startValues, const ROOT::Math::IMultiGradFunction *func=0)
ROOT::Math::GSLMultiFit * fGSLMultiFit
STL namespace.
virtual double MinValue() const
return minimum function value
int Set(const std::vector< Func > &funcVec, const double *x)
set the solver parameters
Definition: GSLMultiFit.h:123
std::vector< double > fCovMatrix
virtual std::string VariableName(unsigned int ivar) const
get name of variables (override if minimizer support storing of variable names)
int PrintLevel() const
minimizer configuration parameters
Definition: Minimizer.h:411
#define MATH_ERROR_MSGVAL(loc, str, x)
Definition: Error.h:69
double sqrt(double)
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)
set the function to minimize
Double_t x[n]
Definition: legend1.C:17
const double * CovarMatrix() const
return covariance matrix of the parameters
Definition: GSLMultiFit.h:181
void SetMinValue(double val)
virtual unsigned int NPar() const
total number of parameter defined
unsigned int MaxIterations() const
max iterations
Definition: Minimizer.h:417
Definition: GSLMultiFit.h:170
virtual unsigned int NDim() const
number of dimensions
std::string Name() const
Definition: GSLMultiFit.h:152
const ROOT::Math::IMultiGenFunction * ObjFunction() const
return pointer to used objective function
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:50
double Tolerance() const
absolute tolerance
Definition: Minimizer.h:420
double Edm() const
Definition: GSLMultiFit.h:209
const int NPoints
Definition: testNdimFit.cxx:35
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
GSLNLSMinimizer(int type=0)
Default constructor.
void SetFinalValues(const double *x)
virtual bool Minimize()
method to perform the minimization
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)
set the function to minimize
LSResidualFunc class description.
double f(double x)
FitMethodFunction class Interface for objective functions (like chi2 and likelihood used in the fit) ...
Definition: Fitter.h:40
const double * X() const
parameter values at the minimum
Definition: GSLMultiFit.h:163
int type
Definition: TGX11.cxx:120
double func(double *x, double *p)
Definition: stressTF1.cxx:213
IMultiGenFunctionTempl< double > IMultiGenFunction
Definition: IFunctionfwd.h:30
std::vector< LSResidualFunc > fResiduals
Namespace for new Math classes and functions.
virtual unsigned int NCalls() const
return the total number of function calls (overrided if needed)
Binding & operator=(OUT(*fun)(void))
virtual const double * MinGradient() const
return pointer to gradient values at the minimum
void SetPrintLevel(int level)
set print level
Definition: Minimizer.h:445
virtual double DataElement(const double *x, unsigned int i, double *g=0) const =0
method returning the data i-th contribution to the fit objective function For example the residual fo...
BasicFitMethodFunction< ROOT::Math::IMultiGenFunction > FitMethodFunction
Definition: Fitter.h:40