]> git.treefish.org Git - phys/latlib.git/commitdiff
...
authorAlex Schmidt <alex@treefish.org>
Thu, 26 Jan 2012 18:55:33 +0000 (19:55 +0100)
committerAlex Schmidt <alex@treefish.org>
Thu, 26 Jan 2012 18:55:33 +0000 (19:55 +0100)
.gitignore [new file with mode: 0644]
CMakeLists.txt [new file with mode: 0644]
obs.cpp [new file with mode: 0644]
obs.h [new file with mode: 0644]
test [deleted file]

diff --git a/.gitignore b/.gitignore
new file mode 100644 (file)
index 0000000..cdc8bbb
--- /dev/null
@@ -0,0 +1,4 @@
+CMakeCache.txt
+CMakeFiles
+cmake_install.cmake
+Makefile
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644 (file)
index 0000000..7ef5aae
--- /dev/null
@@ -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 (file)
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 (file)
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 (file)
index e69de29..0000000