implemented partial virtual equilibration.

[phys/latlib.git] / obs.hpp

diff --git a/obs.hpp b/obs.hpp
index bba6687e344ed0f20fe74e691a5ca104a5fcb3ce..4d20bd6d901bff88f113d7a654a8fd87dca2ca4e 100644 (file)
--- a/obs.hpp
+++ b/obs.hpp
@@ -1,7 +1,7 @@
 #ifndef OBS_HPP
 #define OBS_HPP
 
-#include <list>
+#include <vector>
 #include <string>
 #include <map>
 #include <math.h>
@@ -17,24 +17,35 @@ public:
 
   void computeEasy(const string& compid="one", const int& ival=0) { mean(compid, &measurements, ival); }
 
-  void computeJack(const string& compid, double (*func)(list< vector <obstype> > *vals, void *para)=NULL, void *para=NULL);
-  void computeJack(double (*func)(list< vector <obstype> > *vals, void *para)=NULL, void *para=NULL) { computeJack("one",func,para); }
+  void computeJack(const string& compid, double (*func)(vector< vector <obstype> > *vals, void *para)=NULL, void *para=NULL);
+  void computeJack(double (*func)(vector< vector <obstype> > *vals, void *para)=NULL, void *para=NULL) { computeJack("one",func,para); }
 
   double getMean(const string& compid="one") { return computations[compid].val; }
   double getErr(const string& compid="one") { return computations[compid].err; }
 
+  void reset();
+
 private:
   struct result{
     double val;
     double err;
   };
 
-  list< vector<obstype> > measurements;
+  vector< vector<obstype> > measurements;
   map<string,result> computations;
   
-  void mean(const string& compid, list< vector <double> > *meas, const int& ival);
+  void mean(const string& compid, vector< vector <double> > *meas, const int& ival);
+  void mean(const string& compid, vector< vector <int> > *meas, const int& ival);
 };
 
+
+template <typename obstype> 
+void obs<obstype>::reset()
+{
+  computations.clear();
+  measurements.clear();
+}
+
 template <typename obstype> 
 void obs<obstype>::addMeas(const obstype& val)
 {
@@ -50,25 +61,41 @@ void obs<obstype>::addMeas(obstype val[], int valsize)
 }
 
 template <typename obstype>
-void obs<obstype>::mean(const string& compid, list< vector<double> > *meas, const int& ival)
+void obs<obstype>::mean(const string& compid, vector< vector<double> > *meas, const int& ival)
+{
+  computations[compid].val = 0;
+  computations[compid].err = 0;
+  
+  int nmeas = meas->size();
+
+  for(vector< vector<double> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
+    computations[compid].val += (*measIt)[ival];
+  computations[compid].val /= nmeas;
+  
+  for(vector< vector<double> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
+    computations[compid].err += pow( (*measIt)[ival] - computations[compid].val, 2 );
+  computations[compid].err = sqrt( computations[compid].err ) / nmeas;  
+}
+
+template <typename obstype>
+void obs<obstype>::mean(const string& compid, vector< vector<int> > *meas, const int& ival)
 {
   computations[compid].val = 0;
   computations[compid].err = 0;
   
   int nmeas = meas->size();
 
-  for(list< vector<double> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
+  for(vector< vector<int> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
     computations[compid].val += (*measIt)[ival];
   computations[compid].val /= nmeas;
   
-  for(list< vector<double> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
+  for(vector< vector<int> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
     computations[compid].err += pow( (*measIt)[ival] - computations[compid].val, 2 );
-  computations[compid].err /= nmeas*(nmeas-1);
-  computations[compid].err = sqrt( computations[compid].err );  
+  computations[compid].err = sqrt( computations[compid].err ) / nmeas;
 }                               
     
 template <typename obstype>
-void obs<obstype>::computeJack(const string& compid, double (*func)(list< vector<obstype> > *vals, void *para), void *para)
+void obs<obstype>::computeJack(const string& compid, double (*func)(vector< vector<obstype> > *vals, void *para), void *para)
 {
   int nmeas=measurements.size();
   double manymeans[nmeas];
@@ -77,26 +104,23 @@ void obs<obstype>::computeJack(const string& compid, double (*func)(list< vector
   computations[compid].err = 0;
 
   int imeas=0;
-  for(typename list< vector<obstype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); imeas++)
+  for(typename vector< vector<obstype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); ++removedIt, imeas++)
     {
       vector<obstype> removed = *removedIt;
       
-      typename list< vector<obstype> >::iterator nextAfterIt = removedIt;
-      ++nextAfterIt;
+      *removedIt = measurements.back();
+      measurements.pop_back();
 
-      measurements.erase(removedIt);
-     
       manymeans[imeas] = func(&measurements, para);
       computations[compid].val += manymeans[imeas];
-
-      measurements.insert(nextAfterIt, removed);
-
-      removedIt = nextAfterIt;
+      
+      measurements.push_back( *removedIt );
+      *removedIt = removed;
     }
   computations[compid].val /= nmeas;
 
   for(int imean=0; imean<nmeas; imean++)
-    computations[compid].err += pow( computations[compid].val - manymeans[imean], 2 );
+    computations[compid].err += pow( manymeans[imean] - computations[compid].val, 2 );
   computations[compid].err *= (double)(nmeas-1)/nmeas;
   computations[compid].err = sqrt(computations[compid].err);
 }