Debugged smearing and adding HMC functions for hirep

lib/qcd/QCD.h

@@ -496,11 +496,12 @@ namespace QCD {
 
 #include <qcd/action/Actions.h>
 
+#include <qcd/smearing/Smearing.h>
+
+#include <qcd/representations/hmc_types.h>
+
 #include <qcd/hmc/integrators/Integrator.h>
 #include <qcd/hmc/integrators/Integrator_algorithm.h>
 #include <qcd/hmc/HMC.h>
 
-#include <qcd/smearing/Smearing.h>
-
-#include <qcd/representations/hmc_types.h>
 #endif

lib/qcd/action/ActionBase.h

@@ -66,13 +66,13 @@ template<class GaugeField> struct ActionLevel{
 public:
    
   typedef Action<GaugeField>*  ActPtr; // now force the same colours as the rest of the code
-
-  int multiplier;
+    
+  unsigned int multiplier;
 
   std::vector<ActPtr> actions;
 
-  ActionLevel(int mul = 1) : multiplier(mul) {
-    assert (mul > 0);
+  ActionLevel(unsigned int mul = 1) : actions(0), multiplier(mul) {
+    assert (mul >= 1);
   };
    
   void push_back(ActPtr ptr){
@@ -83,5 +83,6 @@ public:
 template<class GaugeField> using ActionSet = std::vector<ActionLevel< GaugeField > >;
 
 
-}}
+}
+}
 #endif

lib/qcd/hmc/integrators/Integrator.h

@@ -64,7 +64,7 @@ struct IntegratorParameters {
 };
 
 /*! @brief Class for Molecular Dynamics management */
