Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 <cctype> // need to use c version of tolower defined here
44#include <limits>
45
46namespace ROOT {
47
48 namespace Math {
49
50
52{
53 // Constructor implementation : create GSLMultiMin wrapper object
54 //std::cout << "create GSL Minimizer of type " << type << std::endl;
55
57
58 fLSTolerance = 0.1; // line search tolerance (use fixed)
60 if (niter <=0 ) niter = 1000;
61 SetMaxIterations(niter);
63}
64
66{
67 // Constructor implementation from a string
68 std::string algoname(type);
69 std::transform(algoname.begin(), algoname.end(), algoname.begin(), (int(*)(int)) tolower );
70
71 ROOT::Math::EGSLMinimizerType algo = kVectorBFGS2; // default value
72
73 if (algoname == "conjugatefr") algo = kConjugateFR;
74 if (algoname == "conjugatepr") algo = kConjugatePR;
75 if (algoname == "bfgs") algo = kVectorBFGS;
76 if (algoname == "bfgs2") algo = kVectorBFGS2;
77 if (algoname == "steepestdescent") algo = kSteepestDescent;
78
79
80 //std::cout << "create GSL Minimizer of type " << algo << std::endl;
81
83
84 fLSTolerance = 0.1; // use 10**-4
86 if (niter <=0 ) niter = 1000;
87 SetMaxIterations(niter);
89}
90
91
93 assert(fGSLMultiMin != nullptr);
94 delete fGSLMultiMin;
95}
96
97
98
100 // set the function to minimizer
101 // need to calculate numerically the derivatives: do via class MultiNumGradFunction
102 // no need to clone the passed function
104 // function is cloned inside so can be delete afterwards
105 // called base class method setfunction
106 // (note: write explicitly otherwise it will call back itself)
108}
109
110
111unsigned int GSLMinimizer::NCalls() const {
112 // return number of function calls
113 // if original function does not support gradient it is wrapped in MultiNumGradFunction
114 // and we have NCalls available
115 const ROOT::Math::MultiNumGradFunction * fnumgrad = dynamic_cast<const ROOT::Math::MultiNumGradFunction *>(ObjFunction());
116 if (fnumgrad) return fnumgrad->NCalls();
117 // if original function implement gradient, we can get NumCalls a=only if it is a
118 // FitMethodGradFunction
119 const ROOT::Math::FitMethodGradFunction * ffitmethod = dynamic_cast<const ROOT::Math::FitMethodGradFunction *>(ObjFunction());
120 if (ffitmethod) return ffitmethod->NCalls();
121 // not supported in the other case
122 return 0;
123}
124
126 // set initial parameters of the minimizer
127
128 if (fGSLMultiMin == nullptr) return false;
130 if (function == nullptr) {
131 MATH_ERROR_MSG("GSLMinimizer::Minimize","Function has not been set");
132 return false;
133 }
134
135 unsigned int npar = NPar();
136 unsigned int ndim = NDim();
137 if (npar == 0 || npar < NDim() ) {
138 MATH_ERROR_MSGVAL("GSLMinimizer::Minimize","Wrong number of parameters",npar);
139 return false;
140 }
141 if (npar > ndim ) {
142 MATH_WARN_MSGVAL("GSLMinimizer::Minimize","number of parameters larger than function dimension - ignore extra parameters",npar);
143 }
144
145 const double eps = std::numeric_limits<double>::epsilon();
146
147 std::vector<double> startValues;
148 std::vector<double> steps(StepSizes(), StepSizes()+npar);
149
150 std::unique_ptr<MinimTransformFunction> trFunc( CreateTransformation(startValues));
151 if (trFunc) {
152 // need to transform also the steps
153 trFunc->InvStepTransformation(X(), StepSizes(), &steps[0]);
154 steps.resize(trFunc->NDim());
155 }
156
157 // in case all parameters are free - just evaluate the function
158 if (NFree() == 0) {
159 MATH_INFO_MSG("GSLMinimizer::Minimize","There are no free parameter - just compute the function value");
160 double fval = (*function)((double*)nullptr); // no need to pass parameters or used transformed function
161 SetFinalValues(&startValues[0]);
162 SetMinValue(fval);
163 fStatus = 0;
164 return true;
165 }
166
167 // use a global step size = modules of step vectors
168 double stepSize = 0;
169 for (unsigned int i = 0; i < steps.size(); ++i)
170 stepSize += steps[i]*steps[i];
171 stepSize = std::sqrt(stepSize);
172 if (stepSize < eps) {
173 MATH_ERROR_MSGVAL("GSLMinimizer::Minimize","Step size is too small",stepSize);
174 return false;
175 }
176
177 // set parameters in internal GSL minimization class
178 fGSLMultiMin->Set( (trFunc) ? *trFunc : *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) {
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 == nullptr) return false;
238 SetFinalValues(x, trFunc.get());
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 > 0 ) {
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 << "Niterations = " << iter << std::endl;
279 }
280 }
281 }
282 return false;
283 }
284 return false;
285}
286
287const double * GSLMinimizer::MinGradient() const {
288 // return gradient (internal values)
289 return fGSLMultiMin->Gradient();
290}
291
292
293 } // end namespace Math
294
295} // end namespace ROOT
296
#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
#define g(i)
Definition RSha256.hxx: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 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)
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 SetFinalValues(const double *x, const MinimTransformFunction *func=nullptr)
double MinValue() const override
return minimum function value
MinimTransformFunction * CreateTransformation(std::vector< double > &startValues, const ROOT::Math::IMultiGradFunction *func=nullptr)
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
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
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:61
Interface (abstract class) for multi-dimensional functions providing a gradient calculation.
Definition IFunction.h:168
double Tolerance() const
absolute tolerance
Definition Minimizer.h:295
void SetMaxIterations(unsigned int maxiter)
set maximum iterations (one iteration can have many function calls)
Definition Minimizer.h:329
int fStatus
status of minimizer
Definition Minimizer.h:366
unsigned int MaxIterations() const
max iterations
Definition Minimizer.h:292
void SetPrintLevel(int level)
set print level
Definition Minimizer.h:323
int PrintLevel() const
minimizer configuration parameters
Definition Minimizer.h:286
MultiNumGradFunction class to wrap a normal function in a gradient function using numerical gradient ...
EGSLMinimizerType
enumeration specifying the types of GSL minimizers
Double_t x[n]
Definition legend1.C:17
Namespace for new Math classes and functions.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...