Logo ROOT  
Reference Guide
GSLMultiRootSolver.h
Go to the documentation of this file.
1// @(#)root/mathmore:$Id$
2// Author: L. Moneta Wed Dec 20 17:26:06 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// Header file for the class GSLMultiRootBaseSolver,
26// GSLMultiRootSolver and GSLMultiRootDerivSolver
27
28#ifndef ROOT_Math_GSLMultiRootSolver
29#define ROOT_Math_GSLMultiRootSolver
30
31#include "gsl/gsl_vector.h"
32#include "gsl/gsl_matrix.h"
33#include "gsl/gsl_multiroots.h"
34#include "gsl/gsl_blas.h"
36
37#include "Math/IFunction.h"
38#include "Math/Error.h"
39
40#include <vector>
41#include <string>
42#include <cassert>
43
44
45namespace ROOT {
46
47 namespace Math {
48
49
50/**
51 GSLMultiRootBaseSolver, internal class for implementing GSL multi-root finders
52 This is the base class for GSLMultiRootSolver (solver not using derivatives) and
53 GSLMUltiRootDerivSolver (solver using derivatives)
54
55 @ingroup MultiRoot
56*/
58
59public:
60
61 /**
62 virtual Destructor
63 */
65
66
67public:
68
69
70 /// init the solver with function list and initial values
71 bool InitSolver(const std::vector<ROOT::Math::IMultiGenFunction*> & funcVec, const double * x) {
72 // create a vector of the fit contributions
73 // create function wrapper from an iterator of functions
74 unsigned int n = funcVec.size();
75 if (n == 0) return false;
76
77 unsigned int ndim = funcVec[0]->NDim(); // should also be = n
78
79 if (ndim != n) {
80 MATH_ERROR_MSGVAL("GSLMultiRootSolver::InitSolver","Wrong function dimension",ndim);
81 MATH_ERROR_MSGVAL("GSLMultiRootSolver::InitSolver","Number of functions",n);
82 return false;
83 }
84
85
86 // set function list and initial values in solver
87 int iret = SetSolver(funcVec,x);
88 return (iret == 0);
89 }
90
91 /// return name
92 virtual const std::string & Name() const = 0;
93
94 /// perform an iteration
95 virtual int Iterate() = 0;
96
97 /// solution values at the current iteration
98 const double * X() const {
99 gsl_vector * x = GetRoot();
100 return x->data;
101 }
102
103 /// return function values
104 const double * FVal() const {
105 gsl_vector * f = GetF();
106 return f->data;
107 }
108
109 /// return function steps
110 const double * Dx() const {
111 gsl_vector * dx = GetDx();
112 return dx->data;
113 }
114
115 /// test using abs and relative tolerance
116 /// |dx| < absTol + relTol*|x| for every component
117 int TestDelta(double absTol, double relTol) const {
118 gsl_vector * x = GetRoot();
119 gsl_vector * dx = GetDx();
120 if (x == 0 || dx == 0) return -1;
121 return gsl_multiroot_test_delta(dx, x, absTol, relTol);
122 }
123
124 /// test using abs tolerance
125 /// Sum |f|_i < absTol
126 int TestResidual(double absTol) const {
127 gsl_vector * f = GetF();
128 if (f == 0) return -1;
129 return gsl_multiroot_test_residual(f, absTol);
130 }
131
132
133private:
134
135 // accessor to be implemented by the derived classes
136
137 virtual int SetSolver(const std::vector<ROOT::Math::IMultiGenFunction*> & funcVec, const double * x) = 0;
138
139 virtual gsl_vector * GetRoot() const = 0;
140
141 virtual gsl_vector * GetF() const = 0;
142
143 virtual gsl_vector * GetDx() const = 0;
144
145
146};
147
148
149/**
150 GSLMultiRootSolver, internal class for implementing GSL multi-root finders
151 not using derivatives
152
153 @ingroup MultiRoot
154*/
156
157public:
158
159 /**
160 Constructor from type and simension of system (number of functions)
161 */
162 GSLMultiRootSolver (const gsl_multiroot_fsolver_type * type, int n ) :
163 fSolver(0),
164 fVec(0),
165 fName(std::string("undefined"))
166 {
168 }
169
170 /**
171 Destructor (no operations)
172 */
174 if (fSolver) gsl_multiroot_fsolver_free(fSolver);
175 if (fVec != 0) gsl_vector_free(fVec);
176 }
177
178private:
179 // usually copying is non trivial, so we make this unaccessible
180
181 /**
182 Copy constructor
183 */
185
186 /**
187 Assignment operator
188 */
190 if (this == &rhs) return *this; // time saving self-test
191 return *this;
192 }
193
194
195public:
196
197
198 void CreateSolver(const gsl_multiroot_fsolver_type * type, unsigned int n) {
199
200 /// create the solver from the type and size of number of fitting points and number of parameters
201 if (fSolver) gsl_multiroot_fsolver_free(fSolver);
202 fSolver = gsl_multiroot_fsolver_alloc(type, n);
203 fName = std::string(gsl_multiroot_fsolver_name(fSolver) );
204 }
205
206
207 /// set the solver parameters
208 virtual int SetSolver(const std::vector<ROOT::Math::IMultiGenFunction*> & funcVec, const double * x) {
209 // create a vector of the fit contributions
210 // create function wrapper from an iterator of functions
211 assert(fSolver !=0);
212 unsigned int n = funcVec.size();
213
214 fFunctions.SetFunctions(funcVec, funcVec.size() );
215 // set initial values and create cached vector
216 if (fVec != 0) gsl_vector_free(fVec);
217 fVec = gsl_vector_alloc( n);
218 std::copy(x,x+n, fVec->data);
219 // solver should have been already created at this point
220 assert(fSolver != 0);
221 return gsl_multiroot_fsolver_set(fSolver, fFunctions.GetFunctions(), fVec);
222 }
223
224 virtual const std::string & Name() const {
225 return fName;
226 }
227
228 virtual int Iterate() {
229 if (fSolver == 0) return -1;
230 return gsl_multiroot_fsolver_iterate(fSolver);
231 }
232
233 /// solution values at the current iteration
234 virtual gsl_vector * GetRoot() const {
235 if (fSolver == 0) return 0;
236 return gsl_multiroot_fsolver_root(fSolver);
237 }
238
239 /// return function values
240 virtual gsl_vector * GetF() const {
241 if (fSolver == 0) return 0;
242 return gsl_multiroot_fsolver_f(fSolver);
243 }
244
245 /// return function steps
246 virtual gsl_vector * GetDx() const {
247 if (fSolver == 0) return 0;
248 return gsl_multiroot_fsolver_dx(fSolver);
249 }
250
251
252private:
253
255 gsl_multiroot_fsolver * fSolver;
256 // cached vector to avoid re-allocating every time a new one
257 mutable gsl_vector * fVec;
258 std::string fName; // solver nane
259
260};
261
262/**
263 GSLMultiRootDerivSolver, internal class for implementing GSL multi-root finders
264 using derivatives
265
266 @ingroup MultiRoot
267*/
269
270public:
271
272 /**
273 Constructor
274 */
275 GSLMultiRootDerivSolver (const gsl_multiroot_fdfsolver_type * type, int n ) :
276 fDerivSolver(0),
277 fVec(0),
278 fName(std::string("undefined"))
279 {
281 }
282
283 /**
284 Destructor (no operations)
285 */
287 if (fDerivSolver) gsl_multiroot_fdfsolver_free(fDerivSolver);
288 if (fVec != 0) gsl_vector_free(fVec);
289 }
290
291private:
292 // usually copying is non trivial, so we make this unaccessible
293
294 /**
295 Copy constructor
296 */
298
299 /**
300 Assignment operator
301 */
303 if (this == &rhs) return *this; // time saving self-test
304 return *this;
305 }
306
307
308public:
309
310
311 /// create the solver from the type and size of number of fitting points and number of parameters
312 void CreateSolver(const gsl_multiroot_fdfsolver_type * type, unsigned int n) {
313
314 /// create the solver from the type and size of number of fitting points and number of parameters
315 if (fDerivSolver) gsl_multiroot_fdfsolver_free(fDerivSolver);
316 fDerivSolver = gsl_multiroot_fdfsolver_alloc(type, n);
317 fName = std::string(gsl_multiroot_fdfsolver_name(fDerivSolver) );
318 }
319
320
321
322 /// set the solver parameters for the case of derivative
323 virtual int SetSolver(const std::vector<ROOT::Math::IMultiGenFunction*> & funcVec, const double * x) {
324 // create a vector of the fit contributions
325 // need to create a vecctor of gradient functions, convert and store in the class
326 // the new vector pointer
327 assert(fDerivSolver !=0);
328 unsigned int n = funcVec.size();
329 fGradFuncVec.reserve( n );
330 for (unsigned int i = 0; i < n; ++i) {
331 ROOT::Math::IMultiGradFunction * func = dynamic_cast<ROOT::Math::IMultiGradFunction *>(funcVec[i] );
332 if (func == 0) {
333 MATH_ERROR_MSG("GSLMultiRootSolver::SetSolver","Function does not provide gradient interface");
334 return -1;
335 }
336 fGradFuncVec.push_back( func);
337 }
338
339 fDerivFunctions.SetFunctions(fGradFuncVec, funcVec.size() );
340 // set initial values
341 if (fVec != 0) gsl_vector_free(fVec);
342 fVec = gsl_vector_alloc( n);
343 std::copy(x,x+n, fVec->data);
344
345 return gsl_multiroot_fdfsolver_set(fDerivSolver, fDerivFunctions.GetFunctions(), fVec);
346 }
347
348 virtual const std::string & Name() const {
349 return fName;
350 }
351
352 virtual int Iterate() {
353 if (fDerivSolver == 0) return -1;
354 return gsl_multiroot_fdfsolver_iterate(fDerivSolver);
355 }
356
357 /// solution values at the current iteration
358 virtual gsl_vector * GetRoot() const {
359 if (fDerivSolver == 0) return 0;
360 return gsl_multiroot_fdfsolver_root(fDerivSolver);
361 }
362
363 /// return function values
364 virtual gsl_vector * GetF() const {
365 if (fDerivSolver == 0) return 0;
366 return gsl_multiroot_fdfsolver_f(fDerivSolver);
367 }
368
369 /// return function steps
370 virtual gsl_vector * GetDx() const {
371 if (fDerivSolver == 0) return 0;
372 return gsl_multiroot_fdfsolver_dx(fDerivSolver);
373 }
374
375
376
377private:
378
380 gsl_multiroot_fdfsolver * fDerivSolver;
381 // cached vector to avoid re-allocating every time a new one
382 mutable gsl_vector * fVec;
383 std::vector<ROOT::Math::IMultiGradFunction*> fGradFuncVec;
384 std::string fName; // solver nane
385
386};
387
388 } // end namespace Math
389
390} // end namespace ROOT
391
392
393#endif /* ROOT_Math_GSLMultiRootSolver */
#define MATH_ERROR_MSGVAL(loc, txt, x)
Definition: Error.h:108
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:82
#define f(i)
Definition: RSha256.hxx:104
int type
Definition: TGX11.cxx:120
GSLMultiRootBaseSolver, internal class for implementing GSL multi-root finders This is the base class...
int TestResidual(double absTol) const
test using abs tolerance Sum |f|_i < absTol
virtual gsl_vector * GetDx() const =0
virtual int SetSolver(const std::vector< ROOT::Math::IMultiGenFunction * > &funcVec, const double *x)=0
virtual gsl_vector * GetRoot() const =0
virtual gsl_vector * GetF() const =0
const double * FVal() const
return function values
const double * X() const
solution values at the current iteration
virtual const std::string & Name() const =0
return name
virtual ~GSLMultiRootBaseSolver()
virtual Destructor
const double * Dx() const
return function steps
bool InitSolver(const std::vector< ROOT::Math::IMultiGenFunction * > &funcVec, const double *x)
init the solver with function list and initial values
virtual int Iterate()=0
perform an iteration
int TestDelta(double absTol, double relTol) const
test using abs and relative tolerance |dx| < absTol + relTol*|x| for every component
wrapper to a multi-dim function with derivatives for multi roots algorithm
void SetFunctions(const FuncVector &f, unsigned int n)
Fill gsl function structure from a C++ function iterator and size and number of residuals.
GSLMultiRootDerivSolver, internal class for implementing GSL multi-root finders using derivatives.
gsl_multiroot_fdfsolver * fDerivSolver
virtual const std::string & Name() const
return name
virtual gsl_vector * GetF() const
return function values
virtual int Iterate()
perform an iteration
virtual ~GSLMultiRootDerivSolver()
Destructor (no operations)
virtual gsl_vector * GetRoot() const
solution values at the current iteration
GSLMultiRootDerivSolver & operator=(const GSLMultiRootDerivSolver &rhs)
Assignment operator.
GSLMultiRootDerivSolver(const gsl_multiroot_fdfsolver_type *type, int n)
Constructor.
std::vector< ROOT::Math::IMultiGradFunction * > fGradFuncVec
void CreateSolver(const gsl_multiroot_fdfsolver_type *type, unsigned int n)
create the solver from the type and size of number of fitting points and number of parameters
GSLMultiRootDerivFunctionWrapper fDerivFunctions
virtual gsl_vector * GetDx() const
return function steps
GSLMultiRootDerivSolver(const GSLMultiRootDerivSolver &)
Copy constructor.
virtual int SetSolver(const std::vector< ROOT::Math::IMultiGenFunction * > &funcVec, const double *x)
set the solver parameters for the case of derivative
wrapper to a multi-dim function without derivatives for multi roots algorithm
void SetFunctions(const FuncVector &f, unsigned int n)
Fill gsl function structure from a C++ function iterator and size and number of residuals.
GSLMultiRootSolver, internal class for implementing GSL multi-root finders not using derivatives.
virtual gsl_vector * GetF() const
return function values
virtual gsl_vector * GetRoot() const
solution values at the current iteration
virtual ~GSLMultiRootSolver()
Destructor (no operations)
void CreateSolver(const gsl_multiroot_fsolver_type *type, unsigned int n)
GSLMultiRootFunctionWrapper fFunctions
gsl_multiroot_fsolver * fSolver
GSLMultiRootSolver(const gsl_multiroot_fsolver_type *type, int n)
Constructor from type and simension of system (number of functions)
virtual int SetSolver(const std::vector< ROOT::Math::IMultiGenFunction * > &funcVec, const double *x)
set the solver parameters
virtual const std::string & Name() const
return name
virtual int Iterate()
perform an iteration
GSLMultiRootSolver & operator=(const GSLMultiRootSolver &rhs)
Assignment operator.
virtual gsl_vector * GetDx() const
return function steps
GSLMultiRootSolver(const GSLMultiRootSolver &)
Copy constructor.
Interface (abstract class) for multi-dimensional functions providing a gradient calculation.
Definition: IFunction.h:327
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
VSD Structures.
Definition: StringConv.hxx:21