Jpp
Public Types | Public Member Functions | Public Attributes | Static Public Attributes | Private Member Functions | List of all members
JFIT::JSimplex< JModel_t > Class Template Reference

Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++, W.H. More...

#include <JSimplex.hh>

Inheritance diagram for JFIT::JSimplex< JModel_t >:
JEEP::JMessage< JSimplex< JModel_t > >

Public Types

typedef double result_type
 

Public Member Functions

 JSimplex ()
 Default constructor. More...
 
template<class JFunction_t , class T >
double operator() (const JFunction_t &fit, T __begin, T __end)
 Multi-dimensional fit. More...
 
template<class JFunction_t , class T >
double operator() (const JFunction_t &fit, T __begin, T __end, const JModel_t &step)
 1D fit. More...
 

Public Attributes

JModel_t value
 
std::vector< JModel_t > step
 
int numberOfIterations
 

Static Public Attributes

static int MAXIMUM_ITERATIONS = 1000
 maximal number of iterations More...
 
static double EPSILON = 1.0e-4
 maximal distance to minimum More...
 
static int debug
 debug level (default is off). More...
 

Private Member Functions

template<class JFunction_t , class T >
double evaluate (const JFunction_t &fit, T __begin, T __end) const
 Evaluate chi2 for given data set. More...
 

Detailed Description

template<class JModel_t>
class JFIT::JSimplex< JModel_t >

Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++, W.H.

Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Cambridge University Press.

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

The data member JSimplex::value corresponds to the start, current or final value of the model of the fit procedure. The data member JSimplex::step corresponds to the step directions. Note that the step directions may change during the fit. The template fit function should return the data type JGandalf::result_type which is composed of the values of the chi2 and gradient of a hit. The function operator returns the chi2 of the fit.

Definition at line 42 of file JSimplex.hh.

Member Typedef Documentation

◆ result_type

template<class JModel_t>
typedef double JFIT::JSimplex< JModel_t >::result_type

Definition at line 47 of file JSimplex.hh.

Constructor & Destructor Documentation

◆ JSimplex()

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

Default constructor.

Definition at line 55 of file JSimplex.hh.

56  {}

Member Function Documentation

◆ operator()() [1/2]

template<class JModel_t>
template<class JFunction_t , class T >
double JFIT::JSimplex< JModel_t >::operator() ( const JFunction_t &  fit,
__begin,
__end 
)
inline

Multi-dimensional fit.

The given fit function should return the equivalent of chi2 for the current value of the given model and a given data point.

Parameters
fitfit function
__beginbegin of data
__endend of data
Returns
chi2

Definition at line 71 of file JSimplex.hh.

72  {
73  using namespace std;
74  using namespace JPP;
75 
76  double chi2_old = evaluate(fit, __begin, __end);
77 
78  const int N = step.size();
79 
80  if (N != 0) {
81 
82  double chi2[N];
83 
85 
87 
88  DEBUG("old: " << FIXED(12,5) << chi2_old << endl);
89 
90  const JModel_t p0(value);
91 
92  for (int i = 0; i != N; ++i) {
93 
94  DEBUG("step: " << i << ' ' << setw(5) << numberOfIterations << endl);
95 
96  chi2[i] = (*this)(fit, __begin, __end, step[i]);
97  }
98 
99  // overall step direction of last iteration
100 
101  JModel_t wall = value - p0;
102 
103  const double chi2_new = (*this)(fit, __begin, __end, wall);
104 
105  DEBUG("new: " << FIXED(12,5) << chi2_new << endl);
106 
107 
108  if (fabs(chi2_old - chi2_new) < EPSILON*fabs(chi2_old)) {
109  return chi2_new;
110  }
111 
112  // double overall step
113 
114  wall = value - p0;
115 
116  value += wall;
117 
118  const double fe = evaluate(fit, __begin, __end);
119 
120  value -= wall;
121 
122 
123  for (int i = N-1; i != 0; --i) {
124  chi2[i] = chi2[i-1] - chi2[i];
125  }
126 
127  chi2[0] = chi2_old - chi2[0];
128 
129 
130  double df = 0.0;
131 
132  for (int i = 0; i != N; ++i) {
133  if (chi2[i] > df) {
134  df = chi2[i];
135  }
136  }
137 
138  const double fn = chi2_new;
139  const double f0 = chi2_old;
140  const double ff = f0 - fn - df;
141 
142  // shift step directions
143 
144  if (fe < f0 && 2.0*(f0 - 2.0*fn + fe)*ff*ff < (f0-fe)*(f0-fe)*df) {
145 
146  for (int i = 0; i != N - 1; ++i) {
147  step[i] = step[i+1];
148  }
149 
150  step[N-1] = wall;
151  }
152 
153  chi2_old = chi2_new;
154  }
155  }
156 
157  return chi2_old;
158  }

