Logo ROOT   6.13/01
Reference Guide
GSLMultiFit.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 class GSLMultiFit
26 
27 #ifndef ROOT_Math_GSLMultiFit
28 #define ROOT_Math_GSLMultiFit
29 
30 #include "gsl/gsl_vector.h"
31 #include "gsl/gsl_matrix.h"
32 #include "gsl/gsl_multifit_nlin.h"
33 #include "gsl/gsl_blas.h"
34 #include "gsl/gsl_version.h"
36 
37 #include "Math/IFunction.h"
38 #include <string>
39 #include <cassert>
40 
41 
42 namespace ROOT {
43 
44  namespace Math {
45 
46 
47 /**
48  GSLMultiFit, internal class for implementing GSL non linear least square GSL fitting
49 
50  @ingroup MultiMin
51 */
52 class GSLMultiFit {
53 
54 public:
55 
56  /**
57  Default constructor
58  No need to specify the type so far since only one solver exists so far
59  */
60  GSLMultiFit (const gsl_multifit_fdfsolver_type * type = 0) :
61  fSolver(0),
62  fVec(0),
63  fTmp(0),
64  fCov(0),
65 #if GSL_MAJOR_VERSION > 1
66  fJac(0),
67 #endif
68  fType(type)
69  {
70  if (fType == 0) fType = gsl_multifit_fdfsolver_lmsder; // default value
71  }
72 
73  /**
74  Destructor (no operations)
75  */
77  if (fSolver) gsl_multifit_fdfsolver_free(fSolver);
78  if (fVec != 0) gsl_vector_free(fVec);
79  if (fTmp != 0) gsl_vector_free(fTmp);
80  if (fCov != 0) gsl_matrix_free(fCov);
81 #if GSL_MAJOR_VERSION > 1
82  if (fJac != 0) gsl_matrix_free(fJac);
83 #endif
84  }
85 
86 private:
87  // usually copying is non trivial, so we make this unaccessible
88 
89  /**
90  Copy constructor
91  */
93 
94  /**
95  Assignment operator
96  */
98  if (this == &rhs) return *this; // time saving self-test
99  return *this;
100  }
101 
102 
103 public:
104 
105  /// create the minimizer from the type and size of number of fitting points and number of parameters
106  void CreateSolver(unsigned int npoints, unsigned int npar) {
107  if (fSolver) gsl_multifit_fdfsolver_free(fSolver);
108  fSolver = gsl_multifit_fdfsolver_alloc(fType, npoints, npar);
109  if (fVec != 0) gsl_vector_free(fVec);
110  fVec = gsl_vector_alloc( npar );
111  if (fTmp != 0) gsl_vector_free(fTmp);
112  fTmp = gsl_vector_alloc( npar );
113  if (fCov != 0) gsl_matrix_free(fCov);
114  fCov = gsl_matrix_alloc( npar, npar );
115 #if GSL_MAJOR_VERSION > 1
116  if (fJac != 0) gsl_matrix_free(fJac);
117  fJac = gsl_matrix_alloc( npoints, npar );
118 #endif
119  }
120 
121  /// set the solver parameters
122  template<class Func>
123  int Set(const std::vector<Func> & funcVec, const double * x) {
124  // create a vector of the fit contributions
125  // create function wrapper from an iterator of functions
126  unsigned int npts = funcVec.size();
127  if (npts == 0) return -1;
128 
129  unsigned int npar = funcVec[0].NDim();
130  // Remove unused typedef to remove warning in GCC48
131  // http://gcc.gnu.org/gcc-4.8/porting_to.html
132  // typedef typename std::vector<Func> FuncVec;
133  //FuncIt funcIter = funcVec.begin();
134  fFunc.SetFunction(funcVec, npts, npar);
135  // create solver object
136  CreateSolver( npts, npar );
137  assert(fSolver != 0);
138  // set initial values
139  assert(fVec != 0);
140  std::copy(x,x+npar, fVec->data);
141  assert(fTmp != 0);
142  gsl_vector_set_zero(fTmp);
143  assert(fCov != 0);
144  gsl_matrix_set_zero(fCov);
145 #if GSL_MAJOR_VERSION > 1
146  assert(fJac != 0);
147  gsl_matrix_set_zero(fJac);
148 #endif
149  return gsl_multifit_fdfsolver_set(fSolver, fFunc.GetFunc(), fVec);
150  }
151 
152  std::string Name() const {
153  if (fSolver == 0) return "undefined";
154  return std::string(gsl_multifit_fdfsolver_name(fSolver) );
155  }
156 
157  int Iterate() {
158  if (fSolver == 0) return -1;
159  return gsl_multifit_fdfsolver_iterate(fSolver);
160  }
161 
162  /// parameter values at the minimum
163  const double * X() const {
164  if (fSolver == 0) return 0;
165  gsl_vector * x = gsl_multifit_fdfsolver_position(fSolver);
166  return x->data;
167  }
168 
169  /// gradient value at the minimum
170  const double * Gradient() const {
171  if (fSolver == 0) return 0;
172 #if GSL_MAJOR_VERSION > 1
173  fType->gradient(fSolver->state, fVec);
174 #else
175  gsl_multifit_gradient(fSolver->J, fSolver->f,fVec);
176 #endif
177  return fVec->data;
178  }
179 
180  /// return covariance matrix of the parameters
181  const double * CovarMatrix() const {
182  if (fSolver == 0) return 0;
183  static double kEpsrel = 0.0001;
184 #if GSL_MAJOR_VERSION > 1
185  gsl_multifit_fdfsolver_jac (fSolver, fJac);
186  int ret = gsl_multifit_covar(fJac, kEpsrel, fCov);
187 #else
188  int ret = gsl_multifit_covar(fSolver->J, kEpsrel, fCov);
189 #endif
190  if (ret != GSL_SUCCESS) return 0;
191  return fCov->data;
192  }
193 
194  /// test gradient (ask from solver gradient vector)
195  int TestGradient(double absTol) const {
196  if (fSolver == 0) return -1;
197  Gradient();
198  return gsl_multifit_test_gradient( fVec, absTol);
199  }
200 
201  /// test using abs and relative tolerance
202  /// |dx| < absTol + relTol*|x| for every component
203  int TestDelta(double absTol, double relTol) const {
204  if (fSolver == 0) return -1;
205  return gsl_multifit_test_delta(fSolver->dx, fSolver->x, absTol, relTol);
206  }
207 
208  // calculate edm 1/2 * ( grad^T * Cov * grad )
209  double Edm() const {
210  // product C * g
211  double edm = -1;
212  const double * g = Gradient();
213  if (g == 0) return edm;
214  const double * c = CovarMatrix();
215  if (c == 0) return edm;
216  if (fTmp == 0) return edm;
217  int status = gsl_blas_dgemv(CblasNoTrans, 1.0, fCov, fVec, 0.,fTmp);
218  if (status == 0) status |= gsl_blas_ddot(fVec, fTmp, &edm);
219  if (status != 0) return -1;
220  // need to divide by 2 ??
221  return 0.5*edm;
222 
223  }
224 
225 
226 private:
227 
229  gsl_multifit_fdfsolver * fSolver;
230  // cached vector to avoid re-allocating every time a new one
231  mutable gsl_vector * fVec;
232  mutable gsl_vector * fTmp;
233  mutable gsl_matrix * fCov;
234 #if GSL_MAJOR_VERSION > 1
235  mutable gsl_matrix * fJac;
236 #endif
237  const gsl_multifit_fdfsolver_type * fType;
238 
239 };
240 
241  } // end namespace Math
242 
243 } // end namespace ROOT
244 
245 
246 #endif /* ROOT_Math_GSLMultiFit */
const gsl_multifit_fdfsolver_type * fType
Definition: GSLMultiFit.h:237
Namespace for new ROOT classes and functions.
Definition: TFoamSampler.h:19
GSLMultiFit, internal class for implementing GSL non linear least square GSL fitting.
Definition: GSLMultiFit.h:52
GSLMultiFit(const gsl_multifit_fdfsolver_type *type=0)
Default constructor No need to specify the type so far since only one solver exists so far...
Definition: GSLMultiFit.h:60
int TestDelta(double absTol, double relTol) const
test using abs and relative tolerance |dx| < absTol + relTol*|x| for every component ...
Definition: GSLMultiFit.h:203
int Set(const std::vector< Func > &funcVec, const double *x)
set the solver parameters
Definition: GSLMultiFit.h:123
gsl_multifit_fdfsolver * fSolver
Definition: GSLMultiFit.h:229
GSLMultiFit(const GSLMultiFit &)
Copy constructor.
Definition: GSLMultiFit.h:92
const double * CovarMatrix() const
return covariance matrix of the parameters
Definition: GSLMultiFit.h:181
const double * Gradient() const
gradient value at the minimum
Definition: GSLMultiFit.h:170
GSLMultiFit & operator=(const GSLMultiFit &rhs)
Assignment operator.
Definition: GSLMultiFit.h:97
std::string Name() const
Definition: GSLMultiFit.h:152
double Edm() const
Definition: GSLMultiFit.h:209
* x
Deprecated and error prone model selection interface.
Definition: TRolke.cxx:630
GSLMultiFitFunctionWrapper fFunc
Definition: GSLMultiFit.h:228
void SetFunction(const FuncVector &f, unsigned int nres, unsigned int npar)
Fill gsl function structure from a C++ function iterator and size and number of residuals.
const double * X() const
parameter values at the minimum
Definition: GSLMultiFit.h:163
Namespace for new Math classes and functions.
wrapper to a multi-dim function withtout derivatives for multi-dimensional minimization algorithm ...
void CreateSolver(unsigned int npoints, unsigned int npar)
create the minimizer from the type and size of number of fitting points and number of parameters ...
Definition: GSLMultiFit.h:106
~GSLMultiFit()
Destructor (no operations)
Definition: GSLMultiFit.h:76
int TestGradient(double absTol) const
test gradient (ask from solver gradient vector)
Definition: GSLMultiFit.h:195