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

#include <JGandalf.hh>

Inheritance diagram for JFIT::JGandalf< JModel_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 multiple data sets. More...

Public Attributes
std::vector< parameter_type >	parameters
	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 = 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 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

struct {
result_type result

JMATH::JMatrixNS V

}	current

struct {
JModel_t value

JModel_t error

result_type result

JMATH::JMatrixNS V

}	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.
Alternatively, the type definition can be size_t or int.
In that case, the model class 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 type comprises the values of the chi2 and the gradient for a given data point.
The function operator returns the minimal chi2 and combined gradient of the data points.

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 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

fit	fit function
__begin	begin of data
__end	end of data
args	optional data

Returns: chi2

Definition at line 128 of file JGandalf.hh.

     {
       using namespace std;
       using namespace JPP;
 
       // note that all model values should be assigned to the start value of the model before use
       // because the actual list of model parameters can vary from fit to fit
       // (e.g. if model consists of a container).
 
       const size_t N = parameters.size();
 
       V.resize(N);
       h.resize(N);
 
       previous.result.chi2 = numeric_limits<double>::max();
 
       previous.V.resize(N);
 
       current.result.chi2     = numeric_limits<double>::max();
       current.result.gradient = value;
       current.result.gradient = zero;
 
       current.V.resize(N);
 
       error     = value;
       error     = zero;
 
       lambda    = LAMBDA_MIN;
 
 
       for (numberOfIterations = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations) {
 
         DEBUG("step:     " << numberOfIterations << endl);
 
         reset();
 
         update(fit, __begin, __end, args...);
 
         DEBUG("lambda:   " << FIXED(12,5) << lambda              << endl);
         DEBUG("chi2:     " << FIXED(12,5) << current.result.chi2 << endl);
 
         if (current.result.chi2 < previous.result.chi2) {
 
           if (numberOfIterations != 0) {
 
             if (fabs(previous.result.chi2 - current.result.chi2) < EPSILON*fabs(previous.result.chi2)) {
               return current.result;                                // normal end
             }
 
             if (lambda > LAMBDA_MIN) {
               lambda /= LAMBDA_DOWN;
             }
 
             // store current values
 
             previous.value  = value;
             previous.error  = error;
             previous.result = current.result;
             
             previous.V.swap(V);
           }
 
         } else {
 
           value   = previous.value;                                 // restore value
 
           lambda *= LAMBDA_UP;
           
           if (lambda > LAMBDA_MAX) {
             break;
           }
 
           reset();
 
           update(fit, __begin, __end, args...);
         }
 
 
         V.swap(current.V);
 
         DEBUG("Hesse matrix:" << endl << V << endl);
 
 
         // force definite positiveness
 
         for (size_t i = 0; i != N; ++i) {
 
           if (V(i,i) < PIVOT) {
             V(i,i) = PIVOT;
           }
 
           h[i] = 1.0 / sqrt(V(i,i));
         }
 
 
         // normalisation
         
         for (size_t row = 0; row != N; ++row) {
           for (size_t col = 0; col != row; ++col) {
             V(row,col) *= h[row] * h[col];
             V(col,row)  = V(row,col);
           }
         }
         
         for (size_t i = 0; i != N; ++i) {
           V(i,i) = 1.0 + lambda;
         }
 
 
         // invert Hesse matrix
 
         try {
           V.invert();
         }
         catch (const JException& error) {
 
           ERROR("JGandalf: " << error.what() << endl << V << endl);
 
           break;
         }
 
 
         // update value and error       
         
         for (size_t row = 0; row != N; ++row) {
 
           DEBUG("u[" << noshowpos << row << "] = " << showpos << FIXED(15,5) << get(value, parameters[row]));
 
           for (size_t col = 0; col != N; ++col) {
 
             V(row,col) *= h[row] * h[col];
 
             get(value, parameters[row]) -= V(row,col) * get(current.result.gradient, parameters[col]);
           }
           
           DEBUG(" -> " << FIXED(15,5) << get(value, parameters[row]) << noshowpos << endl);
 
           get(error, parameters[row]) = h[row];
         }
       }
 
       // abnormal end
 
       value   = previous.value;                                     // restore value
       error   = previous.error;                                     // restore error
 
       V.swap(previous.V);                                           // restore Hesse matrix
 
       return previous.result;
     }

