Moving parameters outside of the HMCrunner

lib/qcd/hmc/GenericHMCrunner.h

@@ -7,6 +7,7 @@ Source file: ./lib/qcd/hmc/GenericHmcRunner.h
 Copyright (C) 2015
 
 Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -49,17 +50,39 @@ class BinaryHmcRunnerTemplate {
   // that can be injected from outside  
   std::vector< HmcObservable<typename Implementation::Field>* > ObservablesList;
 
+  IntegratorParameters MDparameters;
 
   GridCartesian *UGrid;
   GridCartesian *FGrid;
   GridRedBlackCartesian *UrbGrid;
   GridRedBlackCartesian *FrbGrid;
 
+  std::vector<int> SerialSeed;
+  std::vector<int> ParallelSeed;
+
+  void RNGSeeds(std::vector<int> S, std::vector<int> P){
+    SerialSeed = S;
+    ParallelSeed = P;
+  }
+
   virtual void BuildTheAction(int argc, char **argv) = 0;  // necessary?
 
-// add here the smearing implementation?
+  // A couple of wrapper classes
-template <class IOCheckpointer = BinaryHmcCheckpointer<Implementation> >
+  template <class IOCheckpointer>
   void Run(int argc, char **argv, IOCheckpointer &Checkpoint) {
+    NoSmearing<Implementation> S;
+    Runner(argc, argv, Checkpoint, S);
+  }
+
+  template <class IOCheckpointer, class SmearingPolicy>
+  void Run(int argc, char **argv, IOCheckpointer &CP, SmearingPolicy &S) {
+    Runner(argc, argv, CP, S);
+  }
+  //////////////////////////////
+
+  template <class SmearingPolicy, class IOCheckpointer >
+  void Runner(int argc, char **argv, IOCheckpointer &Checkpoint,
+           SmearingPolicy &Smearing) {
     StartType_t StartType = HotStart;
 
     std::string arg;
@@ -111,17 +134,10 @@ template <class IOCheckpointer = BinaryHmcCheckpointer<Implementation> >
     GridParallelRNG pRNG(UGrid);
     Field U(UGrid);
 
-    // This outside
+    typedef MinimumNorm2<Implementation, SmearingPolicy,
-    std::vector<int> SerSeed({1, 2, 3, 4, 5});
-    std::vector<int> ParSeed({6, 7, 8, 9, 10});
-
-    // these decisions outside
-    NoSmearing<Implementation> SmearingPolicy;
-    typedef MinimumNorm2<Implementation, NoSmearing<Implementation>,
                          RepresentationsPolicy>
         IntegratorType;  // change here to change the algorithm
-    IntegratorParameters MDpar(20, 1.0);
+    IntegratorType MDynamics(UGrid, MDparameters, TheAction, Smearing);
-    IntegratorType MDynamics(UGrid, MDpar, TheAction, SmearingPolicy);
 
     HMCparameters HMCpar;
     HMCpar.StartTrajectory = StartTraj;
@@ -131,20 +147,20 @@ template <class IOCheckpointer = BinaryHmcCheckpointer<Implementation> >
     if (StartType == HotStart) {
       // Hot start
       HMCpar.MetropolisTest = true;
-      sRNG.SeedFixedIntegers(SerSeed);
+      sRNG.SeedFixedIntegers(SerialSeed);
-      pRNG.SeedFixedIntegers(ParSeed);
+      pRNG.SeedFixedIntegers(ParallelSeed);
       Implementation::HotConfiguration(pRNG, U);
     } else if (StartType == ColdStart) {
       // Cold start
       HMCpar.MetropolisTest = true;
-      sRNG.SeedFixedIntegers(SerSeed);
+      sRNG.SeedFixedIntegers(SerialSeed);
-      pRNG.SeedFixedIntegers(ParSeed);
+      pRNG.SeedFixedIntegers(ParallelSeed);
       Implementation::ColdConfiguration(pRNG, U);
     } else if (StartType == TepidStart) {
       // Tepid start
       HMCpar.MetropolisTest = true;
-      sRNG.SeedFixedIntegers(SerSeed);
+      sRNG.SeedFixedIntegers(SerialSeed);
-      pRNG.SeedFixedIntegers(ParSeed);
+      pRNG.SeedFixedIntegers(ParallelSeed);
       Implementation::TepidConfiguration(pRNG, U);
     } else if (StartType == CheckpointStart) {
       HMCpar.MetropolisTest = true;
@@ -152,10 +168,9 @@ template <class IOCheckpointer = BinaryHmcCheckpointer<Implementation> >
       Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG);
     }
 
-    SmearingPolicy.set_Field(U);
+    Smearing.set_Field(U);
 
     HybridMonteCarlo<IntegratorType> HMC(HMCpar, MDynamics, sRNG, pRNG, U);
-    //HMC.AddObservable(&Checkpoint);
 
     for (int obs = 0; obs < ObservablesList.size(); obs++)
       HMC.AddObservable(ObservablesList[obs]); 
@@ -170,16 +185,12 @@ typedef BinaryHmcRunnerTemplate<PeriodicGimplR> BinaryHmcRunner;
 typedef BinaryHmcRunnerTemplate<PeriodicGimplF> BinaryHmcRunnerF;
 typedef BinaryHmcRunnerTemplate<PeriodicGimplD> BinaryHmcRunnerD;
 
-//typedef BinaryHmcRunnerTemplate<PeriodicGimplR, NerscHmcCheckpointer<PeriodicGimplR> > NerscTestHmcRunner;
-
-
 template <class RepresentationsPolicy>
 using BinaryHmcRunnerTemplateHirep =
     BinaryHmcRunnerTemplate<PeriodicGimplR, RepresentationsPolicy>;
 
 
 
-//typedef BinaryHmcRunnerTemplate<ScalarImplR, BinaryHmcCheckpointer<ScalarImplR>, ScalarFields> ScalarBinaryHmcRunner;
     typedef BinaryHmcRunnerTemplate<ScalarImplR, ScalarFields> ScalarBinaryHmcRunner;
 }
 }

