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