Jpp - the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JQuantile.hh
Go to the documentation of this file.
1 #ifndef __JTOOLS__JQUANTILE__
2 #define __JTOOLS__JQUANTILE__
3 
4 #include <ostream>
5 #include <iomanip>
6 #include <cmath>
7 #include <limits>
8 #include <algorithm>
9 #include <map>
10 
11 #include "JLang/JException.hh"
12 #include "JLang/JTitle.hh"
13 #include "JLang/JVectorize.hh"
14 #include "JLang/JManip.hh"
15 
16 #include "JMath/JMath.hh"
17 
18 
19 /**
20  * \author mdejong
21  */
22 
23 namespace JTOOLS {}
24 namespace JPP { using namespace JTOOLS; }
25 
26 namespace JTOOLS {
27 
28  using JMATH::JMath;
29  using JLANG::JTitle;
31  using JLANG::JNoValue;
32  using JLANG::array_type;
33  using JLANG::make_array;
34 
35 
36  /**
37  * Auxiliary data structure for running average, standard deviation and quantiles.
38  *
39  * See Knuth TAOCP vol 2, 3rd edition, page 232.\n
40  * This class acts as a zero-dimensional histogram.\n
41  * Note that if a weight is used, it should strictly be positive.
42  */
43  struct JQuantile :
44  public JTitle,
45  public JMath<JQuantile>
46  {
47  /**
48  * Constructor.
49  *
50  * \param title title
51  * \param quantiles quantiles
52  */
53  JQuantile(const JTitle& title = "",
54  const bool quantiles = false) :
55  JTitle(title),
57  {
58  reset();
59  }
60 
61 
62  /**
63  * Constructor.
64  *
65  * \param title title
66  * \param buffer input data
67  * \param quantiles quantiles
68  * \param w weight
69  */
70  template<class JElement_t, class JAllocator_t>
73  const bool quantiles = false,
74  const double w = 1.0) :
75  JTitle(title),
77  {
78  reset();
79 
80  put(buffer, w);
81  }
82 
83 
84  /**
85  * Reset.
86  */
87  void reset()
88  {
89  mean = 0.0;
90  sigma = 0.0;
91  total = 0.0;
92  count = 0;
93  xmin = std::numeric_limits<double>::max();
94  xmax = std::numeric_limits<double>::lowest();
95 
96  buffer.clear();
97  }
98 
99 
100  /**
101  * Add quantile.
102  *
103  * \param Q quantile
104  * \return this quantile
105  */
107  {
108  mean += Q.mean;
109  sigma += Q.sigma;
110  total += Q.total;
111  count += Q.count;
112  xmin = std::min(xmin, Q.xmin);
113  xmax = std::max(xmax, Q.xmax);
114 
115  if (quantiles) {
116  std::copy(Q.buffer.begin(), Q.buffer.end(), std::inserter(buffer, buffer.end()));
117  }
118 
119  return *this;
120  }
121 
122 
123  /**
124  * Put value.
125  *
126  * \param x value
127  * \param w weight
128  */
129  void put(const double x, const double w = 1.0)
130  {
131  total += w;
132  count += 1;
133 
134  if (count == 1) {
135 
136  mean = x;
137  sigma = 0.0;
138 
139  } else {
140 
141  const double new_mean = mean + w * (x - mean) / total;
142  const double new_sigma = sigma + w * (x - mean) * (x - new_mean);
143 
144  // set up for next iteration
145 
146  mean = new_mean;
147  sigma = new_sigma;
148  }
149 
150  xmin = std::min(xmin, x);
151  xmax = std::max(xmax, x);
152 
153  if (quantiles) {
154  buffer.insert(std::make_pair(x,w));
155  }
156  }
157 
158 
159  /**
160  * Put data.
161  *
162  * \param buffer input data
163  * \param w weight
164  */
165  template<class JElement_t, class JAllocator_t>
167  const double w = 1.0)
168  {
169  for (typename array_type<JElement_t, JAllocator_t>::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
170  put(*i, w);
171  }
172  }
173 
174 
175  /**
176  * Get total count.
177  *
178  * \return count
179  */
180  long long int getCount() const
181  {
182  return count;
183  }
184 
185 
186  /**
187  * Get total weight.
188  *
189  * \return weight
190  */
191  double getTotal() const
192  {
193  return total;
194  }
195 
196 
197  /**
198  * Get minimum.
199  *
200  * \return minimum
201  */
202  double getMin() const
203  {
204  return xmin;
205  }
206 
207 
208  /**
209  * Get maximum.
210  *
211  * \return maximum
212  */
213  double getMax() const
214  {
215  return xmax;
216  }
217 
218 
219  /**
220  * Get mean value.
221  *
222  * \return mean value
223  */
224  double getMean() const
225  {
226  if (count != 0.0)
227  return mean;
228  else
229  THROW(JDivisionByZero, "JQuantile::getMean()");
230  }
231 
232 
233  /**
234  * Get standard deviation
235  *
236  * \return standard deviation
237  */
238  double getSTDev() const
239  {
240  if (count > 1)
241  return sqrt(count * sigma/(total * (count - 1)));
242  else
243  THROW(JDivisionByZero, "JQuantile::getSTDev()");
244  }
245 
246 
247  /**
248  * Get maximal deviation from average.
249  *
250  * \param relative if true, relative to average, else absolute
251  * \return deviation
252  */
253  double getDeviation(const bool relative = true) const
254  {
255  if (relative)
256  return std::max(getMax() - getMean(), getMean() - getMin());
257  else
258  return getMax() - getMin();
259  }
260 
261 
262  /**
263  * Test relative accuracy.
264  *
265  * \param precision relative precision
266  * \return true if reached accuracy; else false
267  */
268  bool hasAccuracy(const double precision) const
269  {
270  return getCount() > 1 && getSTDev() < precision * getMean();
271  }
272 
273 
274  /**
275  * Get quantile.
276  *
277  * \param Q quantile
278  * \param reverse reverse
279  * \return value
280  */
281  double getQuantile(const double Q, const bool reverse = false) const
282  {
283  if (quantiles) {
284 
285  double W = 0.0;
286 
287  for (std::map<double, double>::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
288  W += i->second;
289  }
290 
291  if (reverse)
292  return getQuantile(buffer.rbegin(), buffer.rend(), Q*W);
293  else
294  return getQuantile(buffer. begin(), buffer. end(), Q*W);
295  }
296 
297  THROW(JNoValue, "Option 'quantiles' at JQuantile() incompatible with method getQuantile().");
298  }
299 
300 
301  /**
302  * Print quantile.
303  *
304  * \param out output stream
305  * \param lpr long print
306  */
307  std::ostream& print(std::ostream& out, bool lpr = true) const
308  {
309  using namespace std;
310 
311  const int nc = getTitle().size();
312 
313  if (lpr) {
314  out << setw(nc) << left << " " << ' '
315  << setw(10) << left << " mean" << ' '
316  << setw(10) << left << " STD" << ' '
317  << setw(10) << left << " deviation" << endl;
318  }
319 
320  out << setw(nc) << left << getTitle() << ' '
321  << SCIENTIFIC(10,2) << getMean() << ' '
322  << SCIENTIFIC(10,2) << getSTDev() << ' '
323  << SCIENTIFIC(10,2) << getDeviation(false) << endl;
324 
325  return out;
326  }
327 
328 
329  /**
330  * Print quantile.
331  *
332  * \param out output stream
333  * \param quantile quantile
334  * \return output stream
335  */
336  friend inline std::ostream& operator<<(std::ostream& out, const JQuantile& quantile)
337  {
338  return quantile.print(out, getLongprint(out));
339  }
340 
341  protected:
342  /**
343  * Get quantile.
344  *
345  * \param __begin begin of data
346  * \param __end end of data
347  * \param W weight
348  * \return value
349  */
350  template<class T>
351  static double getQuantile(T __begin, T __end, const double W)
352  {
353  double w = 0.0;
354 
355  for (T i = __begin; i != __end; ++i) {
356 
357  w += i->second;
358 
359  if (w >= W) {
360  return i->first;
361  }
362  }
363 
364  THROW(JNoValue, "Invalid weight " << W);
365  }
366 
367 
368  bool quantiles;
369  double mean;
370  double sigma;
371  double total;
372  long long int count;
373  double xmin;
374  double xmax;
376  };
377 }
378 
379 #endif
data_type w[N+1][M+1]
Definition: JPolint.hh:741
Exceptions.
bool hasAccuracy(const double precision) const
Test relative accuracy.
Definition: JQuantile.hh:268
Q(UTCMax_s-UTCMin_s)-livetime_s
Auxiliary base class for aritmetic operations of derived class types.
Definition: JMath.hh:110
double getTotal() const
Get total weight.
Definition: JQuantile.hh:191
void reset()
Reset.
Definition: JQuantile.hh:87
double getSTDev() const
Get standard deviation.
Definition: JQuantile.hh:238
Auxiliary data structure for running average, standard deviation and quantiles.
Definition: JQuantile.hh:43
friend std::ostream & operator<<(std::ostream &out, const JQuantile &quantile)
Print quantile.
Definition: JQuantile.hh:336
double getMax() const
Get maximum.
Definition: JQuantile.hh:213
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
long long int count
Definition: JQuantile.hh:372
const std::string & getTitle() const
Get title.
Definition: JTitle.hh:55
double getMean() const
Get mean value.
Definition: JQuantile.hh:224
double getQuantile(const double Q, const bool reverse=false) const
Get quantile.
Definition: JQuantile.hh:281
Exception for missing value.
Definition: JException.hh:198
JQuantile(const JTitle &title="", const bool quantiles=false)
Constructor.
Definition: JQuantile.hh:53
double getMin() const
Get minimum.
Definition: JQuantile.hh:202
std::multimap< double, double > buffer
Definition: JQuantile.hh:375
std::ostream & print(std::ostream &out, bool lpr=true) const
Print quantile.
Definition: JQuantile.hh:307
Auxiliary class for title.
Definition: JTitle.hh:19
const array_type< JValue_t > & make_array(const JValue_t(&array)[N])
Method to create array of values.
Definition: JVectorize.hh:54
do set_variable OUTPUT_DIRECTORY $WORKDIR T
bool getLongprint(std::ostream &out)
Get long print option.
Definition: JManip.hh:121
void put(const double x, const double w=1.0)
Put value.
Definition: JQuantile.hh:129
void put(const array_type< JElement_t, JAllocator_t > &buffer, const double w=1.0)
Put data.
Definition: JQuantile.hh:166
long long int getCount() const
Get total count.
Definition: JQuantile.hh:180
Auxiliary methods to convert data members or return values of member methods of a set of objects to a...
I/O manipulators.
JQuantile & add(const JQuantile &Q)
Add quantile.
Definition: JQuantile.hh:106
double getDeviation(const bool relative=true) const
Get maximal deviation from average.
Definition: JQuantile.hh:253
std::string title
Definition: JTitle.hh:73
Exception for division by zero.
Definition: JException.hh:270
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:139
JQuantile(const JTitle &title, const array_type< JElement_t, JAllocator_t > &buffer, const bool quantiles=false, const double w=1.0)
Constructor.
Definition: JQuantile.hh:71
Base class for data structures with artithmetic capabilities.
static double getQuantile(T __begin, T __end, const double W)
Get quantile.
Definition: JQuantile.hh:351
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:484
Auxiliary data structure for return type of make methods.
Definition: JVectorize.hh:26