obs.cpp

   1 #include "obs.h"
   2
   3 #include <math.h>
   4
   5 using namespace std;
   6
   7 template <class obstype>
   8 void obs<obstype>::addMeas(obstype val[], int valsize)
   9 {
  10   for(int i=0; i<valsize; i++) measurements.push_back(val[i]);
  11 }
  12
  13 template <class obstype>
  14 void obs<obstype>::addMeas(const obstype& val)
  15 {
  16   measurements.push_back(val);
  17 }
  18
  19 template <class obstype>
  20 void obs<obstype>::mean(const string& compid, const list<double>& meas)
  21 {
  22   computations[compid][0] = 0;
  23   computations[compid][1] = 0;
  24   int nmeas = meas.size();
  25
  26   for(list<double>::iterator measIt = meas.begin(); measIt != meas.end(); ++measIt)
  27     computations[compid][0] += *measIt;
  28   computations[compid][0] /= nmeas;
  29
  30   for(list<double>::iterator measIt = meas.begin(); measIt != meas.end(); ++measIt)
  31     computations[compid][1] += pow( *measIt - computations[compid][0], 2 );
  32   computations[compid][1] /= nmeas*(nmeas-1);
  33   computations[compid][1] = sqrt(computations[compid][1]);
  34 }
  35
  36 template <class obstype>
  37 void obs<obstype>::computeJack(const string& compid, double (*func)(const list<obstype>& vals, void *para), void *para)
  38 {
  39   int nmeas=measurements.size();
  40   double manymeans[nmeas];
  41
  42   computations[compid][0] = 0;
  43   computations[compid][1] = 0;
  44
  45   int imeas=0;
  46   for(typename list<obstype>::iterator removedIt = measurements.begin(); removedIt != measurements.end(); imeas++)
  47     {
  48       obstype removed = *removedIt;
  49
  50       typename list<obstype>::iterator nextAfterIt = removedIt;
  51       ++nextAfterIt;
  52
  53       measurements.erase(removedIt);
  54
  55       manymeans[imeas] = func(measurements, para);
  56       computations[compid][0] += manymeans[imeas];
  57
  58       measurements.insert(nextAfterIt, removed);
  59     }
  60   computations[compid][0] /= nmeas;
  61
  62   for(int imean=0; imean<nmeas; imean++)
  63     computations[compid][1] += pow( computations[compid][0] - manymeans[imean], 2 );
  64   computations[compid][1] *= (double)(nmeas-1)/nmeas;
  65   computations[compid][1] = sqrt(computations[compid][1]);
  66 }