From 85cb7253b975501063e89285d46824d0988b82ee Mon Sep 17 00:00:00 2001
From: Alexander Schmidt <alex@treefish.org>
Date: Fri, 21 Feb 2014 00:02:18 +0100
Subject: [PATCH] Using precalc jack functions.

---
 latlib                       |  2 +-
 u1casc-ordinary/obs_phi2.hpp | 32 +++++++++++++++++++-------------
 u1casc-ordinary/obs_plaq.hpp | 32 +++++++++++++++++++-------------
 3 files changed, 39 insertions(+), 27 deletions(-)

diff --git a/latlib b/latlib
index 33ee1f8..e9c0e0b 160000
--- a/latlib
+++ b/latlib
@@ -1 +1 @@
-Subproject commit 33ee1f81e21dbb384b48c8083aff790ca313558e
+Subproject commit e9c0e0b98886738917300348a522759c25d4a9c7
diff --git a/u1casc-ordinary/obs_phi2.hpp b/u1casc-ordinary/obs_phi2.hpp
index 9afc390..3e4fe48 100644
--- a/u1casc-ordinary/obs_phi2.hpp
+++ b/u1casc-ordinary/obs_phi2.hpp
@@ -27,7 +27,9 @@ private:
   obsmem* OM;
 
   void phi2Compute();
-  static double phi2Sus(vector< vector<double> > *vals, void *para);
+
+  static void prePhi2Sus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para);
+  static double phi2Sus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para);
 
   sim *Sim;
   obstat<double,double> oPhi2[2];
@@ -61,7 +63,7 @@ void obs_phi2::_finish() {
 
   for (int iflav=0; iflav<2; iflav++) {
     compid_phi2[iflav] = oPhi2[iflav].computeEasy();
-    compid_phi2sus[iflav] = oPhi2[iflav].computeJack(obs_phi2::phi2Sus, &(Sim->lsize4));
+    compid_phi2sus[iflav] = oPhi2[iflav].computeJack(obs_phi2::prePhi2Sus, obs_phi2::phi2Sus, &(Sim->lsize4));
 
     *out << "\t" << oPhi2[iflav].getMean(compid_phi2[iflav]) << "\t" << oPhi2[iflav].getErr(compid_phi2[iflav]);
     *out << "\t" << oPhi2[iflav].getMean(compid_phi2sus[iflav]) << "\t" << oPhi2[iflav].getErr(compid_phi2sus[iflav]);
@@ -84,17 +86,21 @@ void obs_phi2::phi2Compute()
   }
 }
 
-double obs_phi2::phi2Sus(vector< vector<double> > *vals, void *para) {
-  double mean=0, mean2=0;
-  for(vector< vector<double> >::iterator valIt = vals->begin(); valIt != vals->end(); ++valIt)  
-    {
-      mean += (*valIt)[0];
-      mean2 += pow((*valIt)[0],2);
-    }
-  mean /= vals->size();
-  mean2 /= vals->size();
-  
-  return ( mean2 - pow(mean,2) ) * *(int*)para;
+void obs_phi2::prePhi2Sus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para) {
+  preCalculated->push_back(0);
+  preCalculated->push_back(0);
+
+  for(vector< vector<double> >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) {
+    (*preCalculated)[0] += (*valIt)[0];
+    (*preCalculated)[1] += pow((*valIt)[0],2);
+  }
+}
+
+double obs_phi2::phi2Sus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para) {
+  return (( (*preCalculated)[1] - pow((*excludedMeas)[0], 2) )
+	  - pow( (*preCalculated)[0] - (*excludedMeas)[0], 2 ))
+    * *(int*)para
+    / (nmeas-1);
 }
 
 #endif
diff --git a/u1casc-ordinary/obs_plaq.hpp b/u1casc-ordinary/obs_plaq.hpp
index 5dacf68..2f09096 100644
--- a/u1casc-ordinary/obs_plaq.hpp
+++ b/u1casc-ordinary/obs_plaq.hpp
@@ -27,7 +27,9 @@ private:
   obsmem* OM;
 
   void plaqCompute();
-  static double plaqSus(vector< vector<double> > *vals, void *para);
+
+  static void prePlaqSus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para);
+  static double plaqSus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para);
 
   sim *Sim;
   obstat<double,double> oPlaq;
@@ -57,7 +59,7 @@ void obs_plaq::_finish() {
   int compid_plaq, compid_plaqsus;
 
   compid_plaq = oPlaq.computeEasy();
-  compid_plaqsus = oPlaq.computeJack(obs_plaq::plaqSus, &(Sim->lsize4));
+  compid_plaqsus = oPlaq.computeJack(obs_plaq::prePlaqSus, obs_plaq::plaqSus, &(Sim->lsize4));
 
   *out << "\t" << oPlaq.getMean(compid_plaq) << "\t" << oPlaq.getErr(compid_plaq);
   *out << "\t" << oPlaq.getMean(compid_plaqsus) << "\t" << oPlaq.getErr(compid_plaqsus);
@@ -81,17 +83,21 @@ void obs_plaq::plaqCompute()
   OM->plaq /= Sim->lsize4*6;
 }
 
-double obs_plaq::plaqSus(vector< vector<double> > *vals, void *para) {
-  double mean=0, mean2=0;
-  for(vector< vector<double> >::iterator valIt = vals->begin(); valIt != vals->end(); ++valIt)  
-    {
-      mean += (*valIt)[0];
-      mean2 += pow((*valIt)[0],2);
-    }
-  mean /= vals->size();
-  mean2 /= vals->size();
-  
-  return ( mean2 - pow(mean,2) ) * *(int*)para * 6;
+void obs_plaq::prePlaqSus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para) {
+  preCalculated->push_back(0);
+  preCalculated->push_back(0);
+
+  for(vector< vector<double> >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) {
+    (*preCalculated)[0] += (*valIt)[0];
+    (*preCalculated)[1] += pow((*valIt)[0],2);
+  }
+}
+
+double obs_plaq::plaqSus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para) {
+  return (( (*preCalculated)[1] - pow((*excludedMeas)[0], 2) )
+	  - pow( (*preCalculated)[0] - (*excludedMeas)[0], 2 ))
+    * *(int*)para * 6
+    / (nmeas-1);
 }
 
 #endif
-- 
2.39.5