Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Classes | Public Types | Public Member Functions | Public Attributes | Static Public Attributes | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
JFIT::JGandalf< JModel_t > Class Template Reference

Fit method based on the Levenberg-Marquardt method. More...

#include <JGandalf.hh>

Inheritance diagram for JFIT::JGandalf< JModel_t >:
JEEP::JMessage< T >

Classes

struct  result_type
 Data structure for return value of fit function. More...
 

Public Types

typedef JModel_t::parameter_type parameter_type
 Data type of fit parameter. More...
 

Public Member Functions

 JGandalf ()
 Default constructor. More...
 
template<class JFunction_t , class T , class... Args>
result_type operator() (const JFunction_t &fit, T __begin, T __end, Args...args)
 Multi-dimensional fit of two data sets. More...
 

Public Attributes

double lambda
 
JModel_t value
 
JModel_t error
 
std::vector< parameter_type > parameters
 
int numberOfIterations
 
JMATH::JMatrixNS H
 

Static Public Attributes

static int MAXIMUM_ITERATIONS = 1000
 maximal number of iterations More...
 
static double EPSILON = 1.0e-3
 maximal distance to minimum More...
 
static double LAMBDA_MIN = 0.01
 minimal value control parameter More...
 
static double LAMBDA_MAX = 100.0
 maximal value control parameter More...
 
static double LAMBDA_UP = 9.0
 multiplication factor control parameter More...
 
static double LAMBDA_DOWN = 11.0
 multiplication factor control parameter More...
 
static double PIVOT = 1.0e-3
 minimal value diagonal element of matrix More...
 
static int debug = 0
 debug level (default is off). More...
 

Private Member Functions

void reset ()
 Reset. More...
 
template<class JFunction_t , class T , class... Args>
void __evaluate__ (const JFunction_t &fit, T __begin, T __end, Args...args)
 Recursive method for evaluation of fit. More...
 
template<class JFunction_t >
void __evaluate__ (const JFunction_t &fit)
 Termination method for evaluation of fit. More...
 

Static Private Member Functions

template<class T >
static double alias (const T &model, const typename T::parameter_type parameter)
 Read/write access to parameter value by data member. More...
 
template<class T >
static double & alias (T &model, const typename T::parameter_type parameter)
 Read/write access to parameter value by data member. More...
 
template<class T >
static double alias (const T &model, const size_t index)
 Read/write access to parameter value by index. More...
 
template<class T >
static double & alias (T &model, const size_t index)
 Read/write access to parameter value by index. More...
 
template<class T >
static double alias (const T &model, const int index)
 Read/write access to parameter value by index. More...
 
template<class T >
static double & alias (T &model, const int index)
 Read/write access to parameter value by index. More...
 

Private Attributes

result_type successor
 
JModel_t previous
 
std::vector< double > h
 

Detailed Description

template<class JModel_t>
class JFIT::JGandalf< JModel_t >

Fit method based on the Levenberg-Marquardt method.

The template argument refers to the model that should be fitted to the data.
This data structure should have arithmetic capabilities.

The data member JGandalf::value corresponds to the start c.q. final value of the model of the fit procedure and JGandalf::error to the uncertainties.
The co-variance matrix is stored in data member JGandalf::H.
The data member JGandalf::parameters constitutes a list of those parameters of the model that should actually be fitted.
For this, the model should contain the type definition for parameter_type.
Normally, this type definition corresponds to a pointer to a data member of the model.
Alternatively, the type definition can be size_t or int.
In that case, the model class should provide for the element access operator[].
The template fit function in JGandalf::operator() should provide for an implementation of the function.
This operator should return the data type JGandalf::result_type which is composed of the value of the chi2 and the gradient of a data point, respectively.
The function operator returns the chi2 of the overall fit.

Definition at line 52 of file JGandalf.hh.

Member Typedef Documentation

template<class JModel_t>
typedef JModel_t::parameter_type JFIT::JGandalf< JModel_t >::parameter_type

Data type of fit parameter.

Definition at line 63 of file JGandalf.hh.

Constructor & Destructor Documentation

template<class JModel_t>
JFIT::JGandalf< JModel_t >::JGandalf ( )
inline

Default constructor.

Definition at line 111 of file JGandalf.hh.

112  {}

Member Function Documentation

template<class JModel_t>
template<class JFunction_t , class T , class... Args>
result_type JFIT::JGandalf< JModel_t >::operator() ( const JFunction_t &  fit,
T  __begin,
T  __end,
Args...  args 
)
inline

Multi-dimensional fit of two data sets.

The fit function should return the equivalent of chi2 for the current value of the model and the given data point as well as the partial derivatives.

Parameters
fitfit function
__beginbegin of data
__endend of data
Returns
chi2

Definition at line 127 of file JGandalf.hh.

128  {
129  using namespace std;
130  using namespace JPP;
131 
132  const int N = parameters.size();
133 
134  H.resize(N);
135  h.resize(N);
136 
137  previous = value;
138 
139  successor.gradient = value;
140  successor.gradient = zero;
141 
142  error = value;
143  error = zero;
144 
145  lambda = LAMBDA_MIN;
146 
147  result_type precessor = result_type(numeric_limits<double>::max(), value);
148 
149  for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
150 
151  DEBUG("step: " << numberOfIterations << endl);
152 
153  reset();
154 
155  __evaluate__(fit, __begin, __end, args...);
156 
157  DEBUG("lambda: " << SCIENTIFIC(12,5) << lambda << endl);
158  DEBUG("chi2: " << FIXED (12,5) << successor.chi2 << endl);
159 
160  if (successor.chi2 < precessor.chi2) {
161 
162  if (numberOfIterations != 0) {
163 
164  if (fabs(precessor.chi2 - successor.chi2) < EPSILON*fabs(precessor.chi2)) {
165  return successor;
166  }
167 
168  if (lambda > LAMBDA_MIN) {
169  lambda /= LAMBDA_DOWN;
170  }
171  }
172 
173  precessor = successor;
174  previous = value;
175 
176  } else {
177 
178  value = previous;
179  lambda *= LAMBDA_UP;
180 
181  if (lambda > LAMBDA_MAX) {
182  return precessor; // no improvement found
183  }
184 
185  reset();
186 
187  __evaluate__(fit, __begin, __end, args...);
188  }
189 
190  DEBUG("Hesse matrix:" << endl);
191  DEBUG(H << endl);
192 
193 
194  // force definite positiveness
195 
196  for (int i = 0; i != N; ++i) {
197 
198  if (H(i,i) < PIVOT) {
199  H(i,i) = PIVOT;
200  }
201 
202  h[i] = 1.0 / sqrt(H(i,i));
203  }
204 
205 
206  // normalisation
207 
208  for (int i = 0; i != N; ++i) {
209  for (int j = 0; j != i; ++j) {
210  H(j,i) *= h[i] * h[j];
211  H(i,j) = H(j,i);
212  }
213  }
214 
215  for (int i = 0; i != N; ++i) {
216  H(i,i) = 1.0 + lambda;
217  }
218 
219 
220  try {
221  H.invert();
222  }
223  catch (const JException& error) {
224  ERROR("JGandalf: " << error.what() << endl);
225  return precessor;
226  }
227 
228 
229  for (int i = 0; i != N; ++i) {
230 
231  DEBUG("u[" << noshowpos << i << "] = " << showpos << FIXED(15,5) << alias(value, parameters[i]));
232 
233  for (int j = 0; j != N; ++j) {
234  alias(value, parameters[i]) -= H(i,j) * alias(successor.gradient, parameters[j]) * h[i] * h[j];
235  }
236 
237  DEBUG(' ' << FIXED(15,5) << alias(value, parameters[i]) << noshowpos << endl);
238 
239  alias(error, parameters[i]) = h[i];
240  }
241  }
242 
243  return precessor;
244  }
JFIT::JGandalf::lambda
double lambda
Definition: JGandalf.hh:255
JFIT::JGandalf::EPSILON
static double EPSILON
maximal distance to minimum
Definition: JGandalf.hh:248
JFIT::JGandalf::value
JModel_t value
Definition: JGandalf.hh:256
JFIT::JGandalf::PIVOT
static double PIVOT
minimal value diagonal element of matrix
Definition: JGandalf.hh:253
JMATH::zero
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
FIXED
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:481
JFIT::JGandalf::parameters
std::vector< parameter_type > parameters
Definition: JGandalf.hh:258
JMATH::JMatrixND::resize
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:359
JFIT::JGandalf::LAMBDA_MIN
static double LAMBDA_MIN
minimal value control parameter
Definition: JGandalf.hh:249
JFIT::JGandalf::LAMBDA_DOWN
static double LAMBDA_DOWN
multiplication factor control parameter
Definition: JGandalf.hh:252
JFIT::JGandalf::previous
JModel_t previous
Definition: JGandalf.hh:403
JFIT::JGandalf::reset
void reset()
Reset.
Definition: JGandalf.hh:266
JFIT::JGandalf::successor
result_type successor
Definition: JGandalf.hh:402
JFIT::JGandalf::result_type::gradient
JModel_t gradient
d(chi2)/d(...)
Definition: JGandalf.hh:104
JFIT::JGandalf::result_type::chi2
double chi2
chi2
Definition: JGandalf.hh:103
JFIT::JGandalf::H
JMATH::JMatrixNS H
Definition: JGandalf.hh:260
JFIT::JGandalf::alias
static double alias(const T &model, const typename T::parameter_type parameter)
Read/write access to parameter value by data member.
Definition: JGandalf.hh:326
JFIT::JGandalf::LAMBDA_UP
static double LAMBDA_UP
multiplication factor control parameter
Definition: JGandalf.hh:251
JMATH::JMatrixNS::invert
void invert()
Invert matrix according LDU decomposition.
Definition: JMatrixNS.hh:80
JFIT::JGandalf::__evaluate__
void __evaluate__(const JFunction_t &fit, T __begin, T __end, Args...args)
Recursive method for evaluation of fit.
Definition: JGandalf.hh:286
JFIT::JGandalf::MAXIMUM_ITERATIONS
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JGandalf.hh:247
JTOOLS::j
int j
Definition: JPolint.hh:634
JFIT::JGandalf::numberOfIterations
int numberOfIterations
Definition: JGandalf.hh:259
JFIT::JGandalf::LAMBDA_MAX
static double LAMBDA_MAX
maximal value control parameter
Definition: JGandalf.hh:250
JFIT::JGandalf::error
JModel_t error
Definition: JGandalf.hh:257
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:518
N
then usage $script[input file[working directory[option]]] nWhere option can be N
Definition: JMuonPostfit.sh:37
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
JFIT::JGandalf::h
std::vector< double > h
Definition: JGandalf.hh:404
template<class JModel_t>
void JFIT::JGandalf< JModel_t >::reset ( )
inlineprivate

Reset.

Definition at line 266 of file JGandalf.hh.

267  {
268  using namespace std;
269  using namespace JPP;
270 
271  successor.chi2 = 0.0;
272  successor.gradient = zero;
273 
274  H.reset();
275  }
JMATH::zero
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
JMATH::JMatrixND::reset
JMatrixND & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:372
JFIT::JGandalf::successor
result_type successor
Definition: JGandalf.hh:402
JFIT::JGandalf::result_type::gradient
JModel_t gradient
d(chi2)/d(...)
Definition: JGandalf.hh:104
JFIT::JGandalf::result_type::chi2
double chi2
chi2
Definition: JGandalf.hh:103
JFIT::JGandalf::H
JMATH::JMatrixNS H
Definition: JGandalf.hh:260
template<class JModel_t>
template<class JFunction_t , class T , class... Args>
void JFIT::JGandalf< JModel_t >::__evaluate__ ( const JFunction_t &  fit,
T  __begin,
T  __end,
Args...  args 
)
inlineprivate

Recursive method for evaluation of fit.

Parameters
fitfit function
__beginbegin of data
__endend of data

Definition at line 286 of file JGandalf.hh.

287  {
288  using namespace std;
289 
290  for (T hit = __begin; hit != __end; ++hit) {
291 
292  const result_type& result = fit(value, *hit);
293 
294  successor.chi2 += result.chi2;
295  successor.gradient += result.gradient;
296 
297  for (unsigned int i = 0; i != parameters.size(); ++i) {
298  for (unsigned int j = i; j != parameters.size(); ++j) {
299  H(i,j) += alias(result.gradient, parameters[i]) * alias(result.gradient, parameters[j]);
300  }
301  }
302  }
303 
304  __evaluate__(fit, args...);
305  }
JFIT::JGandalf::value
JModel_t value
Definition: JGandalf.hh:256
JFIT::JGandalf::parameters
std::vector< parameter_type > parameters
Definition: JGandalf.hh:258
JFIT::JGandalf::successor
result_type successor
Definition: JGandalf.hh:402
JTOOLS::result
return result
Definition: JPolint.hh:695
T
do set_variable OUTPUT_DIRECTORY $WORKDIR T
Definition: JCalibrateHeight.sh:61
JFIT::JGandalf::result_type::gradient
JModel_t gradient
d(chi2)/d(...)
Definition: JGandalf.hh:104
JFIT::JGandalf::result_type::chi2
double chi2
chi2
Definition: JGandalf.hh:103
JFIT::JGandalf::H
JMATH::JMatrixNS H
Definition: JGandalf.hh:260
JFIT::JGandalf::alias
static double alias(const T &model, const typename T::parameter_type parameter)
Read/write access to parameter value by data member.
Definition: JGandalf.hh:326
JFIT::JGandalf::__evaluate__
void __evaluate__(const JFunction_t &fit, T __begin, T __end, Args...args)
Recursive method for evaluation of fit.
Definition: JGandalf.hh:286
JTOOLS::j
int j
Definition: JPolint.hh:634
template<class JModel_t>
template<class JFunction_t >
void JFIT::JGandalf< JModel_t >::__evaluate__ ( const JFunction_t &  fit)
inlineprivate

