summaryrefslogtreecommitdiff
path: root/graphics/asymptote/statistics.h
diff options
context:
space:
mode:
Diffstat (limited to 'graphics/asymptote/statistics.h')
-rw-r--r--graphics/asymptote/statistics.h95
1 files changed, 88 insertions, 7 deletions
diff --git a/graphics/asymptote/statistics.h b/graphics/asymptote/statistics.h
index a8cca39bf4..5588f7a99d 100644
--- a/graphics/asymptote/statistics.h
+++ b/graphics/asymptote/statistics.h
@@ -1,18 +1,46 @@
#ifndef __statistics_h__
#define __statistics_h__ 1
+#include <queue>
+#include <cfloat>
+
namespace utils {
+template <class T, class S, class C>
+void clearpq(std::priority_queue<T, S, C>& q) {
+ struct HackedQueue : private std::priority_queue<T, S, C> {
+ static S& Container(std::priority_queue<T, S, C>& q) {
+ return q.*&HackedQueue::c;
+ }
+ };
+ HackedQueue::Container(q).clear();
+}
+
class statistics {
- unsigned int N;
+ size_t N;
double A;
double varL;
double varH;
+ double m,M;
+ double Median;
+ bool computeMedian;
+
+ // These heap stores are only used when computeMedian=true.
+ // Max heap stores the smaller half elements:
+ std::priority_queue<double> s;
+ // Min heap stores the greater half elements:
+ std::priority_queue<double,std::vector<double>,std::greater<double> > g;
+
public:
- statistics() {clear();}
- void clear() {N=0; A=varL=varH=0.0;}
+ statistics(bool computeMedian=false) : computeMedian(computeMedian) {
+ clear();
+ }
+ void clear() {N=0; A=varL=varH=0.0; m=DBL_MAX; M=-m; clearpq(s); clearpq(g);}
double count() {return N;}
double mean() {return A;}
+ double max() {return M;}
+ double min() {return m;}
+ double sum() {return N*A;}
void add(double t) {
++N;
double diff=t-A;
@@ -22,10 +50,53 @@ public:
varL += v;
else
varH += v;
+
+ if(t < m) m=t;
+ if(t > M) M=t;
+
+ if(computeMedian) {
+ if(N == 1)
+ s.push(Median=t);
+ else {
+ if(s.size() > g.size()) { // left side heap has more elements
+
+ if(t < Median) {
+ g.push(s.top());
+ s.pop();
+ s.push(t);
+ } else
+ g.push(t);
+
+ Median=0.5*(s.top()+g.top());
+ }
+
+ else if(s.size() == g.size()) { // both heaps are balanced
+ if(t < Median) {
+ s.push(t);
+ Median=(double) s.top();
+ } else {
+ g.push(t);
+ Median=(double) g.top();
+ }
+ }
+
+ else { // right side heap has more elements
+ if(t > Median) {
+ s.push(g.top());
+ g.pop();
+ g.push(t);
+ } else
+ s.push(t);
+
+ Median=0.5*(s.top()+g.top());
+ }
+ }
+ }
}
+
double stdev(double var, double f) {
- double factor=N > f ? f/(N-f) : 0.0;
- return sqrt(var*factor);
+ if(N <= f) return DBL_MAX;
+ return sqrt(var*f/(N-f));
}
double stdev() {
return stdev(varL+varH,1.0);
@@ -36,8 +107,18 @@ public:
double stdevH() {
return stdev(varH,2.0);
}
- void output(const char *text, unsigned int m) {
- std::cout << text << ":\n"
+ double stderror() {
+ return stdev()/sqrt(N);
+ }
+ double median() {
+ if(!computeMedian) {
+ std::cerr << "Constructor requires median=true" << std::endl;
+ exit(-1);
+ }
+ return Median;
+ }
+ void output(const char *text, size_t m) {
+ std::cout << text << ": \n"
<< m << "\t"
<< A << "\t"
<< stdevL() << "\t"