From 86f08c6b9ae1d7d1f31a7ca9d5bebba153e3ff62 Mon Sep 17 00:00:00 2001
From: Christopher Kelly <ckelly@bnl.gov>
Date: Tue, 18 May 2021 13:57:44 -0400
Subject: [PATCH] Added a check that the initial EOFA action agrees with
 |eta|^2, thus checking the quality of the rational approximation in the
 heatbath

---
 Grid/qcd/action/ActionParams.h                |  9 ++++---
 .../pseudofermion/ExactOneFlavourRatio.h      | 26 ++++++++++++++++++-
 tests/forces/Test_mobius_force_eofa.cc        | 15 -----------
 3 files changed, 31 insertions(+), 19 deletions(-)

diff --git a/Grid/qcd/action/ActionParams.h b/Grid/qcd/action/ActionParams.h
index 12d2b91b..e4b47664 100644
--- a/Grid/qcd/action/ActionParams.h
+++ b/Grid/qcd/action/ActionParams.h
@@ -66,7 +66,8 @@ struct StaggeredImplParams {
 				    RealD, tolerance, 
 				    int,   degree, 
 				    int,   precision,
-				    int,   BoundsCheckFreq);
+				    int,   BoundsCheckFreq,
+				    RealD, BoundsCheckTol);
     
   // MaxIter and tolerance, vectors??
     
@@ -77,14 +78,16 @@ struct StaggeredImplParams {
 				RealD tol      = 1.0e-8, 
                            	int _degree    = 10,
 				int _precision = 64,
-				int _BoundsCheckFreq=20)
+				int _BoundsCheckFreq=20,
+				double _BoundsCheckTol=1e-6)
       : lo(_lo),
 	hi(_hi),
 	MaxIter(_maxit),
 	tolerance(tol),
 	degree(_degree),
         precision(_precision),
-        BoundsCheckFreq(_BoundsCheckFreq){};
+        BoundsCheckFreq(_BoundsCheckFreq),
+        BoundsCheckTol(_BoundsCheckTol){};
   };
 
 
diff --git a/Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h b/Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h
index bea67415..800866a9 100644
--- a/Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h
+++ b/Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h
@@ -64,6 +64,8 @@ NAMESPACE_BEGIN(Grid);
       SchurRedBlackDiagMooeeSolve<FermionField> DerivativeSolverR;
       FermionField Phi; // the pseudofermion field for this trajectory
 
+      RealD norm2_eta; //|eta|^2 where eta is the random gaussian field used to generate the pseudofermion field
+      bool initial_action; //true for the first call to S after refresh, for which the identity S = |eta|^2 holds provided the rational approx is good
     public:
 
       ExactOneFlavourRatioPseudoFermionAction(AbstractEOFAFermion<Impl>& _Lop, 
@@ -87,7 +89,8 @@ NAMESPACE_BEGIN(Grid);
 	DerivativeSolverL(DerivCGL, false, true), DerivativeSolverR(DerivCGR, false, true), 
 	Phi(_Lop.FermionGrid()), 
 	param(p), 
-        use_heatbath_forecasting(use_fc)
+	use_heatbath_forecasting(use_fc),
+	initial_action(false)
       {
         AlgRemez remez(param.lo, param.hi, param.precision);
 
@@ -216,8 +219,14 @@ NAMESPACE_BEGIN(Grid);
         Lop.RefreshShiftCoefficients(0.0);
         Rop.RefreshShiftCoefficients(-1.0);
 
+	//Mark that the next call to S is the first after refresh
+	initial_action = true;
+
+
 	// Bounds check
 	RealD EtaDagEta = norm2(eta);
+	norm2_eta = EtaDagEta;
+
 	//	RealD PhiDagMPhi= norm2(eta);
 
       };
@@ -290,6 +299,21 @@ NAMESPACE_BEGIN(Grid);
         Rop.Omega(tmp[1], tmp[0], 1, 1);
         action += Rop.k * innerProduct(spProj_Phi, tmp[0]).real();
 
+	if(initial_action){
+	  //For the first call to S after refresh,  S = |eta|^2. We can use this to ensure the rational approx is good
+	  RealD diff = action - norm2_eta;
+
+	  //S_init = eta^dag M^{-1/2} M M^{-1/2} eta
+	  //S_init - eta^dag eta =  eta^dag ( M^{-1/2} M M^{-1/2} - 1 ) eta
+	  RealD test = sqrt(fabs(diff)/norm2_eta); //test the quality of the rational approx
+
+	  std::cout << GridLogMessage << action_name() << " initial action " << action << " expect " << norm2_eta << "; diff " << diff << std::endl;
+	  std::cout << GridLogMessage << action_name() << " sqrt( eta^dag ( M^{-1/2} M M^{-1/2} - 1 ) eta )/sqrt( eta^dag eta ) = " << test << "  expect 0 (tol " << param.BoundsCheckTol << ")" << std::endl;
+
+	  assert( ( test < param.BoundsCheckTol ) && " Initial action check failed" );
+	  initial_action = false;
+	}
+
         return action;
       };
 
diff --git a/tests/forces/Test_mobius_force_eofa.cc b/tests/forces/Test_mobius_force_eofa.cc
index f7c25006..1d25771a 100644
--- a/tests/forces/Test_mobius_force_eofa.cc
+++ b/tests/forces/Test_mobius_force_eofa.cc
@@ -132,21 +132,6 @@ int main (int argc, char** argv)
 
   Meofa.refresh(U, sRNG, RNG5 );
 
-
-    
-
-
-
-
-
-
-
-
-
-
-
-
-
   RealD S = Meofa.S(U); // pdag M p
 
   // get the deriv of phidag M phi with respect to "U"