Jpp  18.6.0-rc.1
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 > JROOT::JRootfit_t< JFs_t > JROOT::JRootfit< JFs_t >

Classes

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

Public Types

typedef JFIT_LOCAL::JTypedef_t
< 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 multiple data sets. More...
 

Public Attributes

std::vector< parameter_typeparameters
 fit parameters More...
 
int numberOfIterations
 number of iterations More...
 
double lambda
 control parameter More...
 
JModel_t value
 value More...
 
JModel_t error
 error More...
 
JMATH::JMatrixNS V
 Hesse matrix. More...
 
result_type result
 

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 = 10.0
 multiplication factor control parameter More...
 
static double LAMBDA_DOWN = 10.0
 multiplication factor control parameter More...
 
static double PIVOT = std::numeric_limits<double>::epsilon()
 minimal value diagonal element of Hesse matrix More...
 
static int debug = 0
 debug level (default is off). More...
 

Private Member Functions

void reset ()
 Reset current parameters. More...
 
template<class JFunction_t , class T , class... Args>
void update (const JFunction_t &fit, T __begin, T __end, Args...args)
 Recursive method to update current parameters. More...
 
template<class JFunction_t >
void update (const JFunction_t &fit)
 Termination method to update current parameters. More...
 

Static Private Member Functions

static double get (const JModel_t &model, double JModel_t::*parameter)
 Read/write access to parameter value by data member. More...
 
static double & get (JModel_t &model, double JModel_t::*parameter)
 Read/write access to parameter value by data member. More...
 
static double get (const JModel_t &model, const size_t index)
 Read/write access to parameter value by index. More...
 
static double & get (JModel_t &model, const size_t index)
 Read/write access to parameter value by index. More...
 
static double get (const JModel_t &model, const int index)
 Read/write access to parameter value by index. More...
 
static double & get (JModel_t &model, const int index)
 Read/write access to parameter value by index. More...
 

Private Attributes

std::vector< double > h
 
JMATH::JVectorND x
 
struct {
   result_type   result
 
current
 
struct {
   JModel_t   value
 
   JModel_t   error
 