Termination method for evaluation of fit.

Parameters
fitfit function

Definition at line 314 of file JGandalf.hh.

315  {}
template<class JModel_t>
template<class T >
static double JFIT::JGandalf< JModel_t >::alias ( const T &  model,
const typename T::parameter_type  parameter 
)
inlinestaticprivate

Read/write access to parameter value by data member.

Parameters
modelmodel
parameterparameter
Returns
value

Definition at line 326 of file JGandalf.hh.

327  {
328  return model.*parameter;
329  }
template<class JModel_t>
template<class T >
static double& JFIT::JGandalf< JModel_t >::alias ( T &  model,
const typename T::parameter_type  parameter 
)
inlinestaticprivate

Read/write access to parameter value by data member.

Parameters
modelmodel
parameterparameter
Returns
value

Definition at line 340 of file JGandalf.hh.

341  {
342  return model.*parameter;
343  }
template<class JModel_t>
template<class T >
static double JFIT::JGandalf< JModel_t >::alias ( const T &  model,
const size_t  index 
)
inlinestaticprivate

Read/write access to parameter value by index.

Parameters
modelmodel
indexindex
Returns
value

Definition at line 354 of file JGandalf.hh.

355  {
356  return model[index];
357  }
template<class JModel_t>
template<class T >
static double& JFIT::JGandalf< JModel_t >::alias ( T &  model,
const size_t  index 
)
inlinestaticprivate

