Logo ROOT   6.18/05
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
42namespace 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*/
53
54public:
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
86private:
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
103public:
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
226private:
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 */
#define c(i)
Definition: RSha256.hxx:101
#define g(i)
Definition: RSha256.hxx:105
int type
Definition: TGX11.cxx:120
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:52
int TestGradient(double absTol) const
test gradient (ask from solver gradient vector)
Definition: GSLMultiFit.h:195
GSLMultiFit(const GSLMultiFit &)
Copy constructor.
Definition: GSLMultiFit.h:92
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 TestDelta(double absTol, double relTol) const
test using abs and relative tolerance |dx| < absTol + relTol*|x| for every component
Definition: GSLMultiFit.h:203
GSLMultiFitFunctionWrapper fFunc
Definition: GSLMultiFit.h:228
gsl_multifit_fdfsolver * fSolver
Definition: GSLMultiFit.h:229
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
GSLMultiFit & operator=(const GSLMultiFit &rhs)
Assignment operator.
Definition: GSLMultiFit.h:97
const double * Gradient() const
gradient value at the minimum
Definition: GSLMultiFit.h:170
std::string Name() const
Definition: GSLMultiFit.h:152
const gsl_multifit_fdfsolver_type * fType
Definition: GSLMultiFit.h:237
double Edm() const
Definition: GSLMultiFit.h:209
const double * CovarMatrix() const
return covariance matrix of the parameters
Definition: GSLMultiFit.h:181
const double * X() const
parameter values at the minimum
Definition: GSLMultiFit.h:163
int Set(const std::vector< Func > &funcVec, const double *x)
set the solver parameters
Definition: GSLMultiFit.h:123
Double_t x[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21