lib/qcd/hmc/integrators/Integrator.h

@@ -46,27 +46,28 @@ namespace Grid {
 namespace QCD {
 
 struct IntegratorParameters {
-  unsigned int
-      Nexp;  // number of terms in the Taylor expansion of the exponential
   unsigned int MDsteps;  // number of outer steps
   RealD trajL;           // trajectory length
   RealD stepsize;        // trajectory stepsize
 
-  IntegratorParameters(int MDsteps_, RealD trajL_ = 1.0,
+  IntegratorParameters(int MDsteps_ = 10, RealD trajL_ = 1.0)
-                       unsigned int Nexp_ = 12)
+      : MDsteps(MDsteps_),
-      : Nexp(Nexp_),
-        MDsteps(MDsteps_),
         trajL(trajL_),
         stepsize(trajL / MDsteps){
             // empty body constructor
   };
 
+  void set(int MDsteps_, RealD trajL_){
+    MDsteps = MDsteps_;
+    trajL = trajL_;
+    stepsize = trajL/MDsteps;
+  }
+
+
   void print_parameters() {
     std::cout << GridLogMessage << "[Integrator] Trajectory length  : " << trajL << std::endl;
     std::cout << GridLogMessage << "[Integrator] Number of MD steps : " << MDsteps << std::endl;
     std::cout << GridLogMessage << "[Integrator] Step size          : " << stepsize << std::endl;
-    std::cout << GridLogMessage << "[Integrator] Exponential approx.: " << Nexp << std::endl;
-
   }
 };
 

lib/qcd/smearing/APEsmearing.h

@@ -92,19 +92,19 @@
 
 	  						temp_Sigma = -rho_numu*staple*iLambda_nu;  //ok
 	        				//-r_numu*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)*Lambda_nu(x)
-	  						Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
+	  						Gimpl::AddLink(SigmaTerm, temp_Sigma, mu);
 
 	    					sh_field = Cshift(iLambda_nu, mu, 1);// general also for Gparity?
 
 	    					temp_Sigma = rho_numu*sh_field*staple; //ok
 	    					//r_numu*Lambda_nu(mu)*U_nu(x+mu)*Udag_mu(x+nu)*Udag_nu(x)
-	    					Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
+	    					Gimpl::AddLink(SigmaTerm, temp_Sigma, mu);
 
 	    					sh_field = Cshift(iLambda_mu, nu, 1);
 
 	    					temp_Sigma = -rho_munu*staple*U_nu*sh_field*adj(U_nu); //ok
 	    					//-r_munu*U_nu(x+mu)*Udag_mu(x+nu)*Lambda_mu(x+nu)*Udag_nu(x)
-	    					Gimpl::AddGaugeLink(SigmaTerm, temp_Sigma, mu);
+	    					Gimpl::AddLink(SigmaTerm, temp_Sigma, mu);
 
 	    					staple = zero;
 	    					sh_field = Cshift(U_nu, mu, 1);
@@ -116,7 +116,7 @@
 	    					sh_field = Cshift(u_tmp, mu, 1);
 	    					temp_Sigma += -rho_numu*sh_field*adj(U_mu)*U_nu;
 	    					sh_field = Cshift(temp_Sigma, nu, -1);
-	    					Gimpl::AddGaugeLink(SigmaTerm, sh_field, mu);
+	    					Gimpl::AddLink(SigmaTerm, sh_field, mu);
 
 	    				}
 	    			}

