Added all elements for Hirep HMC

TODO: Test and debug

lib/qcd/action/ActionBase.h

@@ -97,7 +97,8 @@ struct ActionLevelHirep {
   //std::vector<ActPtr> actions;
   // construct a tuple of vectors of the actions for the corresponding higher
   // representation fields
-  typename AccessTypes<Action, Repr>::VectorCollection actions_hirep;
+  typedef typename AccessTypes<Action, Repr>::VectorCollection action_collection;
+  action_collection actions_hirep;
   typedef typename  AccessTypes<Action, Repr>::ClassCollection actions_hirep_ptrs_type;
 
   std::vector<ActPtr>& actions;
@@ -109,7 +110,7 @@ struct ActionLevelHirep {
 
   ActionLevelHirep(unsigned int mul = 1) : actions(std::get<0>(actions_hirep)), multiplier(mul) {
     // initialize the hirep vectors to zero.
-    //apply(&ActionLevelHirep::resize, actions_hirep, 0); //need a working resize
+    //apply(this->resize, actions_hirep, 0); //need a working resize
     assert(mul >= 1);
   };
 
@@ -128,18 +129,19 @@ struct ActionLevelHirep {
 
   }
 
-  
+  template <std::size_t I>
+  auto getRepresentation(Repr& R)->decltype(std::get<I>(R).U)  {return std::get<I>(R).U;}
 
   // Loop on tuple for a callable function
-  template <std::size_t I = 0, class Tuple, typename Callable, typename ...Args>
+  template <std::size_t I = 1, typename Callable, typename ...Args>
-  inline typename std::enable_if<(I == std::tuple_size<Tuple>::value), void>::type apply(
+  inline typename std::enable_if<I == std::tuple_size<action_collection>::value, void>::type apply(
-      Callable&, Tuple& , Args...) {}
+      Callable, Repr& R,Args...) const {}
 
-  template <std::size_t I = 0, class Tuple, typename Callable, typename ...Args>
+  template <std::size_t I = 1, typename Callable, typename ...Args>
-  inline typename std::enable_if<(I < std::tuple_size<Tuple>::value), void>::type apply(
+  inline typename std::enable_if<I < std::tuple_size<action_collection>::value, void>::type apply(
-      Callable& fn,  Tuple& T, Args... arguments) {
+      Callable fn, Repr& R, Args... arguments) const {
-    fn(std::get<I>(T), arguments...);
+    fn(std::get<I>(actions_hirep), std::get<I>(R.rep), arguments...);
-    apply<I + 1>(T, fn, arguments...);
+    apply<I + 1>(fn, R, arguments...);
   }  
 
 };

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -111,12 +111,15 @@ namespace Grid {
 ///////
 // Single flavour four spinors with colour index
 ///////
-template <class S, int Nrepresentation = Nc>
+template <class S, class Representation = FundamentalRepresentation >
 class WilsonImpl
-    : public PeriodicGaugeImpl<GaugeImplTypes<S, Nrepresentation> > {
+    : public PeriodicGaugeImpl<GaugeImplTypes<S, Representation::Dimension > > {
  public:
-  typedef PeriodicGaugeImpl<GaugeImplTypes<S, Nrepresentation> > Gimpl;
+  static const int Nrepresentation = Representation::Dimension;
-  constexpr bool is_fundamental() const{return Nrepresentation == Nc ? 1 : 0;}
+  typedef PeriodicGaugeImpl<GaugeImplTypes<S, Representation::Dimension > > Gimpl;
+  
+  //Necessary?
+  constexpr bool is_fundamental() const{return Representation::Dimension == Nc ? 1 : 0;}
 
   INHERIT_GIMPL_TYPES(Gimpl);
 
@@ -501,13 +504,13 @@ class GparityWilsonImpl
   }
 };
 
-typedef WilsonImpl<vComplex,  Nc> WilsonImplR;   // Real.. whichever prec
+typedef WilsonImpl<vComplex,  FundamentalRepresentation> WilsonImplR;   // Real.. whichever prec
-typedef WilsonImpl<vComplexF, Nc> WilsonImplF;  // Float
+typedef WilsonImpl<vComplexF, FundamentalRepresentation> WilsonImplF;  // Float
-typedef WilsonImpl<vComplexD, Nc> WilsonImplD;  // Double
+typedef WilsonImpl<vComplexD, FundamentalRepresentation> WilsonImplD;  // Double
 
-typedef WilsonImpl<vComplex,  SU_Adjoint<Nc>::Dimension > WilsonAdjImplR;   // Real.. whichever prec
+typedef WilsonImpl<vComplex,  AdjointRepresentation > WilsonAdjImplR;   // Real.. whichever prec
-typedef WilsonImpl<vComplexF, SU_Adjoint<Nc>::Dimension > WilsonAdjImplF;  // Float
+typedef WilsonImpl<vComplexF, AdjointRepresentation > WilsonAdjImplF;  // Float
-typedef WilsonImpl<vComplexD, SU_Adjoint<Nc>::Dimension > WilsonAdjImplD;  // Double
+typedef WilsonImpl<vComplexD, AdjointRepresentation > WilsonAdjImplD;  // Double
 
 typedef DomainWallRedBlack5dImpl<vComplex, Nc>
     DomainWallRedBlack5dImplR;  // Real.. whichever prec

lib/qcd/hmc/HmcRunner.h

@@ -113,7 +113,7 @@ class NerscHmcRunnerTemplate {
     //////////////
     typedef MinimumNorm2<GaugeField, SmearedConfiguration<Gimpl>, RepresentationsPolicy >
         IntegratorType;  // change here to change the algorithm
-    IntegratorParameters MDpar(20);
+    IntegratorParameters MDpar(20, 1.0);
     IntegratorType MDynamics(UGrid, MDpar, TheAction, SmearingPolicy);
 
     // Checkpoint strategy

lib/qcd/hmc/integrators/Integrator.h

@@ -64,7 +64,7 @@ struct IntegratorParameters {
 };
 
 /*! @brief Class for Molecular Dynamics management */
-template <class GaugeField, class SmearingPolicy ,class RepresentationPolicy >
+template <class GaugeField, class SmearingPolicy, class RepresentationPolicy>
 class Integrator {
  protected:
   typedef IntegratorParameters ParameterType;
@@ -81,7 +81,7 @@ class Integrator {
 
   SmearingPolicy& Smearer;
 
-  RepresentationPolicy Representations; 
+  RepresentationPolicy Representations;
 
   // Should match any legal (SU(n)) gauge field
   // Need to use this template to match Ncol to pass to SU<N> class
@@ -112,26 +112,26 @@ class Integrator {
 
   // to be used by the actionlevel class to iterate
   // over the representations
-  template <class Level>
+  struct _updateP {
-  void update_P_hireps(Level repr_level, GaugeField& Mom, GaugeField& U,
+    template <class FieldType, class GF, class Repr>
-                     double ep) {
+    void operator()(std::vector<Action<FieldType>*> repr_set, Repr& Rep,
-    typedef typename Level::LatticeField FieldType;
+                    GF& Mom, GF& U, double ep) {
-    FieldType Ur = repr_level->getRepresentation();// update U is better
+      for (int a = 0; a < repr_set.size(); ++a) {
-    for (int a = 0; a < repr_level.size(); ++a) {
+        FieldType forceR(U._grid);
-      FieldType forceR(U._grid);
+        // Implement smearing only for the fundamental representation now
-      // Implement smearing only for the fundamental representation now
+        repr_set.at(a)->deriv(Rep.U, forceR);
-      repr_level.at(a)->deriv(Ur, forceR);
+        GF force =
-      GaugeField force = repr_level.at(a)->RtoFundamentalProject(forceR);
+            Rep.RtoFundamentalProject(forceR);  // Ta for the fundamental rep
-      std::cout << GridLogIntegrator
+        std::cout << GridLogIntegrator << "Hirep Force average: "
-                << "Hirep Force average: " << norm2(force) / (U._grid->gSites())
+                  << norm2(force) / (U._grid->gSites()) << std::endl;
-                << std::endl;
+        Mom -= force * ep;
-      Mom -= force * ep;
+      }
     }
-  }
+  } update_P_hireps{};
 
   void update_P(GaugeField& Mom, GaugeField& U, int level, double ep) {
     // input U actually not used in the fundamental case
-  	// Fundamental updates, include smearing
+    // Fundamental updates, include smearing
     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);
@@ -147,8 +147,9 @@ class Integrator {
                 << std::endl;
       Mom -= force * ep;
     }
-    // Add here the other representations
+
-    //apply(update_P_hireps, as[level], Args...)
+    // Force from the other representations
+    as[level].apply(update_P_hireps, Representations, Mom, U, ep);
   }
 
   void update_U(GaugeField& U, double ep) {
@@ -172,15 +173,21 @@ 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
+    Representations.update(U);  // void functions if fundamental representation
   }
 
   virtual void step(GaugeField& U, int level, int first, int last) = 0;
 
  public:
   Integrator(GridBase* grid, IntegratorParameters Par,
-             ActionSetHirep<GaugeField, RepresentationPolicy>& Aset, SmearingPolicy& Sm)
+             ActionSetHirep<GaugeField, RepresentationPolicy>& Aset,
-      : Params(Par), as(Aset), P(grid), levels(Aset.size()), Smearer(Sm), Representations(grid) {
+             SmearingPolicy& 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
@@ -188,13 +195,24 @@ class Integrator {
 
   virtual ~Integrator() {}
 
+  // to be used by the actionlevel class to iterate
+  // over the representations
+  struct _refresh {
+    template <class FieldType, class Repr>
+    void operator()(std::vector<Action<FieldType>*> repr_set, Repr& Rep,
+                    GridParallelRNG& pRNG) {
+      for (int a = 0; a < repr_set.size(); ++a)
+        repr_set.at(a)->refresh(Rep.U, pRNG);
+    }
+  } refresh_hireps{};
+
   // Initialization of momenta and actions
   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 
+    // necessary to keep the fields updated even after a reject
     // of the Metropolis
     Smearer.set_GaugeField(U);
     // Set the (eventual) representations gauge fields
@@ -211,12 +229,27 @@ class Integrator {
             Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
         as[level].actions.at(actionID)->refresh(Us, pRNG);
       }
+
+      as[level].apply(refresh_hireps, Representations, pRNG);
     }
- 
-
-
   }
 
+  // to be used by the actionlevel class to iterate
+  // over the representations
+  struct _S {
+    template <class FieldType, class Repr>
+    void operator()(std::vector<Action<FieldType>*> repr_set, Repr& Rep,
+                    int level, RealD& H) {
+      RealD H_hirep = 0.0;
+      for (int a = 0; a < repr_set.size(); ++a) {
+        RealD Hterm = repr_set.at(a)->S(Rep.U);
+        std::cout << GridLogMessage << "S Level " << level << " term " << a
+                  << " H Hirep = " << Hterm << std::endl;
+        H += Hterm;
+      }
+    }
+  } S_hireps{};
+
   // Calculate action
   RealD S(GaugeField& U) {  // here also U not used
 
@@ -245,6 +278,7 @@ class Integrator {
                   << actionID << " H = " << Hterm << std::endl;
         H += Hterm;
       }
+      as[level].apply(S_hireps, Representations, level, H);
     }
 
     return H;
@@ -257,8 +291,7 @@ 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/representations/adjoint.h

@@ -1,12 +1,11 @@
 /*
  *	Policy classes for the HMC
  *	Author: Guido Cossu
-*/	
+*/
 
 #ifndef ADJOINT_H
 #define ADJOINT_H
 
-
 namespace Grid {
 namespace QCD {
 
@@ -19,17 +18,16 @@ namespace QCD {
 template <int ncolour>
 class AdjointRep {
  public:
- 	// typdef to be used by the Representations class in HMC to get the
+  // typdef to be used by the Representations class in HMC to get the
- 	// types for the higher representation fields
+  // types for the higher representation fields
   typedef typename SU_Adjoint<ncolour>::LatticeAdjMatrix LatticeMatrix;
   typedef typename SU_Adjoint<ncolour>::LatticeAdjField LatticeField;
-  const int Dimension = ncolour * ncolour - 1;
+  static const int Dimension = ncolour * ncolour - 1;
 
   LatticeField U;
-  
 
-  explicit AdjointRep(GridBase* grid) : U(grid) {}
+  explicit AdjointRep(GridBase *grid) : U(grid) {}
-  void update_representation(const LatticeGaugeField& Uin) {
+  void 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]
@@ -38,31 +36,59 @@ class AdjointRep {
     // T_F is 1/2 for the fundamental representation
     conformable(U, Uin);
     U = zero;
-    LatticeGaugeField tmp(Uin._grid);
+    LatticeColourMatrix tmp(Uin._grid);
 
-    Vector<typename SU<ncolour>::Matrix> ta(ncolour * ncolour - 1);
+    Vector<typename SU<ncolour>::Matrix> ta(Dimension);
 
     // 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++) {
+    for (int mu = 0; mu < Nd; mu++) {
-      tmp = 2.0 * adj(Uin) * ta[a] * Uin;
+      auto Uin_mu = peekLorentz(Uin, mu);
-      for (int b = 0; b < (ncolour * ncolour - 1); b++) {
+      auto U_mu = peekLorentz(U, mu);
-        auto Tr = TensorRemove(trace(tmp * ta[b]));
+      for (int a = 0; a < Dimension; a++) {
-        pokeColour(U, Tr, a, b);
+        tmp = 2.0 * adj(Uin_mu) * ta[a] * Uin_mu;
+        for (int b = 0; b < (ncolour * ncolour - 1); b++)
+          pokeColour(U_mu, trace(tmp * ta[b]), a, b);
       }
-    }  	 
+      pokeLorentz(U, U_mu, mu);
-    
+    }
+  }
+
+  LatticeGaugeField RtoFundamentalProject(const LatticeField &in,
+                                          Real scale = 1.0) const {
+    LatticeGaugeField out(in._grid);
+
+    for (int mu = 0; mu < Nd; mu++) {
+      LatticeColourMatrix out_mu(in._grid);  // fundamental representation
+      LatticeMatrix in_mu = peekLorentz(in, mu);
+
+      out_mu = zero;
+
+      typename SU<ncolour>::LatticeAlgebraVector h(in._grid);
+      projectOnAlgebra(h, in_mu, scale);
+      FundamentalLieAlgebraMatrix(h, out_mu, 1.0);  // apply scale only once
+      pokeLorentz(out, out_mu, mu);
+    }
+    return out;
+  }
+
+ private:
+  void projectOnAlgebra(typename SU<ncolour>::LatticeAlgebraVector &h_out,
+                        const LatticeMatrix &in, Real scale = 1.0) const {
+    SU_Adjoint<ncolour>::projectOnAlgebra(h_out, in, scale);
+  }
+
+  void FundamentalLieAlgebraMatrix(
+      typename SU<ncolour>::LatticeAlgebraVector &h,
+      typename SU<ncolour>::LatticeMatrix &out, Real scale = 1.0) const {
+    SU<ncolour>::FundamentalLieAlgebraMatrix(h, out, scale);
   }
 };
 
-typedef	 AdjointRep<Nc> AdjointRepresentation;
+typedef AdjointRep<Nc> AdjointRepresentation;
-
 }
 }
 
-
-
-
 #endif

lib/qcd/representations/fundamental.h

@@ -18,7 +18,7 @@ namespace QCD {
 template <int ncolour>
 class FundamentalRep {
  public:
-  const int Dimension = ncolour;
+  static const int Dimension = ncolour;
 
   // typdef to be used by the Representations class in HMC to get the
   // types for the higher representation fields
@@ -27,6 +27,11 @@ class FundamentalRep {
   
   explicit FundamentalRep(GridBase* grid) {} //do nothing
   void update_representation(const LatticeGaugeField& Uin) {} // do nothing
+
+  LatticeField RtoFundamentalProject(const LatticeField& in, Real scale = 1.0) const{
+    return (scale * in);
+  }
+
 };
 
 typedef	 FundamentalRep<Nc> FundamentalRepresentation;

lib/qcd/representations/hmc_types.h

@@ -28,7 +28,7 @@ class Representations {
   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
+  // type repr_type<I>::LatticeField
 
   Representation_type rep;
 

lib/qcd/utils/SUnAdjoint.h

@@ -111,7 +111,7 @@ class SU_Adjoint : public SU<ncolour> {
   }
 
   // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 )
-  static void projectOnAlgebra(typename SU<ncolour>::LatticeAlgebraVector &h_out, LatticeAdjMatrix &in, Real scale = 1.0) {
+  static void projectOnAlgebra(typename SU<ncolour>::LatticeAlgebraVector &h_out, const LatticeAdjMatrix &in, Real scale = 1.0) {
     conformable(h_out, in);
     h_out = zero;
     AMatrix iTa;
@@ -124,7 +124,7 @@ class SU_Adjoint : public SU<ncolour> {
   }
 
   // a projector that keeps the generators stored to avoid the overhead of recomputing. 
-  static void projector(typename SU<ncolour>::LatticeAlgebraVector &h_out, LatticeAdjMatrix &in, Real scale = 1.0) {
+  static void projector(typename SU<ncolour>::LatticeAlgebraVector &h_out, const LatticeAdjMatrix &in, Real scale = 1.0) {
     conformable(h_out, in);
     static std::vector<AMatrix> iTa(Dimension);  // to store the generators
     h_out = zero;

tests/Test_hmc_WilsonAdjointFermionGauge.cc

@@ -40,7 +40,7 @@ namespace Grid {
 namespace QCD {
 
 // Here change the allowed (higher) representations
-typedef Representations< FundamentalRepresentation, FundamentalRepresentation > TheRepresentations;
+typedef Representations< FundamentalRepresentation, AdjointRepresentation > TheRepresentations;
 
 
 class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > {