Jpp  18.0.0-rc.2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JGradient.hh
Go to the documentation of this file.
1 #ifndef __JFIT__JGRADIENT__
2 #define __JFIT__JGRADIENT__
3 
4 #include <limits>
5 #include <vector>
6 #include <cmath>
7 #include <memory>
8 #include <ostream>
9 #include <iomanip>
10 
11 #include "JLang/JManip.hh"
12 #include "JLang/JException.hh"
13 #include "Jeep/JMessage.hh"
14 
15 
16 /**
17  * \author mdejong
18  */
19 
20 namespace JFIT {}
21 namespace JPP { using namespace JFIT; }
22 
23 namespace JFIT {
24 
25  /**
26  * Auxiliary data structure for fit parameter.
27  */
28  struct JParameter_t {
29  /**
30  * Virtual destructor.
31  */
32  virtual ~JParameter_t()
33  {}
34 
35  /**
36  * Apply step.
37  *
38  * \param step step
39  */
40  virtual void apply(const double step) = 0;
41  };
42 
43 
44  /**
45  * Auxiliary data structure for editable parameter.
46  */
47  struct JEditor_t :
48  public std::shared_ptr<JParameter_t>
49  {
50  /**
51  * Constructor.
52  *
53  * \param name name
54  * \param parameter parameter
55  * \param value value
56  */
58  JParameter_t* parameter,
59  const double value) :
60  std::shared_ptr<JParameter_t>(parameter),
61  name (name),
62  value(value)
63  {}
64 
66  double value;
67  };
68 
69 
70  /**
71  * Conjugate gradient fit.
72  */
73  struct JGradient :
74  public std::vector<JEditor_t>
75  {
76  /**
77  * Constructor.
78  *
79  * \param Nmax maximum number of iterations
80  * \param Nextra maximum number of extra steps
81  * \param epsilon epsilon
82  * \param debug debug
83  */
84  JGradient(const size_t Nmax = std::numeric_limits<size_t>::max(),
85  const size_t Nextra = 0,
86  const double epsilon = 1.0e-4,
87  const double debug = 3) :
88  Nmax (Nmax),
89  Nextra (Nextra),
91  debug (debug)
92  {}
93 
94 
95  /**
96  * Fit.
97  *
98  * The template parameter should provide for the following syntax.
99  * <pre>
100  * double getChi2(int option);
101  * </pre>
102  * The value of the option corresponds to the following cases.
103  * - 0 => step wise improvement of the chi2;
104  * - 1 => evaluation of the chi2 before the determination of the gradient of the chi2; and
105  * - 2 => evaluation of the derivative of the chi2 to each fit parameter.
106  *
107  * \param getChi2 chi2 function
108  */
109  template<class T>
110  double operator()(const T& getChi2)
111  {
112  using namespace std;
113  using namespace JPP;
114 
115  if (this->empty()) {
116  return numeric_limits<double>::max();
117  }
118 
119  this->evaluate(getChi2);
120 
121  const size_t N = this->size();
122 
123  vector<double> H(N);
124  vector<double> G(N);
125 
126  for (size_t i = 0; i != N; ++i) {
127  G[i] = -1.0 * gradient[i];
128  H[i] = G[i];
129  gradient[i] = H[i];
130  }
131 
132  size_t number_of_iterations = 0;
133 
134  for ( ; number_of_iterations != Nmax; ++number_of_iterations) {
135 
136  DEBUG("chi2[0] " << setw(4) << number_of_iterations << ' ' << FIXED(12,5) << chi2[0] << endl);
137 
138  // minimise chi2 in direction of gradient
139 
140  chi2[1] = chi2[0];
141  chi2[2] = chi2[1];
142 
143  size_t m = 0;
144 
145  for (double ds = 1.0; ds > 1.0e-3; ) {
146 
147  this->move(+1.0 * ds);
148 
149  chi2[3] = getChi2(0);
150 
151  DEBUG("chi2[3] " << setw(4) << m << ' ' << FIXED(12,5) << chi2[3] << ' ' << FIXED(12,5) << ds << endl);
152 
153  if (chi2[3] < chi2[2]) {
154 
155  chi2[1] = chi2[2];
156  chi2[2] = chi2[3];
157 
158  m = 0;
159 
160  continue;
161  }
162 
163  if (ds == 1.0) {
164 
165  if (m == 0) {
166  chi2[4] = chi2[3];
167  }
168 
169  if (m != Nextra) {
170 
171  ++m;
172 
173  continue;
174 
175  } else {
176 
177  for ( ; m != 0; --m) {
178  this->move(-1.0 * ds);
179  }
180 
181  chi2[3] = chi2[4];
182  }
183  }
184 
185  this->move(-1.0 * ds);
186 
187  if (chi2[2] < chi2[3]) {
188 
189  // final step based on parabolic interpolation through following points
190  //
191  // x1 = -1 * ds -> chi2[1]
192  // x2 = 0 * ds -> chi2[2]
193  // x3 = +1 * ds -> chi2[3]
194 
195  const double f21 = chi2[2] - chi2[1]; // f(x2) - f(x1)
196  const double f23 = chi2[2] - chi2[3]; // f(x2) - f(x3)
197 
198  const double xs = 0.5 * (f21 - f23) / (f23 + f21);
199 
200  this->move(+1.0 * xs * ds);
201 
202  chi2[3] = getChi2(0);
203 
204  if (chi2[3] < chi2[2]) {
205 
206  chi2[2] = chi2[3];
207 
208  } else {
209 
210  this->move(-1.0 * xs * ds);
211 
212  chi2[2] = getChi2(0);
213  }
214 
215  DEBUG("chi2[2] " << setw(4) << number_of_iterations << ' ' << FIXED(12,5) << chi2[2] << ' ' << SCIENTIFIC(12,5) << ds << endl);
216 
217  break;
218 
219  } else {
220 
221  ds *= 0.5;
222  }
223  }
224 
225  if (fabs(chi2[2] - chi2[0]) < epsilon * 0.5 * (fabs(chi2[0]) + fabs(chi2[2]))) {
226 
227  chi2[0] = chi2[2];
228 
229  break;
230  }
231 
232  this->evaluate(getChi2);
233 
234  double gg = 0.0;
235  double dgg = 0.0;
236 
237  for (size_t i = 0; i != N; ++i){
238  gg += G[i]*G[i];
239  dgg += (gradient[i] + G[i]) * gradient[i];
240  }
241 
242  if (gg == 0.0) {
243  break;
244  }
245 
246  dgg /= gg;
247 
248  for (size_t i = 0; i != N; ++i){
249  G[i] = -1.0 * gradient[i];
250  H[i] = G[i] + dgg * H[i];
251  gradient[i] = H[i];
252  }
253  }
254 
255  DEBUG("chi2[0] " << setw(4) << number_of_iterations << ' ' << FIXED(12,5) << chi2[0] << endl);
256 
257  return chi2[0];
258  }
259 
260 
261  size_t Nmax; //!< maximum number of iterations
262  size_t Nextra; //!< maximum number of extra steps
263  double epsilon; //!< epsilon
264  int debug; //!< debug
265 
266  private:
267  /**
268  * Evaluate gradient.
269  */
270  template<class T>
271  void evaluate(const T& getChi2)
272  {
273  using namespace std;
274  using namespace JPP;
275 
276  const size_t N = this->size();
277 
278  gradient.resize(N);
279 
280  for (std::vector<double>::iterator i = gradient.begin(); i != gradient.end(); ++i) {
281  *i = 0.0;
282  }
283 
284  chi2[0] = getChi2(0);
285 
286  size_t width = 1;
287 
288  for (size_t i = 0; i != N; ++i) {
289  if ((*this)[i].name.size() > width) {
290  width = (*this)[i].name.size();
291  }
292  }
293 
294  double V = 0.0;
295 
296  for (size_t i = 0; i != N; ++i) {
297 
298  if ((*this)[i].value != 0.0) {
299 
300  (*this)[i]->apply(+0.5 * (*this)[i].value);
301 
302  chi2[1] = getChi2(2);
303 
304  (*this)[i]->apply(-0.5 * (*this)[i].value);
305  (*this)[i]->apply(-0.5 * (*this)[i].value);
306 
307  chi2[2] = getChi2(2);
308 
309  gradient[i] = chi2[1] - chi2[2];
310 
311  (*this)[i]->apply(+0.5 * (*this)[i].value);
312 
313  DEBUG(setw(width) << left << (*this)[i].name << right << ' ' << FIXED(12,5) << (*this)[i].value << ' ' << FIXED(12,5) << gradient[i] << endl);
314 
315  } else {
316 
317  gradient[i] = 0.0;
318  }
319 
320  V += gradient[i] * gradient[i];
321  }
322 
323  DEBUG(setw(width) << left << "|gradient|" << right << ' ' << FIXED(12,5) << sqrt(V) << endl);
324  }
325 
326 
327  /**
328  * Move.
329  *
330  * \param factor factor
331  */
332  void move(const double factor)
333  {
334  if (factor > 0.0) {
335  for (size_t i = 0; i != this->size(); ++i) {
336  (*this)[ i ]->apply((*this)[ i ].value * gradient[ i ] * factor);
337  }
338  } else if (factor < 0.0) {
339  for (size_t i = this->size(); i != 0; --i) {
340  (*this)[i-1]->apply((*this)[i-1].value * gradient[i-1] * factor);
341  }
342  }
343  }
344 
346  double chi2[5];
347  };
348 }
349 
350 #endif
void move(const double factor)
Move.
Definition: JGradient.hh:332
virtual ~JParameter_t()
Virtual destructor.
Definition: JGradient.hh:32
Exceptions.
JGradient(const size_t Nmax=std::numeric_limits< size_t >::max(), const size_t Nextra=0, const double epsilon=1.0e-4, const double debug=3)
Constructor.
Definition: JGradient.hh:84
then echo Enter input within $TIMEOUT_S seconds echo n User name
Definition: JCookie.sh:42
double operator()(const T &getChi2)
Fit.
Definition: JGradient.hh:110
static const double H
Planck constant [eV s].
void evaluate(const T &getChi2)
Evaluate gradient.
Definition: JGradient.hh:271
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
size_t Nextra
maximum number of extra steps
Definition: JGradient.hh:262
double epsilon
epsilon
Definition: JGradient.hh:263
do set_variable OUTPUT_DIRECTORY $WORKDIR T
then awk string
std::vector< double > gradient
Definition: JGradient.hh:345
virtual void apply(const double step)=0
Apply step.
size_t Nmax
maximum number of iterations
Definition: JGradient.hh:261
Auxiliary data structure for fit parameter.
Definition: JGradient.hh:28
std::string name
Definition: JGradient.hh:65
General purpose messaging.
I/O manipulators.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:36
int debug
debug
Definition: JGradient.hh:264
Auxiliary data structure for editable parameter.
Definition: JGradient.hh:47
double chi2[5]
Definition: JGradient.hh:346
double getChi2(const double P)
Get chi2 corresponding to given probability.
Definition: JFitToolkit.hh:56
Conjugate gradient fit.
Definition: JGradient.hh:73
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:484
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
JEditor_t(const std::string &name, JParameter_t *parameter, const double value)
Constructor.
Definition: JGradient.hh:57