tests/hmc/Test_hmc_WilsonFermionGauge_Binary.cc

@@ -49,6 +49,8 @@ class HMCRunnerParameters : Serializable {
                                   std::string, format,
                                   std::string, conf_prefix,
                                   std::string, rng_prefix,
+                                  double, rho,
+                                  int, SmearingLevels,
                                   );
 
   HMCRunnerParameters() {}
@@ -115,11 +117,21 @@ class HmcRunner : public BinaryHmcRunner {
     // Then force all checkpoint to have few common functions
     // return an object that is then passed to the Run function
 
-    PlaquetteLogger<BinaryHmcRunner::ImplPolicy> PlaqLog(std::string("Plaquette"));
+    PlaquetteLogger<BinaryHmcRunner::ImplPolicy> PlaqLog(
+        std::string("Plaquette"));
     ObservablesList.push_back(&PlaqLog);
     ObservablesList.push_back(&Checkpoint);
 
-    Run(argc, argv, Checkpoint);
+    // Smearing section, omit if not needed
+    double rho = 0.1;  // smearing parameter
+    int Nsmear = 2;    // number of smearing levels
+    Smear_Stout<BinaryHmcRunner::ImplPolicy> Stout(rho);
+    SmearedConfiguration<BinaryHmcRunner::ImplPolicy> SmearingPolicy(
+        UGrid, Nsmear, Stout);
+    ///////////////////
+
+    // Run(argc, argv, Checkpoint, SmearingPolicy);
+    Run(argc, argv, Checkpoint);  // no smearing
   };
 };
 }
@@ -134,5 +146,12 @@ int main(int argc, char **argv) {
 
   HmcRunner TheHMC;
 
+  // Seeds for the random number generators
+  std::vector<int> SerSeed({1, 2, 3, 4, 5});
+  std::vector<int> ParSeed({6, 7, 8, 9, 10});
+  TheHMC.RNGSeeds(SerSeed, ParSeed);
+
+  TheHMC.MDparameters.set(20, 1.0);// MDsteps, traj length
+
   TheHMC.BuildTheAction(argc, argv);
 }