ROOT  6.06/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 #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) ),
49  fGrad( std::vector<double>(f.NDim() ) )
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),
61  fGrad( std::vector<double>(f.NDim() ) )
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 );
78  // use gradient
79  double val = fFunc.DataElement( xExt, i, &fGrad[0]);
80  // transform gradient
81  fTransform->GradientTransformation( x, &fGrad.front(), g);
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
128  mutable std::vector<double> fGrad; // cached vector of gradient values
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
223  MultiNumGradFunction * gradFunction = new MultiNumGradFunction(*fChi2Func);
224  MinimTransformFunction * trFuncRaw = CreateTransformation(startValues, gradFunction);
225  // need to transform in a FitTransformFunction which is set in the residual functions
226  std::auto_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)
368  return fGSLMultiFit->Gradient();
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
void SetMaxIterations(unsigned int maxiter)
set maximum iterations (one iteration can have many function calls)
Definition: Minimizer.h:459
~GSLNLSMinimizer()
Destructor (no operations)
virtual const double * MinGradient() const
return pointer to gradient values at the minimum
int TestDelta(double absTol, double relTol) const
test using abs and relative tolerance |dx| < absTol + relTol*|x| for every component ...
Definition: GSLMultiFit.h:183
std::vector< double > fErrors
Namespace for new ROOT classes and functions.
Definition: ROOT.py:1
const ROOT::Math::FitMethodFunction * fChi2Func
GSLMultiFit, internal class for implementing GSL non linear least square GSL fitting.
Definition: GSLMultiFit.h:52
const double * Gradient() const
gradient value at the minimum
Definition: GSLMultiFit.h:145
#define assert(cond)
Definition: unittest.h:542
unsigned int MaxIterations() const
max iterations
Definition: Minimizer.h:425
virtual unsigned int NPoints() const
return the number of data points used in evaluating the function
MinimTransformFunction class to perform a transformations on the variables to deal with fixed or limi...
TAlienJobStatus * status
Definition: TAlienJob.cxx:51
MinimTransformFunction * CreateTransformation(std::vector< double > &startValues, const ROOT::Math::IMultiGradFunction *func=0)
ROOT::Math::GSLMultiFit * fGSLMultiFit
double Tolerance() const
absolute tolerance
Definition: Minimizer.h:428
STL namespace.
int Set(const std::vector< Func > &funcVec, const double *x)
set the solver parameters
Definition: GSLMultiFit.h:105
std::vector< double > fCovMatrix
ClassImp(TIterator) Bool_t TIterator return false
Compare two iterator objects.
Definition: TIterator.cxx:20
#define MATH_ERROR_MSGVAL(loc, str, x)
Definition: Error.h:69
std::string Name() const
Definition: GSLMultiFit.h:127
const double * CovarMatrix() const
return covariance matrix of the parameters
Definition: GSLMultiFit.h:156
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 * X() const
parameter values at the minimum
Definition: GSLMultiFit.h:138
std::map< std::string, std::string >::const_iterator iter
Definition: TAlienJob.cxx:54
void SetMinValue(double val)
virtual int CovMatrixStatus() const
return covariance matrix status
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:50
double Edm() const
Definition: GSLMultiFit.h:189
virtual unsigned int NDim() const
number of dimensions
const int NPoints
Definition: testNdimFit.cxx:35
GSLNLSMinimizer(int type=0)
Default constructor.
void SetFinalValues(const double *x)
virtual double MinValue() const
return minimum function value
virtual std::string VariableName(unsigned int ivar) const
get name of variables (override if minimizer support storing of variable names)
virtual const double * X() const
return pointer to X values at the minimum
virtual unsigned int NDim() const =0
Retrieve the dimension of the function.
int PrintLevel() const
minimizer configuration parameters
Definition: Minimizer.h:419
virtual bool Minimize()
method to perform the minimization
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)
set the function to minimize
MultiNumGradFunction class to wrap a normal function in a gradient function using numerical gradient ...
virtual unsigned int NPar() const
total number of parameter defined
int TestGradient(double absTol) const
test gradient (ask from solver gradient vector)
Definition: GSLMultiFit.h:175
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 ROOT::Math::IMultiGenFunction * ObjFunction() const
return pointer to used objective function
int type
Definition: TGX11.cxx:120
double func(double *x, double *p)
Definition: stressTF1.cxx:213
IBaseFunctionMultiDim IMultiGenFunction
Definition: IFunctionfwd.h:28
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...
std::vector< LSResidualFunc > fResiduals
Namespace for new Math classes and functions.
Binding & operator=(OUT(*fun)(void))
virtual unsigned int NCalls() const
return the total number of function calls (overrided if needed)
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