From 3276bdd1e9796fec71e169e6c41d77da72b3a4fb Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 25 Feb 2021 15:28:11 +0100 Subject: Many changes. --- complex_normal.hpp | 4 +- dynamics.hpp | 28 +++++----- langevin.cpp | 56 ++++++++++--------- p-spin.hpp | 11 ++-- stokes.hpp | 154 ++++++++++++++++++++++++++++++++++++++++++----------- stokes_test.cpp | 4 +- types.hpp | 3 ++ 7 files changed, 176 insertions(+), 84 deletions(-) create mode 100644 types.hpp diff --git a/complex_normal.hpp b/complex_normal.hpp index 8ed75d1..2ffcea1 100644 --- a/complex_normal.hpp +++ b/complex_normal.hpp @@ -6,8 +6,8 @@ template class complex_normal_distribution { private: std::normal_distribution d; - double δx; - double δy; + T δx; + T δy; std::complex μ; std::complex eθ; diff --git a/dynamics.hpp b/dynamics.hpp index d421d13..10b1be2 100644 --- a/dynamics.hpp +++ b/dynamics.hpp @@ -13,10 +13,10 @@ class gradientDescentStallException: public std::exception { } }; -template -std::tuple> gradientDescent(const Tensor& J, const Vector& z0, double ε, double γ0 = 1, double δγ = 2) { +template +std::tuple> gradientDescent(const Tensor& J, const Vector& z0, Real ε, Real γ0 = 1, Real δγ = 2) { Vector z = z0; - double γ = γ0; + Real γ = γ0; auto [W, dW] = WdW(J, z); @@ -42,12 +42,12 @@ std::tuple> gradientDescent(const Tensor& J, c return {W, z}; } -template -Vector findSaddle(const Tensor& J, const Vector& z0, double ε, double δW = 2, double γ0 = 1, double δγ = 2) { +template +Vector findSaddle(const Tensor& J, const Vector& z0, Real ε, Real δW = 2, Real γ0 = 1, Real δγ = 2) { Vector z = z0; Vector ζ = euclideanToStereographic(z); - double W; + Real W; std::tie(W, std::ignore) = WdW(J, z); Vector dH; @@ -61,7 +61,7 @@ Vector findSaddle(const Tensor& J, const Vector& z0, Vector ζNew = ζ - dζ; Vector zNew = stereographicToEuclidean(ζNew); - double WNew; + Real WNew; std::tie(WNew, std::ignore) = WdW(J, zNew); if (WNew < W) { // If Newton's step lowered the objective, accept it! @@ -90,20 +90,20 @@ Vector randomVector(unsigned N, Distribution d, Generator& r) { return z; } -template -std::tuple> metropolis(const Tensor& J, const Vector& z0, - std::function&, const Vector&)>& energy, - double T, double γ, unsigned N, Distribution d, Generator& r) { +template +std::tuple> metropolis(const Tensor& J, const Vector& z0, + std::function&, const Vector&)>& energy, + Real T, Real γ, unsigned N, Distribution d, Generator& r) { Vector z = z0; - double E = energy(J, z); + Real E = energy(J, z); - std::uniform_real_distribution D(0, 1); + std::uniform_real_distribution D(0, 1); for (unsigned i = 0; i < N; i++) { Vector zNew = normalize(z + γ * randomVector(z.size(), d, r)); - double ENew = energy(J, zNew); + Real ENew = energy(J, zNew); if (E - ENew > T * log(D(r))) { z = zNew; diff --git a/langevin.cpp b/langevin.cpp index 51cb2a0..8c191f3 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -14,14 +14,14 @@ #define PSPIN_P 3 const int p = PSPIN_P; // polynomial degree of Hamiltonian -using Complex = std::complex; +using Complex = std::complex; using ComplexVector = Vector; using ComplexMatrix = Matrix; using ComplexTensor = Tensor; using Rng = randutils::random_generator; -std::list> saddlesCloserThan(const std::unordered_map& saddles, double δ) { +std::list> saddlesCloserThan(const std::unordered_map& saddles, Real δ) { std::list> pairs; for (auto it1 = saddles.begin(); it1 != saddles.end(); it1++) { @@ -36,17 +36,17 @@ std::list> saddlesCloserThan(const std::unordered_m } template -std::tuple matchImaginaryEnergies(const ComplexTensor& J0, const ComplexVector& z10, const ComplexVector& z20, double ε, double Δ, Generator& r) { - double σ = sqrt(factorial(p) / (2.0 * pow(z10.size(), p - 1))); - complex_normal_distribution<> dJ(0, σ, 0); +std::tuple matchImaginaryEnergies(const ComplexTensor& J0, const ComplexVector& z10, const ComplexVector& z20, Real ε, Real Δ, Generator& r) { + Real σ = sqrt(factorial(p) / (2.0 * pow(z10.size(), p - 1))); + complex_normal_distribution dJ(0, σ, 0); ComplexTensor J = J0; Complex H1, H2; ComplexVector z1, z2; std::tie(H1, std::ignore, std::ignore) = hamGradHess(J, z10); std::tie(H2, std::ignore, std::ignore) = hamGradHess(J, z20); - double prevdist = abs(imag(H1-H2)); - double γ = 0.1 * prevdist; + Real prevdist = abs(imag(H1-H2)); + Real γ = 0.1 * prevdist; std::function)> perturbJ = [&γ, &dJ, &r] (ComplexTensor& JJ, std::array ind) { @@ -62,7 +62,7 @@ std::tuple matchImaginaryEnergies(c z1 = findSaddle(Jp, z10, Δ); z2 = findSaddle(Jp, z20, Δ); - double dist = (z1 - z2).norm(); + Real dist = (z1 - z2).norm(); std::tie(H1, std::ignore, std::ignore) = hamGradHess(Jp, z1); std::tie(H2, std::ignore, std::ignore) = hamGradHess(Jp, z2); @@ -87,16 +87,16 @@ std::tuple matchImaginaryEnergies(c 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 + Real T = 1; // temperature + Real Rκ = 1; // real part of distribution parameter + Real Iκ = 0; // imaginary part of distribution parameter // simulation parameters - double ε = 1e-12; - double εJ = 1e-5; - double δ = 1e-2; // threshold for determining saddle - double Δ = 1e-3; - double γ = 1e-2; // step size + Real ε = 1e-12; + Real εJ = 1e-5; + Real δ = 1e-2; // threshold for determining saddle + Real Δ = 1e-3; + Real γ = 1e-2; // step size unsigned t = 1000; // number of Langevin steps unsigned M = 100; unsigned n = 100; @@ -140,18 +140,18 @@ int main(int argc, char* argv[]) { } Complex κ(Rκ, Iκ); - double σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1))); + Real σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1))); Rng r; - complex_normal_distribution<> d(0, 1, 0); + complex_normal_distribution d(0, 1, 0); - ComplexTensor J = generateCouplings(N, complex_normal_distribution<>(0, σ, κ), r.engine()); + ComplexTensor J = generateCouplings(N, complex_normal_distribution(0, σ, κ), r.engine()); ComplexVector z = normalize(randomVector(N, d, r.engine())); - std::function energyNormGrad = [] + std::function energyNormGrad = [] (const ComplexTensor& J, const ComplexVector& z) { - double W; + Real W; std::tie(W, std::ignore) = WdW(J, z); return W; }; @@ -185,18 +185,16 @@ int main(int argc, char* argv[]) { ComplexVector z1 = nearbySaddles.front()[0]; ComplexVector z2 = nearbySaddles.front()[1]; - for (unsigned j = 2; j < 8; j++) { - std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, pow(10, -2.0 * j), ε, r); + std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, 1e-14, ε, r); - Rope stokes(10, z1, z2, J); + Rope stokes(10, z1, z2, J); - for (unsigned i = 0; i < 8; i++) { - stokes.relax(J, 1e5, 0.1); + for (unsigned i = 0; i < 9; i++) { + stokes.relaxDiscreteGradient(J, 1e6, 0.1, 0); - std::cout << stokes.n() << " " << stokes.cost(J) << " " << stokes.length() << " " << stokes.error(J) << std::endl; + std::cout << stokes.n() << " " << stokes.cost(J) << std::endl; - stokes = stokes.interpolate(); - } + stokes = stokes.interpolate(); } return 0; diff --git a/p-spin.hpp b/p-spin.hpp index bd3cacc..f1dc07f 100644 --- a/p-spin.hpp +++ b/p-spin.hpp @@ -2,12 +2,13 @@ #include +#include "types.hpp" #include "tensor.hpp" #include "factorial.hpp" template Vector normalize(const Eigen::MatrixBase& z) { - return z * sqrt((double)z.size() / (typename Derived::Scalar)(z.transpose() * z)); + return z * sqrt((Real)z.size() / (typename Derived::Scalar)(z.transpose() * z)); } template @@ -16,7 +17,7 @@ std::tuple, Matrix> hamGradHess(const Tensor Jzz = Jz * z; Scalar Jzzz = Jzz.transpose() * z; - double pBang = factorial(p); + Real pBang = factorial(p); Matrix hessian = ((p - 1) * p / pBang) * Jz; Vector gradient = (p / pBang) * Jzz; @@ -31,18 +32,18 @@ Vector zDot(const Vector& z, const Vector& dH) { } template -std::tuple> WdW(const Tensor& J, const Vector& z) { +std::tuple> WdW(const Tensor& J, const Vector& z) { Vector dH; Matrix ddH; std::tie(std::ignore, dH, ddH) = hamGradHess(J, z); Vector dzdt = zDot(z, dH); - double a = z.squaredNorm(); + Real a = z.squaredNorm(); Scalar A = (Scalar)(z.transpose() * dzdt) / a; Scalar B = dH.dot(z) / a; - double W = dzdt.squaredNorm(); + Real W = dzdt.squaredNorm(); Vector dW = ddH * (dzdt - A * z.conjugate()) + 2 * (conj(A) * B * z).real() - conj(B) * dzdt - conj(A) * dH.conjugate(); diff --git a/stokes.hpp b/stokes.hpp index 6cc5c4a..cbf5ebb 100644 --- a/stokes.hpp +++ b/stokes.hpp @@ -1,4 +1,6 @@ #include "p-spin.hpp" +#include "complex_normal.hpp" +#include "dynamics.hpp" #include class ropeRelaxationStallException: public std::exception { @@ -9,13 +11,13 @@ class ropeRelaxationStallException: public std::exception { template Vector variation(const Vector& z, const Vector& z´, const Vector& z´´, const Vector& dH, const Matrix& ddH) { - double z² = z.squaredNorm(); - double z´² = z´.squaredNorm(); + Real z² = z.squaredNorm(); + Real z´² = z´.squaredNorm(); Vector ż = zDot(z, dH); - double ż² = ż.squaredNorm(); + Real ż² = ż.squaredNorm(); - double Reż·z´ = real(ż.dot(z´)); + Real Reż·z´ = real(ż.dot(z´)); Matrix dż = (dH.conjugate() - (dH.dot(z) / z²) * z.conjugate()) * z.adjoint() / z²; Matrix dżc = -ddH + (ddH * z.conjugate()) * z.transpose() / z² @@ -33,7 +35,7 @@ Vector variation(const Vector& z, const Vector& z´, con ) * z.conjugate() ) / z²; - double dReż·z´ = real(ż.dot(z´´) + ż´.dot(z´)); + Real dReż·z´ = real(ż.dot(z´´) + ż´.dot(z´)); Vector ddtdLdz´ = - ( ( @@ -72,7 +74,7 @@ class Rope { } for (unsigned i = 0; i < N + 2; i++) { - z[i] = normalize(z.front() + (z.back() - z.front()) * ((double)i / (N + 1.0))); + z[i] = normalize(z.front() + (z.back() - z.front()) * ((Real)i / (N + 1.0))); } } @@ -80,8 +82,8 @@ class Rope { return z.size() - 2; } - double length() const { - double l = 0; + Real length() const { + Real l = 0; for (unsigned i = 0; i < z.size() - 1; i++) { l += (z[i + 1] - z[i]).norm(); @@ -91,14 +93,14 @@ class Rope { } template - double error(const Tensor& J) const { + Real error(const Tensor& J) const { Scalar H0, HN; std::tie(H0, std::ignore, std::ignore) = hamGradHess(J, z.front()); std::tie(HN, std::ignore, std::ignore) = hamGradHess(J, z.back()); - double ImH = imag((H0 + HN) / 2.0); + Real ImH = imag((H0 + HN) / 2.0); - double err = 0; + Real err = 0; for (unsigned i = 1; i < z.size() - 1; i++) { Scalar Hi; @@ -119,7 +121,7 @@ class Rope { } template - std::vector> δz(const Tensor& J) const { + std::vector> generateGradientδz(const Tensor& J) const { std::vector> δz(z.size()); #pragma omp parallel for @@ -137,7 +139,7 @@ class Rope { } // We return a δz with average norm of one. - double mag = 0; + Real mag = 0; for (unsigned i = 1; i < z.size() - 1; i++) { mag += δz[i].norm(); } @@ -150,20 +152,87 @@ class Rope { } template - double perturb(const Tensor& J, double δ0) { - std::vector> Δz = δz(J); + std::vector> generateDiscreteGradientδz(const Tensor& J, Real γ) const { + std::vector> δz(z.size()); + + std::vector> ż(z.size()); + std::vector> dH(z.size()); + std::vector> ddH(z.size()); + +#pragma omp parallel for + for (unsigned i = 1; i < z.size() - 1; i++) { + std::tie(std::ignore, dH[i], ddH[i]) = hamGradHess(J, z[i]); + ż[i] = zDot(z[i], dH[i]); + } + +#pragma omp parallel for + for (unsigned i = 1; i < z.size() - 1; i++) { + Real z² = z[i].squaredNorm(); + Matrix dż = (dH[i].conjugate() - (dH[i].dot(z[i]) / z²) * z[i].conjugate()) * z[i].adjoint() / z²; + Matrix dżc = -ddH[i] + (ddH[i] * z[i].conjugate()) * z[i].transpose() / z² + + (z[i].dot(dH[i]) / z²) * ( + Matrix::Identity(ddH[i].rows(), ddH[i].cols()) - z[i].conjugate() * z[i].transpose() / z² + ); + + Vector dC = -0.5 * (dżc * dz(i) + dż * dz(i).conjugate() - (dżc * ż[i] + dż * ż[i].conjugate()) * real(ż[i].dot(dz(i))) / ż[i].squaredNorm()) / (ż[i].norm() * dz(i).norm()); + + if (i > 1) { + dC += -0.5 * (ż[i - 1].conjugate() - dz(i - 1).conjugate() * real(ż[i - 1].dot(dz(i - 1))) / dz(i - 1).squaredNorm()) / (ż[i - 1].norm() * dz(i - 1).norm()); + } + + if (i < z.size() - 2) { + dC += 0.5 * (ż[i + 1].conjugate() - dz(i + 1).conjugate() * real(ż[i + 1].dot(dz(i + 1))) / dz(i + 1).squaredNorm()) / (ż[i + 1].norm() * dz(i + 1).norm()); + } + + dC += - γ * (z[i - 1] + z[i + 1]).conjugate(); + δz[i] = dC; + } + + for (unsigned i = 1; i < z.size() - 1; i++) { + δz[i] = δz[i].conjugate() - (δz[i].dot(z[i]) / z[i].squaredNorm()) * z[i].conjugate(); + } + + // We return a δz with average norm of one. + Real mag = 0; + for (unsigned i = 1; i < z.size() - 1; i++) { + mag += δz[i].norm(); + } + + for (unsigned i = 1; i < z.size() - 1; i++) { + δz[i] /= mag / n(); + } + + return δz; + } + + template + std::vector> generateRandomδz(Gen& r) const { + std::vector> δz(z.size()); + + complex_normal_distribution<> d(0, 1, 0); + for (unsigned i = 1; i < z.size() - 1; i++) { + δz[i] = randomVector(z[0].size(), d, r); + } + + return δz; + } + + template + Real perturb(const Tensor& J, Real δ₀, const std::vector>& δz, Real γ = 0) { Rope rNew = *this; - double δ = δ0; + Real δ = δ₀; // We rescale the step size by the distance between points. - double Δl = length() / (n() + 1); + Real Δl = length() / (n() + 1); while (true) { for (unsigned i = 1; i < z.size() - 1; i++) { - rNew.z[i] = normalize(z[i] - (δ * Δl) * Δz[i]); + rNew.z[i] = normalize(z[i] - (δ * Δl) * δz[i]); } - if (rNew.cost(J) < cost(J)) { + rNew.spread(); + + if (rNew.cost(J, γ) < cost(J, γ)) { break; } else { δ /= 2; @@ -180,15 +249,15 @@ class Rope { } void spread() { - double l = length(); + Real l = length(); - double a = 0; + Real a = 0; unsigned pos = 0; std::vector> zNew = z; for (unsigned i = 1; i < z.size() - 1; i++) { - double b = i * l / (z.size() - 1); + Real b = i * l / (z.size() - 1); while (b > a) { pos++; @@ -204,27 +273,44 @@ class Rope { } template - double step(const Tensor& J, double δ0) { - double δ = perturb(J, δ0); - spread(); - - return δ; + void relaxGradient(const Tensor& J, unsigned N, Real δ0) { + Real δ = δ0; + try { + for (unsigned i = 0; i < N; i++) { + std::vector> δz = generateGradientδz(J); + δ = 1.1 * perturb(J, δ, δz); + } + } catch (std::exception& e) { + } } template - void relax(const Tensor& J, unsigned N, double δ0) { - double δ = δ0; + void relaxDiscreteGradient(const Tensor& J, unsigned N, Real δ0, Real γ) { + Real δ = δ0; try { for (unsigned i = 0; i < N; i++) { - δ = 1.1 * step(J, δ); + std::vector> δz = generateDiscreteGradientδz(J, γ); + δ = 1.1 * perturb(J, δ, δz, γ); } } catch (std::exception& e) { } } + template + void relaxRandom(const Tensor& J, unsigned N, Real δ0, Gen& r) { + Real δ = δ0; + for (unsigned i = 0; i < N; i++) { + try { + std::vector> δz = generateRandomδz(r); + δ = 1.1 * perturb(J, δ, δz); + } catch (std::exception& e) { + } + } + } + template - double cost(const Tensor& J) const { - double c = 0; + Real cost(const Tensor& J, Real γ = 0) const { + Real c = 0; for (unsigned i = 1; i < z.size() - 1; i++) { Vector dH; @@ -235,6 +321,10 @@ class Rope { c += 1.0 - real(zD.dot(dz(i))) / zD.norm() / dz(i).norm(); } + for (unsigned i = 0; i < z.size() - 1; i++) { + c += γ * (z[i + 1] - z[i]).squaredNorm(); + } + return c; } diff --git a/stokes_test.cpp b/stokes_test.cpp index 7d96c8c..a9347a7 100644 --- a/stokes_test.cpp +++ b/stokes_test.cpp @@ -4,8 +4,8 @@ using namespace std::complex_literals; int main() { - Vector> z(2), dz(2), ddz(2), dH(2); - Matrix> ddH(2, 2); + Vector> z(2), dz(2), ddz(2), dH(2); + Matrix> ddH(2, 2); z << -0.75067 - 0.190132 * 1i, -0.625994 + 0.665987 * 1i; dz << -0.0636149 - 0.469166 * 1i, 0.820037 - 0.449064 * 1i; diff --git a/types.hpp b/types.hpp new file mode 100644 index 0000000..86ccb8f --- /dev/null +++ b/types.hpp @@ -0,0 +1,3 @@ +#pragma once + +using Real = double; -- cgit v1.2.3-70-g09d2