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 meastype, typename restype>
  12 class obstat
  13 {
  14 public:
  15   void addMeas(const meastype& val);
  16   void addMeas(meastype val[], int valsize);
  17
  18   int computeEasy(const int& ival=0) { return mean(&measurements, ival); }
  19
  20   int computeJack(restype (*func)(vector< vector <meastype> > *vals, void *para), void *para=NULL);
  21   int computeJack(void (*preMeasFunc)(vector< vector <meastype> > *allVals, vector <meastype> *preCalculated, void *para),
  22                   restype (*measFunc)(vector <meastype> *preCalculated, vector <meastype> *excludedmeas, int nmeas, void *para), void *para=NULL);
  23
  24   restype getMean(int compid) { return computations[compid].val; }
  25   restype getErr(int compid) { return computations[compid].err; }
  26
  27   void reset();
  28
  29 private:
  30   struct result{
  31     restype val;
  32     restype err;
  33   };
  34
  35   vector< vector<meastype> > measurements;
  36   vector<result> computations;
  37
  38   int mean(vector< vector <meastype> > *meas, const int& ival);
  39 };
  40
  41 template <typename meastype, typename restype>
  42 void obstat<meastype,restype>::reset()
  43 {
  44   computations.clear();
  45   measurements.clear();
  46 }
  47
  48 template <typename meastype, typename restype>
  49 void obstat<meastype,restype>::addMeas(const meastype& val)
  50 {
  51   measurements.push_back( vector<meastype>(1,val) );
  52 }
  53
  54 template <typename meastype, typename restype>
  55 void obstat<meastype,restype>::addMeas(meastype val[], int valsize)
  56 {
  57   vector<meastype> tmpvec;
  58   for (int i=0; i<valsize; i++)
  59     tmpvec.push_back(val[i]);
  60   measurements.push_back(tmpvec);
  61 }
  62
  63 template <typename meastype, typename restype>
  64 int obstat<meastype,restype>::mean(vector< vector<meastype> > *meas, const int& ival)
  65 {
  66   result meanres;
  67
  68   meanres.val = 0;
  69   meanres.err = 0;
  70
  71   int nmeas = meas->size();
  72
  73   for (typename vector< vector<meastype> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
  74     meanres.val += (*measIt)[ival];
  75   meanres.val /= nmeas;
  76
  77   for (typename vector< vector<meastype> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
  78     meanres.err += pow( (*measIt)[ival] - meanres.val, 2 );
  79   meanres.err = sqrt( meanres.err ) / (restype)nmeas;
  80
  81   computations.push_back(meanres);
  82
  83   return computations.size()-1;
  84 }
  85
  86 template <typename meastype, typename restype>
  87 int obstat<meastype,restype>::computeJack(restype (*func)(vector< vector<meastype> > *vals, void *para), void *para)
  88 {
  89   int nmeas=measurements.size();
  90   restype *manymeans = new restype[nmeas];
  91   result jackres;
  92
  93   jackres.val = 0;
  94   jackres.err = 0;
  95
  96   int imeas=0;
  97   for(typename vector< vector<meastype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); ++removedIt, imeas++)
  98     {
  99       vector<meastype> removed = *removedIt;
 100
 101       *removedIt = measurements.back();
 102       measurements.pop_back();
 103
 104       manymeans[imeas] = func(&measurements, para);
 105       jackres.val += manymeans[imeas];
 106
 107       measurements.push_back( *removedIt );
 108       *removedIt = removed;
 109     }
 110   jackres.val /= nmeas;
 111
 112   for(int imean=0; imean<nmeas; imean++)
 113     jackres.err += pow( manymeans[imean] - jackres.val, 2 );
 114   jackres.err *= (double)(nmeas-1)/nmeas;
 115   jackres.err = sqrt(jackres.err);
 116
 117   computations.push_back(jackres);
 118
 119   delete [] manymeans;
 120
 121   return computations.size()-1;
 122 }
 123
 124 template <typename meastype, typename restype>
 125 int obstat<meastype,restype>::computeJack(void (*preMeasFunc)(vector< vector <meastype> > *allVals, vector<meastype> *preCalculated, void *para),
 126                                           restype (*measFunc)(vector<meastype> *preCalculated, vector<meastype> *excludedmeas, int nmeas, void *para),
 127                                           void *para) {
 128   int nmeas=measurements.size();
 129   restype *manymeans = new restype[nmeas];
 130   result jackres;
 131   vector<meastype> preCalculated;
 132
 133   jackres.val = 0;
 134   jackres.err = 0;
 135
 136   preMeasFunc(&measurements, &preCalculated, para);
 137
 138   int imeas=0;
 139   for(typename vector< vector<meastype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); ++removedIt, imeas++) {
 140     manymeans[imeas] = measFunc(&preCalculated, &(*removedIt), measurements.size(), para);
 141     jackres.val += manymeans[imeas];
 142   }
 143   jackres.val /= nmeas;
 144
 145   for(int imean=0; imean<nmeas; imean++)
 146     jackres.err += pow( manymeans[imean] - jackres.val, 2 );
 147   jackres.err *= (double)(nmeas-1)/nmeas;
 148   jackres.err = sqrt(jackres.err);
 149
 150   computations.push_back(jackres);
 151
 152   delete [] manymeans;
 153
 154   return computations.size()-1;
 155 }
 156
 157 #endif