Logo ROOT   6.10/09
Reference Guide
GSLMinimizer.cxx
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Author: L. Moneta Tue Dec 19 15:41:39 2006
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
7  * *
8  * This library is free software; you can redistribute it and/or *
9  * modify it under the terms of the GNU General Public License *
10  * as published by the Free Software Foundation; either version 2 *
11  * of the License, or (at your option) any later version. *
12  * *
13  * This library is distributed in the hope that it will be useful, *
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
16  * General Public License for more details. *
17  * *
18  * You should have received a copy of the GNU General Public License *
19  * along with this library (see file COPYING); if not, write *
20  * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
21  * 330, Boston, MA 02111-1307 USA, or contact the author. *
22  * *
23  **********************************************************************/
24 
25 // Implementation file for class GSLMinimizer
26 
27 #include "Math/GSLMinimizer.h"
28 
29 #include "GSLMultiMinimizer.h"
30 
32 #include "Math/FitMethodFunction.h"
33 
35 
36 #include <cassert>
37 
38 #include <iostream>
39 #include <iomanip>
40 #include <cmath>
41 #include <algorithm>
42 #include <functional>
43 #include <ctype.h> // need to use c version of tolower defined here
44 #include <limits>
45 
46 namespace ROOT {
47 
48  namespace Math {
49 
50 
53 {
54  // Constructor implementation : create GSLMultiMin wrapper object
55  //std::cout << "create GSL Minimizer of type " << type << std::endl;
56 
58 
59  fLSTolerance = 0.1; // line search tolerance (use fixed)
61  if (niter <=0 ) niter = 1000;
62  SetMaxIterations(niter);
64 }
65 
67 {
68  // Constructor implementation from a string
69  std::string algoname(type);
70  std::transform(algoname.begin(), algoname.end(), algoname.begin(), (int(*)(int)) tolower );
71 
72  ROOT::Math::EGSLMinimizerType algo = kVectorBFGS2; // default value
73 
74  if (algoname == "conjugatefr") algo = kConjugateFR;
75  if (algoname == "conjugatepr") algo = kConjugatePR;
76  if (algoname == "bfgs") algo = kVectorBFGS;
77  if (algoname == "bfgs2") algo = kVectorBFGS2;
78  if (algoname == "steepestdescent") algo = kSteepestDescent;
79 
80 
81  //std::cout << "create GSL Minimizer of type " << algo << std::endl;
82 
83  fGSLMultiMin = new GSLMultiMinimizer(algo);
84 
85  fLSTolerance = 0.1; // use 10**-4
87  if (niter <=0 ) niter = 1000;
88  SetMaxIterations(niter);
90 }
91 
92 
94  assert(fGSLMultiMin != 0);
95  delete fGSLMultiMin;
96 }
97 
98 
99 
101  // set the function to minimizer
102  // need to calculate numerically the derivatives: do via class MultiNumGradFunction
103  // no need to clone the passed function
104  ROOT::Math::MultiNumGradFunction gradFunc(func);
105  // function is cloned inside so can be delete afterwards
106  // called base class method setfunction
107  // (note: write explicitly otherwise it will call back itself)
108  BasicMinimizer::SetFunction(gradFunc);
109 }
110 
111 
112 unsigned int GSLMinimizer::NCalls() const {
113  // return numbr of function calls
114  // if method support
115  const ROOT::Math::MultiNumGradFunction * fnumgrad = dynamic_cast<const ROOT::Math::MultiNumGradFunction *>(ObjFunction());
116  if (fnumgrad) return fnumgrad->NCalls();
117  const ROOT::Math::FitMethodGradFunction * ffitmethod = dynamic_cast<const ROOT::Math::FitMethodGradFunction *>(ObjFunction());
118  if (ffitmethod) return ffitmethod->NCalls();
119  // not supported in the other case
120  return 0;
121 }
122 
124  // set initial parameters of the minimizer
125 
126  if (fGSLMultiMin == 0) return false;
127  const ROOT::Math::IMultiGradFunction * function = GradObjFunction();
128  if (function == 0) {
129  MATH_ERROR_MSG("GSLMinimizer::Minimize","Function has not been set");
130  return false;
131  }
132 
133  unsigned int npar = NPar();
134  unsigned int ndim = NDim();
135  if (npar == 0 || npar < NDim() ) {
136  MATH_ERROR_MSGVAL("GSLMinimizer::Minimize","Wrong number of parameters",npar);
137  return false;
138  }
139  if (npar > ndim ) {
140  MATH_WARN_MSGVAL("GSLMinimizer::Minimize","number of parameters larger than function dimension - ignore extra parameters",npar);
141  }
142 
143  const double eps = std::numeric_limits<double>::epsilon();
144 
145  std::vector<double> startValues;
146  std::vector<double> steps(StepSizes(), StepSizes()+npar);
147 
148  MinimTransformFunction * trFunc = CreateTransformation(startValues);
149  if (trFunc) {
150  function = trFunc;
151  // need to transform also the steps
152  trFunc->InvStepTransformation(X(), StepSizes(), &steps[0]);
153  steps.resize(trFunc->NDim());
154  }
155 
156  // in case all parameters are free - just evaluate the function
157  if (NFree() == 0) {
158  MATH_INFO_MSG("GSLMinimizer::Minimize","There are no free parameter - just compute the function value");
159  double fval = (*function)((double*)0); // no need to pass parameters
160  SetFinalValues(&startValues[0]);
161  SetMinValue(fval);
162  fStatus = 0;
163  return true;
164  }
165 
166  // use a global step size = modules of step vectors
167  double stepSize = 0;
168  for (unsigned int i = 0; i < steps.size(); ++i)
169  stepSize += steps[i]*steps[i];
170  stepSize = std::sqrt(stepSize);
171  if (stepSize < eps) {
172  MATH_ERROR_MSGVAL("GSLMinimizer::Minimize","Step size is too small",stepSize);
173  return false;
174  }
175 
176 
177  // set parameters in internal GSL minimization class
178  fGSLMultiMin->Set(*function, &startValues.front(), stepSize, fLSTolerance );
179 
180 
181  int debugLevel = PrintLevel();
182 
183  if (debugLevel >=1 ) std::cout <<"Minimize using GSLMinimizer " << fGSLMultiMin->Name() << std::endl;
184 
185 
186  //std::cout <<"print Level " << debugLevel << std::endl;
187  //debugLevel = 3;
188 
189  // start iteration
190  unsigned int iter = 0;
191  int status;
192  bool minFound = false;
193  bool iterFailed = false;
194  do {
195  status = fGSLMultiMin->Iterate();
196  if (status) {
197  iterFailed = true;
198  break;
199  }
200 
201  status = fGSLMultiMin->TestGradient( Tolerance() );
202  if (status == GSL_SUCCESS) {
203  minFound = true;
204  }
205 
206  if (debugLevel >=2) {
207  std::cout << "----------> Iteration " << std::setw(4) << iter;
208  int pr = std::cout.precision(18);
209  std::cout << " FVAL = " << fGSLMultiMin->Minimum() << std::endl;
210  std::cout.precision(pr);
211  if (debugLevel >=3) {
212  std::cout << " Parameter Values : ";
213  const double * xtmp = fGSLMultiMin->X();
214  std::cout << std::endl;
215  if (trFunc != 0 ) {
216  xtmp = trFunc->Transformation(xtmp);
217  }
218  for (unsigned int i = 0; i < NDim(); ++i) {
219  std::cout << " " << VariableName(i) << " = " << xtmp[i];
220  // avoid nan
221  // if (std::isnan(xtmp[i])) status = -11;
222  }
223  std::cout << std::endl;
224  }
225  }
226 
227 
228  iter++;
229 
230 
231  }
232  while (status == GSL_CONTINUE && iter < MaxIterations() );
233 
234 
235  // save state with values and function value
236  double * x = fGSLMultiMin->X();
237  if (x == 0) return false;
238  SetFinalValues(x);
239 
240  double minVal = fGSLMultiMin->Minimum();
241  SetMinValue(minVal);
242 
243  fStatus = status;
244 
245 
246  if (minFound) {
247  if (debugLevel >=1 ) {
248  std::cout << "GSLMinimizer: Minimum Found" << std::endl;
249  int pr = std::cout.precision(18);
250  std::cout << "FVAL = " << MinValue() << std::endl;
251  std::cout.precision(pr);
252 // std::cout << "Edm = " << fState.Edm() << std::endl;
253  std::cout << "Niterations = " << iter << std::endl;
254  unsigned int ncalls = NCalls();
255  if (ncalls) std::cout << "NCalls = " << ncalls << std::endl;
256  for (unsigned int i = 0; i < NDim(); ++i)
257  std::cout << VariableName(i) << "\t = " << X()[i] << std::endl;
258  }
259  return true;
260  }
261  else {
262  if (debugLevel >= -1 ) {
263  std::cout << "GSLMinimizer: Minimization did not converge" << std::endl;
264  if (iterFailed) {
265  if (status == GSL_ENOPROG) // case status 27
266  std::cout << "\t Iteration is not making progress towards solution" << std::endl;
267  else
268  std::cout << "\t Iteration failed with status " << status << std::endl;
269 
270  if (debugLevel >= 1) {
271  double * g = fGSLMultiMin->Gradient();
272  double dg2 = 0;
273  for (unsigned int i = 0; i < NDim(); ++i) dg2 += g[i] * g[1];
274  std::cout << "Grad module is " << std::sqrt(dg2) << std::endl;
275  for (unsigned int i = 0; i < NDim(); ++i)
276  std::cout << VariableName(i) << "\t = " << X()[i] << std::endl;
277  std::cout << "FVAL = " << MinValue() << std::endl;
278 // std::cout << "Edm = " << fState.Edm() << std::endl;
279  std::cout << "Niterations = " << iter << std::endl;
280  }
281  }
282  }
283  return false;
284  }
285  return false;
286 }
287 
288 const double * GSLMinimizer::MinGradient() const {
289  // return gradient (internal values)
290  return fGSLMultiMin->Gradient();
291 }
292 
293 
294  } // end namespace Math
295 
296 } // end namespace ROOT
297 
Interface (abstract class) for multi-dimensional functions providing a gradient calculation.
Definition: IFunction.h:330
ROOT::Math::GSLMultiMinimizer * fGSLMultiMin
Definition: GSLMinimizer.h:166
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
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
double * Gradient() const
gradient value at the minimum
const double * Transformation(const double *x) const
transform from internal to external result is cached also inside the class
MinimTransformFunction class to perform a transformations on the variables to deal with fixed or limi...
EGSLMinimizerType
enumeration specifying the types of GSL minimizers
Definition: GSLMinimizer.h:56
Base Minimizer class, which defines the basic funcionality of various minimizer implementations (apar...
MinimTransformFunction * CreateTransformation(std::vector< double > &startValues, const ROOT::Math::IMultiGradFunction *func=0)
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)
int PrintLevel() const
minimizer configuration parameters
Definition: Minimizer.h:411
double * X() const
x values at the minimum
virtual const double * StepSizes() const
accessor methods
#define MATH_ERROR_MSGVAL(loc, str, x)
Definition: Error.h:69
virtual bool Minimize()
method to perform the minimization
double sqrt(double)
Double_t x[n]
Definition: legend1.C:17
void SetMinValue(double val)
virtual unsigned int NPar() const
total number of parameter defined
unsigned int MaxIterations() const
max iterations
Definition: Minimizer.h:417
unsigned int NDim() const
Retrieve the dimension of the function.
virtual unsigned int NDim() const
number of dimensions
const ROOT::Math::IMultiGenFunction * ObjFunction() const
return pointer to used objective function
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)
set the function to minimize
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:50
double Tolerance() const
absolute tolerance
Definition: Minimizer.h:420
GSLMinimizer(ROOT::Math::EGSLMinimizerType type=ROOT::Math::kConjugateFR)
Default constructor.
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
#define MATH_INFO_MSG(loc, str)
Definition: Error.h:44
void SetFinalValues(const double *x)
virtual ~GSLMinimizer()
Destructor.
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)
set the function to minimize
int Set(const ROOT::Math::IMultiGradFunction &func, const double *x, double stepSize, double tol)
set the function to be minimize the initial minimizer parameters, step size and tolerance in the line...
int TestGradient(double absTol) const
test gradient (ask from minimizer gradient vector)
REAL epsilon
Definition: triangle.c:617
MultiNumGradFunction class to wrap a normal function in a gradient function using numerical gradient ...
const ROOT::Math::IMultiGradFunction * GradObjFunction() const
return pointer to used gradient object function (NULL if gradient is not supported) ...
virtual unsigned int NCalls() const
number of function calls to reach the minimum
virtual unsigned int NFree() const
number of free variables (real dimension of the problem)
FitMethodFunction class Interface for objective functions (like chi2 and likelihood used in the fit) ...
Definition: Fitter.h:40
void InvStepTransformation(const double *x, const double *sext, double *sint) const
inverse transformation for steps (external -> internal) at external point x
int type
Definition: TGX11.cxx:120
double func(double *x, double *p)
Definition: stressTF1.cxx:213
#define MATH_WARN_MSGVAL(loc, str, x)
Definition: Error.h:67
GSLMultiMinimizer class , for minimizing multi-dimensional function using derivatives.
Namespace for new Math classes and functions.
virtual unsigned int NCalls() const
return the total number of function calls (overrided if needed)
double Minimum() const
function value at the minimum
void SetPrintLevel(int level)
set print level
Definition: Minimizer.h:445
virtual const double * MinGradient() const
return pointer to gradient values at the minimum