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 number of function calls
114 // if original function does not support gradient it is wrapped in MultiNumGradFuction
115 // and we have NCalls available
116 const ROOT::Math::MultiNumGradFunction * fnumgrad = dynamic_cast<const ROOT::Math::MultiNumGradFunction *>(ObjFunction());
117 if (fnumgrad) return fnumgrad->NCalls();
118 // if original function implement gradient, we can get NumCalls a=only if it is a
119 // FitMethodGradFunction
120 const ROOT::Math::FitMethodGradFunction * ffitmethod = dynamic_cast<const ROOT::Math::FitMethodGradFunction *>(ObjFunction());
121 if (ffitmethod) return ffitmethod->NCalls();
122 // not supported in the other case
123 return 0;
124}
125
127 // set initial parameters of the minimizer
128
129 if (fGSLMultiMin == 0) return false;
131 if (function == 0) {
132 MATH_ERROR_MSG("GSLMinimizer::Minimize","Function has not been set");
133 return false;
134 }
135
136 unsigned int npar = NPar();
137 unsigned int ndim = NDim();
138 if (npar == 0 || npar < NDim() ) {
139 MATH_ERROR_MSGVAL("GSLMinimizer::Minimize","Wrong number of parameters",npar);
140 return false;
141 }
142 if (npar > ndim ) {
143 MATH_WARN_MSGVAL("GSLMinimizer::Minimize","number of parameters larger than function dimension - ignore extra parameters",npar);
144 }
145
146 const double eps = std::numeric_limits<double>::epsilon();
147
148 std::vector<double> startValues;
149 std::vector<double> steps(StepSizes(), StepSizes()+npar);
150
151 std::unique_ptr<MinimTransformFunction> trFunc( CreateTransformation(startValues));
152 if (trFunc) {
153 // need to transform also the steps
154 trFunc->InvStepTransformation(X(), StepSizes(), &steps[0]);
155 steps.resize(trFunc->NDim());
156 }
157
158 // in case all parameters are free - just evaluate the function
159 if (NFree() == 0) {
160 MATH_INFO_MSG("GSLMinimizer::Minimize","There are no free parameter - just compute the function value");
161 double fval = (*function)((double*)0); // no need to pass parameters or used transformed function
162 SetFinalValues(&startValues[0]);
163 SetMinValue(fval);
164 fStatus = 0;
165 return true;
166 }
167
168 // use a global step size = modules of step vectors
169 double stepSize = 0;
170 for (unsigned int i = 0; i < steps.size(); ++i)
171 stepSize += steps[i]*steps[i];
172 stepSize = std::sqrt(stepSize);
173 if (stepSize < eps) {
174 MATH_ERROR_MSGVAL("GSLMinimizer::Minimize","Step size is too small",stepSize);
175 return false;
176 }
177
178 // set parameters in internal GSL minimization class
179 fGSLMultiMin->Set( (trFunc) ? *trFunc : *function, &startValues.front(), stepSize, fLSTolerance );
180
181
182 int debugLevel = PrintLevel();
183
184 if (debugLevel >=1 ) std::cout <<"Minimize using GSLMinimizer " << fGSLMultiMin->Name() << std::endl;
185
186
187 //std::cout <<"print Level " << debugLevel << std::endl;
188 //debugLevel = 3;
189
190 // start iteration
191 unsigned int iter = 0;
192 int status;
193 bool minFound = false;
194 bool iterFailed = false;
195 do {
196 status = fGSLMultiMin->Iterate();
197 if (status) {
198 iterFailed = true;
199 break;
200 }
201
202 status = fGSLMultiMin->TestGradient( Tolerance() );
203 if (status == GSL_SUCCESS) {
204 minFound = true;
205 }
206
207 if (debugLevel >=2) {
208 std::cout << "----------> Iteration " << std::setw(4) << iter;
209 int pr = std::cout.precision(18);
210 std::cout << " FVAL = " << fGSLMultiMin->Minimum() << std::endl;
211 std::cout.precision(pr);
212 if (debugLevel >=3) {
213 std::cout << " Parameter Values : ";
214 const double * xtmp = fGSLMultiMin->X();
215 std::cout << std::endl;
216 if (trFunc) {
217 xtmp = trFunc->Transformation(xtmp);
218 }
219 for (unsigned int i = 0; i < NDim(); ++i) {
220 std::cout << " " << VariableName(i) << " = " << xtmp[i];
221 // avoid nan
222 // if (std::isnan(xtmp[i])) status = -11;
223 }
224 std::cout << std::endl;
225 }
226 }
227
228
229 iter++;
230
231
232 }
233 while (status == GSL_CONTINUE && iter < MaxIterations() );
234
235
236 // save state with values and function value
237 double * x = fGSLMultiMin->X();
238 if (x == 0) return false;
239 SetFinalValues(x, trFunc.get());
240
241 double minVal = fGSLMultiMin->Minimum();
242 SetMinValue(minVal);
243
244 fStatus = status;
245
246
247 if (minFound) {
248 if (debugLevel >=1 ) {
249 std::cout << "GSLMinimizer: Minimum Found" << std::endl;
250 int pr = std::cout.precision(18);
251 std::cout << "FVAL = " << MinValue() << std::endl;
252 std::cout.precision(pr);
253// std::cout << "Edm = " << fState.Edm() << std::endl;
254 std::cout << "Niterations = " << iter << std::endl;
255 unsigned int ncalls = NCalls();
256 if (ncalls) std::cout << "NCalls = " << ncalls << std::endl;
257 for (unsigned int i = 0; i < NDim(); ++i)
258 std::cout << VariableName(i) << "\t = " << X()[i] << std::endl;
259 }
260 return true;
261 }
262 else {
263 if (debugLevel >= -1 ) {
264 std::cout << "GSLMinimizer: Minimization did not converge" << std::endl;
265 if (iterFailed) {
266 if (status == GSL_ENOPROG) // case status 27
267 std::cout << "\t Iteration is not making progress towards solution" << std::endl;
268 else
269 std::cout << "\t Iteration failed with status " << status << std::endl;
270
271 if (debugLevel >= 1) {
272 double * g = fGSLMultiMin->Gradient();
273 double dg2 = 0;
274 for (unsigned int i = 0; i < NDim(); ++i) dg2 += g[i] * g[1];
275 std::cout << "Grad module is " << std::sqrt(dg2) << std::endl;
276 for (unsigned int i = 0; i < NDim(); ++i)
277 std::cout << VariableName(i) << "\t = " << X()[i] << std::endl;
278 std::cout << "FVAL = " << MinValue() << 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:77
#define MATH_ERROR_MSGVAL(loc, txt, x)
Definition: Error.h:109
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:83
#define MATH_WARN_MSGVAL(loc, txt, x)
Definition: Error.h:105
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t g
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
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 (override if needed)
Base Minimizer class, which defines the basic functionality of various minimizer implementations (apa...
virtual unsigned int NPar() const
total number of parameter defined
unsigned int NFree() const override
number of free variables (real dimension of the problem)
unsigned int NDim() const override
number of dimensions
void SetMinValue(double val)
void SetFinalValues(const double *x, const MinimTransformFunction *func=nullptr)
double MinValue() const override
return minimum function value
virtual const double * StepSizes() const
accessor methods
void SetFunction(const ROOT::Math::IMultiGenFunction &func) override
set the function to minimize
const ROOT::Math::IMultiGenFunction * ObjFunction() const
return pointer to used objective function
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)
const double * X() const override
return pointer to X values at the minimum
std::string VariableName(unsigned int ivar) const override
get name of variables (override if minimizer support storing of variable names)
const double * MinGradient() const override
return pointer to gradient values at the minimum
ROOT::Math::GSLMultiMinimizer * fGSLMultiMin
Definition: GSLMinimizer.h:159
void SetFunction(const ROOT::Math::IMultiGenFunction &func) override
set the function to minimize
~GSLMinimizer() override
Destructor.
unsigned int NCalls() const override
number of function calls to reach the minimum
GSLMinimizer(ROOT::Math::EGSLMinimizerType type=ROOT::Math::kConjugateFR)
Default constructor.
bool Minimize() override
method to perform the minimization
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:343
double Tolerance() const
absolute tolerance
Definition: Minimizer.h:416
void SetMaxIterations(unsigned int maxiter)
set maximum iterations (one iteration can have many function calls)
Definition: Minimizer.h:450
int fStatus
status of minimizer
Definition: Minimizer.h:492
unsigned int MaxIterations() const
max iterations
Definition: Minimizer.h:413
void SetPrintLevel(int level)
set print level
Definition: Minimizer.h:444
int PrintLevel() const
minimizer configuration parameters
Definition: Minimizer.h:407
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:49
@ kSteepestDescent
Definition: GSLMinimizer.h:54
Double_t x[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:167
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
double epsilon
Definition: triangle.c:618