Jpp  master_rocky-43-ge265d140c
the software that should make you happy
JQuantiles.hh
Go to the documentation of this file.
1 #ifndef __JTOOLS__JQUANTILES__
2 #define __JTOOLS__JQUANTILES__
3 
4 #include <limits>
5 #include <cmath>
6 
7 #include "JLang/JException.hh"
8 #include "JTools/JRange.hh"
9 #include "JTools/JElement.hh"
10 #include "JTools/JResult.hh"
12 #include "JTools/JToolsToolkit.hh"
13 
14 
15 /**
16  * \author mdejong
17  */
18 
19 namespace JTOOLS {}
20 namespace JPP { using namespace JTOOLS; }
21 
22 namespace JTOOLS {
23 
24  using JLANG::JException;
26 
27 
28  /**
29  * Locate maximum or minimun of function.
30  *
31  * Golden section search code is taken from reference:
32  * Numerical Recipes in C++, W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
33  * Cambridge University Press.
34  *
35  * xa < xb < xc
36  *
37  * is = +1 -> There is a minimum, i.e: f(xb) < min(f(xa),f(xc))
38  * is = -1 -> There is a maximum, i.e: f(xb) > max(f(xa),f(xc))
39  *
40  * \param xa
41  * \param xb
42  * \param xc
43  * \param f function
44  * \param is sign (+1 -> minimim, -1 -> maximum)
45  * \param eps relative precision
46  */
47  template<class JFunction1D_t>
48  double search(const double xa,
49  const double xb,
50  const double xc,
51  const JFunction1D_t& f,
52  const int is,
53  const double eps = 1.0e-6)
54  {
55  static const double R = 0.61803399;
56  static const double C = 1.0 - R;
57 
58  double x0 = xa;
59  double x3 = xc;
60  double x1, x2;
61 
62  if (fabs(xc-xb) > fabs(xb-xa)) {
63  x1 = xb;
64  x2 = xb + C*(xc-xb);
65  } else {
66  x2 = xb;
67  x1 = xb - C*(xb-xa);
68  }
69 
70  double f1 = is * get_value(f(x1));
71  double f2 = is * get_value(f(x2));
72 
73  while (fabs(x3-x0) > eps*(fabs(x1)+fabs(x2))) {
74 
75  if (f2 < f1) {
76 
77  x0 = x1;
78  x1 = x2;
79  x2 = R*x2 + C*x3;
80 
81  f1 = f2;
82  f2 = is * get_value(f(x2));
83 
84  } else {
85 
86  x3 = x2;
87  x2 = x1;
88  x1 = R*x1 + C*x0;
89 
90  f2 = f1;
91  f1 = is * get_value(f(x1));
92  }
93  }
94 
95  if (f1 < f2)
96  return x1;
97  else
98  return x2;
99  }
100 
101 
102  /**
103  * Quantile calculator for a given function.
104  * It is assumed that the function has a single maximum.
105  */
106  class JQuantiles :
107  public JRange<double>
108  {
109  public:
110 
112 
113  /**
114  * Default constructor.
115  */
117  Xmax(0.0),
118  Ymax(0.0),
119  fwhm(0.0),
120  sum (0.0)
121  {}
122 
123 
124  /**
125  * Constructor.
126  *
127  * \param f1 functional collection
128  * \param Q quantile
129  * \param eps relative precision
130  */
131  template<class JFunction1D_t>
132  JQuantiles(const JFunction1D_t& f1,
133  const double Q = 1.0,
134  const double eps = 1.0e-6) :
135  Xmax(0.0),
136  Ymax(0.0),
137  fwhm(0.0),
138  sum (0.0)
139  {
140  set(f1, Q, eps);
141  }
142 
143 
144  /**
145  * Constructor.
146  *
147  * \param abscissa abscissa
148  * \param f1 function
149  * \param Q quantile
150  * \param eps relative precision
151  */
152  template<class JFunction1D_t>
153  JQuantiles(const JAbscissa_t& abscissa,
154  const JFunction1D_t& f1,
155  const double Q = 1.0,
156  const double eps = 1.0e-6) :
157  Xmax(0.0),
158  Ymax(0.0),
159  fwhm(0.0),
160  sum (0.0)
161  {
162  set(abscissa, f1, Q, eps);
163  }
164 
165 
166  /**
167  * Set quantiles.
168  *
169  * \param f1 functional collection
170  * \param Q quantile
171  * \param eps relative precision
172  */
173  template<class JFunction1D_t>
174  void set(const JFunction1D_t& f1,
175  const double Q = 1.0,
176  const double eps = 1.0e-6)
177  {
178  typedef typename JFunction1D_t::const_iterator const_iterator;
179 
180  if (f1.empty()) {
181  throw JEmptyCollection("JQuantiles() no data.");
182  }
183 
184 
185  // maximum
186 
187  const_iterator p = f1.begin();
188 
189  for (const_iterator i = f1.begin(); i != f1.end(); ++i) {
190  if (i->getY() > p->getY()) {
191  p = i;
192  }
193  }
194 
195 
196  // x at maximum
197 
198  Xmax = p->getX();
199 
200  if (p != f1.begin()) {
201 
202  const double xa = (--p)->getX();
203  const double xb = (++p)->getX();
204 
205  if (++p != f1.end()) {
206 
207  const double xc = p->getX();
208 
209  Xmax = search(xa, xb, xc, f1, -1, eps);
210  }
211  }
212 
213  Ymax = get_value(f1(Xmax));
214 
215 
216  // integral & quantile
217 
218  if (Q > 0.0 && Q <= 1.0) {
219 
221 
222  try {
223 
224  sum = makeCDF(f1, buffer);
225 
226  setLowerLimit(buffer(0.5 * (1.0 - Q)));
227  setUpperLimit(buffer(0.5 * (1.0 + Q)));
228  }
229  catch(const JException& error) {
230  sum = 0.0;
231  }
232 
233  } else {
234 
236 
237  if (Q > 1.0) {
238  setLowerLimit(f1. begin()->getX());
239  setUpperLimit(f1.rbegin()->getX());
240  } else if (Q <= 0.0) {
243  }
244  }
245 
246 
247  // FWHM
248 
249  fwhm = 0.0;
250 
251  for (double xmin = f1.begin()->getX(), xmax = Xmax, v = 0.5*Ymax; ; ) {
252 
253  const double x = 0.5 * (xmin + xmax);
254  const double y = get_value(f1(x));
255 
256  if (fabs(y - v) < eps*v || xmax - xmin < eps) {
257  fwhm -= x;
258  break;
259  }
260 
261  if (y > v)
262  xmax = x;
263  else
264  xmin = x;
265  }
266 
267  for (double xmin = Xmax, xmax = f1.rbegin()->getX(), v = 0.5*Ymax; ; ) {
268 
269  const double x = 0.5 * (xmin + xmax);
270  const double y = get_value(f1(x));
271 
272  if (fabs(y - v) < eps*v || xmax - xmin < eps) {
273  fwhm += x;
274  break;
275  }
276 
277  if (y > v)
278  xmin = x;
279  else
280  xmax = x;
281  }
282  }
283 
284 
285  /**
286  * Set quantiles.
287  *
288  * \param abscissa abscissa
289  * \param f1 function
290  * \param Q quantile
291  * \param eps relative precision
292  */
293  template<class JFunction1D_t>
294  void set(const JAbscissa_t& abscissa,
295  const JFunction1D_t& f1,
296  const double Q = 1.0,
297  const double eps = 1.0e-6)
298  {
300 
301  buffer.configure(abscissa, f1);
302  buffer.compile();
303 
304  set(buffer, Q, eps);
305  }
306 
307 
308  /**
309  * Get position of maximum.
310  *
311  * \return x value at maximum
312  */
313  double getX() const
314  {
315  return Xmax;
316  }
317 
318 
319  /**
320  * Get value of maximum.
321  *
322  * \return y value at maximum
323  */
324  double getY() const
325  {
326  return Ymax;
327  }
328 
329 
330  /**
331  * Get Full Width at Half Maximum.
332  *
333  * \return FWHM
334  */
335  double getFWHM() const
336  {
337  return fwhm;
338  }
339 
340 
341  /**
342  * Get integral of function.
343  *
344  * \return integral
345  */
346  double getIntegral() const
347  {
348  return sum;
349  }
350 
351 
352  protected:
353  double Xmax;
354  double Ymax;
355  double fwhm;
356  double sum;
357  };
358 }
359 
360 #endif
The elements in a collection are sorted according to their abscissa values and a given distance opera...
Exceptions.
Auxiliary class to define a range between two values.
This include file containes various data structures that can be used as specific return types for the...
This include file contains various recursive methods to operate on multi-dimensional collections.
Exception for an empty collection.
Definition: JException.hh:162
General exception.
Definition: JException.hh:24
General purpose class for collection of elements, see: <a href="JTools.PDF";>Collection of elements...
Definition: JCollection.hh:79
Quantile calculator for a given function.
Definition: JQuantiles.hh:108
double getIntegral() const
Get integral of function.
Definition: JQuantiles.hh:346
void set(const JFunction1D_t &f1, const double Q=1.0, const double eps=1.0e-6)
Set quantiles.
Definition: JQuantiles.hh:174
double getY() const
Get value of maximum.
Definition: JQuantiles.hh:324
void set(const JAbscissa_t &abscissa, const JFunction1D_t &f1, const double Q=1.0, const double eps=1.0e-6)
Set quantiles.
Definition: JQuantiles.hh:294
JQuantiles(const JAbscissa_t &abscissa, const JFunction1D_t &f1, const double Q=1.0, const double eps=1.0e-6)
Constructor.
Definition: JQuantiles.hh:153
JQuantiles()
Default constructor.
Definition: JQuantiles.hh:116
double getX() const
Get position of maximum.
Definition: JQuantiles.hh:313
JAbstractCollection< double > JAbscissa_t
Definition: JQuantiles.hh:111
JQuantiles(const JFunction1D_t &f1, const double Q=1.0, const double eps=1.0e-6)
Constructor.
Definition: JQuantiles.hh:132
double getFWHM() const
Get Full Width at Half Maximum.
Definition: JQuantiles.hh:335
Range of values.
Definition: JRange.hh:42
void setUpperLimit(argument_type y)
Set upper limit.
Definition: JRange.hh:235
void setLowerLimit(argument_type x)
Set lower limit.
Definition: JRange.hh:224
Template class for spline interpolation in 1D.
Definition: JSpline.hh:738
const JPolynome f1(1.0, 2.0, 3.0)
Function.
const double xmax
Definition: JQuadrature.cc:24
const double xmin
Definition: JQuadrature.cc:23
static const double C
Physics constants.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary classes and methods for multi-dimensional interpolations and histograms.
double search(const double xa, const double xb, const double xc, const JFunction1D_t &f, const int is, const double eps=1.0e-6)
Locate maximum or minimun of function.
Definition: JQuantiles.hh:48
JContainer_t::ordinate_type makeCDF(const JContainer_t &input, JMappableCollection< JKey_t, JValue_t > &output, const typename JContainer_t::ordinate_type eps=JMATH::zero)
Conversion of data points to cumulative probability distribition (CDF).
JContainer_t::ordinate_type getIntegral(const JContainer_t &input)
Get integral of input data points.
JResultEvaluator< JResult_t >::result_type get_value(const JResult_t &value)
Helper method to recursively evaluate a to function value.
Definition: JResult.hh:998
data_type v[N+1][M+1]
Definition: JPolint.hh:866
Abstract interface for abscissa values of a collection of elements.