summaryrefslogtreecommitdiff
path: root/graphics/asymptote/statistics.h
blob: 5588f7a99d302394e6b3dee425bfb4f8cccef347 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
#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 {
  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(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;
    A += diff/N;
    double v=diff*(t-A);
    if(diff < 0.0)
      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) {
    if(N <= f) return DBL_MAX;
    return sqrt(var*f/(N-f));
  }
  double stdev() {
    return stdev(varL+varH,1.0);
  }
  double stdevL() {
    return stdev(varL,2.0);
  }
  double stdevH() {
    return stdev(varH,2.0);
  }
  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"
              << stdevH() << std::endl;
  }
};

}

#endif