Logo ROOT   6.10/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 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 
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  fName(std::string("undefined"))
166  {
167  CreateSolver(type, n);
168  }
169 
170  /**
171  Destructor (no operations)
172  */
173  virtual ~GSLMultiRootSolver () {
174  if (fSolver) gsl_multiroot_fsolver_free(fSolver);
175  if (fVec != 0) gsl_vector_free(fVec);
176  }
177 
178 private:
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 
195 public:
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 
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  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 
270 public:
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  {
280  CreateSolver(type, n);
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 
291 private:
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 
308 public:
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 
377 private:
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 */
Interface (abstract class) for multi-dimensional functions providing a gradient calculation.
Definition: IFunction.h:330
const double absTol
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
GSLMultiRootBaseSolver, internal class for implementing GSL multi-root finders This is the base class...
virtual gsl_vector * GetRoot() const
solution values at the current iteration
GSLMultiRootSolver(const GSLMultiRootSolver &)
Copy constructor.
virtual gsl_vector * GetRoot() const
solution values at the current iteration
virtual int SetSolver(const std::vector< ROOT::Math::IMultiGenFunction *> &funcVec, const double *x)=0
virtual const std::string & Name() const
return name
GSLMultiRootDerivSolver(const GSLMultiRootDerivSolver &)
Copy constructor.
gsl_multiroot_fsolver * fSolver
STL namespace.
GSLMultiRootDerivSolver(const gsl_multiroot_fdfsolver_type *type, int n)
Constructor.
std::vector< ROOT::Math::IMultiGradFunction * > fGradFuncVec
int TestResidual(double absTol) const
test using abs tolerance Sum |f|_i < absTol
#define MATH_ERROR_MSGVAL(loc, str, x)
Definition: Error.h:69
virtual ~GSLMultiRootSolver()
Destructor (no operations)
virtual gsl_vector * GetDx() const =0
bool InitSolver(const std::vector< ROOT::Math::IMultiGenFunction *> &funcVec, const double *x)
init the solver with function list and initial values
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 without derivatives for multi roots algorithm
Double_t x[n]
Definition: legend1.C:17
virtual ~GSLMultiRootBaseSolver()
virtual Destructor
const double * FVal() const
return function values
virtual int Iterate()=0
perform an iteration
virtual gsl_vector * GetRoot() const =0
virtual ~GSLMultiRootDerivSolver()
Destructor (no operations)
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:50
virtual gsl_vector * GetF() const
return function values
virtual const std::string & Name() const =0
return name
virtual int Iterate()
perform an iteration
const double * Dx() const
return function steps
virtual gsl_vector * GetF() const
return function values
GSLMultiRootSolver, internal class for implementing GSL multi-root finders not using derivatives...
virtual int SetSolver(const std::vector< ROOT::Math::IMultiGenFunction *> &funcVec, const double *x)
set the solver parameters for the case of derivative
virtual int SetSolver(const std::vector< ROOT::Math::IMultiGenFunction *> &funcVec, const double *x)
set the solver parameters
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 ...
virtual gsl_vector * GetDx() const
return function steps
double f(double x)
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
Namespace for new Math classes and functions.
wrapper to a multi-dim function with derivatives for multi roots algorithm
Binding & operator=(OUT(*fun)(void))
gsl_multiroot_fdfsolver * fDerivSolver
virtual gsl_vector * GetDx() const
return function steps
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)
const double * X() const
solution values at the current iteration
virtual gsl_vector * GetF() const =0
virtual const std::string & Name() const
return name
const Int_t n
Definition: legend1.C:16