From 459d89d7614842d51d04a7cfa6941199d08bcaf9 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 29 Dec 2020 13:41:46 +0100 Subject: Initial commit. --- .gitignore | 3 + compile_flags.txt | 3 + langevin.cpp | 245 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 251 insertions(+) create mode 100644 .gitignore create mode 100644 compile_flags.txt create mode 100644 langevin.cpp diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0c0a063 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +* +!/**/ +!*.* diff --git a/compile_flags.txt b/compile_flags.txt new file mode 100644 index 0000000..02d9f07 --- /dev/null +++ b/compile_flags.txt @@ -0,0 +1,3 @@ +-std=c++20 +-O3 +-I/usr/include/eigen3 diff --git a/langevin.cpp b/langevin.cpp new file mode 100644 index 0000000..73cf4c0 --- /dev/null +++ b/langevin.cpp @@ -0,0 +1,245 @@ +#include +#include +#include +#include +#include +#include + +#include + +#include "pcg-cpp/include/pcg_random.hpp" +#include "randutils/randutils.hpp" + +#define P_SPIN_P 3 +const unsigned p = P_SPIN_P; + +using Tensor = Eigen::Tensor, p>; +using Scalar = Eigen::Tensor, 0>; +using Vector = Eigen::Tensor, 1>; +using Matrix = Eigen::Tensor, 2>; + +const std::array, 1> ip = {Eigen::IndexPair(0,0)}; + +using Rng = randutils::random_generator; + +template +class complex_normal_distribution { + private: + std::complex μ; + T σ; + std::complex κ; + std::normal_distribution dist; + + public: + complex_normal_distribution(std::complex μ, T σ, std::complex κ) : + μ(μ), σ(σ), κ(κ), dist(0, σ / sqrt(2)) {} + + template + std::complex operator()(Generator& g) { + std::complex x(dist(g) * sqrt(1 + std::abs(κ)), dist(g) * sqrt(1 - std::abs(κ))); + return μ + std::polar((T)1, std::arg(κ)) * x; + } +}; + +long unsigned factorial(unsigned p) { + if (p == 0) { + return 1; + } else { + return p * factorial(p - 1); + } +} + +Tensor generate_coupling_tensor(unsigned N, std::complex κ, Rng& r) { + double σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1))); + complex_normal_distribution<> dist(0, σ, κ); + + Tensor J(N, N, N); + + for (unsigned i1 = 0; i1 < N; i1++) { + for (unsigned i2 = i1; i2 < N; i2++) { + for (unsigned i3 = i2; i3 < N; i3++) { + std::complex x = dist(r.engine()); + + J(i1, i2, i3) = x; + J(i1, i3, i2) = x; + J(i2, i1, i3) = x; + J(i2, i3, i1) = x; + J(i3, i1, i2) = x; + J(i3, i2, i1) = x; + } + } + } + + return J; +} + +template +Eigen::Tensor, r> conj(const Eigen::Tensor, r>& t) { + return t.unaryExpr(static_cast (*)(const std::complex&)>(&std::conj)); +} + +template +double norm(const Eigen::Tensor, r>& t) { + Eigen::Tensor t2 = t.unaryExpr(static_cast&)>(&std::norm)).sum(); + return t2(0); +} + +Vector generate_initial_vector(unsigned N) { + Vector z(N); + z.setConstant(1); + return z; +} + +std::complex hamiltonian(const Tensor& J, const Vector& z) { + Eigen::Tensor, 0> t = z.contract(z.contract(z.contract(J, ip), ip), ip); + return t(0) / (double)factorial(p); +} + +Matrix identity(unsigned N) { + Matrix I(N, N); + for (unsigned i = 0; i < N; i++) { + I(i) = 1; + } + return I; +} + +std::tuple gradHess(const Tensor& J, const Vector& z, std::complex ε) { + Matrix J1 = z.contract(J, ip); + Matrix H = ((p - 1) * (double)p / factorial(p)) * J1 - ((double)p * ε) * identity(z.size()); + + Vector J2 = z.contract(J1, ip); + Vector g = ((double)p / factorial(p)) * J2 - ((double)p * ε) * z; + + return {g, H}; +} + +std::tuple, Vector> WdW(const Tensor& J, const Vector& z, std::complex ε) { + Vector grad; + Matrix hess; + + std::tie(grad, hess) = gradHess(J, z, ε); + + double W = norm(grad); + Matrix conjHess = conj(hess); + Scalar zz = z.pow(2).sum(); + Vector zc = conj(z); + + Vector dWdz = grad.contract(conjHess, ip) - (pow(p, 2) / 2.0 * ((double)z.size() - zz(0))) * zc; + Scalar dWdε = (-(double)p) * grad.contract(zc, ip); + Vector dεdz = (1 / (double)z.size()) * ((1 - 1/(double)p) * grad - ((double)p * ε) * z - (1 /(double)p) * z.contract(hess, ip)); + + return {W, dWdz, dWdε(0), dεdz}; + +} + +void gradientDescent(const Tensor& J, Vector& z, std::complex& ε, double δ, double γ) { + double W; + Vector dz; + std::complex dε; + + std::tie(W, dz, dε, std::ignore) = WdW(J, z, ε); + + while (W > δ * z.size()) { + std::cout << W << std::endl; + z = z - (γ / z.size()) * dz; + ε = ε - (γ / z.size()) * dε; + + std::tie(W, dz, dε, std::ignore) = WdW(J, z, ε); + } +} + +Vector generateKick(const Vector& z, Rng& r) { + Vector k(z.size()); + + for (unsigned i = 0; i < k.size(); i++) { + k(i) = std::complex(r.variate(0, 1), r.variate(0, 1)) / sqrt(2); + } + + Scalar kz = z.contract(k, ip); + k = k - kz(0) * z; + + return k; +} + +void langevin(const Tensor& J, Vector& z, std::complex& ε, double T, unsigned t, double γ, Rng& r) { + double W; + Vector dz; + Vector dεdz; + + std::tie(W, dz, std::ignore, dεdz) = WdW(J, z, ε); + + for (unsigned i = 0; i < t; i++) { + std::cout << W << std::endl; + + Vector η = generateKick(z, r); + Vector δz = (γ / z.size()) * (-dz + T * η); + Scalar δε = dεdz.contract(δz, ip); + + z = z + δz; + ε = ε + δε(0); + + std::tie(W, dz, std::ignore, dεdz) = WdW(J, z, ε); + } +} + +int main(int argc, char* argv[]) { + // model parameters + unsigned N = 10; // number of spins + double T = 1; // temperature + double Rκ = 1; // real part of distribution parameter + double Iκ = 0; // imaginary part of distribution parameter + + // simulation parameters + double δ = 1e-2; // threshold for determining saddle + double γ = 1e-2; // step size + + int opt; + + while ((opt = getopt(argc, argv, "N:T:E:e:r:i:g:")) != -1) { + switch (opt) { + case 'N': + N = (unsigned)atof(optarg); + break; + case 'T': + T = atof(optarg); + break; + case 'e': + δ = atof(optarg); + break; + case 'g': + γ = atof(optarg); + case 'r': + Rκ = atof(optarg); + break; + case 'i': + Iκ = atof(optarg); + break; + default: + exit(1); + } + } + + std::complex κ(Rκ, Iκ); + + Rng r; + + Tensor J = generate_coupling_tensor(N, κ, r); + Vector z = generate_initial_vector(N); + std::complex ε = hamiltonian(J, z) / (double)N; + + Scalar zz = z.pow(2).sum(); + + std::cout << zz(0) << std::endl; + + gradientDescent(J, z, ε, δ, γ); + langevin(J, z, ε, T, 1000, γ, r); + + zz = z.pow(2).sum(); + double a = norm(z); + + std::cout << zz(0) / (double)N << std::endl; + std::cout << ε << std::endl; + std::cout << hamiltonian(J, z) / (double)N << std::endl; + std::cout << a / N << std::endl; +} + -- cgit v1.2.3-70-g09d2