Fixed MPI cmake config.

[phys/latlib.git] / obstat.hpp

diff --git a/obstat.hpp b/obstat.hpp
index d66467b48a51bec64133938d850397d1d7e384f9..57a3431e6faa36078ecad0935256a6b8bc3e45d7 100644 (file)
--- a/obstat.hpp
+++ b/obstat.hpp
@@ -8,121 +8,150 @@
 
 using namespace std;
 
 
 using namespace std;
 

-template <typename obstype>
+template <typename meastype, typename restype>

 class obstat
 { 
 public:
 class obstat
 { 
 public:

-  void addMeas(const obstype& val);
-  void addMeas(obstype val[], int valsize);
+  void addMeas(const meastype& val);
+  void addMeas(meastype val[], int valsize);

 
 

-  void computeEasy(const string& compid="one", const int& ival=0) { mean(compid, &measurements, ival); }
+  int computeEasy(const int& ival=0) { return mean(&measurements, ival); }

 
 

-  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); }
+  int computeJack(restype (*func)(vector< vector <meastype> > *vals, void *para), void *para=NULL);
+  int computeJack(void (*preMeasFunc)(vector< vector <meastype> > *allVals, vector <meastype> *preCalculated, void *para), 
+                 restype (*measFunc)(vector <meastype> *preCalculated, vector <meastype> *excludedmeas, int nmeas, void *para), void *para=NULL);

 
 

-  double getMean(const string& compid="one") { return computations[compid].val; }
-  double getErr(const string& compid="one") { return computations[compid].err; }
+  restype getMean(int compid) { return computations[compid].val; }
+  restype getErr(int compid) { return computations[compid].err; }

 
   void reset();
 
 private:
   struct result{
 
   void reset();
 
 private:
   struct result{

-    double val;
-    double err;
+    restype val;
+    restype err;

   };
 
   };
 

-  vector< vector<obstype> > measurements;
-  map<string,result> computations;
+  vector< vector<meastype> > measurements;
+  vector<result> computations;

   
   

-  void mean(const string& compid, vector< vector <double> > *meas, const int& ival);
-  void mean(const string& compid, vector< vector <int> > *meas, const int& ival);
+  int mean(vector< vector <meastype> > *meas, const int& ival);

 };
 
 };
 

-
-template <typename obstype> 
-void obstat<obstype>::reset()
+template <typename meastype, typename restype>
+void obstat<meastype,restype>::reset()

 {
   computations.clear();
   measurements.clear();
 }
 
 {
   computations.clear();
   measurements.clear();
 }
 

-template <typename obstype> 
-void obstat<obstype>::addMeas(const obstype& val)
+template <typename meastype, typename restype>
+void obstat<meastype,restype>::addMeas(const meastype& val)

 {
 {

-  measurements.push_back( vector<obstype>(1,val) );
+  measurements.push_back( vector<meastype>(1,val) );

 }
 
 }
 

-template <typename obstype>
-void obstat<obstype>::addMeas(obstype val[], int valsize)
+template <typename meastype, typename restype>
+void obstat<meastype,restype>::addMeas(meastype val[], int valsize)

 {
 {

-  vector<obstype> tmpvec;
-  for(int i=0; i<valsize; i++) tmpvec.push_back(val[i]);
+  vector<meastype> tmpvec;
+  for (int i=0; i<valsize; i++) 
+    tmpvec.push_back(val[i]);

   measurements.push_back(tmpvec);
 }
 
   measurements.push_back(tmpvec);
 }
 

