Logo ROOT  
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) ),
48  fGrad( std::vector<double>(f.NDim() ) )
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),
60  fGrad( std::vector<double>(f.NDim() ) )
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 );
77  // use gradient
78  double val = fFunc.DataElement( xExt, i, &fGrad[0]);
79  // transform gradient
80  fTransform->GradientTransformation( x, &fGrad.front(), g);
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
127  mutable std::vector<double> fGrad; // cached vector of gradient values
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
222  MultiNumGradFunction * gradFunction = new MultiNumGradFunction(*fChi2Func);
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)
367  return fGSLMultiFit->Gradient();
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 
ROOT::Math::GSLMultiFit::TestDelta
int TestDelta(double absTol, double relTol) const
test using abs and relative tolerance |dx| < absTol + relTol*|x| for every component
Definition: GSLMultiFit.h:261
ROOT::Math::BasicMinimizer::SetMinValue
void SetMinValue(double val)
Definition: BasicMinimizer.h:185
ROOT::Math::GSLMultiFit::Iterate
int Iterate()
Definition: GSLMultiFit.h:215
ROOT::Math::BasicMinimizer::VariableName
virtual std::string VariableName(unsigned int ivar) const
get name of variables (override if minimizer support storing of variable names)
Definition: BasicMinimizer.cxx:237
GSLMultiFit.h
f
#define f(i)
Definition: RSha256.hxx:122
ROOT::Math::GSLNLSMinimizer::Minimize
virtual bool Minimize()
method to perform the minimization
Definition: GSLNLSMinimizer.cxx:210
ROOT::Math::GSLMultiFit::X
const double * X() const
parameter values at the minimum
Definition: GSLMultiFit.h:221
ROOT::Math::Minimizer::PrintLevel
int PrintLevel() const
minimizer configuration parameters
Definition: Minimizer.h:412
GSLNLSMinimizer.h
ROOT::Math::BasicMinimizer::NPar
virtual unsigned int NPar() const
total number of parameter defined
Definition: BasicMinimizer.h:158
ROOT::Math::BasicMinimizer::SetFunction
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)
set the function to minimize
Definition: BasicMinimizer.cxx:250
ROOT::Math::GSLMultiFit
GSLMultiFit, internal class for implementing GSL non linear least square GSL fitting.
Definition: GSLMultiFit.h:91
ROOT::Math::GSLNLSMinimizer::GSLNLSMinimizer
GSLNLSMinimizer(int type=0)
Default constructor.
Definition: GSLNLSMinimizer.cxx:146
ROOT::Math::Minimizer::MaxIterations
unsigned int MaxIterations() const
max iterations
Definition: Minimizer.h:418
GSL_SUCCESS
#define GSL_SUCCESS
Definition: RooAdaptiveGaussKronrodIntegrator1D.cxx:400
operator=
Binding & operator=(OUT(*fun)(void))
Definition: TRInterface_Binding.h:11
ROOT::Math::GSLNLSMinimizer::MinGradient
virtual const double * MinGradient() const
return pointer to gradient values at the minimum
Definition: GSLNLSMinimizer.cxx:375
ROOT::Math::MinimTransformFunction
MinimTransformFunction class to perform a transformations on the variables to deal with fixed or limi...
Definition: MinimTransformFunction.h:49
ROOT::Math::BasicFitMethodFunction::NCalls
virtual unsigned int NCalls() const
return the total number of function calls (overrided if needed)
Definition: FitMethodFunction.h:100
x
Double_t x[n]
Definition: legend1.C:17
ROOT::Math::IGradientFunctionMultiDimTempl
Interface (abstract class) for multi-dimensional functions providing a gradient calculation.
Definition: IFunction.h:326
FitMethodFunction.h
MATH_ERROR_MSG
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:83
ROOT::Math::GSLNLSMinimizer::fResiduals
std::vector< LSResidualFunc > fResiduals
Definition: GSLNLSMinimizer.h:278
ROOT::Math::GSLNLSMinimizer::~GSLNLSMinimizer
~GSLNLSMinimizer()
Destructor (no operations)
Definition: GSLNLSMinimizer.cxx:171
ROOT::Math::GSLNLSMinimizer::CovMatrixStatus
virtual int CovMatrixStatus() const
return covariance matrix status
Definition: GSLNLSMinimizer.cxx:389
ROOT::Math::GSLNLSMinimizer::fEdm
double fEdm
Definition: GSLNLSMinimizer.h:274
ROOT::Math::BasicMinimizer::MinValue
virtual double MinValue() const
return minimum function value
Definition: BasicMinimizer.h:146
ROOT::Math::GSLNLSMinimizer::fGSLMultiFit
ROOT::Math::GSLMultiFit * fGSLMultiFit
Definition: GSLNLSMinimizer.h:271
MinimTransformFunction.h
ROOT::Math::IMultiGenFunction
IMultiGenFunctionTempl< double > IMultiGenFunction
Definition: IFunctionfwd.h:42
ROOT::Math::GSLMultiFit::TestGradient
int TestGradient(double absTol) const
test gradient (ask from solver gradient vector)
Definition: GSLMultiFit.h:253
Error.h
ROOT::Math::BasicMinimizer::SetFinalValues
void SetFinalValues(const double *x)
Definition: BasicMinimizer.cxx:332
ROOT::Math::GSLNLSMinimizer::SetFunction
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)
set the function to minimize
Definition: GSLNLSMinimizer.cxx:178
ROOT::Math::LSResidualFunc
LSResidualFunc class description.
Definition: GSLNLSMinimizer.h:101
MATH_ERROR_MSGVAL
#define MATH_ERROR_MSGVAL(loc, txt, x)
Definition: Error.h:109
ROOT::Math::GSLNLSMinimizer::fNFree
unsigned int fNFree
Definition: GSLNLSMinimizer.h:268
ROOT::Math::Minimizer::Tolerance
double Tolerance() const
absolute tolerance
Definition: Minimizer.h:421
ROOT::Math::Minimizer::fStatus
int fStatus
Definition: Minimizer.h:494
double
double
Definition: Converters.cxx:921
sqrt
double sqrt(double)
ROOT::Math::GSLMultiFit::Set
int Set(const std::vector< Func > &funcVec, const double *x)
set the solver parameters
Definition: GSLMultiFit.h:181
ROOT::Math::BasicMinimizer::X
virtual const double * X() const
return pointer to X values at the minimum
Definition: BasicMinimizer.h:149
ROOT::Math::BasicFitMethodFunction::NPoints
virtual unsigned int NPoints() const
return the number of data points used in evaluating the function
Definition: FitMethodFunction.h:90
ROOT::Math::GSLMultiFit::Edm
double Edm() const
Definition: GSLMultiFit.h:267
ROOT::Math::MinimizerOptions::DefaultMaxIterations
static int DefaultMaxIterations()
Definition: MinimizerOptions.cxx:97
ROOT::Math::GSLNLSMinimizer::fCovMatrix
std::vector< double > fCovMatrix
Definition: GSLNLSMinimizer.h:277
ROOT::Math::MinimizerOptions::DefaultTolerance
static double DefaultTolerance()
Definition: MinimizerOptions.cxx:94
ROOT::Math::GSLMultiFit::CovarMatrix
const double * CovarMatrix() const
return covariance matrix of the parameters
Definition: GSLMultiFit.h:239
ROOT::Math::BasicMinimizer::NDim
virtual unsigned int NDim() const
number of dimensions
Definition: BasicMinimizer.h:152
ROOT::Math::MinimizerOptions::DefaultPrintLevel
static int DefaultPrintLevel()
Definition: MinimizerOptions.cxx:99
ROOT::Math::GSLMultiFit::Name
std::string Name() const
Definition: GSLMultiFit.h:210
ROOT::Math::MultiNumGradFunction
MultiNumGradFunction class to wrap a normal function in a gradient function using numerical gradient ...
Definition: MultiNumGradFunction.h:87
ROOT::Math::BasicFitMethodFunction::DataElement
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...
fSize
size_t fSize
Definition: DeclareConverters.h:342
ROOT::Math::BasicMinimizer::CreateTransformation
MinimTransformFunction * CreateTransformation(std::vector< double > &startValues, const ROOT::Math::IMultiGradFunction *func=0)
Definition: BasicMinimizer.cxx:282
ROOT::Math::GSLNLSMinimizer::fChi2Func
const ROOT::Math::FitMethodFunction * fChi2Func
Definition: GSLNLSMinimizer.h:272
ROOT::Math::BasicFitMethodFunction
FitMethodFunction class Interface for objective functions (like chi2 and likelihood used in the fit) ...
Definition: Fitter.h:40
type
int type
Definition: TGX11.cxx:121
ROOT::Math::IBaseFunctionMultiDimTempl
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
MultiNumGradFunction.h
ROOT::Math::GSLMultiFit::Gradient
const double * Gradient() const
gradient value at the minimum
Definition: GSLMultiFit.h:228
ROOT::Math::FitMethodFunction
BasicFitMethodFunction< ROOT::Math::IMultiGenFunction > FitMethodFunction
Definition: Fitter.h:40
ROOT::Math::GSLNLSMinimizer::fErrors
std::vector< double > fErrors
Definition: GSLNLSMinimizer.h:276
ROOT
VSD Structures.
Definition: StringConv.hxx:21
ROOT::Math::GSLNLSMinimizer::CovMatrix
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...
Definition: GSLNLSMinimizer.cxx:381
Math
ROOT::Math::BasicMinimizer::ObjFunction
const ROOT::Math::IMultiGenFunction * ObjFunction() const
return pointer to used objective function
Definition: BasicMinimizer.h:161
ROOT::Math::GSLNLSMinimizer::fSize
unsigned int fSize
Definition: GSLNLSMinimizer.h:269
g
#define g(i)
Definition: RSha256.hxx:123