Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

some refactoring to get tempered trajectories working #5

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

some refactoring to get tempered trajectories working #5

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion adjointfield.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -222,7 +222,7 @@ template<typename Float> inline adjointu1<Float> operator*(const Float &x, const
}


template<typename Float, class Group> adjointfield<Float, su2> operator*(const Float &x, const adjointfield<Float, Group> &A) {
template<typename Float, class Group> adjointfield<Float, Group> operator*(const Float &x, const adjointfield<Float, Group> &A) {
adjointfield<Float, Group> res(A.getLx(), A.getLy(), A.getLz(), A.getLt());
for(size_t i = 0; i < A.getSize(); i++) {
res[i] = x * A[i];
Expand Down
9 changes: 5 additions & 4 deletions hmc-u1.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ int main(int ac, char* av[]) {
("nsteps", po::value<size_t>(&n_steps)->default_value(1000), "n_steps")
("tau", po::value<double>(&tau)->default_value(1.), "trajectory length tau")
("exponent", po::value<size_t>(&exponent)->default_value(0), "exponent for rounding")
("integrator", po::value<size_t>(&integs)->default_value(0), "itegration scheme to be used: 0=leapfrog, 1=lp_leapfrog, 2=omf4, 3=lp_omf4, 4=Euler, 5=RUTH, 6=omf2")
("integrator", po::value<size_t>(&integs)->default_value(0), "itegration scheme to be used: 0=leapfrog, 1=lp_leapfrog, 2=omf4, 3=lp_omf4, 4=Euler, 5=RUTH, 6=omf2, 7=temperedLeapfrog")
;

int err = parse_commandline(ac, av, desc, gparams);
Expand Down Expand Up @@ -71,7 +71,7 @@ int main(int ac, char* av[]) {
cout << "## Plaquette after rnd trafo: " << plaquette*normalisation << endl;

// Molecular Dynamics parameters
md_params mdparams(n_steps, tau);
md_params<std::mt19937> mdparams(n_steps, tau);

// generate list of monomials
gaugemonomial<double, _u1> gm(0);
Expand All @@ -82,7 +82,7 @@ int main(int ac, char* av[]) {
monomial_list.push_back(&gm);
monomial_list.push_back(&km);

integrator<double, _u1> * md_integ = set_integrator<double, _u1>(integs, exponent);
integrator<double, _u1, std::mt19937> * md_integ = set_integrator<double, _u1, std::mt19937>(integs, exponent);

std::ofstream os;
if(gparams.icounter == 0)
Expand All @@ -98,8 +98,9 @@ int main(int ac, char* av[]) {
}
// PRNG engine
std::mt19937 engine(gparams.seed+i);
mdparams.setengine(&engine);
// perform the MD update
md_update(U, engine, mdparams, monomial_list, *md_integ);
md_update(U, mdparams, monomial_list, *md_integ);

double energy = gauge_energy(U);
double E = 0., Q = 0.;
Expand Down
9 changes: 5 additions & 4 deletions hmc.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ int main(int ac, char* av[]) {
("nsteps", po::value<size_t>(&n_steps)->default_value(1000), "n_steps")
("tau", po::value<double>(&tau)->default_value(1.), "trajectory length tau")
("exponent", po::value<size_t>(&exponent)->default_value(0), "exponent for rounding")
("integrator", po::value<size_t>(&integs)->default_value(0), "itegration scheme to be used: 0=leapfrog, 1=lp_leapfrog, 2=omf4, 3=lp_omf4, 4=Euler, 5=RUTH, 6=omf2")
("integrator", po::value<size_t>(&integs)->default_value(0), "itegration scheme to be used: 0=leapfrog, 1=lp_leapfrog, 2=omf4, 3=lp_omf4, 4=Euler, 5=RUTH, 6=omf2, 7=temperedLeapfrog")
;

int err = parse_commandline(ac, av, desc, gparams);
Expand Down Expand Up @@ -69,7 +69,7 @@ int main(int ac, char* av[]) {
cout << "## Plaquette after rnd trafo: " << plaquette*normalisation << endl;

// Molecular Dynamics parameters
md_params mdparams(n_steps, tau);
md_params<std::mt19937> mdparams(n_steps, tau);

// generate list of monomials
gaugemonomial<double, su2> gm(0);
Expand All @@ -80,7 +80,7 @@ int main(int ac, char* av[]) {
monomial_list.push_back(&gm);
monomial_list.push_back(&km);

integrator<double, su2> * md_integ = set_integrator<double, su2>(integs, exponent);
integrator<double, su2, std::mt19937> * md_integ = set_integrator<double, su2, std::mt19937>(integs, exponent);

std::ofstream os;
if(gparams.icounter == 0)
Expand All @@ -96,8 +96,9 @@ int main(int ac, char* av[]) {
}
// PRNG engine
std::mt19937 engine(gparams.seed+i);
mdparams.setengine(&engine);
// perform the MD update
md_update(U, engine, mdparams, monomial_list, *md_integ);
md_update(U, mdparams, monomial_list, *md_integ);

double energy = gauge_energy(U);
rate += mdparams.getaccept();
Expand Down
109 changes: 80 additions & 29 deletions integrator.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,23 +10,26 @@
#include<list>
#include<iostream>
#include<cmath>
#include<random>

enum integrators { LEAPFROG = 0, LP_LEAPFROG = 1, OMF4 = 2, LP_OMF4 = 3, EULER = 4, RUTH = 5, OMF2 = 6};
enum integrators { LEAPFROG = 0, LP_LEAPFROG = 1, OMF4 = 2, LP_OMF4 = 3, EULER = 4, RUTH = 5, OMF2 = 6, T_LEAPFROG=7};

// virtual integrator class
template<typename Float, class Group> class integrator{
template<typename Float, class Group, class URNG> class integrator{
public:
integrator() {}
virtual void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h,
md_params const &params, const bool restore = true) = 0;
virtual void integrate(std::list<monomial<Float, Group>*> &monomial_list,
hamiltonian_field<Float, Group> &h,
md_params<URNG> &params, const bool restore = true) = 0;
};

// leapfrog integration scheme
template<typename Float, class Group> class leapfrog : public integrator<Float, Group> {
template<typename Float, class Group, class URNG> class leapfrog : public integrator<Float, Group, URNG> {
public:
leapfrog() {}
void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h,
md_params const &params, const bool restore=true) {
void integrate(std::list<monomial<Float, Group>*> &monomial_list,
hamiltonian_field<Float, Group> &h,
md_params<URNG> &params, const bool restore=true) {

adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());

Expand All @@ -48,13 +51,57 @@ public:
}
};

// tempered leapfrog integration scheme
template<typename Float, class Group, class URNG> class tempered_leapfrog : public integrator<Float, Group, URNG> {
public:
tempered_leapfrog() {}
void integrate(std::list<monomial<Float, Group>*> &monomial_list,
hamiltonian_field<Float, Group> &h,
md_params<URNG> &params, const bool restore=true) {

const size_t nsteps = params.getnsteps();
// std::normal_distribution<Float> normal(1., params.gettemperedwidth());
std::normal_distribution<Float> normal(1., params.gettemperedwidth());
Float a = Float(normal(*params.getengine()));
std::cout << "tempered a " << a << std::endl;
adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());

bool is_even = (params.getnsteps()%2 == 0);
Float dtau = params.gettau()/Float(nsteps);
// initial half-step for the momenta
deriv = sqrt(a) * deriv;
update_momenta(monomial_list, deriv, h, dtau/2.);
// first full step for gauge
update_gauge(h, dtau);
// nsteps-1 full steps
for(size_t i = 0; i < nsteps-1; i++) {
update_momenta(monomial_list, deriv, h, dtau/2);
if ((is_even && i < nsteps/2-1) || (!is_even && i < nsteps/2)) {
deriv = a * deriv;
}
if ((i >= nsteps/2)) {
deriv = (1./a) * deriv;
}
update_momenta(monomial_list, deriv, h, dtau/2);
update_gauge(h, dtau);
}

// final half-step for the momenta
update_momenta(monomial_list, deriv, h, dtau/2.);
deriv = (1./sqrt(a)) * deriv;
// restore SU
if(restore) h.U->restoreSU();
}
};


// leapfrog with rounding to allow for lower precision during MD update
template<typename Float, class Group> class lp_leapfrog : public integrator<Float, Group> {
template<typename Float, class Group, class URNG> class lp_leapfrog : public integrator<Float, Group, URNG> {
public:
lp_leapfrog() : n_prec(16) {}
lp_leapfrog(size_t n) : n_prec(n) {}
void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h,
md_params const &params, const bool restore = true) {
md_params<URNG> &params, const bool restore = true) {
adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
const size_t N = pow(10, n_prec);

Expand All @@ -79,11 +126,11 @@ private:
};

// OMF2 integration scheme (also called 2MN, second order minimal norm)
template<typename Float, class Group> class omf2 : public integrator<Float, Group> {
template<typename Float, class Group, class URNG> class omf2 : public integrator<Float, Group, URNG> {
public:
omf2() : lambda(0.1938) {}
void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h,
md_params const &params, const bool restore = true) {
md_params<URNG> &params, const bool restore = true) {

adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
Float dtau = params.gettau()/Float(params.getnsteps());
Expand Down Expand Up @@ -115,11 +162,11 @@ private:


// OMF4 integration scheme
template<typename Float, class Group> class omf4 : public integrator<Float, Group> {
template<typename Float, class Group, class URNG> class omf4 : public integrator<Float, Group, URNG> {
public:
omf4() : rho(0.2539785108410595), theta(-0.03230286765269967), vartheta(0.08398315262876693), lambda(0.6822365335719091) {}
void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h,
md_params const &params, const bool restore = true) {
md_params<URNG> &params, const bool restore = true) {

adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
Float dtau = params.gettau()/Float(params.getnsteps());
Expand Down Expand Up @@ -154,12 +201,12 @@ private:
};

// OMF4 integration scheme in low precision
template<typename Float, class Group> class lp_omf4 : public integrator<Float, Group> {
template<typename Float, class Group, class URNG> class lp_omf4 : public integrator<Float, Group, URNG> {
public:
lp_omf4() : rho(0.2539785108410595), theta(-0.03230286765269967), vartheta(0.08398315262876693), lambda(0.6822365335719091), n_prec(16) {}
lp_omf4(size_t n) : rho(0.2539785108410595), theta(-0.03230286765269967), vartheta(0.08398315262876693), lambda(0.6822365335719091), n_prec(n) {}
void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h,
md_params const &params, const bool restore = true) {
md_params<URNG> &params, const bool restore = true) {

adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
const size_t N = pow(10, n_prec);
Expand Down Expand Up @@ -197,11 +244,11 @@ private:
};

// semi-implicit symplectic Euler integration scheme
template<typename Float, class Group> class euler : public integrator<Float, Group> {
template<typename Float, class Group, class URNG> class euler : public integrator<Float, Group, URNG> {
public:
euler() {}
void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h,
md_params const &params, const bool restore = true) {
md_params<URNG> &params, const bool restore = true) {

adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());
Float dtau = params.gettau()/Float(params.getnsteps());
Expand All @@ -216,11 +263,11 @@ public:
};

// third order symplectic, but non-reversible integrator (Ruth)
template<typename Float, class Group> class ruth : public integrator<Float, Group> {
template<typename Float, class Group, class URNG> class ruth : public integrator<Float, Group, URNG> {
public:
ruth() {}
void integrate(std::list<monomial<Float, Group>*> &monomial_list, hamiltonian_field<Float, Group> &h,
md_params const &params, const bool restore = true) {
md_params<URNG> &params, const bool restore = true) {

adjointfield<Float, Group> deriv(h.U->getLx(), h.U->getLy(), h.U->getLz(), h.U->getLt(), h.U->getndims());

Expand All @@ -240,39 +287,43 @@ public:
};


template<typename Float, class Group> integrator<Float, Group>* set_integrator(const size_t integs, const size_t exponent) {
integrator<Float, Group> * integ;
template<typename Float, class Group, class URNG> integrator<Float, Group, URNG>* set_integrator(const size_t integs, const size_t exponent) {
integrator<Float, Group, URNG> * integ;
if(static_cast<integrators>(integs) == LEAPFROG) {
integ = new leapfrog<Float, Group>();
integ = new leapfrog<Float, Group, URNG>();
std::cerr << "leapfrog" << std::endl;
}
else if(static_cast<integrators>(integs) == LP_LEAPFROG) {
integ = new lp_leapfrog<Float, Group>(exponent);
integ = new lp_leapfrog<Float, Group, URNG>(exponent);
std::cerr << "lp_leapfrog" << std::endl;
}
else if(static_cast<integrators>(integs) == OMF4) {
integ = new omf4<Float, Group>();
integ = new omf4<Float, Group, URNG>();
std::cerr << "omf4" << std::endl;
}
else if(static_cast<integrators>(integs) == LP_OMF4) {
integ = new lp_omf4<Float, Group>(exponent);
integ = new lp_omf4<Float, Group, URNG>(exponent);
std::cerr << "lp_omf4" << std::endl;
}
else if(static_cast<integrators>(integs) == EULER) {
integ = new euler<Float, Group>();
integ = new euler<Float, Group, URNG>();
std::cerr << "euler" << std::endl;
}
else if(static_cast<integrators>(integs) == RUTH) {
integ = new ruth<Float, Group>();
integ = new ruth<Float, Group, URNG>();
std::cerr << "ruth" << std::endl;
}
else if(static_cast<integrators>(integs) == OMF2) {
integ = new omf2<Float, Group>();
integ = new omf2<Float, Group, URNG>();
std::cerr << "omf2" << std::endl;
}
else if(static_cast<integrators>(integs) == T_LEAPFROG) {
integ = new tempered_leapfrog<Float, Group, URNG>();
std::cerr << "tempered leapfrog" << std::endl;
}
else {
std::cerr << "Integrator does not match, using default" << std::endl;
integ = new leapfrog<Float, Group>();
integ = new leapfrog<Float, Group, URNG>();
}
return integ;
}
7 changes: 4 additions & 3 deletions kramers.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ int main(int ac, char* av[]) {
hotstart(U, gparams.seed, gparams.heat);
}
// Molecular Dynamics parameters
md_params mdparams(n_steps, tau);
md_params<std::mt19937> mdparams(n_steps, tau);
mdparams.setkmax(k_max);
mdparams.setgamma(gamma);

Expand All @@ -81,7 +81,7 @@ int main(int ac, char* av[]) {
monomial_list.push_back(&gm);
monomial_list.push_back(&km);

integrator<double, su2> * md_integ = set_integrator<double, su2>(integs, exponent);
integrator<double, su2, std::mt19937> * md_integ = set_integrator<double, su2, std::mt19937>(integs, exponent);

mdparams.setkmax(5);
if(!gparams.acceptreject) mdparams.disableacceptreject();
Expand All @@ -98,8 +98,9 @@ int main(int ac, char* av[]) {

// PRNG engine
std::mt19937 engine(gparams.seed + i);
mdparams.setengine(&engine);
// perform the MD update
kramers_md_update(U, engine, mdparams, monomial_list, *md_integ);
kramers_md_update(U, mdparams, monomial_list, *md_integ);

double energy = gauge_energy(U);
rate += mdparams.getaccept();
Expand Down
12 changes: 6 additions & 6 deletions kramers_md_update.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,14 +18,14 @@ using std::vector;


template<class URNG, typename Float, class Group> void kramers_md_update(gaugeconfig<su2> &U,
URNG &engine,
md_params &params,
md_params<URNG> &params,
std::list<monomial<Float, Group>*> &monomial_list,
integrator<Float, Group> &md_integ) {
integrator<Float, Group, URNG> &md_integ) {
adjointfield<Float, Group> momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims());
URNG * engine = params.getengine();
// generate standard normal distributed random momenta
// normal distribution checked!
initnormal(engine, momenta);
initnormal(*engine, momenta);

// for the accept reject step
std::uniform_real_distribution<Float> uniform(0., 1.);
Expand All @@ -37,7 +37,7 @@ template<class URNG, typename Float, class Group> void kramers_md_update(gaugeco

for(size_t k = 0; k < params.getkmax(); k++) {
// first momenta update
initnormal<URNG, Float>(engine, eta);
initnormal<URNG, Float>(*engine, eta);
Float epsilon = params.gettau()/Float(params.getnsteps());
for(size_t i = 0; i < momenta.getSize(); i++) {
momenta[i].seta(momenta[i].geta()*exp(-gamma*epsilon) + sqrt(1 - exp(-2*gamma*epsilon))*eta[i].geta());
Expand Down Expand Up @@ -77,7 +77,7 @@ template<class URNG, typename Float, class Group> void kramers_md_update(gaugeco
params.setaccept(true);
if(params.getacceptreject()) {
if(delta_H > 0) {
if(uniform(engine) > exp(-delta_H)) {
if(uniform(*engine) > exp(-delta_H)) {
params.setaccept(false);
}
}
Expand Down
4 changes: 2 additions & 2 deletions lyapunov.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,9 +22,9 @@ using std::vector;

template<typename Float, class URNG, class Group> void compute_lyapunov(gaugeconfig<Group> &U,
URNG &engine,
md_params params,
md_params<URNG> params,
std::list<monomial<Float, Group>*> &monomial_list,
integrator<Float, Group> &md_integ,
integrator<Float, Group, URNG> &md_integ,
std::string const &path,
const size_t d = 12) {

Expand Down
Loading