-template <typename obstype>
-void obstat<obstype>::mean(const string& compid, vector< vector<double> > *meas, const int& ival)
+template <typename meastype, typename restype>
+int obstat<meastype,restype>::mean(vector< vector<meastype> > *meas, const int& ival)

 {
 {

-  computations[compid].val = 0;
-  computations[compid].err = 0;
+  result meanres;
+
+  meanres.val = 0;
+  meanres.err = 0;

   
   int nmeas = meas->size();
 
   
   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 (typename vector< vector<meastype> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
+    meanres.val += (*measIt)[ival];
+  meanres.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;  
-}
+  for (typename vector< vector<meastype> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
+    meanres.err += pow( (*measIt)[ival] - meanres.val, 2 );
+  meanres.err = sqrt( meanres.err ) / (restype)nmeas;   

 
 

-template <typename obstype>
-void obstat<obstype>::mean(const string& compid, vector< vector<int> > *meas, const int& ival)
-{
-  computations[compid].val = 0;
-  computations[compid].err = 0;
-  
-  int nmeas = meas->size();
+  computations.push_back(meanres);

 
 

-  for(vector< vector<int> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
-    computations[compid].val += (*measIt)[ival];
-  computations[compid].val /= nmeas;
-  
-  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 = sqrt( computations[compid].err ) / nmeas;
-}                               
-    
-template <typename obstype>
-void obstat<obstype>::computeJack(const string& compid, double (*func)(vector< vector<obstype> > *vals, void *para), void *para)
+  return computations.size()-1;
+}
+
+template <typename meastype, typename restype>
+int obstat<meastype,restype>::computeJack(restype (*func)(vector< vector<meastype> > *vals, void *para), void *para)

 {
   int nmeas=measurements.size();
 {
   int nmeas=measurements.size();

-  double manymeans[nmeas];
+  restype *manymeans = new restype[nmeas];
+  result jackres;

 
 

-  computations[compid].val = 0;
-  computations[compid].err = 0;
+  jackres.val = 0;
+  jackres.err = 0;

 
   int imeas=0;
 
   int imeas=0;

-  for(typename vector< vector<obstype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); ++removedIt, imeas++)
+  for(typename vector< vector<meastype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); ++removedIt, imeas++)

     {
     {

-      vector<obstype> removed = *removedIt;
+      vector<meastype> removed = *removedIt;

       
       *removedIt = measurements.back();
       measurements.pop_back();
 
       manymeans[imeas] = func(&measurements, para);
       
       *removedIt = measurements.back();
       measurements.pop_back();
 
       manymeans[imeas] = func(&measurements, para);

-      computations[compid].val += manymeans[imeas];
+      jackres.val += manymeans[imeas];

       
       measurements.push_back( *removedIt );
       *removedIt = removed;
     }
       
       measurements.push_back( *removedIt );
       *removedIt = removed;
     }

-  computations[compid].val /= nmeas;
+  jackres.val /= nmeas;
+
+  for(int imean=0; imean<nmeas; imean++)
+    jackres.err += pow( manymeans[imean] - jackres.val, 2 );
+  jackres.err *= (double)(nmeas-1)/nmeas;
+  jackres.err = sqrt(jackres.err);
+
+  computations.push_back(jackres);

 
 

+  delete [] manymeans;
+
+  return computations.size()-1;
+}
+
+template <typename meastype, typename restype>
+int obstat<meastype,restype>::computeJack(void (*preMeasFunc)(vector< vector <meastype> > *allVals, vector<meastype> *preCalculated, void *para), 
+                                         restype (*measFunc)(vector<meastype> *preCalculated, vector<meastype> *excludedmeas, int nmeas, void *para),
+                                         void *para) {
+  int nmeas=measurements.size();
+  restype *manymeans = new restype[nmeas];
+  result jackres;
+  vector<meastype> preCalculated;
+
+  jackres.val = 0;
+  jackres.err = 0;
+
+  preMeasFunc(&measurements, &preCalculated, para);
+
+  int imeas=0;
+  for(typename vector< vector<meastype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); ++removedIt, imeas++) {
+    manymeans[imeas] = measFunc(&preCalculated, &(*removedIt), measurements.size(), para);
+    jackres.val += manymeans[imeas]; 
+  }
+  jackres.val /= nmeas;
+  

   for(int imean=0; imean<nmeas; imean++)
   for(int imean=0; imean<nmeas; imean++)

-    computations[compid].err += pow( manymeans[imean] - computations[compid].val, 2 );
-  computations[compid].err *= (double)(nmeas-1)/nmeas;
-  computations[compid].err = sqrt(computations[compid].err);
+    jackres.err += pow( manymeans[imean] - jackres.val, 2 );
+  jackres.err *= (double)(nmeas-1)/nmeas;
+  jackres.err = sqrt(jackres.err);
+
+  computations.push_back(jackres);
+
+  delete [] manymeans;
+
+  return computations.size()-1;

 }
 
 #endif
 }
 
 #endif