Cleaning up and aligning variable naming between action deriv versions

Grid/qcd/action/gauge/PlaqPlusRectangleAction.h

@@ -71,30 +71,30 @@ public:
     return action;
   };
 
-  virtual void deriv(const GaugeField &Umu,GaugeField & dSdU) {
+  virtual void deriv(const GaugeField &U, GaugeField &dSdU) {
-    //extend Ta to include Lorentz indexes
+    // extend Ta to include Lorentz indexes
-    RealD factor_p = c_plaq/RealD(Nc)*0.5;
+    RealD factor_p = c_plaq / RealD(Nc) * 0.5;
-    RealD factor_r = c_rect/RealD(Nc)*0.5;
+    RealD factor_r = c_rect / RealD(Nc) * 0.5;
 
-    GridBase *grid = Umu.Grid();
+    GridBase *grid = U.Grid();
 
-    std::vector<GaugeLinkField> U (Nd,grid);
+    std::vector<GaugeLinkField> Umu(Nd, grid);
-    for(int mu=0;mu<Nd;mu++){
+    for (int mu = 0; mu < Nd; mu++) {
-      U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
+      Umu[mu] = PeekIndex<LorentzIndex>(U, mu);
     }
-    std::vector<GaugeLinkField> RectStaple(Nd,grid), Staple(Nd,grid);
+    std::vector<GaugeLinkField> RectStaple(Nd, grid), Staple(Nd, grid);
-    WilsonLoops<Gimpl>::StapleAndRectStapleAll(Staple, RectStaple, U, workspace);
+    WilsonLoops<Gimpl>::StapleAndRectStapleAll(Staple, RectStaple, Umu,
+                                               workspace);
 
     GaugeLinkField dSdU_mu(grid);
     GaugeLinkField staple(grid);
 
-    for (int mu=0; mu < Nd; mu++){
+    for (int mu = 0; mu < Nd; mu++) {
-      dSdU_mu = Ta(U[mu]*Staple[mu])*factor_p;
+      dSdU_mu = Ta(Umu[mu] * Staple[mu]) * factor_p;
-      dSdU_mu = dSdU_mu + Ta(U[mu]*RectStaple[mu])*factor_r;
+      dSdU_mu = dSdU_mu + Ta(Umu[mu] * RectStaple[mu]) * factor_r;
 
       PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
     }
-
   };
 
 };

Grid/qcd/action/gauge/WilsonGaugeAction.h

@@ -68,27 +68,23 @@ public:
     // extend Ta to include Lorentz indexes
 
     RealD factor = 0.5 * beta / RealD(Nc);
+    GridBase *grid = U.Grid();
 
-    // GaugeLinkField Umu(U.Grid());
+    GaugeLinkField dSdU_mu(grid);
-    std::vector<GaugeLinkField> Uvec(Nd, U.Grid());
+    std::vector<GaugeLinkField> Umu(Nd, grid);
-    for (int d = 0; d < Nd; d++) {
+    for (int mu = 0; mu < Nd; mu++) {
-      Uvec[d] = PeekIndex<LorentzIndex>(U, d);
+      Umu[mu] = PeekIndex<LorentzIndex>(U, mu);
     }
 
-    GaugeLinkField dSdU_mu(U.Grid());
     for (int mu = 0; mu < Nd; mu++) {
-
-      // Umu = PeekIndex<LorentzIndex>(U, mu);
-
       // Staple in direction mu
-      // WilsonLoops<Gimpl>::Staple(dSdU_mu, U, mu);
+      WilsonLoops<Gimpl>::Staple(dSdU_mu, Umu, mu);
-      // dSdU_mu = Ta(Umu * dSdU_mu) * factor;
+      dSdU_mu = Ta(Umu[mu] * dSdU_mu) * factor;
-      WilsonLoops<Gimpl>::Staple(dSdU_mu, Uvec, mu);
-      dSdU_mu = Ta(Uvec[mu] * dSdU_mu) * factor;
 
       PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
     }
   }
+
 private:
   RealD beta;  
  };

Grid/qcd/utils/WilsonLoops.h

@@ -292,18 +292,16 @@ public:
   //////////////////////////////////////////////////
   // the sum over all nu-oriented staples for nu != mu on each site
   //////////////////////////////////////////////////
-  static void Staple(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
+  static void Staple(GaugeMat &staple, const GaugeLorentz &U, int mu) {
 
-    GridBase *grid = Umu.Grid();
+    std::vector<GaugeMat> Umu(Nd, U.grid());
-
-    std::vector<GaugeMat> U(Nd, grid);
     for (int d = 0; d < Nd; d++) {
-      U[d] = PeekIndex<LorentzIndex>(Umu, d);
+      Umu[d] = PeekIndex<LorentzIndex>(U, d);
     }
-    Staple(staple, U, mu);
+    Staple(staple, Umu, mu);
   }
 
-  static void Staple(GaugeMat &staple, const std::vector<GaugeMat> &U, int mu) {
+  static void Staple(GaugeMat &staple, const std::vector<GaugeMat> &Umu, int mu) {
     staple = Zero();
 
     for (int nu = 0; nu < Nd; nu++) {
@@ -320,11 +318,8 @@ public:
         //
 
         staple += Gimpl::ShiftStaple(
-				     Gimpl::CovShiftForward(
+          Gimpl::CovShiftForward(Umu[nu], nu,
-							    U[nu], nu,
+            Gimpl::CovShiftBackward(Umu[mu], mu, Gimpl::CovShiftIdentityBackward(Umu[nu], nu))), mu);
-							    Gimpl::CovShiftBackward(
-										    U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))),
-				     mu);
 
         //  __
         // |
@@ -333,8 +328,8 @@ public:
         //
 
         staple += Gimpl::ShiftStaple(
-				     Gimpl::CovShiftBackward(U[nu], nu,
+          Gimpl::CovShiftBackward(Umu[nu], nu,
-							     Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu);
+            Gimpl::CovShiftBackward(Umu[mu], mu, Umu[nu])), mu);
       }
     }
   }

