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