From 9f6e969663a4126813ae8a37d87acdb698ef8e62 Mon Sep 17 00:00:00 2001
From: Alex Schmidt <alex@treefish.org>
Date: Thu, 26 Jan 2012 19:55:33 +0100
Subject: [PATCH] ...

---
 .gitignore     |   4 ++
 CMakeLists.txt |   3 +
 obs.cpp        | 165 +++++++++++++++++++++++++++++++++++++++++++++++++
 obs.h          |  31 ++++++++++
 test           |   0
 5 files changed, 203 insertions(+)
 create mode 100644 .gitignore
 create mode 100644 CMakeLists.txt
 create mode 100644 obs.cpp
 create mode 100644 obs.h
 delete mode 100644 test

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..cdc8bbb
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,4 @@
+CMakeCache.txt
+CMakeFiles
+cmake_install.cmake
+Makefile
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000..7ef5aae
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,3 @@
+project(latlib)
+
+add_library(obs obs.cpp)
\ No newline at end of file
diff --git a/obs.cpp b/obs.cpp
new file mode 100644
index 0000000..f9a7dbf
--- /dev/null
+++ b/obs.cpp
@@ -0,0 +1,165 @@
+#include "obs.h"
+
+#include <math.h>
+
+using namespace std;
+
+namespace obs
+{
+  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 );
+	      }
+	  }
+
+	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;
+  }
+
+  void reset()
+  {
+    mymeas.clear();
+  }
+}
diff --git a/obs.h b/obs.h
new file mode 100644
index 0000000..fa40a40
--- /dev/null
+++ b/obs.h
@@ -0,0 +1,31 @@
+#ifndef OBS_H
+#define OBS_H
+
+#include <vector>
+
+using namespace std;
+
+namespace obs
+{
+  struct obsa
+  {
+    double mean;
+    double err;
+  };
+
+  struct meas
+  {
+    int id;
+    vector<double> val;
+  };
+
+  double addMeas(int id, double val);
+  double* addMeas(int id, double *val, int valsize);
+
+  obsa jackObs(int id, double (*func)(vector< vector<double> > vals));
+  obsa normObs(int id);
+
+  void reset();
+}
+
+#endif
diff --git a/test b/test
deleted file mode 100644
index e69de29..0000000
-- 
2.39.5