obs.hpp

   1 #ifndef OBS_HPP
   2 #define OBS_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 obs
  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 obs<obstype>::reset()
  44 {
  45   computations.clear();
  46   measurements.clear();
  47 }
  48
  49 template <typename obstype>
  50 void obs<obstype>::addMeas(const obstype& val)
  51 {
  52   measurements.push_back( vector<obstype>(1,val) );
  53 }
  54
  55 template <typename obstype>
  56 void obs<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 obs<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 /= nmeas-1;
  78   computations[compid].err = sqrt( computations[compid].err );
  79 }
  80
  81 template <typename obstype>
  82 void obs<obstype>::mean(const string& compid, vector< vector<int> > *meas, const int& ival)
  83 {
  84   computations[compid].val = 0;
  85   computations[compid].err = 0;
  86
  87   int nmeas = meas->size();
  88
  89   for(vector< vector<int> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
  90     computations[compid].val += (*measIt)[ival];
  91   computations[compid].val /= nmeas;
  92
  93   for(vector< vector<int> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
  94     computations[compid].err += pow( (*measIt)[ival] - computations[compid].val, 2 );
  95   computations[compid].err /= nmeas-1;
  96   computations[compid].err = sqrt( computations[compid].err );
  97 }
  98
  99 template <typename obstype>
 100 void obs<obstype>::computeJack(const string& compid, double (*func)(vector< vector<obstype> > *vals, void *para), void *para)
 101 {
 102   int nmeas=measurements.size();
 103   double manymeans[nmeas];
 104
 105   computations[compid].val = 0;
 106   computations[compid].err = 0;
 107
 108   int imeas=0;
 109   for(typename vector< vector<obstype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); ++removedIt, imeas++)
 110     {
 111       vector<obstype> removed = *removedIt;
 112
 113       *removedIt = measurements.back();
 114       measurements.pop_back();
 115
 116       manymeans[imeas] = func(&measurements, para);
 117       computations[compid].val += manymeans[imeas];
 118
 119       measurements.push_back( *removedIt );
 120       *removedIt = removed;
 121     }
 122   computations[compid].val /= nmeas;
 123
 124   for(int imean=0; imean<nmeas; imean++)
 125     computations[compid].err += pow( manymeans[imean] - computations[compid].val, 2 );
 126   computations[compid].err *= (double)(nmeas-1)/nmeas;
 127   computations[compid].err = sqrt(computations[compid].err);
 128 }
 129
 130 #endif