Read/write access to parameter value by index.

Parameters
modelmodel
indexindex
Returns
value

Definition at line 368 of file JGandalf.hh.

369  {
370  return model[index];
371  }
template<class JModel_t>
template<class T >
static double JFIT::JGandalf< JModel_t >::alias ( const T &  model,
const int  index 
)
inlinestaticprivate

Read/write access to parameter value by index.

Parameters
modelmodel
indexindex
Returns
value

Definition at line 382 of file JGandalf.hh.

383  {
384  return model[index];
385  }
template<class JModel_t>
template<class T >
static double& JFIT::JGandalf< JModel_t >::alias ( T &  model,
const int  index 
)
inlinestaticprivate

Read/write access to parameter value by index.

Parameters
modelmodel
indexindex
Returns
value

Definition at line 396 of file JGandalf.hh.

397  {
398  return model[index];
399  }

Member Data Documentation

template<class JModel_t>
int JFIT::JGandalf< JModel_t >::MAXIMUM_ITERATIONS = 1000
static

maximal number of iterations

maximal number of iterations.

Definition at line 247 of file JGandalf.hh.

template<class JModel_t>
double JFIT::JGandalf< JModel_t >::EPSILON = 1.0e-3
static

maximal distance to minimum

maximal distance to minimum.

Definition at line 248 of file JGandalf.hh.

template<class JModel_t>
double JFIT::JGandalf< JModel_t >::LAMBDA_MIN = 0.01
static

minimal value control parameter

Definition at line 249 of file JGandalf.hh.

template<class JModel_t>
double JFIT::JGandalf< JModel_t >::LAMBDA_MAX = 100.0
static

maximal value control parameter

Definition at line 250 of file JGandalf.hh.

template<class JModel_t>
double JFIT::JGandalf< JModel_t >::LAMBDA_UP = 9.0
static

multiplication factor control parameter

Definition at line 251 of file JGandalf.hh.

template<class JModel_t>
double JFIT::JGandalf< JModel_t >::LAMBDA_DOWN = 11.0
static

multiplication factor control parameter

Definition at line 252 of file JGandalf.hh.

template<class JModel_t>
double JFIT::JGandalf< JModel_t >::PIVOT = 1.0e-3
static

minimal value diagonal element of matrix

Definition at line 253 of file JGandalf.hh.

template<class JModel_t>
double JFIT::JGandalf< JModel_t >::lambda

Definition at line 255 of file JGandalf.hh.

template<class JModel_t>
JModel_t JFIT::JGandalf< JModel_t >::value

Definition at line 256 of file JGandalf.hh.

template<class JModel_t>
JModel_t JFIT::JGandalf< JModel_t >::error

Definition at line 257 of file JGandalf.hh.

template<class JModel_t>
std::vector<parameter_type> JFIT::JGandalf< JModel_t >::parameters

Definition at line 258 of file JGandalf.hh.

template<class JModel_t>
int JFIT::JGandalf< JModel_t >::numberOfIterations

Definition at line 259 of file JGandalf.hh.

template<class JModel_t>
JMATH::JMatrixNS JFIT::JGandalf< JModel_t >::H

Definition at line 260 of file JGandalf.hh.

template<class JModel_t>
result_type JFIT::JGandalf< JModel_t >::successor
private

Definition at line 402 of file JGandalf.hh.

template<class JModel_t>
JModel_t JFIT::JGandalf< JModel_t >::previous
private

Definition at line 403 of file JGandalf.hh.

template<class JModel_t>
std::vector<double> JFIT::JGandalf< JModel_t >::h
private

Definition at line 404 of file JGandalf.hh.

template<class T>
int JEEP::JMessage< T >::debug = 0
staticinherited

debug level (default is off).

Definition at line 45 of file JMessage.hh.

The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.5