From 68fe0769a11a15e0f6ac19ab79c53d2a79cd9b36 Mon Sep 17 00:00:00 2001
From: neo <cossu@post.kek.jp>
Date: Mon, 6 Jul 2015 16:17:32 +0900
Subject: [PATCH] Added minimum norm integrator

Little rearrangement of HMC and integrator classes
---
 lib/qcd/hmc/HMC.h                         |  15 +--
 lib/qcd/hmc/integrators/Integrator_base.h | 111 +++++++++++++++++-----
 tests/Test_hmc_WilsonGauge.cc             |   7 +-
 3 files changed, 94 insertions(+), 39 deletions(-)

diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h
index 01fb2f9c..f2bb4bcb 100644
--- a/lib/qcd/hmc/HMC.h
+++ b/lib/qcd/hmc/HMC.h
@@ -10,7 +10,6 @@
 
 #include <string>
 #include <memory>
-#include <unistd.h>
 
 namespace Grid{
   namespace QCD{
@@ -30,10 +29,8 @@ namespace Grid{
     class HybridMonteCarlo{
       const HMCparameters Params;
       GridSerialRNG sRNG;
-      GridParallelRNG pRNG;
       Integrator<IntegType>& MD;
       
-
       bool metropolis_test(const RealD DeltaH){
 	RealD rn_test;
 	RealD prob = std::exp(-DeltaH);
@@ -54,7 +51,7 @@ namespace Grid{
 
       RealD evolve_step(LatticeLorentzColourMatrix& U){
 
-	MD.init(U,pRNG);     // set U and initialize P and phi's 
+	MD.init(U);     // set U and initialize P and phi's 
 	RealD H0 = MD.S(U);     // current state            
 	std::cout<<"Total H before = "<< H0 << "\n";
       
@@ -70,14 +67,12 @@ namespace Grid{
       
     public:
     HybridMonteCarlo(HMCparameters Pms, 
-		     Integrator<IntegType>& MolDyn, 
-		     GridBase* grid):
-      Params(Pms),MD(MolDyn),pRNG(grid){
+		     Integrator<IntegType>& MolDyn):
+      Params(Pms),MD(MolDyn){
 	//FIXME
 
-	// initialize RNGs
+	// initialize RNGs also with seed
 	sRNG.SeedRandomDevice();
-	pRNG.SeedRandomDevice();
       }
       ~HybridMonteCarlo(){};
 
@@ -85,7 +80,7 @@ namespace Grid{
 
       void evolve(LatticeLorentzColourMatrix& Uin){
 	Real DeltaH;
-      
+	
 	// Thermalizations
 	for(int iter=1; iter <= Params.ThermalizationSteps; ++iter){
 	  std::cout << "-- # Thermalization step = "<< iter <<  "\n";
diff --git a/lib/qcd/hmc/integrators/Integrator_base.h b/lib/qcd/hmc/integrators/Integrator_base.h
index 4891eb5e..9b47a19d 100644
--- a/lib/qcd/hmc/integrators/Integrator_base.h
+++ b/lib/qcd/hmc/integrators/Integrator_base.h
@@ -24,9 +24,6 @@ namespace Grid{
     typedef std::vector<ActionLevel> ActionSet;
     typedef std::vector<Observer*> ObserverList;
     
-    class LeapFrog;
-   
-
     struct IntegratorParameters{
       int Nexp;
       int MDsteps;  // number of outer steps
@@ -46,20 +43,20 @@ namespace Grid{
     }
 
     
-    
-
-
     template< class IntegratorPolicy >
     class Integrator{
     private:
       IntegratorParameters Params;
       const ActionSet as;
       const std::vector<int> Nrel; //relative step size per level
-      //ObserverList observers; // not yet
       std::unique_ptr<LatticeLorentzColourMatrix> P;
-      
-      IntegratorPolicy TheIntegrator;// contains parameters too
+      GridParallelRNG pRNG;
+      //ObserverList observers; // not yet
+     
+      IntegratorPolicy TheIntegrator;
 
+      void register_observers();
+      void notify_observers();
 
       void update_P(LatticeLorentzColourMatrix&U, int level,double ep){
 	for(int a=0; a<as[level].size(); ++a){
@@ -77,36 +74,32 @@ namespace Grid{
 	for (int mu = 0; mu < Nd; mu++){
 	  Umu=peekLorentz(U, mu);
 	  Pmu=peekLorentz(*P, mu);
-	  Umu = expMat(Pmu, ep)*Umu;
+	  Umu = expMat(Pmu, ep, Params.Nexp)*Umu;
 	  pokeLorentz(U, Umu, mu);
 	}
 
       }
       
-      void register_observers();
-      void notify_observers();
+
       
       friend void IntegratorPolicy::step (LatticeLorentzColourMatrix& U, 
 				     int level, std::vector<int>& clock,
-				     Integrator<LeapFrog>* Integ);
+				     Integrator<IntegratorPolicy>* Integ);
     public:
-      Integrator(IntegratorParameters Par,
+    Integrator(GridBase* grid, IntegratorParameters Par,
 		 ActionSet& Aset, std::vector<int> Nrel_):
-	Params(Par),as(Aset),Nrel(Nrel_){
+      Params(Par),as(Aset),Nrel(Nrel_),P(new LatticeLorentzColourMatrix(grid)),pRNG(grid){
 	assert(as.size() == Nrel.size());
+	pRNG.SeedRandomDevice();
       };
       
       ~Integrator(){}
 
 
       //Initialization of momenta and actions
-      void init(LatticeLorentzColourMatrix& U,
-		GridParallelRNG& pRNG){
+      void init(LatticeLorentzColourMatrix& U){
 	std::cout<< "Integrator init\n";
-	if (!P){
-	  std::unique_ptr<LatticeLorentzColourMatrix> Pnew(new LatticeLorentzColourMatrix(U._grid));	  
-	  P = std::move(Pnew);
-	}
+
 	MDutils::generate_momenta(*P,pRNG);
 	for(int level=0; level< as.size(); ++level){
 	  for(int actionID=0; actionID<as.at(level).size(); ++actionID){
@@ -116,7 +109,7 @@ namespace Grid{
       }
 
       
-
+      // Calculate action
       RealD S(LatticeLorentzColourMatrix& U){
 	LatticeComplex Hloc(U._grid);
 	Hloc = zero;
@@ -153,7 +146,74 @@ namespace Grid{
     class MinimumNorm2{
       const double lambda = 0.1931833275037836;
     public:
-      void step (LatticeLorentzColourMatrix& U, int level, std::vector<int>& clock);
+      void step (LatticeLorentzColourMatrix& U, 
+		 int level, std::vector<int>& clock,
+		 Integrator<MinimumNorm2>* Integ){
+	// level  : current level
+	// fl  : final level
+	// eps : current step size
+
+	int fl = Integ->as.size() -1;
+	double eps = Integ->Params.stepsize;
+	
+	for(int l=0; l<=level; ++l) eps/= 2.0*Integ->Nrel[l];
+	
+	int fin = Integ->Nrel[0];
+	for(int l=1; l<=level; ++l) fin*= 2.0*Integ->Nrel[l];
+	fin = 3*Integ->Params.MDsteps*fin -1;
+	
+	
+	for(int e=0; e<Integ->Nrel[level]; ++e){
+	  
+	  if(clock[level] == 0){    // initial half step 
+	    Integ->update_P(U,level,lambda*eps);
+	    ++clock[level];
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< clock[level] <<std::endl;
+	  }
+	  
+	  if(level == fl){          // lowest level 
+	    Integ->update_U(U,0.5*eps);
+	    
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"U "<< (clock[level]+1) <<std::endl;
+	  }else{                 // recursive function call 
+	    step(U,level+1,clock, Integ);
+	  }
+	  
+	  Integ->update_P(U,level,(1.0-2.0*lambda)*eps);
+	  ++clock[level];
+	  for(int l=0; l<level;++l) std::cout<<"   ";
+	  std::cout<<"P "<< (clock[level]) <<std::endl;
+	  
+	  if(level == fl){          // lowest level 
+	    Integ->update_U(U,0.5*eps);
+	    
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"U "<< (clock[level]+1) <<std::endl;
+	  }else{                 // recursive function call 
+	    step(U,level+1,clock, Integ);
+	  }    
+	  
+	  
+	  if(clock[level] == fin){  // final half step
+	    Integ->update_P(U,level,lambda*eps);
+	    
+	    ++clock[level];
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< clock[level] <<std::endl;
+	  }else{                  // bulk step
+	    Integ->update_P(U,level,lambda*2.0*eps);
+	    
+	    clock[level]+=2;
+	    for(int l=0; l<level;++l) std::cout<<"   ";
+	    std::cout<<"P "<< clock[level] <<std::endl;
+	  }
+	}
+	
+	
+	
+      }
       
     };
     
@@ -162,7 +222,6 @@ namespace Grid{
       void step (LatticeLorentzColourMatrix& U, 
 		 int level, std::vector<int>& clock,
 		 Integrator<LeapFrog>* Integ){
-	// cl  : current level
 	// fl  : final level
 	// eps : current step size
 	
@@ -179,7 +238,7 @@ namespace Grid{
 	for(int e=0; e<Integ->Nrel[level]; ++e){
 	  
 	  if(clock[level] == 0){    // initial half step
-	    Integ->update_P(U, level,eps/2);
+	    Integ->update_P(U, level,eps/2.0);
 	    ++clock[level];
 	    for(int l=0; l<level;++l) std::cout<<"   ";
 	    std::cout<<"P "<< 0.5*clock[level] <<std::endl;
@@ -192,7 +251,7 @@ namespace Grid{
 	    step(U, level+1,clock, Integ);
 	  }
 	  if(clock[level] == fin){  // final half step
-	    Integ->update_P(U, level,eps/2);
+	    Integ->update_P(U, level,eps/2.0);
 	    
 	    ++clock[level];
 	    for(int l=0; l<level;++l) std::cout<<"   ";
diff --git a/tests/Test_hmc_WilsonGauge.cc b/tests/Test_hmc_WilsonGauge.cc
index 4fc468c9..98f75b3f 100644
--- a/tests/Test_hmc_WilsonGauge.cc
+++ b/tests/Test_hmc_WilsonGauge.cc
@@ -29,7 +29,7 @@ int main (int argc, char ** argv)
   SU3::HotConfiguration(pRNG, U);
  
 
-  // simplify template?
+  // simplify template declaration? Strip the lorentz from the second template
   WilsonGaugeAction<LatticeLorentzColourMatrix, LatticeColourMatrix> Waction(6.0);
 
   //Collect actions
@@ -39,13 +39,14 @@ int main (int argc, char ** argv)
   FullSet.push_back(Level1);
 
   // Create integrator
+  typedef MinimumNorm2  IntegratorAlgorithm;// change here to change the algorithm
   IntegratorParameters MDpar(12,30,1.0);
   std::vector<int> rel ={1};
-  Integrator<LeapFrog> MDleapfrog(MDpar, FullSet,rel);
+  Integrator<IntegratorAlgorithm> MDynamics(&Fine,MDpar, FullSet,rel);
 
   // Create HMC
   HMCparameters HMCpar;
-  HybridMonteCarlo<LeapFrog>  HMC(HMCpar, MDleapfrog, &Fine);
+  HybridMonteCarlo<IntegratorAlgorithm>  HMC(HMCpar, MDynamics);
 
   HMC.evolve(U);