Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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< T >

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 = 0
 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.

Constructor & Destructor Documentation

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

Default constructor.

Definition at line 53 of file JSimplex.hh.

54  {}

Member Function Documentation

template<class JModel_t>
template<class JFunction_t , class T >
double JFIT::JSimplex< JModel_t >::operator() ( const JFunction_t &  fit,
T  __begin,
T  __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 69 of file JSimplex.hh.

70  {
71  using namespace std;
72  using namespace JPP;
73 
74  double chi2_old = evaluate(fit, __begin, __end);
75 
76  const int N = step.size();
77 
78  if (N != 0) {
79 
80  double chi2[N];
81 
82  numberOfIterations = 0;
83 
84  while (numberOfIterations != MAXIMUM_ITERATIONS) {
85 
86  DEBUG("old: " << FIXED(12,5) << chi2_old << endl);
87 
88  const JModel_t p0(value);
89 
90  for (int i = 0; i != N; ++i) {
91 
92  DEBUG("step: " << i << ' ' << setw(5) << numberOfIterations << endl);
93 
94  chi2[i] = (*this)(fit, __begin, __end, step[i]);
95  }
96 
97  // overall step direction of last iteration
98 
99  JModel_t wall = value - p0;
100 
101  const double chi2_new = (*this)(fit, __begin, __end, wall);
102 
103  DEBUG("new: " << FIXED(12,5) << chi2_new << endl);
104 
105 
106  if (fabs(chi2_old - chi2_new) < EPSILON*fabs(chi2_old)) {
107  return chi2_new;
108  }
109 
110  // double overall step
111 
112  wall = value - p0;
113 
114  value += wall;
115 
116  const double fe = evaluate(fit, __begin, __end);
117 
118  value -= wall;
119 
120 
121  for (int i = N-1; i != 0; --i) {
122  chi2[i] = chi2[i-1] - chi2[i];
123  }
124 
125  chi2[0] = chi2_old - chi2[0];
126 
127 
128  double df = 0.0;
129 
130  for (int i = 0; i != N; ++i) {
131  if (chi2[i] > df) {
132  df = chi2[i];
133  }
134  }
135 
136  const double fn = chi2_new;
137  const double f0 = chi2_old;
138  const double ff = f0 - fn - df;
139 
140  // shift step directions
141 
142  if (fe < f0 && 2.0*(f0 - 2.0*fn + fe)*ff*ff < (f0-fe)*(f0-fe)*df) {
143 
144  for (int i = 0; i != N - 1; ++i) {
145  step[i] = step[i+1];
146  }
147 
148  step[N-1] = wall;
149  }
150 
151  chi2_old = chi2_new;
152  }
153  }
154 
155  return chi2_old;
156  }
FIXED
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:461
JFIT::JSimplex::numberOfIterations
int numberOfIterations
Definition: JSimplex.hh:229
JFIT::JSimplex::value
JModel_t value
Definition: JSimplex.hh:227
JFIT::JSimplex::evaluate
double evaluate(const JFunction_t &fit, T __begin, T __end) const
Evaluate chi2 for given data set.
Definition: JSimplex.hh:240
JFIT::JSimplex::MAXIMUM_ITERATIONS
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JSimplex.hh:224
JFIT::JSimplex::EPSILON
static double EPSILON
maximal distance to minimum
Definition: JSimplex.hh:225
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
JFIT::JSimplex::step
std::vector< JModel_t > step
Definition: JSimplex.hh:228
template<class JModel_t>
template<class JFunction_t , class T >
double JFIT::JSimplex< JModel_t >::operator() ( const JFunction_t &  fit,
T  __begin,
T  __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 171 of file JSimplex.hh.

172  {
173  using namespace std;
174  using namespace JPP;
175 
176  double lambda = 0.5; // control parameter
177  double factor = 1.0; // multiplication factor
178 
179  double chi2_old = evaluate(fit, __begin, __end);
180 
181  for (int i = 0; numberOfIterations != MAXIMUM_ITERATIONS; ++numberOfIterations, ++i) {
182 
183  value += lambda * step;
184 
185  const double chi2_new = evaluate(fit, __begin, __end);
186 
187  DEBUG("step: " << setw(3) << i << ' ' << FIXED(12,5) << chi2_old << ' ' << FIXED(12,5) << chi2_new << ' ' << FIXED(5,2) << lambda << endl);
188 
189  if (fabs(chi2_old - chi2_new) < EPSILON*fabs(chi2_old)) {
190 
191  if (chi2_new > chi2_old) {
192 
193  value -= lambda * step; // undo last step
194 
195  return chi2_old;
196 
197  } else {
198 
199  return chi2_new;
200  }
201  }
202 
203  if (chi2_new < chi2_old) {
204 
205  chi2_old = chi2_new;
206 
207  } else {
208 
209  value -= lambda * step; // step back
210  lambda = -lambda; // change direction
211 
212  if (i != 0) {
213  factor = 0.5; // reduce step size
214  }
215  }
216 
217  lambda = factor * lambda;
218  }
219 
220  return chi2_old;
221  }
FIXED
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:461
JFIT::JSimplex::numberOfIterations
int numberOfIterations
Definition: JSimplex.hh:229
JFIT::JSimplex::value
JModel_t value
Definition: JSimplex.hh:227
JFIT::JSimplex::evaluate
double evaluate(const JFunction_t &fit, T __begin, T __end) const
Evaluate chi2 for given data set.
Definition: JSimplex.hh:240
JFIT::JSimplex::MAXIMUM_ITERATIONS
static int MAXIMUM_ITERATIONS
maximal number of iterations
Definition: JSimplex.hh:224
JFIT::JSimplex::EPSILON
static double EPSILON
maximal distance to minimum
Definition: JSimplex.hh:225
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:60
JFIT::JSimplex::step
std::vector< JModel_t > step
Definition: JSimplex.hh:228
template<class JModel_t>
template<class JFunction_t , class T >
double JFIT::JSimplex< JModel_t >::evaluate ( const JFunction_t &  fit,
T  __begin,
T  __end 
) const
inlineprivate

Evaluate chi2 for given data set.

Parameters
fitfit function
__beginbegin of data
__endend of data

Definition at line 240 of file JSimplex.hh.

241  {
242  double chi2 = 0.0;
243 
244  for (T hit = __begin; hit != __end; ++hit) {
245  chi2 += fit(value, *hit);
246  }
247 
248  return chi2;
249  }
JFIT::JSimplex::value
JModel_t value
Definition: JSimplex.hh:227

Member Data Documentation

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

maximal number of iterations

maximal number of iterations.

Definition at line 224 of file JSimplex.hh.

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 225 of file JSimplex.hh.

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

Definition at line 227 of file JSimplex.hh.

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

Definition at line 228 of file JSimplex.hh.

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

Definition at line 229 of file JSimplex.hh.

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

debug level (default is off).

Definition at line 43 of file JMessage.hh.

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