Added HMC utitities for the higher representations

TODO: Inherit types for the pseudofermions, Debugging, testing

lib/qcd/QCD.h

@@ -484,6 +484,8 @@ namespace QCD {
 }   //namespace QCD
 } // Grid
 
+
+
 #include <qcd/spin/Dirac.h>
 
 #include <qcd/spin/TwoSpinor.h>
@@ -494,12 +496,12 @@ namespace QCD {
 #include <qcd/utils/SUn.h>
 #include <qcd/utils/SUnAdjoint.h>
 
+#include <qcd/representations/hmc_types.h>
+
 #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>

lib/qcd/action/ActionBase.h

@@ -1,49 +1,52 @@
-    /*************************************************************************************
+/*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid 
+Grid physics library, www.github.com/paboyle/Grid
 
-    Source file: ./lib/qcd/action/ActionBase.h
+Source file: ./lib/qcd/action/ActionBase.h
 
-    Copyright (C) 2015
+Copyright (C) 2015
 
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: neo <cossu@post.kek.jp>
 
-    This program is free software; you can redistribute it and/or modify
+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
+it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
+the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
+(at your option) any later version.
 
-    This program is distributed in the hope that it will be useful,
+This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
+but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
+GNU General Public License for more details.
 
-    You should have received a copy of the GNU General Public License along
+You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
+with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
-    See the full license in the file "LICENSE" in the top level distribution directory
+See the full license in the file "LICENSE" in the top level distribution
-    *************************************************************************************/
+directory
-    /*  END LEGAL */
+*************************************************************************************/
+/*  END LEGAL */
 #ifndef QCD_ACTION_BASE
 #define QCD_ACTION_BASE
 namespace Grid {
-namespace QCD{
+namespace QCD {
-
-template<class GaugeField>
-class Action { 
 
+template <class GaugeField>
+class Action {
  public:
   bool is_smeared = false;
   // Boundary conditions? // Heatbath?
-  virtual void  refresh(const GaugeField &U, GridParallelRNG& pRNG) = 0;// refresh pseudofermions
+  virtual void refresh(const GaugeField& U,
-  virtual RealD S    (const GaugeField &U)                        = 0;  // evaluate the action
+                       GridParallelRNG& pRNG) = 0;  // refresh pseudofermions
-  virtual void  deriv(const GaugeField &U,GaugeField & dSdU )     = 0;  // evaluate the action derivative
+  virtual RealD S(const GaugeField& U) = 0;         // evaluate the action
-  virtual ~Action() {};
+  virtual void deriv(const GaugeField& U,
+                     GaugeField& dSdU) = 0;  // evaluate the action derivative
+  virtual ~Action(){};
 };
 
-// Could derive PseudoFermion action with a PF field, FermionField, and a Grid; implement refresh
+// Could derive PseudoFermion action with a PF field, FermionField, and a Grid;
+// implement refresh
 /*
 template<class GaugeField, class FermionField>
 class PseudoFermionAction : public Action<GaugeField> {
@@ -52,7 +55,8 @@ class PseudoFermionAction : public Action<GaugeField> {
   GridParallelRNG &pRNG;
   GridBase &Grid;
 
-  PseudoFermionAction(GridBase &_Grid,GridParallelRNG &_pRNG) : Grid(_Grid), Phi(&_Grid), pRNG(_pRNG) {
+  PseudoFermionAction(GridBase &_Grid,GridParallelRNG &_pRNG) : Grid(_Grid),
+Phi(&_Grid), pRNG(_pRNG) {
   };
 
   virtual void refresh(const GaugeField &gauge) {
@@ -62,26 +66,90 @@ class PseudoFermionAction : public Action<GaugeField> {
 };
 */
 
-template<class GaugeField> struct ActionLevel{
+template <class GaugeField>
-public:
+struct ActionLevel {
-   
+ public:
-  typedef Action<GaugeField>*  ActPtr; // now force the same colours as the rest of the code
+  typedef Action<GaugeField>*
-    
+      ActPtr;  // now force the same colours as the rest of the code
+
+  //Add supported representations here
+
+
   unsigned int multiplier;
 
   std::vector<ActPtr> actions;
 
   ActionLevel(unsigned int mul = 1) : actions(0), multiplier(mul) {
-    assert (mul >= 1);
+    assert(mul >= 1);
   };
-   
+
-  void push_back(ActPtr ptr){
+  void push_back(ActPtr ptr) { actions.push_back(ptr); }
-    actions.push_back(ptr);
-  }
 };
 
-template<class GaugeField> using ActionSet = std::vector<ActionLevel< GaugeField > >;
 
+template <class GaugeField, class Repr>
+struct ActionLevelHirep {
+ public:
+  unsigned int multiplier; 
+
+  // Fundamental repr actions separated because of the smearing
+  typedef Action<GaugeField>* ActPtr;
+  //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>::ClassCollection actions_hirep_ptrs_type;
+
+  std::vector<ActPtr>& actions;
+
+  // Temporary conversion between ActionLevel and ActionLevelHirep
+  ActionLevelHirep(ActionLevel<GaugeField>& AL ):actions(AL.actions), multiplier(AL.multiplier){}
+
+
+
+  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
+    assert(mul >= 1);
+  };
+
+  void push_back(ActPtr ptr) { actions.push_back(ptr); }
+
+// SFINAE construct, check
+  template <class actionpointer, size_t N>
+  void push_back(actionpointer ptr, decltype(std::tuple_element<N, actions_hirep_ptrs_type>::value)* = 0) {
+    //insert only in the correct vector
+    std::get<N>(actions_hirep).push_back(ptr);
+  };
+
+  template < class ActPtr>
+  static void resize(ActPtr ap, unsigned int n){
+    ap->resize(n);
+
+  }
+
+  
+
+  // Loop on tuple for a callable function
+  template <std::size_t I = 0, class Tuple, typename Callable, typename ...Args>
+  inline typename std::enable_if<(I == std::tuple_size<Tuple>::value), void>::type apply(
+      Callable&, Tuple& , Args...) {}
+
+  template <std::size_t I = 0, class Tuple, typename Callable, typename ...Args>
+  inline typename std::enable_if<(I < std::tuple_size<Tuple>::value), void>::type apply(
+      Callable& fn,  Tuple& T, Args... arguments) {
+    fn(std::get<I>(T), arguments...);
+    apply<I + 1>(T, fn, arguments...);
+  }  
+
+};
+
+
+template <class GaugeField>
+using ActionSet = std::vector<ActionLevel<GaugeField> >;
+
+template <class GaugeField, class R>
+using ActionSetHirep = std::vector<ActionLevelHirep<GaugeField, R> >;
 
 }
 }

lib/qcd/hmc/HmcRunner.h

@@ -32,14 +32,14 @@ directory
 namespace Grid {
 namespace QCD {
 
-template <class Gimpl>
+template <class Gimpl, class RepresentationsPolicy = NoHirep >
 class NerscHmcRunnerTemplate {
  public:
   INHERIT_GIMPL_TYPES(Gimpl);
 
   enum StartType_t { ColdStart, HotStart, TepidStart, CheckpointStart };
 
-  ActionSet<GaugeField> TheAction;
+  ActionSetHirep<GaugeField, RepresentationsPolicy> TheAction;
 
   GridCartesian *UGrid;
   GridCartesian *FGrid;
@@ -111,7 +111,7 @@ class NerscHmcRunnerTemplate {
     SmearedConfiguration<Gimpl> SmearingPolicy(UGrid, Nsmear, Stout);
     std::cout << GridLogDebug << " done\n";
     //////////////
-    typedef MinimumNorm2<GaugeField, SmearedConfiguration<Gimpl> >
+    typedef MinimumNorm2<GaugeField, SmearedConfiguration<Gimpl>, RepresentationsPolicy >
         IntegratorType;  // change here to change the algorithm
     IntegratorParameters MDpar(20);
     IntegratorType MDynamics(UGrid, MDpar, TheAction, SmearingPolicy);
@@ -177,6 +177,12 @@ typedef NerscHmcRunnerTemplate<PeriodicGimplD> PeriodicNerscHmcRunnerD;
 typedef NerscHmcRunnerTemplate<ConjugateGimplR> ConjugateNerscHmcRunner;
 typedef NerscHmcRunnerTemplate<ConjugateGimplF> ConjugateNerscHmcRunnerF;
 typedef NerscHmcRunnerTemplate<ConjugateGimplD> ConjugateNerscHmcRunnerD;
+
+template <class RepresentationsPolicy>
+using NerscHmcRunnerHirep = NerscHmcRunnerTemplate<PeriodicGimplR, RepresentationsPolicy>;
+
+
+
 }
 }
 #endif

lib/qcd/hmc/integrators/Integrator.h

@@ -71,7 +71,7 @@ class Integrator {
 
   IntegratorParameters Params;
 
-  const ActionSet<GaugeField> as;
+  const ActionSetHirep<GaugeField, RepresentationPolicy> as;
 
   int levels;  //
   double t_U;  // Track time passing on each level and for U and for P
@@ -113,7 +113,7 @@ class Integrator {
   // 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,
+  void update_P_hireps(Level repr_level, GaugeField& Mom, GaugeField& U,
                      double ep) {
     typedef typename Level::LatticeField FieldType;
     FieldType Ur = repr_level->getRepresentation();// update U is better
@@ -128,10 +128,10 @@ class Integrator {
       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 in the fundamental case
+  	// 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);
@@ -148,7 +148,7 @@ class Integrator {
       Mom -= force * ep;
     }
     // Add here the other representations
-    // as[level].apply(update_P_hireps, Args...)
+    //apply(update_P_hireps, as[level], Args...)
   }
 
   void update_U(GaugeField& U, double ep) {
@@ -179,7 +179,7 @@ class Integrator {
 
  public:
   Integrator(GridBase* grid, IntegratorParameters Par,
-             ActionSet<GaugeField>& Aset, SmearingPolicy& Sm)
+             ActionSetHirep<GaugeField, RepresentationPolicy>& Aset, 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;

lib/qcd/hmc/integrators/Integrator_algorithm.h

@@ -151,7 +151,7 @@ namespace Grid{
 
       MinimumNorm2(GridBase* grid, 
 		   IntegratorParameters Par,
-		   ActionSet<GaugeField> & Aset,
+		   ActionSetHirep<GaugeField, RepresentationPolicy> & Aset,
 		   SmearingPolicy& Sm):
 	Integrator<GaugeField, SmearingPolicy, RepresentationPolicy>(grid,Par,Aset,Sm) {};
 

lib/qcd/representations/adjoint.h

@@ -21,14 +21,15 @@ class AdjointRep {
  public:
  	// 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;
+  typedef typename SU_Adjoint<ncolour>::LatticeAdjMatrix LatticeMatrix;
+  typedef typename SU_Adjoint<ncolour>::LatticeAdjField LatticeField;
   const int Dimension = ncolour * ncolour - 1;
 
   LatticeField U;
   
 
   explicit AdjointRep(GridBase* grid) : U(grid) {}
-  LatticeField 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]

lib/qcd/representations/fundamental.h

@@ -22,7 +22,8 @@ class FundamentalRep {
 
   // 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;
+  typedef typename SU<ncolour>::LatticeMatrix LatticeMatrix;
+  typedef LatticeGaugeField LatticeField;
   
   explicit FundamentalRep(GridBase* grid) {} //do nothing
   void update_representation(const LatticeGaugeField& Uin) {} // do nothing

lib/qcd/representations/hmc_types.h

@@ -1,61 +1,90 @@
 #ifndef HMC_TYPES_H
 #define HMC_TYPES_H
 
-#include <tuple>
-#include <utility>
 #include <qcd/representations/adjoint.h>
 #include <qcd/representations/fundamental.h>
+#include <tuple>
+#include <utility>
 
 namespace Grid {
 namespace QCD {
 
-
 // Supported types
-//enum {Fundamental, Adjoint} repr_type;
+// enum {Fundamental, Adjoint} repr_type;
 
-// Utility to add support to the HMC for representations other than the fundamental
+// Utility to add support to the HMC for representations other than the
-template<class... Reptypes>
+// fundamental
-class Representations{
+template <class... Reptypes>
-public:
+class Representations {
+ public:
   typedef std::tuple<Reptypes...> Representation_type;
 
+  // Size of the tuple, known at compile time
+  static const int tuple_size = sizeof...(Reptypes);
+  // The collection of types for the gauge fields
+  typedef std::tuple<typename Reptypes::LatticeField...> Representation_Fields;
+
   // To access the Reptypes (FundamentalRepresentation, AdjointRepresentation)
   template <std::size_t N>
-  using repr_type = typename std::tuple_element<N, Representation_type >::type;
+  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
-  explicit Representations(GridBase *grid):rep(Reptypes(grid)...){};
+  explicit Representations(GridBase* grid) : rep(Reptypes(grid)...){};
 
-  int size(){
+  int size() { return tuple_size; }
-    return std::tuple_size< Representation_type >::value;
-  }
 
   // update the fields
   template <std::size_t I = 0>
-  inline typename std::enable_if< I == sizeof...(Reptypes), void >::type update(LatticeGaugeField& U) {}
+  inline typename std::enable_if<(I == tuple_size), void>::type update(
+      LatticeGaugeField& U) {}
 
   template <std::size_t I = 0>
-      inline typename std::enable_if <
+  inline typename std::enable_if<(I < tuple_size), void>::type update(
-      I<sizeof...(Reptypes), void>::type update(LatticeGaugeField& U) {
+      LatticeGaugeField& U) {
     std::get<I>(rep).update_representation(U);
     update<I + 1>(U);
-  }  
+  }
+
+
+
 };
 
+typedef Representations<FundamentalRepresentation> NoHirep;
 
-typedef Representations<FundamentalRepresentation> JustTheFundamental;
+// Helper classes to access the elements
+// Strips the first N parameters from the tuple
+// sequence of classes to obtain the S sequence
+// Creates a type that is a tuple of vectors of the template type A
+template <template <typename> class A, class TupleClass,
+          size_t N = TupleClass::tuple_size, size_t... S>
+struct AccessTypes : AccessTypes<A, TupleClass, N - 1, N - 1, S...> {};
 
+template <template <typename> class A, class TupleClass, size_t... S>
+struct AccessTypes<A, TupleClass, 0, S...> {
+ public:
+  typedef typename TupleClass::Representation_Fields Rfields;
 
+  template <std::size_t N>
+  using elem = typename std::tuple_element<N, Rfields>::type;  // fields types
+
+  typedef std::tuple<std::vector< A< elem<S> >* > ... > VectorCollection;
+  typedef std::tuple< A< elem<S> >* ... > ClassCollection;
+
+  // Debug
+  void return_size() {
+    std::cout << GridLogMessage
+              << "Access:" << std::tuple_size<std::tuple<elem<S>...> >::value
+              << "\n";
+    std::cout << GridLogMessage
+              << "Access vectors:" << std::tuple_size<VectorCollection>::value
+              << "\n";
+  }
+};
 }
 }
 
-
-
 #endif
-
-
-

lib/qcd/utils/SUnAdjoint.h

@@ -48,6 +48,15 @@ class SU_Adjoint : public SU<ncolour> {
   typedef Lattice<vAMatrixF> LatticeAdjMatrixF;
   typedef Lattice<vAMatrixD> LatticeAdjMatrixD;
 
+  typedef Lattice<iVector<iScalar<iMatrix<vComplex, Dimension> >, Nd> >
+      LatticeAdjField;
+  typedef Lattice<iVector<iScalar<iMatrix<vComplexF, Dimension> >, Nd> >
+      LatticeAdjFieldF;
+  typedef Lattice<iVector<iScalar<iMatrix<vComplexD, Dimension> >, Nd> >
+      LatticeAdjFieldD;
+
+
+
 
   template <class cplx>
   static void generator(int Index, iSUnAdjointMatrix<cplx> &iAdjTa) {

tests/Test_hmc_WilsonAdjointFermionGauge.cc

@@ -30,6 +30,7 @@ directory
 *************************************************************************************/
 /*  END LEGAL */
 #include "Grid.h"
+//#include "qcd/hmc/HmcRunner.h"
 
 using namespace std;
 using namespace Grid;
@@ -38,7 +39,11 @@ using namespace Grid::QCD;
 namespace Grid {
 namespace QCD {
 
-class HmcRunner : public NerscHmcRunner {
+// Here change the allowed (higher) representations
+typedef Representations< FundamentalRepresentation, FundamentalRepresentation > TheRepresentations;
+
+
+class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > {
  public:
   void BuildTheAction(int argc, char **argv)
 
@@ -69,13 +74,13 @@ class HmcRunner : public NerscHmcRunner {
     TwoFlavourPseudoFermionAction<ImplPolicy> Nf2(FermOp, CG, CG);
 
     // Set smearing (true/false), default: false
-    Nf2.is_smeared = true;
+    Nf2.is_smeared = false;
 
     // Collect actions
-    ActionLevel<LatticeGaugeField> Level1(1);
+    ActionLevelHirep<LatticeGaugeField, TheRepresentations > Level1(1);
     Level1.push_back(&Nf2);
 
-    ActionLevel<LatticeGaugeField> Level2(4);
+    ActionLevelHirep<LatticeGaugeField, TheRepresentations > Level2(4);
     Level2.push_back(&Waction);
 
     TheAction.push_back(Level1);