From fe5b23e144564e1dcce3e574222b12f80108d4bd Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Sun, 19 Jul 2020 01:42:31 -0400 Subject: [PATCH] RMHMC implementation, originaly from Guido Cossu --- Grid/qcd/hmc/GenericHMCrunner.h | 12 +- Grid/qcd/hmc/HMC.h | 1 + Grid/qcd/hmc/integrators/Integrator.h | 123 ++++++++++++++ .../hmc/integrators/Integrator_algorithm.h | 139 +++++++++++++++- tests/hmc/Test_hmc_WilsonGauge_Implicit.cc | 154 ++++++++++++++++++ 5 files changed, 425 insertions(+), 4 deletions(-) create mode 100644 tests/hmc/Test_hmc_WilsonGauge_Implicit.cc diff --git a/Grid/qcd/hmc/GenericHMCrunner.h b/Grid/qcd/hmc/GenericHMCrunner.h index c2443dd0..4284004c 100644 --- a/Grid/qcd/hmc/GenericHMCrunner.h +++ b/Grid/qcd/hmc/GenericHMCrunner.h @@ -140,7 +140,17 @@ private: // Can move this outside? typedef IntegratorType TheIntegrator; - TheIntegrator MDynamics(UGrid, Parameters.MD, TheAction, Smearing); + // Metric + //TrivialMetric Mtr; + ConjugateGradient CG(1.0e-8,10000); + LaplacianParams LapPar(0.0001, 1.0, 10000, 1e-8, 12, 64); +// RealD Kappa = 1.2; + RealD Kappa = Parameters.Kappa; + std::cout << GridLogMessage << "Kappa = " << Kappa << std::endl; + + // Better to pass the generalised momenta to the integrator + LaplacianAdjointField Laplacian(UGrid, CG, LapPar, Kappa); + TheIntegrator MDynamics(UGrid, Parameters.MD, TheAction, Smearing, Laplacian); if (Parameters.StartingType == "HotStart") { // Hot start diff --git a/Grid/qcd/hmc/HMC.h b/Grid/qcd/hmc/HMC.h index 0f933204..20f50ee6 100644 --- a/Grid/qcd/hmc/HMC.h +++ b/Grid/qcd/hmc/HMC.h @@ -53,6 +53,7 @@ struct HMCparameters: Serializable { bool, MetropolisTest, Integer, NoMetropolisUntil, std::string, StartingType, + RealD, Kappa, IntegratorParameters, MD) HMCparameters() { diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h index 91af7372..7c38d970 100644 --- a/Grid/qcd/hmc/integrators/Integrator.h +++ b/Grid/qcd/hmc/integrators/Integrator.h @@ -137,6 +137,80 @@ protected: as[level].apply(update_P_hireps, Representations, Mom, U, ep); } + void implicit_update_P(Field& U, int level, double ep, bool intermediate = false) { + t_P[level] += ep; + + std::cout << GridLogIntegrator << "[" << level << "] P " + << " dt " << ep << " : t_P " << t_P[level] << std::endl; + // Fundamental updates, include smearing + MomentaField Msum(P.Mom._grid); + Msum = Zero(); + for (int a = 0; a < as[level].actions.size(); ++a) { + // Compute the force terms for the lagrangian part + // We need to compute the derivative of the actions + // only once + Field force(U._grid); + conformable(U._grid, P.Mom._grid); + Field& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared); + as[level].actions.at(a)->deriv(Us, force); // deriv should NOT include Ta + + std::cout << GridLogIntegrator << "Smearing (on/off): " << as[level].actions.at(a)->is_smeared << std::endl; + if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force); + force = FieldImplementation::projectForce(force); // Ta for gauge fields + Real force_abs = std::sqrt(norm2(force) / U._grid->gSites()); + std::cout << GridLogIntegrator << "|Force| site average: " << force_abs + << std::endl; + Msum += force; + } + + MomentaField NewMom = P.Mom; + MomentaField OldMom = P.Mom; + double threshold = 1e-8; + P.M.ImportGauge(U); + MomentaField MomDer(P.Mom._grid); + MomentaField MomDer1(P.Mom._grid); + MomentaField AuxDer(P.Mom._grid); + MomDer1 = Zero(); + MomentaField diff(P.Mom._grid); + double factor = 2.0; + if (intermediate){ + P.DerivativeU(P.Mom, MomDer1); + factor = 1.0; + } + + // Auxiliary fields + P.update_auxiliary_momenta(ep*0.5); + P.AuxiliaryFieldsDerivative(AuxDer); + Msum += AuxDer; + + + // Here run recursively + int counter = 1; + RealD RelativeError; + do { + std::cout << GridLogIntegrator << "UpdateP implicit step "<< counter << std::endl; + + // Compute the derivative of the kinetic term + // with respect to the gauge field + P.DerivativeU(NewMom, MomDer); + Real force_abs = std::sqrt(norm2(MomDer) / U._grid->gSites()); + std::cout << GridLogIntegrator << "|Force| laplacian site average: " << force_abs + << std::endl; + + NewMom = P.Mom - ep* 0.5 * (2.0*Msum + factor*MomDer + MomDer1);// simplify + diff = NewMom - OldMom; + counter++; + RelativeError = std::sqrt(norm2(diff))/std::sqrt(norm2(NewMom)); + std::cout << GridLogIntegrator << "UpdateP RelativeError: " << RelativeError << std::endl; + OldMom = NewMom; + } while (RelativeError > threshold); + + P.Mom = NewMom; + + // update the auxiliary fields momenta + P.update_auxiliary_momenta(ep*0.5); + } + void update_U(Field& U, double ep) { update_U(P, U, ep); @@ -158,6 +232,55 @@ protected: Representations.update(U); // void functions if fundamental representation } + void implicit_update_U(Field&U, double ep){ + t_U += ep; + int fl = levels - 1; + std::cout << GridLogIntegrator << " " << "[" << fl << "] U " << " dt " << ep << " : t_U " << t_U << std::endl; + + MomentaField Mom1(P.Mom._grid); + MomentaField Mom2(P.Mom._grid); + RealD RelativeError; + Field diff(U._grid); + Real threshold = 1e-8; + int counter = 1; + int MaxCounter = 100; + + Field OldU = U; + Field NewU = U; + + P.M.ImportGauge(U); + P.DerivativeP(Mom1); // first term in the derivative + + P.update_auxiliary_fields(ep*0.5); + + + MomentaField sum=Mom1; + do { + std::cout << GridLogIntegrator << "UpdateU implicit step "<< counter << std::endl; + + P.DerivativeP(Mom2); // second term in the derivative, on the updated U + sum = (Mom1 + Mom2); + + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + auto Pmu = PeekIndex(sum, mu); + Umu = expMat(Pmu, ep * 0.5, 12) * Umu; + PokeIndex(NewU, ProjectOnGroup(Umu), mu); + } + + diff = NewU - OldU; + RelativeError = std::sqrt(norm2(diff))/std::sqrt(norm2(NewU)); + std::cout << GridLogIntegrator << "UpdateU RelativeError: " << RelativeError << std::endl; + + P.M.ImportGauge(NewU); + OldU = NewU; // some redundancy to be eliminated + counter++; + } while (RelativeError > threshold && counter < MaxCounter); + + U = NewU; + P.update_auxiliary_fields(ep*0.5); + } + virtual void step(Field& U, int level, int first, int last) = 0; public: diff --git a/Grid/qcd/hmc/integrators/Integrator_algorithm.h b/Grid/qcd/hmc/integrators/Integrator_algorithm.h index b05c4ea8..4a86a755 100644 --- a/Grid/qcd/hmc/integrators/Integrator_algorithm.h +++ b/Grid/qcd/hmc/integrators/Integrator_algorithm.h @@ -101,7 +101,7 @@ public: std::string integrator_name(){return "LeapFrog";} - LeapFrog(GridBase* grid, IntegratorParameters Par, ActionSet& Aset, SmearingPolicy& Sm) + LeapFrog(GridBase* grid, IntegratorParameters Par, ActionSet& Aset, SmearingPolicy& Sm, Metric& M) : Integrator(grid, Par, Aset, Sm){}; void step(Field& U, int level, int _first, int _last) { @@ -144,7 +144,7 @@ private: public: INHERIT_FIELD_TYPES(FieldImplementation); - MinimumNorm2(GridBase* grid, IntegratorParameters Par, ActionSet& Aset, SmearingPolicy& Sm) + MinimumNorm2(GridBase* grid, IntegratorParameters Par, ActionSet& Aset, SmearingPolicy& Sm, Metric& M) : Integrator(grid, Par, Aset, Sm){}; std::string integrator_name(){return "MininumNorm2";} @@ -207,7 +207,7 @@ public: // Looks like dH scales as dt^4. tested wilson/wilson 2 level. ForceGradient(GridBase* grid, IntegratorParameters Par, ActionSet& Aset, - SmearingPolicy& Sm) + SmearingPolicy& Sm, Metric& M) : Integrator( grid, Par, Aset, Sm){}; @@ -271,6 +271,139 @@ public: } }; +//////////////////////////////// +// Riemannian Manifold HMC +// Girolami et al +//////////////////////////////// + + + +// correct +template > +class ImplicitLeapFrog : public Integrator { + public: + typedef ImplicitLeapFrog + Algorithm; + INHERIT_FIELD_TYPES(FieldImplementation); + + // Riemannian manifold metric operator + // Hermitian operator Fisher + + std::string integrator_name(){return "ImplicitLeapFrog";} + + ImplicitLeapFrog(GridBase* grid, IntegratorParameters Par, + ActionSet& Aset, SmearingPolicy& Sm, Metric& M) + : Integrator( + grid, Par, Aset, Sm, M){}; + + void step(Field& U, int level, int _first, int _last) { + int fl = this->as.size() - 1; + // level : current level + // fl : final level + // eps : current step size + + // Get current level step size + RealD eps = this->Params.trajL/this->Params.MDsteps; + for (int l = 0; l <= level; ++l) eps /= this->as[l].multiplier; + + int multiplier = this->as[level].multiplier; + for (int e = 0; e < multiplier; ++e) { + int first_step = _first && (e == 0); + int last_step = _last && (e == multiplier - 1); + + if (first_step) { // initial half step + this->implicit_update_P(U, level, eps / 2.0); + } + + if (level == fl) { // lowest level + this->implicit_update_U(U, eps); + } else { // recursive function call + this->step(U, level + 1, first_step, last_step); + } + + //int mm = last_step ? 1 : 2; + if (last_step){ + this->update_P(U, level, eps / 2.0); + } else { + this->implicit_update_P(U, level, eps, true);// works intermediate step + // this->update_P(U, level, eps); // looks not reversible + } + } + } +}; + + +// This is not completely tested +template > +class ImplicitMinimumNorm2 : public Integrator { + private: + const RealD lambda = 0.1931833275037836; + + public: + INHERIT_FIELD_TYPES(FieldImplementation); + + ImplicitMinimumNorm2(GridBase* grid, IntegratorParameters Par, + ActionSet& Aset, SmearingPolicy& Sm, Metric& M) + : Integrator( + grid, Par, Aset, Sm, M){}; + + std::string integrator_name(){return "ImplicitMininumNorm2";} + + void step(Field& U, int level, int _first, int _last) { + // level : current level + // fl : final level + // eps : current step size + + int fl = this->as.size() - 1; + + RealD eps = this->Params.trajL/this->Params.MDsteps * 2.0; + for (int l = 0; l <= level; ++l) eps /= 2.0 * this->as[l].multiplier; + + // Nesting: 2xupdate_U of size eps/2 + // Next level is eps/2/multiplier + + int multiplier = this->as[level].multiplier; + for (int e = 0; e < multiplier; ++e) { // steps per step + + int first_step = _first && (e == 0); + int last_step = _last && (e == multiplier - 1); + + if (first_step) { // initial half step + this->implicit_update_P(U, level, lambda * eps); + } + + if (level == fl) { // lowest level + this->implicit_update_U(U, 0.5 * eps); + } else { // recursive function call + this->step(U, level + 1, first_step, 0); + } + + this->implicit_update_P(U, level, (1.0 - 2.0 * lambda) * eps, true); + + if (level == fl) { // lowest level + this->implicit_update_U(U, 0.5 * eps); + } else { // recursive function call + this->step(U, level + 1, 0, last_step); + } + + //int mm = (last_step) ? 1 : 2; + //this->update_P(U, level, lambda * eps * mm); + + if (last_step) { + this->update_P(U, level, eps * lambda); + } else { + this->implicit_update_P(U, level, lambda * eps*2.0, true); + } + } + } +}; + NAMESPACE_END(Grid); #endif // INTEGRATOR_INCLUDED diff --git a/tests/hmc/Test_hmc_WilsonGauge_Implicit.cc b/tests/hmc/Test_hmc_WilsonGauge_Implicit.cc new file mode 100644 index 00000000..2abaed68 --- /dev/null +++ b/tests/hmc/Test_hmc_WilsonGauge_Implicit.cc @@ -0,0 +1,154 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc + +Copyright (C) 2015 + +Author: Peter Boyle +Author: neo +Author: Guido Cossu + +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 */ +#include + +namespace Grid{ +struct RMHMCActionParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS(RMHMCActionParameters, + double, gauge_beta) + + template + RMHMCActionParameters(Reader& Reader){ + read(Reader, "Action", *this); + } +}; +} + +int main(int argc, char **argv) { + using namespace Grid; + using namespace Grid::QCD; + + Grid_init(&argc, &argv); + GridLogIntegrator.Active(1); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef GenericHMCRunner HMCWrapper; // Uses the default minimum norm + // Serialiser + typedef Grid::JSONReader Serialiser; + + + HMCWrapper TheHMC; + TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file + + // Reader, file should come from command line + if (TheHMC.ParameterFile.empty()){ + std::cout << "Input file not specified." + << "Use --ParameterFile option in the command line.\nAborting" + << std::endl; + exit(1); + } + Serialiser Reader(TheHMC.ParameterFile); + + RMHMCActionParameters ActionParams(Reader); + + // Grid from the command line + TheHMC.Resources.AddFourDimGrid("gauge"); + + + // Checkpointer definition + CheckpointerParameters CPparams(Reader); +// TheHMC.Resources.LoadBinaryCheckpointer(CPparams); + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar(Reader); + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef PlaquetteMod PlaqObs; + typedef TopologicalChargeMod QObs; + TheHMC.Resources.AddObservable(); + TopologyObsParameters TopParams(Reader); + TheHMC.Resources.AddObservable(TopParams); + ///////////////////////////////////////////////////////////// + // Collect actions + WilsonGaugeActionR Waction(ActionParams.gauge_beta); + + + ActionLevel Level1(1); + Level1.push_back(&Waction); + TheHMC.TheAction.push_back(Level1); + ///////////////////////////////////////////////////////////// + TheHMC.Parameters.initialize(Reader); + TheHMC.Run(); + + Grid_finalize(); + +} // main + +/* Examples for input files + +JSON + +{ + "Checkpointer": { + "config_prefix": "ckpoint_json_lat", + "rng_prefix": "ckpoint_json_rng", + "saveInterval": 10, + "format": "IEEE64BIG" + }, + "RandomNumberGenerator": { + "serial_seeds": "1 2 3 4 6", + "parallel_seeds": "55 7 8 9 11" + }, + "Action":{ + "gauge_beta": 5.8 + }, + "TopologyMeasurement":{ + "interval": 1, + "do_smearing": true, + "Smearing":{ + "steps": 200, + "step_size": 0.01, + "meas_interval": 50, + "maxTau": 2.0 + } + }, + "HMC":{ + "StartTrajectory": 0, + "Trajectories": 10, + "MetropolisTest": true, + "NoMetropolisUntil": 10, + "StartingType": "HotStart", + "MD":{ + "name": "MinimumNorm2", + "MDsteps": 40, + "trajL": 1.0 + } + } +} + + +XML example not provided yet + +*/