From: Alex Schmidt <alex@treefish.org>
Date: Thu, 24 May 2012 11:23:46 +0000 (+0200)
Subject: ...
X-Git-Url: http://git.treefish.org/~alex/phys/latlib.git/commitdiff_plain/d4e3cb0ebd372ab5e931490599fa0c72831bc23f

...
---

diff --git a/obs.cpp b/obs.cpp
index f9a7dbf..75c1d93 100644
--- a/obs.cpp
+++ b/obs.cpp
@@ -4,162 +4,63 @@
 
 using namespace std;
 
-namespace obs
+template <class obstype>
+void obs<obstype>::addMeas(obstype val[], int valsize)
 {
-  vector<meas> mymeas;
-
-  double addMeas(int id, double val)
-  {
-    meas tmpmeas;
-    tmpmeas.id = id;
-    tmpmeas.val.push_back(val);
-    mymeas.push_back(tmpmeas);
-    return val;
-  }
-
-  double* addMeas(int id, double *val, int valsize)
-  {
-    meas tmpmeas;
-    tmpmeas.id = id;
-    for( int ival=0; ival<valsize; ival++ ) tmpmeas.val.push_back( val[ival] );
-    mymeas.push_back(tmpmeas);
-    return val;
-  }
-
-  double meanval(int id)
-  {
-    double meanval = 0;
-    int vecsize = mymeas.size();
-    int ivals=0;
-    for (int i=0; i<vecsize; i++)
-      {
-	if ( mymeas.at(i).id == id )
-	  {
-	    meanval += mymeas.at(i).val[0];
-	    ivals++;
-	  }
-      }
-    return meanval/ivals;
-  }
-
-  double meanerr(int id)
-  {
-    int ivals = 0;
-    int vecsize = mymeas.size();
-    double mymean = meanval(id);
-    double meanerr = 0;
-
-    for (int i=0; i<vecsize; i++)
-      {
-	if ( mymeas.at(i).id == id )
-	  {
-	    meanerr += pow(mymeas.at(i).val[0]-mymean, 2);
-	    ivals++;
-	  }
-      }
-
-    meanerr = meanerr / ( ivals*(ivals-1) );
-    return sqrt(meanerr);
-  }
-
-  obsa normObs(int id)
-  {
-    obsa tmp;
-    tmp.mean = meanval(id);
-    tmp.err = meanerr(id);
-    return tmp;
-  }
-
-  obsa jackObs(int id, double (*func)(vector<double> vals))
-  {
-    int nmeas = 0;
-    vector<double> mymeans;
-    int vecsize = mymeas.size();
-    obsa jackobs;
-
-    jackobs.mean = 0;
-
-    for (int ijack=0; ijack < nmeas || nmeas == 0; ijack++)
-      {
-
-	nmeas = 0;
-	vector<double> myvals;
-
-	for (int i=0; i<vecsize; i++)
-	  {
-	    if ( mymeas.at(i).id == id )
-	      {
-		nmeas++;
-		if ( i == ijack ) continue;
-		myvals.push_back( mymeas.at(i).val[0] );
-	      }
-	  }
-
-	mymeans.push_back( func(myvals) );
-	jackobs.mean += mymeans.at(mymeans.size()-1);
-      }
-
-    jackobs.mean /= nmeas;
-    
-    /* compute jack error */
-    vecsize = mymeans.size();
-    jackobs.err = 0;
-    for (int i=0; i<vecsize; i++)
-      {
-	jackobs.err += pow(mymeans.at(i) - jackobs.mean, 2);
-      }
-    jackobs.err *= (double)(nmeas-1)/nmeas;
-    jackobs.err = sqrt(jackobs.err);
-
-    return jackobs;
-  }
-
-  obsa jackObs(int id, double (*func)(vector< vector<double> > vals))
-  {
-    int nmeas = 0;
-    vector<double> mymeans;
-    int vecsize = mymeas.size();
-    obsa jackobs;
-
-    jackobs.mean = 0;
-
-    for (int ijack=0; ijack < nmeas || nmeas == 0; ijack++)
-      {
-
-	nmeas = 0;
-	vector< vector<double> > myvals;
-
-	for (int i=0; i<vecsize; i++)
-	  {
-	    if ( mymeas.at(i).id == id )
-	      {
-		nmeas++;
-		if ( i == ijack ) continue;
-		myvals.push_back( mymeas[i].val );
-	      }
-	  }
+  for(int i=0; i<valsize; i++) measurements.push_back(val[i]);
+}
 
-	mymeans.push_back( func(myvals) );
-	jackobs.mean += mymeans.at(mymeans.size()-1);
-      }
+template <class obstype>
+void obs<obstype>::addMeas(const obstype& val)
+{
+  measurements.push_back(val);
+}
 
-    jackobs.mean /= nmeas;
+template <class obstype>
+void obs<obstype>::mean(const string& compid, const list<double>& meas)
+{
+  computations[compid][0] = 0;
+  computations[compid][1] = 0;
+  int nmeas = meas.size();
+
+  for(list<double>::iterator measIt = meas.begin(); measIt != meas.end(); ++measIt)
+    computations[compid][0] += *measIt;
+  computations[compid][0] /= nmeas;
+  
+  for(list<double>::iterator measIt = meas.begin(); measIt != meas.end(); ++measIt)
+    computations[compid][1] += pow( *measIt - computations[compid][0], 2 );
+  computations[compid][1] /= nmeas*(nmeas-1);
+  computations[compid][1] = sqrt(computations[compid][1]);	 
+}				 
     
-    /* compute jack error */
-    vecsize = mymeans.size();
-    jackobs.err = 0;
-    for (int i=0; i<vecsize; i++)
-      {
-	jackobs.err += pow(mymeans.at(i) - jackobs.mean, 2);
-      }
-    jackobs.err *= (double)(nmeas-1)/nmeas;
-    jackobs.err = sqrt(jackobs.err);
-
-    return jackobs;
-  }
-
-  void reset()
-  {
-    mymeas.clear();
-  }
+template <class obstype>
+void obs<obstype>::computeJack(const string& compid, double (*func)(const list<obstype>& vals, void *para), void *para)
+{
+  int nmeas=measurements.size();
+  double manymeans[nmeas];
+
+  computations[compid][0] = 0;
+  computations[compid][1] = 0;
+
+  int imeas=0;
+  for(typename list<obstype>::iterator removedIt = measurements.begin(); removedIt != measurements.end(); imeas++)
+    {
+      obstype removed = *removedIt;
+      
+      typename list<obstype>::iterator nextAfterIt = removedIt;
+      ++nextAfterIt;
+
+      measurements.erase(removedIt);
+     
+      manymeans[imeas] = func(measurements, para);
+      computations[compid][0] += manymeans[imeas];
+
+      measurements.insert(nextAfterIt, removed);
+    }
+  computations[compid][0] /= nmeas;
+
+  for(int imean=0; imean<nmeas; imean++)
+    computations[compid][1] += pow( computations[compid][0] - manymeans[imean], 2 );
+  computations[compid][1] *= (double)(nmeas-1)/nmeas;
+  computations[compid][1] = sqrt(computations[compid][1]);
 }
diff --git a/obs.h b/obs.h
index fa40a40..38fce5a 100644
--- a/obs.h
+++ b/obs.h
@@ -1,31 +1,27 @@
 #ifndef OBS_H
 #define OBS_H
 
-#include <vector>
+#include <list>
+#include <string>
+#include <map>
 
 using namespace std;
 
-namespace obs
+template <class obstype>
+class obs
 {
-  struct obsa
-  {
-    double mean;
-    double err;
-  };
+ public:
+  void addMeas(obstype val[], int valsize);
+  void addMeas(const obstype& val);
 
-  struct meas
-  {
-    int id;
-    vector<double> val;
-  };
+  void computeMean(const string& compid) { mean(compid, measurements); }
+  void computeJack(const string& compid, double (*func)(const list<obstype>& vals, void *para), void *para=NULL);
 
-  double addMeas(int id, double val);
-  double* addMeas(int id, double *val, int valsize);
+ private:
+  list<obstype> measurements;
+  map<string,double[2]> computations;
 
-  obsa jackObs(int id, double (*func)(vector< vector<double> > vals));
-  obsa normObs(int id);
-
-  void reset();
-}
+  void mean(const string& compid, const list<double>& meas);
+};
 
 #endif