tests/run_regression_test.py

@@ -0,0 +1,134 @@
+def read_expected(test_name="Test_hmc_Sp_WilsonFundFermionGauge", grid="8.8.8.8", mpi="1.1.1.1"):
+    """
+    Read expected values from file.
+
+    The file contains one or more entries of the following format:
+    <grid> <mpi> <plaquette> <checksum_rng> <checksum_lat>
+    Eg.
+    8.8.8.8 1.1.1.1 0.0256253844 922c392f d1e4cc1c
+    """
+
+    with open(f"{test_name}_expected.txt") as file:
+        for line in file:
+            line_split = line.split()
+            if line_split[0] == grid and line_split[1] == mpi:
+                return float(line_split[2]), line_split[3], line_split[4]
+
+
+def read_output():
+    """
+    Read test output and fish out values of interest.
+    """
+
+    checksum_rng = None
+    checksum_lat = None
+    plaquette = None
+    with open("output.txt", 'r') as file:
+        for line in file:
+            if "Written NERSC" in line:
+                subline = line.split('checksum ')[1]
+                if len(subline.split()) == 1: # this is the rng checksum line
+                    checksum_rng = subline.strip()
+                elif len(subline.split()) == 3: # this is the lat checksum and plaquette value line
+                    checksum_lat = subline.split()[0]
+                    plaquette = float(subline.split()[2])
+                else:
+                    print("Picked wrong line...")
+
+    if (checksum_rng is None) or (checksum_lat is None) or (plaquette is None):
+        print("Error reading values from output file. Make sure you compile the test with CPparams.saveInterval=1 in order to produce the required output.")
+        exit(1)
+        
+    return plaquette, checksum_rng, checksum_lat
+
+
+def compare(actual, expected, what, stop=False):
+    """
+    Compare actual with expected output, and output message if failed.
+    """
+
+    if actual != expected:
+        print(f"{what} comparison failed: actual={actual} , expected={expected}")
+        if stop:
+            exit(1)
+        else:
+            return False
+    return True
+
+
+
+if __name__ == '__main__':
+    import argparse
+    import subprocess
+    import os
+
+    parser = argparse.ArgumentParser(description='Run end-to-end tests and compare results with expectations.')
+    parser.add_argument("test_name", help="File name of the test")
+    parser.add_argument("grid", help="Grid configuration")
+    parser.add_argument("mpi", help="MPI configuration")
+    parser.add_argument("-s", "--stop", action='store_true', help="Flag to stop testing when a test fails.")
+    args = parser.parse_args()
+
+    expected_plaquette, expected_checksum_rng, expected_checksum_lat = read_expected(args.test_name, args.grid, args.mpi)
+
+    result = subprocess.run([f"./{args.test_name} --grid {args.grid} --mpi {args.mpi} --Thermalizations 0 --Trajectories 1 > output.txt"], shell=True, encoding="text")
+    plaquette, checksum_rng, checksum_lat = read_output()
+
+    print(f"Running {args.test_name}")
+    result = compare(plaquette, expected_plaquette, "plaquette", args.stop)
+    result = result and compare(checksum_rng, expected_checksum_rng, "Checksum RNG file ", args.stop)
+    result = result and compare(checksum_lat, expected_checksum_lat, "Checksum LAT file ", args.stop)
+    if result:
+        print("All tests passed!")
+    else:
+        print("Some tests failed...")
+
+    os.remove("output.txt")
+    os.remove("ckpoint_rng.1")
+    os.remove("ckpoint_lat.1")
+
+#result = subprocess.run(["./Test_hmc_Sp_WilsonFundFermionGauge --grid 8.8.8.8 --mpi 1.1.1.1 --Thermalizations 0 --Trajectories 1 > output1.txt"], shell=True, encoding="text")
+
+# expected_value = 0.0256253844
+# checksum_rng = "922c392f"
+# checksum_lat = "d1e4cc1c"
+
+# with open("output1.txt", 'r') as file:
+#     for line in file:
+#         # if "Plaquette" in line:
+#         #     #print(line)
+#         #     plaquette_value = float(line.split('] ')[1])
+#         #     #print(plaquette_value)
+#         #     if plaquette_value == expected_value:
+#         #         print("Success!")
+#         if "Written NERSC" in line:
+#             print(line)
+#             subline = line.split('checksum ')[1]
+#             if len(subline.split()) == 1: # this is the rng checksum line
+#                 print(subline)
+#                 if subline.strip() == checksum_rng:
+#                     print("RNG file checksum success!")
+#                 else:
+#                     print("RNG file checksum failed!")
+#             elif len(subline.split()) == 3: # this is the lat checksum and plaquette value line
+#                 print(subline)
+#                 checksum_value = subline.split()[0]
+#                 plaquette_value = float(subline.split()[2])
+#                 print(checksum_value, plaquette_value)
+#                 if checksum_value == checksum_lat:
+#                     print("LAT file checksum success!")
+#                 else:
+#                     print("LAT file checksum failed!")
+#                 if plaquette_value == expected_value:
+#                     print("Plaquette value success!")
+#                 else:
+#                     print("Plaquette value failed!")
+#             else:
+#                 print("Picked wrong line...")
+
+
+#loc1 = result.find("Plaquette")
+#print(loc1)
+#loc2 = result.find("Smeared")
+#print(loc2)
+#print(result[loc1,loc2])

tests/sp2n/Test_hmc_Sp_WilsonFundFermionGauge_expected.txt

@@ -0,0 +1 @@
+8.8.8.8 1.1.1.1 0.0256253844 922c392f d1e4cc1c