Logo ROOT  
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 
52  BasicMinimizer()
53 {
54  // Constructor implementation : create GSLMultiMin wrapper object
55  //std::cout << "create GSL Minimizer of type " << type << std::endl;
56 
57  fGSLMultiMin = new GSLMultiMinimizer((ROOT::Math::EGSLMinimizerType) type);
58 
59  fLSTolerance = 0.1; // line search tolerance (use fixed)
61  if (niter <=0 ) niter = 1000;
62  SetMaxIterations(niter);
64 }
65 
66 GSLMinimizer::GSLMinimizer( const char * type) : BasicMinimizer()
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
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 
123 bool GSLMinimizer::Minimize() {
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 
ROOT::Math::BasicMinimizer::SetMinValue
void SetMinValue(double val)
Definition: BasicMinimizer.h:185
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
ROOT::Math::GSLMinimizer::SetFunction
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)
set the function to minimize
Definition: GSLMinimizer.cxx:138
ROOT::Math::kConjugateFR
@ kConjugateFR
Definition: GSLMinimizer.h:107
ROOT::Math::GSLMinimizer::fLSTolerance
double fLSTolerance
Definition: GSLMinimizer.h:199
ROOT::Math::GSLMinimizer::NCalls
virtual unsigned int NCalls() const
number of function calls to reach the minimum
Definition: GSLMinimizer.cxx:150
ROOT::Math::Minimizer::PrintLevel
int PrintLevel() const
minimizer configuration parameters
Definition: Minimizer.h:412
ROOT::Math::GSLMinimizer::GSLMinimizer
GSLMinimizer(ROOT::Math::EGSLMinimizerType type=ROOT::Math::kConjugateFR)
Default constructor.
Definition: GSLMinimizer.cxx:89
GSLMultiMinimizer.h
ROOT::Math::BasicMinimizer::NPar
virtual unsigned int NPar() const
total number of parameter defined
Definition: BasicMinimizer.h:158
ROOT::Math::EGSLMinimizerType
EGSLMinimizerType
enumeration specifying the types of GSL minimizers
Definition: GSLMinimizer.h:87
ROOT::Math::BasicMinimizer::SetFunction
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)
set the function to minimize
Definition: BasicMinimizer.cxx:250
ROOT::Math::kConjugatePR
@ kConjugatePR
Definition: GSLMinimizer.h:108
ROOT::Math::Minimizer::MaxIterations
unsigned int MaxIterations() const
max iterations
Definition: Minimizer.h:418
GSL_SUCCESS
#define GSL_SUCCESS
Definition: RooAdaptiveGaussKronrodIntegrator1D.cxx:400
ROOT::Math::BasicMinimizer::NFree
virtual unsigned int NFree() const
number of free variables (real dimension of the problem)
Definition: BasicMinimizer.cxx:366
ROOT::Math::GSLMinimizer::~GSLMinimizer
virtual ~GSLMinimizer()
Destructor.
Definition: GSLMinimizer.cxx:131
ROOT::Math::BasicMinimizer::GradObjFunction
const ROOT::Math::IMultiGradFunction * GradObjFunction() const
return pointer to used gradient object function (NULL if gradient is not supported)
Definition: BasicMinimizer.cxx:358
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
ROOT::Math::MultiNumGradFunction::NCalls
unsigned int NCalls() const
Definition: MultiNumGradFunction.h:148
ROOT::Math::GSLMultiMinimizer::TestGradient
int TestGradient(double absTol) const
test gradient (ask from minimizer gradient vector)
Definition: GSLMultiMinimizer.h:241
FitMethodFunction.h
MATH_ERROR_MSG
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:83
ROOT::Math::GSLMultiMinimizer::Name
std::string Name() const
Definition: GSLMultiMinimizer.h:204
ROOT::Math::BasicMinimizer::MinValue
virtual double MinValue() const
return minimum function value
Definition: BasicMinimizer.h:146
MinimTransformFunction.h
ROOT::Math::GSLMinimizer::MinGradient
virtual const double * MinGradient() const
return pointer to gradient values at the minimum
Definition: GSLMinimizer.cxx:326
ROOT::Math::GSLMultiMinimizer::Gradient
double * Gradient() const
gradient value at the minimum
Definition: GSLMultiMinimizer.h:228
ROOT::Math::BasicMinimizer::SetFinalValues
void SetFinalValues(const double *x)
Definition: BasicMinimizer.cxx:332
epsilon
REAL epsilon
Definition: triangle.c:617
ROOT::Math::kVectorBFGS2
@ kVectorBFGS2
Definition: GSLMinimizer.h:110
MATH_ERROR_MSGVAL
#define MATH_ERROR_MSGVAL(loc, txt, x)
Definition: Error.h:109
ROOT::Math::GSLMultiMinimizer::Set
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...
Definition: GSLMultiMinimizer.h:184
ROOT::Math::Minimizer::Tolerance
double Tolerance() const
absolute tolerance
Definition: Minimizer.h:421
ROOT::Math::Minimizer::fStatus
int fStatus
Definition: Minimizer.h:494
ROOT::Math::BasicMinimizer::StepSizes
virtual const double * StepSizes() const
accessor methods
Definition: BasicMinimizer.h:173
sqrt
double sqrt(double)
ROOT::Math::BasicMinimizer::X
virtual const double * X() const
return pointer to X values at the minimum
Definition: BasicMinimizer.h:149
ROOT::Math::GSLMinimizer::fGSLMultiMin
ROOT::Math::GSLMultiMinimizer * fGSLMultiMin
Definition: GSLMinimizer.h:197
MATH_INFO_MSG
#define MATH_INFO_MSG(loc, str)
Pre-processor macro to report messages which can be configured to use ROOT error or simply an std::io...
Definition: Error.h:77
ROOT::Math::GSLMultiMinimizer::Iterate
int Iterate()
Definition: GSLMultiMinimizer.h:209
ROOT::Math::MinimizerOptions::DefaultMaxIterations
static int DefaultMaxIterations()
Definition: MinimizerOptions.cxx:97
ROOT::Math::GSLMinimizer::Minimize
virtual bool Minimize()
method to perform the minimization
Definition: GSLMinimizer.cxx:161
ROOT::Math::BasicMinimizer::NDim
virtual unsigned int NDim() const
number of dimensions
Definition: BasicMinimizer.h:152
ROOT::Math::GSLMultiMinimizer::Minimum
double Minimum() const
function value at the minimum
Definition: GSLMultiMinimizer.h:222
ROOT::Math::kSteepestDescent
@ kSteepestDescent
Definition: GSLMinimizer.h:111
ROOT::Math::MinimizerOptions::DefaultPrintLevel
static int DefaultPrintLevel()
Definition: MinimizerOptions.cxx:99
ROOT::Math::MultiNumGradFunction
MultiNumGradFunction class to wrap a normal function in a gradient function using numerical gradient ...
Definition: MultiNumGradFunction.h:87
ROOT::Math::BasicMinimizer::CreateTransformation
MinimTransformFunction * CreateTransformation(std::vector< double > &startValues, const ROOT::Math::IMultiGradFunction *func=0)
Definition: BasicMinimizer.cxx:282
MATH_WARN_MSGVAL
#define MATH_WARN_MSGVAL(loc, txt, x)
Definition: Error.h:105
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::GSLMultiMinimizer::X
double * X() const
x values at the minimum
Definition: GSLMultiMinimizer.h:215
ROOT
VSD Structures.
Definition: StringConv.hxx:21
Math
ROOT::Math::BasicMinimizer::ObjFunction
const ROOT::Math::IMultiGenFunction * ObjFunction() const
return pointer to used objective function
Definition: BasicMinimizer.h:161
ROOT::Math::kVectorBFGS
@ kVectorBFGS
Definition: GSLMinimizer.h:109
GSLMinimizer.h
g
#define g(i)
Definition: RSha256.hxx:123