Jpp - the software that should make you happy
 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_typeparameters
 
int numberOfIterations
 
JMATH::JMatrixNS V
 

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

static double alias (const JModel_t &model, double JModel_t::*parameter)
 Read/write access to parameter value by data member. More...
 
static double & alias (JModel_t &model, double JModel_t::*parameter)
 Read/write access to parameter value by data member. More...
 
static double alias (const JModel_t &model, const size_t index)
 Read/write access to parameter value by index. More...
 
static double & alias (JModel_t &model, const size_t index)
 Read/write access to parameter value by index. More...
 
static double alias (const JModel_t &model, const int index)
 Read/write access to parameter value by index. More...
 
static double & alias (JModel_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::V.
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
argsoptional data
Returns
chi2

Definition at line 128 of file JGandalf.hh.

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

Reset.

Definition at line 267 of file JGandalf.hh.

268  {
269  using namespace std;
270  using namespace JPP;
271 
272  successor.chi2 = 0.0;
274 
275  V.reset();
276  }
JMATH::JMatrixNS V
Definition: JGandalf.hh:261
static const JZero zero
Function object to assign zero value.
Definition: JZero.hh:105
JMatrixND & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:375
result_type successor
Definition: JGandalf.hh:398
JModel_t gradient
d(chi2)/d(...)
Definition: JGandalf.hh:104
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
argsoptional data

Definition at line 288 of file JGandalf.hh.

289  {
290  using namespace std;
291 
292  for (T hit = __begin; hit != __end; ++hit) {
293 
294  const result_type& result = fit(value, *hit);
295 
296  successor.chi2 += result.chi2;
297  successor.gradient += result.gradient;
298 
299  for (unsigned int i = 0; i != parameters.size(); ++i) {
300  for (unsigned int j = i; j != parameters.size(); ++j) {
301  V(i,j) += alias(result.gradient, parameters[i]) * alias(result.gradient, parameters[j]);
302  }
303  }
304  }
305 
306  __evaluate__(fit, args...);
307  }
JModel_t value
Definition: JGandalf.hh:257
static double alias(const JModel_t &model, double JModel_t::*parameter)
Read/write access to parameter value by data member.
Definition: JGandalf.hh:327
JMATH::JMatrixNS V
Definition: JGandalf.hh:261
std::vector< parameter_type > parameters
Definition: JGandalf.hh:259
result_type successor
Definition: JGandalf.hh:398
return result
Definition: JPolint.hh:727
do set_variable OUTPUT_DIRECTORY $WORKDIR T
JModel_t gradient
d(chi2)/d(...)
Definition: JGandalf.hh:104
void __evaluate__(const JFunction_t &fit, T __begin, T __end, Args...args)
Recursive method for evaluation of fit.
Definition: JGandalf.hh:288
int j
Definition: JPolint.hh:666
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 316 of file JGandalf.hh.

317  {}
template<class JModel_t>
static double JFIT::JGandalf< JModel_t >::alias ( const JModel_t &  model,
double JModel_t::*  parameter 
)
inlinestaticprivate

Read/write access to parameter value by data member.

Parameters
modelmodel
parameterparameter
Returns
value

Definition at line 327 of file JGandalf.hh.

328  {
329  return model.*parameter;
330  }
template<class JModel_t>
static double& JFIT::JGandalf< JModel_t >::alias ( JModel_t &  model,
double JModel_t::*  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>
static double JFIT::JGandalf< JModel_t >::alias ( const JModel_t &  model,
const size_t  index 
)
inlinestaticprivate

Read/write access to parameter value by index.

Parameters
modelmodel
indexindex
Returns
value

Definition at line 353 of file JGandalf.hh.

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

Read/write access to parameter value by index.

Parameters
modelmodel
indexindex
Returns
value

Definition at line 366 of file JGandalf.hh.

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

Read/write access to parameter value by index.

Parameters
modelmodel
indexindex
Returns
value

Definition at line 379 of file JGandalf.hh.

380  {
381  return model[index];
382  }
template<class JModel_t>
static double& JFIT::JGandalf< JModel_t >::alias ( JModel_t &  model,
const int  index 
)
inlinestaticprivate

Read/write access to parameter value by index.

Parameters
modelmodel
indexindex
Returns
value

Definition at line 392 of file JGandalf.hh.

393  {
394  return model[index];
395  }

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 248 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 249 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 250 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 251 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 252 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 253 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 254 of file JGandalf.hh.

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

Definition at line 256 of file JGandalf.hh.

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

Definition at line 257 of file JGandalf.hh.

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

Definition at line 258 of file JGandalf.hh.

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

Definition at line 259 of file JGandalf.hh.

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

Definition at line 260 of file JGandalf.hh.

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

Definition at line 261 of file JGandalf.hh.

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

Definition at line 398 of file JGandalf.hh.

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

Definition at line 399 of file JGandalf.hh.

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

Definition at line 400 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: