Cleanup

lib/qcd/action/Actions.h

@@ -101,11 +101,10 @@ typedef SymanzikGaugeAction<ConjugateGimplF>        ConjugateSymanzikGaugeAction
 typedef SymanzikGaugeAction<ConjugateGimplD>        ConjugateSymanzikGaugeActionD;
 
 
-  typedef ScalarAction<ScalarImplR>                 ScalarActionR;
-  typedef ScalarAction<ScalarImplF>                 ScalarActionF;
-  typedef ScalarAction<ScalarImplD>                 ScalarActionD;
-  
-  
+typedef ScalarAction<ScalarImplR>                 ScalarActionR;
+typedef ScalarAction<ScalarImplF>                 ScalarActionF;
+typedef ScalarAction<ScalarImplD>                 ScalarActionD;
+
 }}
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////

lib/qcd/action/gauge/GaugeImplTypes.h

@@ -81,7 +81,7 @@ public:
   ///////////////////////////////////////////////////////////
   // Move these to another class
   // HMC auxiliary functions
-  static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){
+  static inline void generate_momenta(Field &P, GridParallelRNG &pRNG) {
     // specific for SU gauge fields
     LinkField Pmu(P._grid);
     Pmu = zero;
@@ -91,9 +91,7 @@ public:
     }
   }
 
-   static inline Field projectForce(Field& P){
-    return Ta(P);
-   }
+  static inline Field projectForce(Field &P) { return Ta(P); }
   
   static inline void update_field(Field& P, Field& U, double ep){
     for (int mu = 0; mu < Nd; mu++) {

lib/qcd/action/scalar/ScalarImpl.h

@@ -45,6 +45,49 @@ namespace Grid {
     }
     
   };
+
+  template <class S, unsigned int N>
+  class ScalarMatrixImplTypes {
+  public:
+    typedef S Simd;
+    
+    template <typename vtype>
+    using iImplField = iScalar<iScalar<iMatrix<vtype, N> > >;
+    
+    typedef iImplField<Simd> SiteField;
+    
+    
+    typedef Lattice<SiteField> Field;
+    
+    static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){
+      gaussian(pRNG, P);
+    }
+    
+    static inline Field projectForce(Field& P){return P;}
+    
+    static inline void update_field(Field& P, Field& U, double ep){
+      U += P*ep;
+    }
+    
+    static inline RealD FieldSquareNorm(Field& U){
+      return (TensorRemove(- sum(trace(U*U))*0.5).real());
+    }
+    
+    static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) {
+      gaussian(pRNG, U);
+    }
+    
+    static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) {
+      gaussian(pRNG, U);
+    }
+    
+    static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) {
+      U = 1.0;
+    }
+    
+  };
+
+
   
   
   typedef ScalarImplTypes<vReal> ScalarImplR;

lib/qcd/action/scalar/ScalarInteractionAction.h

@@ -0,0 +1,84 @@
+/*************************************************************************************
+
+  Grid physics library, www.github.com/paboyle/Grid
+
+  Source file: ./lib/qcd/action/gauge/WilsonGaugeAction.h
+
+  Copyright (C) 2015
+
+Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: neo <cossu@post.kek.jp>
+Author: paboyle <paboyle@ph.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
+  the Free Software Foundation; either version 2 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  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.,
+  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+  See the full license in the file "LICENSE" in the top level distribution
+directory
+  *************************************************************************************/
+/*  END LEGAL */
+
+#ifndef SCALAR_ACTION_H
+#define SCALAR_ACTION_H
+
+namespace Grid {
+  // FIXME drop the QCD namespace everywhere here
+  
+  template <class Impl>
+  class ScalarInteractionAction : public QCD::Action<typename Impl::Field> {
+  public:
+    INHERIT_FIELD_TYPES(Impl);
+    
+  private:
+    RealD mass_square;
+    RealD lambda;
+    
+  public:
+    ScalarAction(RealD ms, RealD l) : mass_square(ms), lambda(l){};
+
+    virtual std::string LogParameters(){
+      std::stringstream sstream;
+      sstream << GridLogMessage << "[ScalarAction] lambda      : " << lambda      << std::endl;
+      sstream << GridLogMessage << "[ScalarAction] mass_square : " << mass_square << std::endl;
+      return sstream.str();
+      
+    }
+    
+    virtual std::string action_name(){return "ScalarAction";}
+    
+    virtual void refresh(const Field &U,
+			 GridParallelRNG &pRNG){};  // noop as no pseudoferms
+    
+    virtual RealD S(const Field &p) {
+      return (mass_square * 0.5 + QCD::Nd) * ScalarObs<Impl>::sumphisquared(p) +
+	(lambda / 24.) * ScalarObs<Impl>::sumphifourth(p) +
+	ScalarObs<Impl>::sumphider(p);
+    };
+    
+    virtual void deriv(const Field &p,
+		       Field &force) {
+      Field tmp(p._grid);
+      Field p2(p._grid);
+      ScalarObs<Impl>::phisquared(p2, p);
+      tmp = -(Cshift(p, 0, -1) + Cshift(p, 0, 1));
+      for (int mu = 1; mu < QCD::Nd; mu++) tmp -= Cshift(p, mu, -1) + Cshift(p, mu, 1);
+      
+      force=+(mass_square + 2. * QCD::Nd) * p + (lambda / 6.) * p2 * p + tmp;
+    };
+  };
+  
+} // Grid
+
+#endif // SCALAR_ACTION_H

lib/qcd/hmc/GenericHMCrunner.h

@@ -94,8 +94,8 @@ class HMCWrapperTemplate: public HMCRunnerBase<ReaderClass> {
       Parameters.StartingType = arg;
     }
 
-    if (GridCmdOptionExists(argv, argv + argc, "--StartTrajectory")) {
-      arg = GridCmdOptionPayload(argv, argv + argc, "--StartTrajectory");
+    if (GridCmdOptionExists(argv, argv + argc, "--StartingTrajectory")) {
+      arg = GridCmdOptionPayload(argv, argv + argc, "--StartingTrajectory");
       std::vector<int> ivec(0);
       GridCmdOptionIntVector(arg, ivec);
       Parameters.StartTrajectory = ivec[0];

lib/qcd/hmc/HMC.h

@@ -40,6 +40,7 @@ directory
 #define HMC_INCLUDED
 
 #include <string>
+#include <list>
 
 namespace Grid {
 namespace QCD {
@@ -158,7 +159,7 @@ class HybridMonteCarlo {
   GridParallelRNG &pRNG; 
 
   Field &Ucur;
-
+  
   IntegratorType &TheIntegrator;
 	ObsListType Observables;
 
@@ -195,7 +196,7 @@ class HybridMonteCarlo {
   /////////////////////////////////////////////////////////
   // Evolution
   /////////////////////////////////////////////////////////
-  RealD evolve_step(Field &U) {
+  RealD evolve_hmc_step(Field &U) {
     TheIntegrator.refresh(U, pRNG);  // set U and initialize P and phi's
 
     RealD H0 = TheIntegrator.S(U);  // initial state action
@@ -211,11 +212,14 @@ class HybridMonteCarlo {
 
     std::cout.precision(17);
     std::cout << GridLogMessage << "Total H after trajectory  = " << H1
-              << "  dH = " << H1 - H0 << "\n";
+	      << "  dH = " << H1 - H0 << "\n";
     std::cout.precision(current_precision);
-
+    
     return (H1 - H0);
   }
+  
+
+  
 
  public:
   /////////////////////////////////////////
@@ -224,7 +228,7 @@ class HybridMonteCarlo {
   HybridMonteCarlo(HMCparameters _Pams, IntegratorType &_Int,
                    GridSerialRNG &_sRNG, GridParallelRNG &_pRNG, 
                    ObsListType _Obs, Field &_U)
-      : Params(_Pams), TheIntegrator(_Int), sRNG(_sRNG), pRNG(_pRNG), Observables(_Obs), Ucur(_U) {}
+    : Params(_Pams), TheIntegrator(_Int), sRNG(_sRNG), pRNG(_pRNG), Observables(_Obs), Ucur(_U) {}
   ~HybridMonteCarlo(){};
 
   void evolve(void) {
@@ -241,13 +245,13 @@ class HybridMonteCarlo {
       std::cout << GridLogMessage << "-- # Trajectory = " << traj << "\n";
       if (traj < Params.StartTrajectory + Params.NoMetropolisUntil) {
       	std::cout << GridLogMessage << "-- Thermalization" << std::endl;
-    	}
-
-    	double t0=usecond();
+      }
+      
+      double t0=usecond();
       Ucopy = Ucur;
 
-      DeltaH = evolve_step(Ucopy);
-
+      DeltaH = evolve_hmc_step(Ucopy);
+      // Metropolis-Hastings test
       bool accept = true;
       if (traj >= Params.StartTrajectory + Params.NoMetropolisUntil) {
         accept = metropolis_test(DeltaH);
@@ -255,9 +259,11 @@ class HybridMonteCarlo {
       	std::cout << GridLogMessage << "Skipping Metropolis test" << std::endl;
       }
 
-      if (accept) {
-        Ucur = Ucopy;
-      }
+      if (accept)
+        Ucur = Ucopy; 
+      
+     
+      
       double t1=usecond();
       std::cout << GridLogMessage << "Total time for trajectory (s): " << (t1-t0)/1e6 << std::endl;
 

lib/qcd/hmc/HMCResourceManager.h

@@ -75,6 +75,13 @@ class HMCResourceManager {
   bool have_RNG;
   bool have_CheckPointer;
 
+  void output_vector_string(const std::vector<std::string> &vs){
+    for (auto &i: vs)
+      std::cout << i << " ";
+    std::cout << std::endl;
+  }
+
+
  public:
   HMCResourceManager() : have_RNG(false), have_CheckPointer(false) {}
 
@@ -93,7 +100,11 @@ class HMCResourceManager {
     std::string cp_type;
     read(Read,"name", cp_type);
     std::cout << "Registered types " << std::endl;
-    std::cout << CPfactory.getBuilderList() << std::endl;
+    // NOTE: operator << is not overloaded for std::vector<string> 
+    // so it complains here
+    //std::cout << CPfactory.getBuilderList() << std::endl;
+    output_vector_string(CPfactory.getBuilderList());
+
 
     CP = CPfactory.create(cp_type, Read);
     CP->print_parameters();
@@ -110,7 +121,7 @@ class HMCResourceManager {
       std::string obs_type;
       read(Read,"name", obs_type);
       std::cout << "Registered types " << std::endl;
-      std::cout << ObsFactory.getBuilderList() << std::endl;
+      output_vector_string(ObsFactory.getBuilderList() );
 
       ObservablesList.emplace_back(ObsFactory.create(obs_type, Read));
       ObservablesList[ObservablesList.size() - 1]->print_parameters();
@@ -271,8 +282,8 @@ private:
     do{
       auto &ActionFactory = HMC_ActionModuleFactory<gauge_string, typename ImplementationPolicy::Field, ReaderClass>::getInstance(); 
       std::string action_type;
-      Read.readDefault("name", action_type);
-      std::cout << ActionFactory.getBuilderList() << std::endl;  // temporary
+      Read.readDefault("name", action_type); 
+      output_vector_string(ActionFactory.getBuilderList() );
       ActionsList.emplace(m, ActionFactory.create(action_type, Read));
     } while (Read.nextElement("Action"));
     ActionsList.find(m)->second->print_parameters();    

lib/qcd/hmc/integrators/Integrator.h

@@ -39,8 +39,8 @@ namespace QCD {
 
 class IntegratorParameters: Serializable {
 public:
-	GRID_SERIALIZABLE_CLASS_MEMBERS(IntegratorParameters,
-		std::string, name,      // name of the integrator
+        GRID_SERIALIZABLE_CLASS_MEMBERS(IntegratorParameters,
+                std::string, name,      // name of the integrator
     unsigned int, MDsteps,  // number of outer steps
     RealD, trajL,           // trajectory length
   )
@@ -55,7 +55,7 @@ public:
   template <class ReaderClass, typename std::enable_if<isReader<ReaderClass>::value, int >::type = 0 >
   IntegratorParameters(ReaderClass & Reader){
     std::cout << "Reading integrator\n";
-  	read(Reader, "Integrator", *this);
+        read(Reader, "Integrator", *this);
   }
 
   void print_parameters() const {
@@ -177,16 +177,16 @@ class Integrator {
   }
 
   void print_actions(){
-  	std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::" << std::endl;
-  	std::cout << GridLogMessage << "[Integrator] Action summary: "<<std::endl;
-  	for (int level = 0; level < as.size(); ++level) {
-  		std::cout << GridLogMessage << "[Integrator] ---- Level: "<< level << std::endl;
-  		for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
-  			std::cout << GridLogMessage << "["<< as[level].actions.at(actionID)->action_name() << "] ID: " << actionID << std::endl;
-  			std::cout << as[level].actions.at(actionID)->LogParameters();
-  		}
-  	}
-  	std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::"<< std::endl;
+        std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::" << std::endl;
+        std::cout << GridLogMessage << "[Integrator] Action summary: "<<std::endl;
+        for (int level = 0; level < as.size(); ++level) {
+                std::cout << GridLogMessage << "[Integrator] ---- Level: "<< level << std::endl;
+                for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
+                        std::cout << GridLogMessage << "["<< as[level].actions.at(actionID)->action_name() << "] ID: " << actionID << std::endl;
+                        std::cout << as[level].actions.at(actionID)->LogParameters();
+                }
+        }
+        std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::"<< std::endl;
 
   }
 
@@ -209,7 +209,7 @@ class Integrator {
     assert(P._grid == U._grid);
     std::cout << GridLogIntegrator << "Integrator refresh\n";
     FieldImplementation::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

lib/qcd/modules/ActionModules.h

@@ -510,28 +510,6 @@ class HMC_ActionModuleFactory
 };
 
 
-/*
-template <char const *str, class ReaderClass >
-class HMC_ScalarActionModuleFactory
-    : public Factory < HMC_ScalarActionModBase , Reader<ReaderClass> > {
- public:
-  typedef Reader<ReaderClass> TheReader; 
-  // use SINGLETON FUNCTOR MACRO HERE
-  HMC_ScalarActionModuleFactory(const HMC_ScalarActionModuleFactory& e) = delete;
-  void operator=(const HMC_ScalarActionModuleFactory& e) = delete;
-  static HMC_ScalarActionModuleFactory& getInstance(void) {
-    static HMC_ScalarActionModuleFactory e;
-    return e;
-  }
-
- private:
-  HMC_ScalarActionModuleFactory(void) = default;
-    std::string obj_type() const {
-        return std::string(str);
-  }
-};
-*/
-
 extern char gauge_string[];
 } // Grid
 

lib/qcd/representations/hmc_types.h

@@ -62,6 +62,7 @@ class Representations {
 
 typedef Representations<FundamentalRepresentation> NoHirep;
 typedef Representations<EmptyRep<typename ScalarImplR::Field> > ScalarFields;
+  //typedef Representations<EmptyRep<typename ScalarMatrixImplR::Field> > ScalarMatrixFields;
 
 // Helper classes to access the elements
 // Strips the first N parameters from the tuple