#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((double) 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