◆ operator()() [2/2]

template<class JModel_t>
template<class JFunction_t , class T >
double JFIT::JSimplex< JModel_t >::operator() ( const JFunction_t &  fit,
__begin,
__end,
const JModel_t &  step 
)
inline

1D fit.

The given fit function should return the equivalent of chi2 for the current value of the given model and a given data point.

Parameters
fitfit function
__beginbegin of data
__endend of data
stepstep direction

Definition at line 173 of file JSimplex.hh.

174  {
175  using namespace std;
176  using namespace JPP;
177 
178  double lambda = 0.5; // control parameter
179  double factor = 1.0; // multiplication factor
180 
181  double chi2_old = evaluate(fit, __begin, __end);
182 
183  for (int i = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations, ++i) {
184 
185  value += lambda * step;
186 
187  const double chi2_new = evaluate(fit, __begin, __end);
188 
189  DEBUG("step: " << setw(3) << i << ' ' << FIXED(12,5) << chi2_old << ' ' << FIXED(12,5) << chi2_new << ' ' << FIXED(5,2) << lambda << endl);
190 
191  if (fabs(chi2_old - chi2_new) < EPSILON*fabs(chi2_old)) {
192 
193  if (chi2_new > chi2_old) {
194 
195  value -= lambda * step; // undo last step
196 
197  return chi2_old;
198 
199  } else {
200 
201  return chi2_new;
202  }
203  }
204 
205  if (chi2_new < chi2_old) {
206 
207  chi2_old = chi2_new;
208 
209  } else {
210 
211  value -= lambda * step; // step back
212  lambda = -lambda; // change direction
213 
214  if (i != 0) {
215  factor = 0.5; // reduce step size
216  }
217  }
218 
219  lambda = factor * lambda;
220  }
221 
222  return chi2_old;
223  }

◆ evaluate()

template<class JModel_t>
template<class JFunction_t , class T >
double JFIT::JSimplex< JModel_t >::evaluate ( const JFunction_t &  fit,
__begin,
__end 
) const
inlineprivate

Evaluate chi2 for given data set.

Parameters
fitfit function
__beginbegin of data
__endend of data

Definition at line 242 of file JSimplex.hh.

243  {
244  double chi2 = 0.0;
245 
246  for (T hit = __begin; hit != __end; ++hit) {
247  chi2 += fit(value, *hit);
248  }
249 
250  return chi2;
251  }

Member Data Documentation

◆ MAXIMUM_ITERATIONS

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

maximal number of iterations

maximal number of iterations.

Definition at line 226 of file JSimplex.hh.

◆ EPSILON

template<class JModel_t>
double JFIT::JSimplex< JModel_t >::EPSILON = 1.0e-4
static

maximal distance to minimum

maximal distance to minimum.

Definition at line 227 of file JSimplex.hh.

◆ value

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

Definition at line 229 of file JSimplex.hh.

◆ step

template<class JModel_t>
std::vector<JModel_t> JFIT::JSimplex< JModel_t >::step

Definition at line 230 of file JSimplex.hh.

◆ numberOfIterations

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

Definition at line 231 of file JSimplex.hh.

◆ debug

int JEEP::JMessage< JSimplex< JModel_t > >::debug
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:
FIXED
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:481
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JFIT::JSimplex::value
JModel_t value
Definition: JSimplex.hh:229
JFIT::JSimplex::evaluate
double evaluate(const JFunction_t &fit, T __begin, T __end) const
Evaluate chi2 for given data set.
Definition: JSimplex.hh:242
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std
Definition: jaanetDictionary.h:36
JFIT::JSimplex::step
std::vector< JModel_t > step
Definition: JSimplex.hh:230
JFIT::JSimplex::EPSILON
static double EPSILON
maximal distance to minimum
Definition: JSimplex.hh:227
JFIT::JSimplex::MAXIMUM_ITERATIONS
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JSimplex.hh:226
JFIT::JSimplex::numberOfIterations
int numberOfIterations
Definition: JSimplex.hh:231