obstat.hpp

   1 #ifndef OBSTAT_HPP
   2 #define OBSTAT_HPP
   3
   4 #include <vector>
   5 #include <string>
   6 #include <map>
   7 #include <math.h>
   8
   9 using namespace std;
  10
  11 template <typename obstype>
  12 class obstat
  13 {
  14 public:
  15   void addMeas(const obstype& val);
  16   void addMeas(obstype val[], int valsize);
  17
  18   void computeEasy(const string& compid="one", const int& ival=0) { mean(compid, &measurements, ival); }
  19
  20   void computeJack(const string& compid, double (*func)(vector< vector <obstype> > *vals, void *para)=NULL, void *para=NULL);
  21   void computeJack(double (*func)(vector< vector <obstype> > *vals, void *para)=NULL, void *para=NULL) { computeJack("one",func,para); }
  22
  23   double getMean(const string& compid="one") { return computations[compid].val; }
  24   double getErr(const string& compid="one") { return computations[compid].err; }
  25
  26   void reset();
  27
  28 private:
  29   struct result{
  30     double val;
  31     double err;
  32   };
  33
  34   vector< vector<obstype> > measurements;
  35   map<string,result> computations;
  36
  37   void mean(const string& compid, vector< vector <double> > *meas, const int& ival);
  38   void mean(const string& compid, vector< vector <int> > *meas, const int& ival);
  39 };
  40
  41
  42 template <typename obstype>
  43 void obstat<obstype>::reset()
  44 {
  45   computations.clear();
  46   measurements.clear();
  47 }
  48
  49 template <typename obstype>
  50 void obstat<obstype>::addMeas(const obstype& val)
  51 {
  52   measurements.push_back( vector<obstype>(1,val) );
  53 }
  54
  55 template <typename obstype>
  56 void obstat<obstype>::addMeas(obstype val[], int valsize)
  57 {
  58   vector<obstype> tmpvec;
  59   for(int i=0; i<valsize; i++) tmpvec.push_back(val[i]);
  60   measurements.push_back(tmpvec);
  61 }
  62
  63 template <typename obstype>
  64 void obstat<obstype>::mean(const string& compid, vector< vector<double> > *meas, const int& ival)
  65 {
  66   computations[compid].val = 0;
  67   computations[compid].err = 0;
  68
  69   int nmeas = meas->size();
  70
  71   for(vector< vector<double> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
  72     computations[compid].val += (*measIt)[ival];
  73   computations[compid].val /= nmeas;
  74
  75   for(vector< vector<double> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
  76     computations[compid].err += pow( (*measIt)[ival] - computations[compid].val, 2 );
  77   computations[compid].err = sqrt( computations[compid].err ) / nmeas;
  78 }
  79
  80 template <typename obstype>
  81 void obstat<obstype>::mean(const string& compid, vector< vector<int> > *meas, const int& ival)
  82 {
  83   computations[compid].val = 0;
  84   computations[compid].err = 0;
  85
  86   int nmeas = meas->size();
  87
  88   for(vector< vector<int> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
  89     computations[compid].val += (*measIt)[ival];
  90   computations[compid].val /= nmeas;
  91
  92   for(vector< vector<int> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
  93     computations[compid].err += pow( (*measIt)[ival] - computations[compid].val, 2 );
  94   computations[compid].err = sqrt( computations[compid].err ) / nmeas;
  95 }
  96
  97 template <typename obstype>
  98 void obstat<obstype>::computeJack(const string& compid, double (*func)(vector< vector<obstype> > *vals, void *para), void *para)
  99 {
 100   int nmeas=measurements.size();
 101   double manymeans[nmeas];
 102
 103   computations[compid].val = 0;
 104   computations[compid].err = 0;
 105
 106   int imeas=0;
 107   for(typename vector< vector<obstype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); ++removedIt, imeas++)
 108     {
 109       vector<obstype> removed = *removedIt;
 110
 111       *removedIt = measurements.back();
 112       measurements.pop_back();
 113
 114       manymeans[imeas] = func(&measurements, para);
 115       computations[compid].val += manymeans[imeas];
 116
 117       measurements.push_back( *removedIt );
 118       *removedIt = removed;
 119     }
 120   computations[compid].val /= nmeas;
 121
 122   for(int imean=0; imean<nmeas; imean++)
 123     computations[compid].err += pow( manymeans[imean] - computations[compid].val, 2 );
 124   computations[compid].err *= (double)(nmeas-1)/nmeas;
 125   computations[compid].err = sqrt(computations[compid].err);
 126 }
 127
 128 #endif