-template <class GaugeField, class SmearingPolicy>
+template <class GaugeField, class SmearingPolicy ,class RepresentationPolicy >
 class Integrator {
  protected:
   typedef IntegratorParameters ParameterType;
@@ -81,6 +81,8 @@ class Integrator {
 
   SmearingPolicy& Smearer;
 
+  RepresentationPolicy Representations; 
+
   // Should match any legal (SU(n)) gauge field
   // Need to use this template to match Ncol to pass to SU<N> class
   template <int Ncol, class vec>
@@ -108,8 +110,28 @@ class Integrator {
               << " dt " << ep << " : t_P " << t_P[level] << std::endl;
   }
 
+  // to be used by the actionlevel class to iterate
+  // over the representations
+  template <class Level>
+  void update_P_core(Level repr_level, GaugeField& Mom, GaugeField& U,
+                     double ep) {
+    typedef typename Level::LatticeField FieldType;
+    FieldType Ur = repr_level->getRepresentation();// update U is better
+    for (int a = 0; a < repr_level.size(); ++a) {
+      FieldType forceR(U._grid);
+      // Implement smearing only for the fundamental representation now
+      repr_level.at(a)->deriv(Ur, forceR);
+      GaugeField force = repr_level.at(a)->RtoFundamentalProject(forceR);
+      std::cout << GridLogIntegrator
+                << "Hirep Force average: " << norm2(force) / (U._grid->gSites())
+                << std::endl;
+      Mom -= force * ep;
+    }
+  }
+  // Add the specialized class for the fundamental case
+
   void update_P(GaugeField& Mom, GaugeField& U, int level, double ep) {
-    // input U actually not used...
+    // input U actually not used in the fundamental case
     for (int a = 0; a < as[level].actions.size(); ++a) {
       GaugeField force(U._grid);
       GaugeField& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
@@ -125,6 +147,8 @@ class Integrator {
                 << std::endl;
       Mom -= force * ep;
     }
+    // Add here the other representations
+    // as[level].apply(update_P_hireps, Args...)
   }
 
   void update_U(GaugeField& U, double ep) {
@@ -147,6 +171,8 @@ class Integrator {
     }
     // Update the smeared fields, can be implemented as observer
     Smearer.set_GaugeField(U);
+    // Update the higher representations fields
+    //Representations.update(U);// void functions if fundamental representation
   }
 
   virtual void step(GaugeField& U, int level, int first, int last) = 0;
@@ -154,7 +180,7 @@ class Integrator {
  public:
   Integrator(GridBase* grid, IntegratorParameters Par,
              ActionSet<GaugeField>& Aset, SmearingPolicy& Sm)
-      : Params(Par), as(Aset), P(grid), levels(Aset.size()), Smearer(Sm) {
+      : Params(Par), as(Aset), P(grid), levels(Aset.size()), Smearer(Sm), Representations(grid) {
     t_P.resize(levels, 0.0);
     t_U = 0.0;
     // initialization of smearer delegated outside of Integrator
@@ -166,8 +192,16 @@ class Integrator {
   void refresh(GaugeField& U, GridParallelRNG& pRNG) {
     std::cout << GridLogIntegrator << "Integrator refresh\n";
     generate_momenta(P, pRNG);
+    
+    // Update the smeared fields, can be implemented as observer
+    // necessary to keep the fields updated even after a reject 
+    // of the Metropolis
+    Smearer.set_GaugeField(U);
+    // Set the (eventual) representations gauge fields
+    // Representations.update(U);
+
     // The Smearer is attached to a pointer of the gauge field
-    // automatically gets the updated field
+    // automatically gets the correct field
     // whether or not has been accepted in the previous sweep
     for (int level = 0; level < as.size(); ++level) {
       for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
@@ -178,6 +212,9 @@ class Integrator {
         as[level].actions.at(actionID)->refresh(Us, pRNG);
       }
     }
+ 
+
+
   }
 
   // Calculate action
@@ -220,7 +257,8 @@ class Integrator {
       t_P[level] = 0;
     }
 
-    for (int step = 0; step < Params.MDsteps; ++step) {  // MD step
+
+     for (int step = 0; step < Params.MDsteps; ++step) {  // MD step
       int first_step = (step == 0);
       int last_step = (step == Params.MDsteps - 1);
       this->step(U, 0, first_step, last_step);

lib/qcd/hmc/integrators/Integrator_algorithm.h

@@ -91,17 +91,19 @@ namespace Grid{
     *  P 1/2                            P 1/2
     */    
 
-    template<class GaugeField, class SmearingPolicy> class LeapFrog :
-      public Integrator<GaugeField, SmearingPolicy> {
+    template<class GaugeField,
+	     class SmearingPolicy,
+	     class RepresentationPolicy = Representations< FundamentalRepresentation > > class LeapFrog :
+      public Integrator<GaugeField, SmearingPolicy, RepresentationPolicy> {
     public:
 
-      typedef LeapFrog<GaugeField, SmearingPolicy> Algorithm;
+      typedef LeapFrog<GaugeField, SmearingPolicy, RepresentationPolicy> Algorithm;
 
       LeapFrog(GridBase* grid, 
 	       IntegratorParameters Par,
 	       ActionSet<GaugeField> & Aset,
 	       SmearingPolicy & Sm):
-	Integrator<GaugeField, SmearingPolicy>(grid,Par,Aset,Sm) {};
+	Integrator<GaugeField, SmearingPolicy, RepresentationPolicy>(grid,Par,Aset,Sm) {};
 
 
       void step (GaugeField& U, int level,int _first, int _last){
@@ -138,8 +140,10 @@ namespace Grid{
       }
     };
 
-    template<class GaugeField, class SmearingPolicy> class MinimumNorm2 :
-      public Integrator<GaugeField, SmearingPolicy> {
+    template<class GaugeField,
+	     class SmearingPolicy,
+	     class RepresentationPolicy = Representations < FundamentalRepresentation > > class MinimumNorm2 :
+      public Integrator<GaugeField, SmearingPolicy, RepresentationPolicy> {
     private:
       const RealD lambda = 0.1931833275037836;
 
@@ -149,7 +153,7 @@ namespace Grid{
 		   IntegratorParameters Par,
 		   ActionSet<GaugeField> & Aset,
 		   SmearingPolicy& Sm):
-	Integrator<GaugeField, SmearingPolicy>(grid,Par,Aset,Sm) {};
+	Integrator<GaugeField, SmearingPolicy, RepresentationPolicy>(grid,Par,Aset,Sm) {};
 
       void step (GaugeField& U, int level, int _first,int _last){
 
@@ -197,8 +201,10 @@ namespace Grid{
     };
 
 
-    template<class GaugeField, class SmearingPolicy> class ForceGradient :
-      public Integrator<GaugeField, SmearingPolicy> {
+    template<class GaugeField,
+	     class SmearingPolicy,
+	     class RepresentationPolicy = Representations< FundamentalRepresentation > > class ForceGradient :
+      public Integrator<GaugeField, SmearingPolicy, RepresentationPolicy> {
     private:
       const RealD lambda = 1.0/6.0;;
       const RealD chi    = 1.0/72.0;
@@ -211,7 +217,7 @@ namespace Grid{
 		  IntegratorParameters Par,
 		  ActionSet<GaugeField> & Aset,
 		  SmearingPolicy &Sm):
-      Integrator<GaugeField, SmearingPolicy>(grid,Par,Aset, Sm) {};
+      Integrator<GaugeField, SmearingPolicy, RepresentationPolicy>(grid,Par,Aset, Sm) {};
 
 
       void FG_update_P(GaugeField&U, int level,double fg_dt,double ep){

lib/qcd/representations/adjoint.h

@@ -19,11 +19,16 @@ namespace QCD {
 template <int ncolour>
 class AdjointRep {
  public:
-  typename SU_Adjoint<ncolour>::LatticeAdjMatrix U;
+ 	// typdef to be used by the Representations class in HMC to get the
+ 	// types for the higher representation fields
+  typedef typename SU_Adjoint<ncolour>::LatticeAdjMatrix LatticeField;
   const int Dimension = ncolour * ncolour - 1;
 
-  explicit AdjointRep(GridBase* grid):U(grid) {}
-  void update_representation(const LatticeGaugeField& Uin) {
+  LatticeField U;
+  
+
+  explicit AdjointRep(GridBase* grid) : U(grid) {}
+  LatticeField update_representation(const LatticeGaugeField& Uin) {
     // Uin is in the fundamental representation
     // get the U in AdjointRep
     // (U_adj)_B = tr[e^a U e^b U^dag]
@@ -34,19 +39,20 @@ class AdjointRep {
     U = zero;
     LatticeGaugeField tmp(Uin._grid);
 
-    Vector<typename SU<ncolour>::Matrix > ta(ncolour * ncolour - 1);
+    Vector<typename SU<ncolour>::Matrix> ta(ncolour * ncolour - 1);
 
-    // FIXME probably not very efficient to get all the generators everytime
+    // FIXME probably not very efficient to get all the generators
+    // everytime
     for (int a = 0; a < Dimension; a++) SU<ncolour>::generator(a, ta[a]);
 
     for (int a = 0; a < Dimension; a++) {
-    	tmp = 2.0 * adj(Uin) * ta[a] * Uin;
+      tmp = 2.0 * adj(Uin) * ta[a] * Uin;
       for (int b = 0; b < (ncolour * ncolour - 1); b++) {
         auto Tr = TensorRemove(trace(tmp * ta[b]));
-        pokeColour(U, Tr, a,b);
+        pokeColour(U, Tr, a, b);
       }
     }  	 
-
+    
   }
 };
 

lib/qcd/representations/fundamental.h

@@ -0,0 +1,39 @@
+/*
+ *	Policy classes for the HMC
+ *	Author: Guido Cossu
+*/	
+
+#ifndef FUNDAMENTAL_H
+#define FUNDAMENTAL_H
+
+
+namespace Grid {
+namespace QCD {
+
+/*
+* This is an helper class for the HMC
+* Empty since HMC updates already the fundamental representation 
+*/
+
+template <int ncolour>
+class FundamentalRep {
+ public:
+  const int Dimension = ncolour;
+
+  // typdef to be used by the Representations class in HMC to get the
+  // types for the higher representation fields
+  typedef typename SU<ncolour>::LatticeMatrix LatticeField;
+  
+  explicit FundamentalRep(GridBase* grid) {} //do nothing
+  void update_representation(const LatticeGaugeField& Uin) {} // do nothing
+};
+
+typedef	 FundamentalRep<Nc> FundamentalRepresentation;
+
+}
+}
+
+
+
+
+#endif

lib/qcd/representations/hmc_types.h

@@ -4,16 +4,27 @@
 #include <tuple>
 #include <utility>
 #include <qcd/representations/adjoint.h>
+#include <qcd/representations/fundamental.h>
 
 namespace Grid {
 namespace QCD {
 
-// Utility to add support for representations other than the fundamental
 
+// Supported types
+//enum {Fundamental, Adjoint} repr_type;
+
+// Utility to add support to the HMC for representations other than the fundamental
 template<class... Reptypes>
 class Representations{
 public:
   typedef std::tuple<Reptypes...> Representation_type;
+
+  // To access the Reptypes (FundamentalRepresentation, AdjointRepresentation)
+  template <std::size_t N>
+  using repr_type = typename std::tuple_element<N, Representation_type >::type;
+  // in order to get the typename of the field use
+  // type repr_type::LatticeField
+
   Representation_type rep;
 
   // Multiple types constructor
@@ -29,12 +40,16 @@ public:
 
   template <std::size_t I = 0>
       inline typename std::enable_if <
-      I<sizeof...(Reptypes), void >::type update(LatticeGaugeField& U) {
+      I<sizeof...(Reptypes), void>::type update(LatticeGaugeField& U) {
     std::get<I>(rep).update_representation(U);
     update<I + 1>(U);
   }  
 };
 
+
+typedef Representations<FundamentalRepresentation> JustTheFundamental;
+
+
 }
 }
 

lib/qcd/utils/SUn.h

@@ -728,6 +728,10 @@ typedef SU<2> SU2;
 typedef SU<3> SU3;
 typedef SU<4> SU4;
 typedef SU<5> SU5;
+
+
+typedef SU<Nc> FundamentalMatrices;
+
 }
 }
 #endif

lib/qcd/utils/SUnAdjoint.h

@@ -137,12 +137,14 @@ class SU_Adjoint : public SU<ncolour> {
 
 
 
-
+// Some useful type names
 
 typedef SU_Adjoint<2> SU2Adjoint;
 typedef SU_Adjoint<3> SU3Adjoint;
 typedef SU_Adjoint<4> SU4Adjoint;
 typedef SU_Adjoint<5> SU5Adjoint;
+
+typedef SU_Adjoint<Nc> AdjointMatrices;
 }
 }