Jpp
JSimplex.hh
Go to the documentation of this file.
1 #ifndef __JFIT__JSIMPLEX__
2 #define __JFIT__JSIMPLEX__
3 
4 #include <cmath>
5 #include <vector>
6 #include <ostream>
7 #include <iomanip>
8 
9 #include "Jeep/JPrint.hh"
10 #include "Jeep/JMessage.hh"
11 
12 
13 /**
14  * \author mdejong
15  */
16 
17 namespace JFIT {}
18 namespace JPP { using namespace JFIT; }
19 
20 namespace JFIT {
21 
22  using JEEP::JMessage;
23 
24 
25  /**
26  * Simple fit method based on Powell's algorithm, see reference:
27  * Numerical Recipes in C++, W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
28  * Cambridge University Press.
29  *
30  * The template argument refers to the model that should be fitted to the data.
31  * This data structure should have arithmetic capabalities.
32  *
33  * The data member JSimplex::value corresponds to the start, current or final value
34  * of the model of the fit procedure.
35  * The data member JSimplex::step corresponds to the step directions.
36  * Note that the step directions may change during the fit.
37  * The template fit function should return the data type JGandalf::result_type
38  * which is composed of the values of the chi2 and gradient of a hit.
39  * The function operator returns the chi2 of the fit.
40  */
41  template<class JModel_t>
42  class JSimplex :
43  public JMessage< JSimplex<JModel_t> >
44  {
45  public:
46 
47  typedef double result_type;
48 
50 
51 
52  /**
53  * Default constructor.
54  */
56  {}
57 
58 
59  /**
60  * Multi-dimensional fit.
61  *
62  * The given fit function should return the equivalent of chi2 for
63  * the current value of the given model and a given data point.
64  *
65  * \param fit fit function
66  * \param __begin begin of data
67  * \param __end end of data
68  * \return chi2
69  */
70  template<class JFunction_t, class T>
71  double operator()(const JFunction_t& fit, T __begin, T __end)
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  }
159 
160 
161  /**
162  * 1D fit.
163  *
164  * The given fit function should return the equivalent of chi2 for
165  * the current value of the given model and a given data point.
166  *
167  * \param fit fit function
168  * \param __begin begin of data
169  * \param __end end of data
170  * \param step step direction
171  */
172  template<class JFunction_t, class T>
173  double operator()(const JFunction_t& fit, T __begin, T __end, const JModel_t& step)
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  }
224 
225 
226  static int MAXIMUM_ITERATIONS; //!< maximal number of iterations
227  static double EPSILON; //!< maximal distance to minimum
228 
229  JModel_t value;
232 
233  private:
234  /**
235  * Evaluate chi2 for given data set.
236  *
237  * \param fit fit function
238  * \param __begin begin of data
239  * \param __end end of data
240  */
241  template<class JFunction_t, class T>
242  inline double evaluate(const JFunction_t& fit, T __begin, T __end) const
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  }
252  };
253 
254 
255  /**
256  * maximal number of iterations.
257  */
258  template<class JModel_t>
260 
261 
262  /**
263  * maximal distance to minimum.
264  */
265  template<class JModel_t>
266  double JSimplex<JModel_t>::EPSILON = 1.0e-4;
267 }
268 
269 #endif
FIXED
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:481
JFIT
Auxiliary classes and methods for linear and iterative data regression.
Definition: JEnergy.hh:15
JMessage.hh
JPrint.hh
JFIT::JSimplex
Simple fit method based on Powell's algorithm, see reference: Numerical Recipes in C++,...
Definition: JSimplex.hh:42
std::vector< JModel_t >
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JFIT::JSimplex::operator()
double operator()(const JFunction_t &fit, T __begin, T __end, const JModel_t &step)
1D fit.
Definition: JSimplex.hh:173
JFIT::JSimplex::value
JModel_t value
Definition: JSimplex.hh:229
JFIT::JSimplex::operator()
double operator()(const JFunction_t &fit, T __begin, T __end)
Multi-dimensional fit.
Definition: JSimplex.hh:71
JFIT::JSimplex::result_type
double result_type
Definition: JSimplex.hh:47
JFIT::JSimplex::evaluate
double evaluate(const JFunction_t &fit, T __begin, T __end) const
Evaluate chi2 for given data set.
Definition: JSimplex.hh:242
JFIT::JSimplex::JSimplex
JSimplex()
Default constructor.
Definition: JSimplex.hh:55
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
std
Definition: jaanetDictionary.h:36
JEEP::JMessage< JSimplex< JModel_t > >::debug
static int debug
debug level (default is off).
Definition: JMessage.hh:45
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
JEEP::JMessage
Auxiliary class for handling debug parameter within a class.
Definition: JMessage.hh:44
JFIT::JSimplex::numberOfIterations
int numberOfIterations
Definition: JSimplex.hh:231