...

[phys/latlib.git] / obs.cpp

diff --git a/obs.cpp b/obs.cpp
index f9a7dbf93f7567214a340f882e93d026a0d0979b..75c1d93f14b26487dc5e653b78abfb78c70bbec4 100644 (file)
--- 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]);
 }