ROOT  6.06/09
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 
45 namespace 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 
59 public:
60 
61  /**
62  virtual Destructor
63  */
65 
66 
67 public:
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 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 
133 private:
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 
157 public:
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  {
166  CreateSolver(type, n);
167  }
168 
169  /**
170  Destructor (no operations)
171  */
172  virtual ~GSLMultiRootSolver () {
173  if (fSolver) gsl_multiroot_fsolver_free(fSolver);
174  if (fVec != 0) gsl_vector_free(fVec);
175  }
176 
177 private:
178  // usually copying is non trivial, so we make this unaccessible
179 
180  /**
181  Copy constructor
182  */
184 
185  /**
186  Assignment operator
187  */
189  if (this == &rhs) return *this; // time saving self-test
190  return *this;
191  }
192 
193 
194 public:
195 
196 
197  void CreateSolver(const gsl_multiroot_fsolver_type * type, unsigned int n) {
198 
199  /// create the solver from the type and size of number of fitting points and number of parameters
200  if (fSolver) gsl_multiroot_fsolver_free(fSolver);
201  fSolver = gsl_multiroot_fsolver_alloc(type, n);
202 
203  }
204 
205 
206  /// set the solver parameters
207  virtual int SetSolver(const std::vector<ROOT::Math::IMultiGenFunction*> & funcVec, const double * x) {
208  // create a vector of the fit contributions
209  // create function wrapper from an iterator of functions
210  assert(fSolver !=0);
211  unsigned int n = funcVec.size();
212 
213  fFunctions.SetFunctions(funcVec, funcVec.size() );
214  // set initial values and create cached vector
215  if (fVec != 0) gsl_vector_free(fVec);
216  fVec = gsl_vector_alloc( n);
217  std::copy(x,x+n, fVec->data);
218  // solver should have been already created at this point
219  assert(fSolver != 0);
220  return gsl_multiroot_fsolver_set(fSolver, fFunctions.GetFunctions(), fVec);
221  }
222 
223  virtual std::string Name() const {
224  if (fSolver == 0 ) return "undefined";
225  return std::string(gsl_multiroot_fsolver_name(fSolver) );
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 
252 private:
253 
255  gsl_multiroot_fsolver * fSolver;
256  // cached vector to avoid re-allocating every time a new one
257  mutable gsl_vector * fVec;
258 
259 };
260 
261 /**
262  GSLMultiRootDerivSolver, internal class for implementing GSL multi-root finders
263  using derivatives
264 
265  @ingroup MultiRoot
266 */
268 
269 public:
270 
271  /**
272  Constructor
273  */
274  GSLMultiRootDerivSolver (const gsl_multiroot_fdfsolver_type * type, int n ) :
275  fDerivSolver(0),
276  fVec(0)
277  {
278  CreateSolver(type, n);
279  }
280 
281  /**
282  Destructor (no operations)
283  */
285  if (fDerivSolver) gsl_multiroot_fdfsolver_free(fDerivSolver);
286  if (fVec != 0) gsl_vector_free(fVec);
287  }
288 
289 private:
290  // usually copying is non trivial, so we make this unaccessible
291 
292  /**
293  Copy constructor
294  */
296 
297  /**
298  Assignment operator
299  */
301  if (this == &rhs) return *this; // time saving self-test
302  return *this;
303  }
304 
305 
306 public:
307 
308 
309  /// create the solver from the type and size of number of fitting points and number of parameters
310  void CreateSolver(const gsl_multiroot_fdfsolver_type * type, unsigned int n) {
311 
312  /// create the solver from the type and size of number of fitting points and number of parameters
313  if (fDerivSolver) gsl_multiroot_fdfsolver_free(fDerivSolver);
314  fDerivSolver = gsl_multiroot_fdfsolver_alloc(type, n);
315  }
316 
317 
318 
319  /// set the solver parameters for the case of derivative
320  virtual int SetSolver(const std::vector<ROOT::Math::IMultiGenFunction*> & funcVec, const double * x) {
321  // create a vector of the fit contributions
322  // need to create a vecctor of gradient functions, convert and store in the class
323  // the new vector pointer
324  assert(fDerivSolver !=0);
325  unsigned int n = funcVec.size();
326  fGradFuncVec.reserve( n );
327  for (unsigned int i = 0; i < n; ++i) {
328  ROOT::Math::IMultiGradFunction * func = dynamic_cast<ROOT::Math::IMultiGradFunction *>(funcVec[i] );
329  if (func == 0) {
330  MATH_ERROR_MSG("GSLMultiRootSolver::SetSolver","Function does not provide gradient interface");
331  return -1;
332  }
333  fGradFuncVec.push_back( func);
334  }
335 
336  fDerivFunctions.SetFunctions(fGradFuncVec, funcVec.size() );
337  // set initial values
338  if (fVec != 0) gsl_vector_free(fVec);
339  fVec = gsl_vector_alloc( n);
340  std::copy(x,x+n, fVec->data);
341 
342  return gsl_multiroot_fdfsolver_set(fDerivSolver, fDerivFunctions.GetFunctions(), fVec);
343  }
344 
345  virtual std::string Name() const {
346  if (fDerivSolver == 0 ) return "undefined";
347  return std::string(gsl_multiroot_fdfsolver_name(fDerivSolver) );
348  }
349 
350  virtual int Iterate() {
351  if (fDerivSolver == 0) return -1;
352  return gsl_multiroot_fdfsolver_iterate(fDerivSolver);
353  }
354 
355  /// solution values at the current iteration
356  virtual gsl_vector * GetRoot() const {
357  if (fDerivSolver == 0) return 0;
358  return gsl_multiroot_fdfsolver_root(fDerivSolver);
359  }
360 
361  /// return function values
362  virtual gsl_vector * GetF() const {
363  if (fDerivSolver == 0) return 0;
364  return gsl_multiroot_fdfsolver_f(fDerivSolver);
365  }
366 
367  /// return function steps
368  virtual gsl_vector * GetDx() const {
369  if (fDerivSolver == 0) return 0;
370  return gsl_multiroot_fdfsolver_dx(fDerivSolver);
371  }
372 
373 
374 
375 private:
376 
378  gsl_multiroot_fdfsolver * fDerivSolver;
379  // cached vector to avoid re-allocating every time a new one
380  mutable gsl_vector * fVec;
381  std::vector<ROOT::Math::IMultiGradFunction*> fGradFuncVec;
382 
383 };
384 
385  } // end namespace Math
386 
387 } // end namespace ROOT
388 
389 
390 #endif /* ROOT_Math_GSLMultiRootSolver */
Interface (abstract class) for multi-dimensional functions providing a gradient calculation.
Definition: IFunction.h:322
const double absTol
Namespace for new ROOT classes and functions.
Definition: ROOT.py:1
int TestDelta(double absTol, double relTol) const
test using abs and relative tolerance |dx| < absTol + relTol*|x| for every component ...
virtual gsl_vector * GetDx() const
return function steps
#define assert(cond)
Definition: unittest.h:542
virtual gsl_vector * GetF() const
return function values
GSLMultiRootBaseSolver, internal class for implementing GSL multi-root finders This is the base class...
virtual int SetSolver(const std::vector< ROOT::Math::IMultiGenFunction * > &funcVec, const double *x)
set the solver parameters for the case of derivative
GSLMultiRootSolver(const GSLMultiRootSolver &)
Copy constructor.
GSLMultiRootDerivSolver(const GSLMultiRootDerivSolver &)
Copy constructor.
gsl_multiroot_fsolver * fSolver
GSLMultiRootDerivSolver(const gsl_multiroot_fdfsolver_type *type, int n)
Constructor.
std::vector< ROOT::Math::IMultiGradFunction * > fGradFuncVec
#define MATH_ERROR_MSGVAL(loc, str, x)
Definition: Error.h:69
int TestResidual(double absTol) const
test using abs tolerance Sum |f|_i < absTol
virtual ~GSLMultiRootSolver()
Destructor (no operations)
virtual gsl_vector * GetDx() const =0
virtual gsl_vector * GetRoot() const
solution values at the current iteration
virtual std::string Name() const
return name
wrapper to a multi-dim function without derivatives for multi roots algorithm
Double_t x[n]
Definition: legend1.C:17
virtual gsl_vector * GetDx() const
return function steps
virtual ~GSLMultiRootBaseSolver()
virtual Destructor
virtual int Iterate()=0
perform an iteration
virtual gsl_vector * GetRoot() const =0
virtual ~GSLMultiRootDerivSolver()
Destructor (no operations)
const double * FVal() const
return function values
void SetFunctions(const FuncVector &f, unsigned int n)
Fill gsl function structure from a C++ function iterator and size and number of residuals.
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:50
virtual std::string Name() const
return name
GSLMultiRootSolver & operator=(const GSLMultiRootSolver &rhs)
Assignment operator.
virtual int Iterate()
perform an iteration
GSLMultiRootDerivSolver & operator=(const GSLMultiRootDerivSolver &rhs)
Assignment operator.
GSLMultiRootSolver, internal class for implementing GSL multi-root finders not using derivatives...
virtual gsl_vector * GetRoot() const
solution values at the current iteration
GSLMultiRootFunctionWrapper fFunctions
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 ...
void SetFunctions(const FuncVector &f, unsigned int n)
Fill gsl function structure from a C++ function iterator and size and number of residuals.
double f(double x)
const double * X() const
solution values at the current iteration
virtual std::string Name() const =0
return name
GSLMultiRootSolver(const gsl_multiroot_fsolver_type *type, int n)
Constructor from type and simension of system (number of functions)
int type
Definition: TGX11.cxx:120
double func(double *x, double *p)
Definition: stressTF1.cxx:213
virtual int SetSolver(const std::vector< ROOT::Math::IMultiGenFunction * > &funcVec, const double *x)
set the solver parameters
const double * Dx() const
return function steps
Namespace for new Math classes and functions.
wrapper to a multi-dim function with derivatives for multi roots algorithm
gsl_multiroot_fdfsolver * fDerivSolver
virtual int Iterate()
perform an iteration
GSLMultiRootDerivFunctionWrapper fDerivFunctions
GSLMultiRootDerivSolver, internal class for implementing GSL multi-root finders using derivatives...
void CreateSolver(const gsl_multiroot_fsolver_type *type, unsigned int n)
virtual gsl_vector * GetF() const =0
const Int_t n
Definition: legend1.C:16
bool InitSolver(const std::vector< ROOT::Math::IMultiGenFunction * > &funcVec, const double *x)
init the solver with function list and initial values
virtual gsl_vector * GetF() const
return function values
virtual int SetSolver(const std::vector< ROOT::Math::IMultiGenFunction * > &funcVec, const double *x)=0