]> git.treefish.org Git - phys/latlib.git/blob - obs.hpp
...
[phys/latlib.git] / 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 = sqrt( computations[compid].err ) / nmeas;   
78 }
79
80 template <typename obstype>
81 void obs<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 obs<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