   result_type   result
 
previous
 

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.
If not defined, the parameters are assumed to be data members of type double.
Alternatively, the type definition can be size_t or int.
In that case, the model should provide for the element access operator[].
The first template parameter in the function operator should provide for an implementation of the actual fit function.
This function should return the data type JGandalf::result_type.
This data structure comprises the values of the chi2 and the gradient for a given data point.
The function operator returns the minimal chi2 and summed gradient of all data points.

Definition at line 84 of file JGandalf.hh.

Member Typedef Documentation

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

Data type of fit parameter.

Definition at line 95 of file JGandalf.hh.

Constructor & Destructor Documentation

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

Default constructor.

Definition at line 143 of file JGandalf.hh.

144  {}

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 multiple data sets.

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

Parameters
fitfit function
__beginbegin of data
__endend of data
argsoptional data
Returns
chi2 and gradient

Definition at line 160 of file JGandalf.hh.

161  {
162  using namespace std;
163  using namespace JPP;
164 
165  // note that all model values should be assigned to the start value of the model before use
166  // because the actual list of model parameters can vary from fit to fit
167  // (e.g. if model consists of a container).
168 
169  const size_t N = parameters.size();
170 
171  V.resize(N);
172  h.resize(N);
173  x.resize(N);
174 
175  previous.result.chi2 = numeric_limits<double>::max();
176 
177  current.result.chi2 = numeric_limits<double>::max();
178  current.result.gradient = value;
179  current.result.gradient = zero;
180 
181  error = value;
182  error = zero;
183 
184  lambda = LAMBDA_MIN;
185 
187 
188  DEBUG("step: " << numberOfIterations << endl);
189 
190  reset();
191 
192  update(fit, __begin, __end, args...);
193 
194  DEBUG("lambda: " << FIXED(12,5) << lambda << endl);
195  DEBUG("chi2: " << FIXED(12,5) << current.result.chi2 << endl);
196 
197  if (current.result.chi2 < previous.result.chi2) {
198 
199  if (numberOfIterations != 0) {
200 
201  if (fabs(previous.result.chi2 - current.result.chi2) < EPSILON*fabs(previous.result.chi2)) {
202 
203  const result_type result = current.result;
204 
205  reset();
206 
207  update(fit, __begin, __end, args...);
208 
209  try {
210  V.invert();
211  }
212  catch (const exception& error) {
213  V.reset();
214  }
215 
216  return result; // normal end
217  }
218 
219  if (lambda > LAMBDA_MIN) {
220  lambda /= LAMBDA_DOWN;
221  }
222  }
223  // store current values
224 
225  previous.value = value;
226  previous.error = error;
227  previous.result = current.result;
228 
229  } else {
230 
231  value = previous.value; // restore value
232 
233  lambda *= LAMBDA_UP;
234 
235  if (lambda > LAMBDA_MAX) {
236  break;
237  }
238 
239  reset();
240 
241  update(fit, __begin, __end, args...);
242  }
243 
244  DEBUG("Hesse matrix:" << endl << V << endl);
245 
246  // force definite positiveness
247 
248  for (size_t i = 0; i != N; ++i) {
249 
250  if (V(i,i) < PIVOT) {
251  V(i,i) = PIVOT;
252  }
253 
254  h[i] = 1.0 / sqrt(V(i,i));
255  }
256 
257  // normalisation
258 
259  for (size_t row = 0; row != N; ++row) {
260  for (size_t col = 0; col != row; ++col) {
261  V(row,col) *= h[row] * h[col];
262  V(col,row) = V(row,col);
263  }
264  }
265 
266  for (size_t i = 0; i != N; ++i) {
267  V(i,i) = 1.0 + lambda;
268  }
269 
270  // solve A x = b
271 
272  for (size_t col = 0; col != N; ++col) {
273  x[col] = h[col] * get(current.result.gradient, parameters[col]);
274  }
275 
276  try {
277  V.solve(x);
278  }
279  catch (const exception& error) {
280 
281  ERROR("JGandalf: " << error.what() << endl << V << endl);
282 
283  break;
284  }
285 
286  // update value and error
287 
288  for (size_t row = 0; row != N; ++row) {
289 
290  DEBUG("u[" << noshowpos << setw(3) << row << "] = " << showpos << FIXED(15,5) << get(value, parameters[row]));
291 
292  get(value, parameters[row]) -= h[row] * x[row];
293  get(error, parameters[row]) = h[row];
294 
295  DEBUG(" -> " << FIXED(15,5) << get(value, parameters[row]) << noshowpos << endl);
296  }
297  model(value);
298  }
299 
300  // abnormal end
301 
302  const result_type result = previous.result;
303 
304  value = previous.value; // restore value
305  error = previous.error; // restore error
306 
307  reset();
308 
309  update(fit, __begin, __end, args...);
310 
311  try {
312  V.invert();
313  }
314  catch (const exception& error) {
315  V.reset();
316  }
317 
318  return result;
319  }
double lambda
control parameter
Definition: JGandalf.hh:332
static double EPSILON
maximal distance to minimum
Definition: JGandalf.hh:323
JMATH::JVectorND x
Definition: JGandalf.hh:475
JModel_t value
value
Definition: JGandalf.hh:333
static double PIVOT
minimal value diagonal element of Hesse matrix
Definition: JGandalf.hh:328
JMATH::JMatrixNS V
Hesse matrix.
Definition: JGandalf.hh:335
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:459
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
std::vector< parameter_type > parameters
fit parameters
Definition: JGandalf.hh:330
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:446
static double LAMBDA_MIN
minimal value control parameter
Definition: JGandalf.hh:324
static double LAMBDA_DOWN
multiplication factor control parameter
Definition: JGandalf.hh:327
void reset()
Reset current parameters.
Definition: JGandalf.hh:341
#define ERROR(A)
Definition: JMessage.hh:66
void update(const JFunction_t &fit, T __begin, T __end, Args...args)
Recursive method to update current parameters.
Definition: JGandalf.hh:361
struct JFIT::JGandalf::@10 current
result_type result
Definition: JGandalf.hh:478
static double LAMBDA_UP
multiplication factor control parameter
Definition: JGandalf.hh:326
void invert()
Invert matrix according LDU decomposition.
Definition: JMatrixNS.hh:75
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
void solve(JVectorND_t &u)
Get solution of equation A x = b.
Definition: JMatrixNS.hh:308
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JGandalf.hh:322
int numberOfIterations
number of iterations
Definition: JGandalf.hh:331
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
Definition: JGandalf.hh:56
static double LAMBDA_MAX
maximal value control parameter
Definition: JGandalf.hh:325
JModel_t error
error
Definition: JGandalf.hh:334
struct JFIT::JGandalf::@11 previous
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std::vector< double > h
Definition: JGandalf.hh:474
template<class JModel_t>
void JFIT::JGandalf< JModel_t >::reset ( )
inlineprivate

Reset current parameters.

Definition at line 341 of file JGandalf.hh.

342  {
343  using namespace JPP;
344 
345  current.result.chi2 = 0.0;
346  current.result.gradient = zero;
347 
348  V.reset();
349  }
JMATH::JMatrixNS V
Hesse matrix.
Definition: JGandalf.hh:335
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:459
struct JFIT::JGandalf::@10 current
template<class JModel_t>
template<class JFunction_t , class T , class... Args>
void JFIT::JGandalf< JModel_t >::update ( const JFunction_t &  fit,
T  __begin,
T  __end,
Args...  args 
)
inlineprivate

Recursive method to update current parameters.

Parameters
fitfit function
__beginbegin of data
__endend of data
argsoptional data

Definition at line 361 of file JGandalf.hh.

362  {
363  for (T i = __begin; i != __end; ++i) {
364 
365  const result_type& result = fit(value, *i);
366 
367  current.result.chi2 += result.chi2;
368  current.result.gradient += result.gradient;
369 
370  for (size_t row = 0; row != parameters.size(); ++row) {
371  for (size_t col = row; col != parameters.size(); ++col) {
372  V(row,col) += get(result.gradient, parameters[row]) * get(result.gradient, parameters[col]);
373  }
374  }
375  }
376 
377  update(fit, args...);
378  }
JModel_t value
value
Definition: JGandalf.hh:333
JMATH::JMatrixNS V
Hesse matrix.
Definition: JGandalf.hh:335
std::vector< parameter_type > parameters
fit parameters
Definition: JGandalf.hh:330
do set_variable OUTPUT_DIRECTORY $WORKDIR T
void update(const JFunction_t &fit, T __begin, T __end, Args...args)
Recursive method to update current parameters.
Definition: JGandalf.hh:361
struct JFIT::JGandalf::@10 current
result_type result
Definition: JGandalf.hh:478
template<class JModel_t>
template<class JFunction_t >
void JFIT::JGandalf< JModel_t >::update ( const JFunction_t &  fit)
inlineprivate

Termination method to update current parameters.

Parameters
fitfit function

Definition at line 387 of file JGandalf.hh.

388  {
389  for (size_t row = 0; row != parameters.size(); ++row) {
390  for (size_t col = 0; col != row; ++col) {
391  V(row,col) = V(col,row);
392  }
393  }
394  }
JMATH::JMatrixNS V
Hesse matrix.
Definition: JGandalf.hh:335
std::vector< parameter_type > parameters
fit parameters
Definition: JGandalf.hh:330
template<class JModel_t>
static double JFIT::JGandalf< JModel_t >::get ( 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 404 of file JGandalf.hh.

405  {
406  return model.*parameter;
407  }
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
Definition: JGandalf.hh:56
template<class JModel_t>
static double& JFIT::JGandalf< JModel_t >::get ( 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 417 of file JGandalf.hh.

418  {
419  return model.*parameter;
420  }
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
Definition: JGandalf.hh:56
template<class JModel_t>
static double JFIT::JGandalf< JModel_t >::get ( 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 430 of file JGandalf.hh.

431  {
432  return model[index];
433  }
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
Definition: JGandalf.hh:56
template<class JModel_t>
static double& JFIT::JGandalf< JModel_t >::get ( JModel_t &  model,
const size_t  index 
)
inlinestaticprivate

Read/write access to parameter value by index.

Parameters
modelmodel
indexindex
Returns
value

Definition at line 443 of file JGandalf.hh.

444  {
445  return model[index];
446  }
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
Definition: JGandalf.hh:56
template<class JModel_t>
static double JFIT::JGandalf< JModel_t >::get ( const JModel_t &  model,
const int  index 
)
inlinestaticprivate

Read/write access to parameter value by index.

Parameters
modelmodel
indexindex
Returns
value

Definition at line 456 of file JGandalf.hh.

457  {
458  return model[index];
459  }
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
Definition: JGandalf.hh:56
template<class JModel_t>
static double& JFIT::JGandalf< JModel_t >::get ( JModel_t &  model,
const int  index 
)
inlinestaticprivate

Read/write access to parameter value by index.

Parameters
modelmodel
indexindex
Returns
value

Definition at line 469 of file JGandalf.hh.

470  {
471  return model[index];
472  }
void model(JModel_t &value)
Auxiliary function to constrain model during fit.
Definition: JGandalf.hh:56

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 322 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 323 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 324 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 325 of file JGandalf.hh.

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

multiplication factor control parameter

Definition at line 326 of file JGandalf.hh.

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

multiplication factor control parameter

Definition at line 327 of file JGandalf.hh.

template<class JModel_t>
double JFIT::JGandalf< JModel_t >::PIVOT = std::numeric_limits<double>::epsilon()
static

minimal value diagonal element of Hesse matrix

minimal value diagonal element of matrix

Definition at line 328 of file JGandalf.hh.

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

fit parameters

Definition at line 330 of file JGandalf.hh.

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

number of iterations

Definition at line 331 of file JGandalf.hh.

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

control parameter

Definition at line 332 of file JGandalf.hh.

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

value

Definition at line 333 of file JGandalf.hh.

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

error

Definition at line 334 of file JGandalf.hh.

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

Hesse matrix.

Definition at line 335 of file JGandalf.hh.

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

Definition at line 474 of file JGandalf.hh.

template<class JModel_t>
JMATH::JVectorND JFIT::JGandalf< JModel_t >::x
private

Definition at line 475 of file JGandalf.hh.

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

Definition at line 478 of file JGandalf.hh.

struct { ... } JFIT::JGandalf< JModel_t >::current
struct { ... } JFIT::JGandalf< JModel_t >::previous
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: