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
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
46namespace 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
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)
109}
110
111
112unsigned 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;
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;
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
288const 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
#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:76
#define MATH_ERROR_MSGVAL(loc, txt, x)
Definition: Error.h:108
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:82
#define MATH_WARN_MSGVAL(loc, txt, x)
Definition: Error.h:104
#define g(i)
Definition: RSha256.hxx:105
int type
Definition: TGX11.cxx:120
double sqrt(double)
FitMethodFunction class Interface for objective functions (like chi2 and likelihood used in the fit) ...
virtual unsigned int NCalls() const
return the total number of function calls (overrided if needed)
Base Minimizer class, which defines the basic funcionality of various minimizer implementations (apar...
virtual double MinValue() const
return minimum function value
virtual unsigned int NPar() const
total number of parameter defined
virtual unsigned int NFree() const
number of free variables (real dimension of the problem)
void SetMinValue(double val)
virtual const double * StepSizes() const
accessor methods
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)
set the function to minimize
const ROOT::Math::IMultiGenFunction * ObjFunction() const
return pointer to used objective function
virtual unsigned int NDim() const
number of dimensions
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
void SetFinalValues(const double *x)
MinimTransformFunction * CreateTransformation(std::vector< double > &startValues, const ROOT::Math::IMultiGradFunction *func=0)
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 const double * MinGradient() const
return pointer to gradient values at the minimum
ROOT::Math::GSLMultiMinimizer * fGSLMultiMin
Definition: GSLMinimizer.h:166
GSLMinimizer(ROOT::Math::EGSLMinimizerType type=ROOT::Math::kConjugateFR)
Default constructor.
virtual bool Minimize()
method to perform the minimization
virtual ~GSLMinimizer()
Destructor.
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)
set the function to minimize
GSLMultiMinimizer class , for minimizing multi-dimensional function using derivatives.
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...
double Minimum() const
function value at the minimum
double * Gradient() const
gradient value at the minimum
double * X() const
x values at the minimum
int TestGradient(double absTol) const
test gradient (ask from minimizer gradient vector)
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
Interface (abstract class) for multi-dimensional functions providing a gradient calculation.
Definition: IFunction.h:327
MinimTransformFunction class to perform a transformations on the variables to deal with fixed or limi...
unsigned int NDim() const
Retrieve the dimension of the function.
const double * Transformation(const double *x) const
transform from internal to external result is cached also inside the class
void InvStepTransformation(const double *x, const double *sext, double *sint) const
inverse transformation for steps (external -> internal) at external point x
double Tolerance() const
absolute tolerance
Definition: Minimizer.h:420
void SetMaxIterations(unsigned int maxiter)
set maximum iterations (one iteration can have many function calls)
Definition: Minimizer.h:451
unsigned int MaxIterations() const
max iterations
Definition: Minimizer.h:417
void SetPrintLevel(int level)
set print level
Definition: Minimizer.h:445
int PrintLevel() const
minimizer configuration parameters
Definition: Minimizer.h:411
MultiNumGradFunction class to wrap a normal function in a gradient function using numerical gradient ...
EGSLMinimizerType
enumeration specifying the types of GSL minimizers
Definition: GSLMinimizer.h:56
@ kSteepestDescent
Definition: GSLMinimizer.h:61
Double_t x[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:151
VSD Structures.
Definition: StringConv.hxx:21
REAL epsilon
Definition: triangle.c:617