template<class JModel_t>

void JFIT::JGandalf< JModel_t >::reset ( )

inlineprivate

Reset current parameters.

Definition at line 299 of file JGandalf.hh.

     {
       using namespace JPP;
 
       current.result.chi2     = 0.0;
       current.result.gradient = zero;
 
       current.V.reset();          
     }

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

fit	fit function
__begin	begin of data
__end	end of data
args	optional data

Definition at line 319 of file JGandalf.hh.

     {
       for (T i = __begin; i != __end; ++i) {
 
         const result_type& result = fit(value, *i);
 
         current.result.chi2     += result.chi2;
         current.result.gradient += result.gradient;
 
         for (size_t row = 0; row != parameters.size(); ++row) {
           for (size_t col = row; col != parameters.size(); ++col) {
             current.V(row,col) += get(result.gradient, parameters[row]) * get(result.gradient, parameters[col]);
           }
         }
       }
 
       update(fit, args...);
     }

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

fit	fit function

Definition at line 345 of file JGandalf.hh.

     {
       for (size_t row = 0; row != parameters.size(); ++row) {
         for (size_t col = 0; col != row; ++col) {
           current.V(row,col) = current.V(col,row);
         }
       }
     }

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

model	model
parameter	parameter

Returns: value

Definition at line 362 of file JGandalf.hh.

     {
       return model.*parameter;
     }

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

model	model
parameter	parameter

Returns: value

Definition at line 375 of file JGandalf.hh.

     {
       return model.*parameter;
     }

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

model	model
index	index

Returns: value

Definition at line 388 of file JGandalf.hh.

     {
       return model[index];
     }

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

model	model
index	index

Returns: value

Definition at line 401 of file JGandalf.hh.

     {
       return model[index];
     }

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

model	model
index	index

Returns: value

Definition at line 414 of file JGandalf.hh.

     {
       return model[index];
     }

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

model	model
index	index

Returns: value

Definition at line 427 of file JGandalf.hh.

     {
       return model[index];
     }

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 280 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 281 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 282 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 283 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 284 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 285 of file JGandalf.hh.

template<class JModel_t>

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

static

minimal value diagonal element of Hesse matrix

minimal value diagonal element of matrix

Definition at line 286 of file JGandalf.hh.

template<class JModel_t>

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

fit parameters

Definition at line 288 of file JGandalf.hh.

template<class JModel_t>

int JFIT::JGandalf< JModel_t >::numberOfIterations

number of iterations

Definition at line 289 of file JGandalf.hh.

template<class JModel_t>

double JFIT::JGandalf< JModel_t >::lambda

control parameter

Definition at line 290 of file JGandalf.hh.

template<class JModel_t>

JModel_t JFIT::JGandalf< JModel_t >::value

value

Definition at line 291 of file JGandalf.hh.

template<class JModel_t>

JModel_t JFIT::JGandalf< JModel_t >::error

error

Definition at line 292 of file JGandalf.hh.

template<class JModel_t>

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

Hesse matrix.

Definition at line 293 of file JGandalf.hh.

template<class JModel_t>

std::vector<double> JFIT::JGandalf< JModel_t >::h

private

Definition at line 432 of file JGandalf.hh.

template<class JModel_t>

result_type JFIT::JGandalf< JModel_t >::result

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

JGandalf.hh

Classes

Public Types

Public Member Functions

Public Attributes

Static Public Attributes

Private Member Functions

Static Private Member Functions

Private Attributes

Detailed Description

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

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation

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