Jpp 19.3.0-rc.3
the software that should make you happy
Loading...
Searching...
No Matches
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
23namespace JTOOLS {}
24namespace JPP { using namespace JTOOLS; }
25
26namespace JTOOLS {
27
28 using JMATH::JMath;
29 using JLANG::JTitle;
31 using JLANG::JNoValue;
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) :
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) :
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 wmin = std::numeric_limits<double>::max();
96 wmax = std::numeric_limits<double>::lowest();
97
98 buffer.clear();
99 }
100
101
102 /**
103 * Add quantile.
104 *
105 * \param Q quantile
106 * \return this quantile
107 */
109 {
110 mean += Q.mean;
111 sigma += Q.sigma;
112 total += Q.total;
113 count += Q.count;
114 xmin = std::min(xmin, Q.xmin);
115 xmax = std::max(xmax, Q.xmax);
116 wmin = std::min(wmin, Q.wmin);
117 wmax = std::max(wmax, Q.wmax);
118
119 if (quantiles) {
120 std::copy(Q.buffer.begin(), Q.buffer.end(), std::inserter(buffer, buffer.end()));
121 }
122
123 return *this;
124 }
125
126
127 /**
128 * Put value.
129 *
130 * \param x value
131 * \param w weight
132 */
133 void put(const double x, const double w = 1.0)
134 {
135 total += w;
136 count += 1;
137
138 if (count == 1) {
139
140 mean = x;
141 sigma = 0.0;
142
143 } else {
144
145 const double new_mean = mean + w * (x - mean) / total;
146 const double new_sigma = sigma + w * (x - mean) * (x - new_mean);
147
148 // set up for next iteration
149
150 mean = new_mean;
151 sigma = new_sigma;
152 }
153
154 xmin = std::min(xmin, x);
155 xmax = std::max(xmax, x);
156 wmin = std::min(wmin, w);
157 wmax = std::max(wmax, w);
158
159 if (quantiles) {
160 buffer.insert(std::make_pair(x,w));
161 }
162 }
163
164
165 /**
166 * Put data.
167 *
168 * \param buffer input data
169 * \param w weight
170 */
171 template<class JElement_t, class JAllocator_t>
173 const double w = 1.0)
174 {
175 for (typename array_type<JElement_t, JAllocator_t>::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
176 put(*i, w);
177 }
178 }
179
180
181 /**
182 * Get total count.
183 *
184 * \return count
185 */
186 long long int getCount() const
187 {
188 return count;
189 }
190
191
192 /**
193 * Get total weight.
194 *
195 * \return weight
196 */
197 double getTotal() const
198 {
199 return total;
200 }
201
202
203 /**
204 * Get minimum value.
205 *
206 * \return minimum value
207 */
208 double getXmin() const
209 {
210 return xmin;
211 }
212
213
214 /**
215 * Get maximum value.
216 *
217 * \return maximum value
218 */
219 double getXmax() const
220 {
221 return xmax;
222 }
223
224
225 /**
226 * Get minimum weight.
227 *
228 * \return minimum weight
229 */
230 double getWmin() const
231 {
232 return wmin;
233 }
234
235
236 /**
237 * Get maximum weight.
238 *
239 * \return maximum weight
240 */
241 double getWmax() const
242 {
243 return wmax;
244 }
245
246
247 /**
248 * Get mean value.
249 *
250 * \return mean value
251 */
252 double getMean() const
253 {
254 if (count != 0.0)
255 return mean;
256 else
257 THROW(JDivisionByZero, "JQuantile::getMean()");
258 }
259
260
261 /**
262 * Get mean value.
263 *
264 * \param value default value
265 * \return mean value
266 */
267 double getMean(const double value) const
268 {
269 if (count != 0.0)
270 return getMean();
271 else
272 return value;
273 }
274
275
276 /**
277 * Get standard deviation
278 *
279 * \return standard deviation
280 */
281 double getSTDev() const
282 {
283 if (count > 1)
284 return sqrt(count * sigma/(total * (count - 1)));
285 else
286 THROW(JDivisionByZero, "JQuantile::getSTDev()");
287 }
288
289
290 /**
291 * Get standard deviation
292 *
293 * \param value default value
294 * \return standard deviation
295 */
296 double getSTDev(const double value) const
297 {
298 if (count > 1)
299 return getSTDev();
300 else
301 return value;
302 }
303
304
305 /**
306 * Get maximal deviation from average.
307 *
308 * \param relative if true, relative to average, else absolute
309 * \return deviation
310 */
311 double getDeviation(const bool relative = true) const
312 {
313 if (relative)
314 return std::max(getXmax() - getMean(), getMean() - getXmin());
315 else
316 return getXmax() - getXmin();
317 }
318
319
320 /**
321 * Test relative accuracy.
322 *
323 * \param precision relative precision
324 * \return true if reached accuracy; else false
325 */
326 bool hasAccuracy(const double precision) const
327 {
328 return getCount() > 1 && getSTDev() < precision * getMean();
329 }
330
331
332 /**
333 * Options for evaluation of quantile.
334 */
336 forward_t = +1, //!< forward
337 symmetric_t = 0, //!< symmatric
338 backward_t = -1 //!< backward
339 };
340
341
342 /**
343 * Get quantile.
344 *
345 * \param Q quantile
346 * \param option option
347 * \return value
348 */
349 double getQuantile(const double Q, const Quantile_t option = forward_t) const
350 {
351 if (quantiles) {
352
353 double W = 0.0;
354
355 for (std::map<double, double>::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
356 W += i->second;
357 }
358
359 switch (option) {
360 case forward_t:
361 return (getQuantile(buffer. begin(), buffer. end(), Q*W));
362 case symmetric_t:
363 return (getQuantile(buffer.rbegin(), buffer.rend(), 0.5*(1.0 - Q)*W) -
364 getQuantile(buffer. begin(), buffer. end(), 0.5*(1.0 - Q)*W));
365 case backward_t:
366 return (getQuantile(buffer.rbegin(), buffer.rend(), Q*W));
367 default:
368 THROW(JNoValue, "Invalid option " << option);
369 }
370 }
371
372 THROW(JNoValue, "Option 'quantiles' at JQuantile() incompatible with method getQuantile().");
373 }
374
375
376 /**
377 * Print quantile.
378 *
379 * \param out output stream
380 * \param lpr long print
381 */
382 std::ostream& print(std::ostream& out, bool lpr = true) const
383 {
384 using namespace std;
385
386 const int nc = getTitle().size();
387
388 if (lpr) {
389 out << setw(nc) << left << " " << ' '
390 << setw(12) << left << " mean" << ' '
391 << setw(12) << left << " STD" << ' '
392 << setw(12) << left << " deviation" << endl;
393 }
394
395 out << setw(nc) << left << getTitle() << ' '
396 << SCIENTIFIC(12,5) << getMean() << ' '
397 << SCIENTIFIC(12,5) << getSTDev() << ' '
398 << SCIENTIFIC(12,5) << getDeviation(false) << endl;
399
400 return out;
401 }
402
403
404 /**
405 * Print quantile.
406 *
407 * \param out output stream
408 * \param quantile quantile
409 * \return output stream
410 */
411 friend inline std::ostream& operator<<(std::ostream& out, const JQuantile& quantile)
412 {
413 return quantile.print(out, getLongprint(out));
414 }
415
416 protected:
417 /**
418 * Get quantile.
419 *
420 * \param __begin begin of data
421 * \param __end end of data
422 * \param W weight
423 * \return value
424 */
425 template<class T>
426 static double getQuantile(T __begin, T __end, const double W)
427 {
428 double w = 0.0;
429
430 for (T i = __begin; i != __end; ++i) {
431
432 w += i->second;
433
434 if (w >= W) {
435 return i->first;
436 }
437 }
438
439 THROW(JNoValue, "Invalid weight " << W);
440 }
441
442
444 double mean;
445 double sigma;
446 double total;
447 long long int count;
448 double xmin;
449 double xmax;
450 double wmin;
451 double wmax;
453 };
454}
455
456#endif
Exceptions.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
I/O manipulators.
bool getLongprint(std::ostream &out)
Get long print option.
Definition JManip.hh:121
Base class for data structures with artithmetic capabilities.
Auxiliary methods to convert data members or return values of member methods of a set of objects to a...
Exception for division by zero.
Exception for missing value.
Auxiliary class for title.
Definition JTitle.hh:19
std::string title
Definition JTitle.hh:73
const std::string & getTitle() const
Get title.
Definition JTitle.hh:55
const array_type< JValue_t > & make_array(const JValue_t(&array)[N])
Method to create array of values.
Definition JVectorize.hh:54
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary classes and methods for multi-dimensional interpolations and histograms.
Auxiliary data structure for return type of make methods.
Definition JVectorize.hh:28
Auxiliary base class for aritmetic operations of derived class types.
Definition JMath.hh:347
Auxiliary data structure for running average, standard deviation and quantiles.
Definition JQuantile.hh:46
double getDeviation(const bool relative=true) const
Get maximal deviation from average.
Definition JQuantile.hh:311
double getQuantile(const double Q, const Quantile_t option=forward_t) const
Get quantile.
Definition JQuantile.hh:349
double getWmin() const
Get minimum weight.
Definition JQuantile.hh:230
void put(const array_type< JElement_t, JAllocator_t > &buffer, const double w=1.0)
Put data.
Definition JQuantile.hh:172
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
Quantile_t
Options for evaluation of quantile.
Definition JQuantile.hh:335
@ backward_t
backward
Definition JQuantile.hh:338
@ symmetric_t
symmatric
Definition JQuantile.hh:337
double getWmax() const
Get maximum weight.
Definition JQuantile.hh:241
double getMean(const double value) const
Get mean value.
Definition JQuantile.hh:267
JQuantile & add(const JQuantile &Q)
Add quantile.
Definition JQuantile.hh:108
std::ostream & print(std::ostream &out, bool lpr=true) const
Print quantile.
Definition JQuantile.hh:382
static double getQuantile(T __begin, T __end, const double W)
Get quantile.
Definition JQuantile.hh:426
double getXmin() const
Get minimum value.
Definition JQuantile.hh:208
friend std::ostream & operator<<(std::ostream &out, const JQuantile &quantile)
Print quantile.
Definition JQuantile.hh:411
double getSTDev() const
Get standard deviation.
Definition JQuantile.hh:281
bool hasAccuracy(const double precision) const
Test relative accuracy.
Definition JQuantile.hh:326
long long int count
Definition JQuantile.hh:447
void reset()
Reset.
Definition JQuantile.hh:87
double getTotal() const
Get total weight.
Definition JQuantile.hh:197
void put(const double x, const double w=1.0)
Put value.
Definition JQuantile.hh:133
std::multimap< double, double > buffer
Definition JQuantile.hh:452
JQuantile(const JTitle &title="", const bool quantiles=false)
Constructor.
Definition JQuantile.hh:53
long long int getCount() const
Get total count.
Definition JQuantile.hh:186
double getMean() const
Get mean value.
Definition JQuantile.hh:252
double getXmax() const
Get maximum value.
Definition JQuantile.hh:219
double getSTDev(const double value) const
Get standard deviation.
Definition JQuantile.hh:296
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488