diff options
| -rw-r--r-- | complex_normal.hpp | 4 | ||||
| -rw-r--r-- | dynamics.hpp | 28 | ||||
| -rw-r--r-- | langevin.cpp | 56 | ||||
| -rw-r--r-- | p-spin.hpp | 11 | ||||
| -rw-r--r-- | stokes.hpp | 154 | ||||
| -rw-r--r-- | stokes_test.cpp | 4 | ||||
| -rw-r--r-- | types.hpp | 3 |
7 files changed, 176 insertions, 84 deletions
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 T = double> class complex_normal_distribution { private: std::normal_distribution<T> d; - double δx; - double δy; + T δx; + T δy; std::complex<T> μ; std::complex<T> 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 <class Scalar, int p> -std::tuple<double, Vector<Scalar>> gradientDescent(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, double ε, double γ0 = 1, double δγ = 2) { +template <class Real, class Scalar, int p> +std::tuple<Real, Vector<Scalar>> gradientDescent(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, Real ε, Real γ0 = 1, Real δγ = 2) { Vector<Scalar> z = z0; - double γ = γ0; + Real γ = γ0; auto [W, dW] = WdW(J, z); @@ -42,12 +42,12 @@ std::tuple<double, Vector<Scalar>> gradientDescent(const Tensor<Scalar, p>& J, c return {W, z}; } -template <class Scalar, int p> -Vector<Scalar> findSaddle(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, double ε, double δW = 2, double γ0 = 1, double δγ = 2) { +template <class Real, class Scalar, int p> +Vector<Scalar> findSaddle(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, Real ε, Real δW = 2, Real γ0 = 1, Real δγ = 2) { Vector<Scalar> z = z0; Vector<Scalar> ζ = euclideanToStereographic(z); - double W; + Real W; std::tie(W, std::ignore) = WdW(J, z); Vector<Scalar> dH; @@ -61,7 +61,7 @@ Vector<Scalar> findSaddle(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, Vector<Scalar> ζNew = ζ - dζ; Vector<Scalar> 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<Scalar> randomVector(unsigned N, Distribution d, Generator& r) { return z; } -template <class Scalar, int p, class Distribution, class Generator> -std::tuple<double, Vector<Scalar>> metropolis(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, - std::function<double(const Tensor<Scalar, p>&, const Vector<Scalar>&)>& energy, - double T, double γ, unsigned N, Distribution d, Generator& r) { +template <class Real, class Scalar, int p, class Distribution, class Generator> +std::tuple<Real, Vector<Scalar>> metropolis(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, + std::function<Real(const Tensor<Scalar, p>&, const Vector<Scalar>&)>& energy, + Real T, Real γ, unsigned N, Distribution d, Generator& r) { Vector<Scalar> z = z0; - double E = energy(J, z); + Real E = energy(J, z); - std::uniform_real_distribution<double> D(0, 1); + std::uniform_real_distribution<Real> D(0, 1); for (unsigned i = 0; i < N; i++) { Vector<Scalar> zNew = normalize(z + γ * randomVector<Scalar>(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<double>; +using Complex = std::complex<Real>; using ComplexVector = Vector<Complex>; using ComplexMatrix = Matrix<Complex>; using ComplexTensor = Tensor<Complex, p>; using Rng = randutils::random_generator<pcg32>; -std::list<std::array<ComplexVector, 2>> saddlesCloserThan(const std::unordered_map<uint64_t, ComplexVector>& saddles, double δ) { +std::list<std::array<ComplexVector, 2>> saddlesCloserThan(const std::unordered_map<uint64_t, ComplexVector>& saddles, Real δ) { std::list<std::array<ComplexVector, 2>> pairs; for (auto it1 = saddles.begin(); it1 != saddles.end(); it1++) { @@ -36,17 +36,17 @@ std::list<std::array<ComplexVector, 2>> saddlesCloserThan(const std::unordered_m } template <class Generator> -std::tuple<ComplexTensor, ComplexVector, ComplexVector> 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<ComplexTensor, ComplexVector, ComplexVector> 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<Real> 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<void(ComplexTensor&, std::array<unsigned, p>)> perturbJ = [&γ, &dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) { @@ -62,7 +62,7 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> 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<ComplexTensor, ComplexVector, ComplexVector> 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<Real> d(0, 1, 0); - ComplexTensor J = generateCouplings<Complex, PSPIN_P>(N, complex_normal_distribution<>(0, σ, κ), r.engine()); + ComplexTensor J = generateCouplings<Complex, PSPIN_P>(N, complex_normal_distribution<Real>(0, σ, κ), r.engine()); ComplexVector z = normalize(randomVector<Complex>(N, d, r.engine())); - std::function<double(const ComplexTensor&, const ComplexVector&)> energyNormGrad = [] + std::function<Real(const ComplexTensor&, const ComplexVector&)> 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<Complex> 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; @@ -2,12 +2,13 @@ #include <eigen3/Eigen/Dense> +#include "types.hpp" #include "tensor.hpp" #include "factorial.hpp" template <typename Derived> Vector<typename Derived::Scalar> normalize(const Eigen::MatrixBase<Derived>& 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 <class Scalar, int p> @@ -16,7 +17,7 @@ std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>> hamGradHess(const Tensor<Scal Vector<Scalar> Jzz = Jz * z; Scalar Jzzz = Jzz.transpose() * z; - double pBang = factorial(p); + Real pBang = factorial(p); Matrix<Scalar> hessian = ((p - 1) * p / pBang) * Jz; Vector<Scalar> gradient = (p / pBang) * Jzz; @@ -31,18 +32,18 @@ Vector<Scalar> zDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) { } template <class Scalar, int p> -std::tuple<double, Vector<Scalar>> WdW(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) { +std::tuple<Real, Vector<Scalar>> WdW(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) { Vector<Scalar> dH; Matrix<Scalar> ddH; std::tie(std::ignore, dH, ddH) = hamGradHess(J, z); Vector<Scalar> 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<Scalar> dW = ddH * (dzdt - A * z.conjugate()) + 2 * (conj(A) * B * z).real() - conj(B) * dzdt - conj(A) * dH.conjugate(); @@ -1,4 +1,6 @@ #include "p-spin.hpp" +#include "complex_normal.hpp" +#include "dynamics.hpp" #include <iostream> class ropeRelaxationStallException: public std::exception { @@ -9,13 +11,13 @@ class ropeRelaxationStallException: public std::exception { template <class Scalar> Vector<Scalar> variation(const Vector<Scalar>& z, const Vector<Scalar>& z´, const Vector<Scalar>& z´´, const Vector<Scalar>& dH, const Matrix<Scalar>& ddH) { - double z² = z.squaredNorm(); - double z´² = z´.squaredNorm(); + Real z² = z.squaredNorm(); + Real z´² = z´.squaredNorm(); Vector<Scalar> ż = zDot(z, dH); - double ż² = ż.squaredNorm(); + Real ż² = ż.squaredNorm(); - double Reż·z´ = real(ż.dot(z´)); + Real Reż·z´ = real(ż.dot(z´)); Matrix<Scalar> dż = (dH.conjugate() - (dH.dot(z) / z²) * z.conjugate()) * z.adjoint() / z²; Matrix<Scalar> dżc = -ddH + (ddH * z.conjugate()) * z.transpose() / z² @@ -33,7 +35,7 @@ Vector<Scalar> variation(const Vector<Scalar>& z, const Vector<Scalar>& z´, con ) * z.conjugate() ) / z²; - double dReż·z´ = real(ż.dot(z´´) + ż´.dot(z´)); + Real dReż·z´ = real(ż.dot(z´´) + ż´.dot(z´)); Vector<Scalar> 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 <int p> - double error(const Tensor<Scalar, p>& J) const { + Real error(const Tensor<Scalar, p>& 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 <int p> - std::vector<Vector<Scalar>> δz(const Tensor<Scalar, p>& J) const { + std::vector<Vector<Scalar>> generateGradientδz(const Tensor<Scalar, p>& J) const { std::vector<Vector<Scalar>> δ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 <int p> - double perturb(const Tensor<Scalar, p>& J, double δ0) { - std::vector<Vector<Scalar>> Δz = δz(J); + std::vector<Vector<Scalar>> generateDiscreteGradientδz(const Tensor<Scalar, p>& J, Real γ) const { + std::vector<Vector<Scalar>> δz(z.size()); + + std::vector<Vector<Scalar>> ż(z.size()); + std::vector<Vector<Scalar>> dH(z.size()); + std::vector<Matrix<Scalar>> 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<Scalar> dż = (dH[i].conjugate() - (dH[i].dot(z[i]) / z²) * z[i].conjugate()) * z[i].adjoint() / z²; + Matrix<Scalar> dżc = -ddH[i] + (ddH[i] * z[i].conjugate()) * z[i].transpose() / z² + + (z[i].dot(dH[i]) / z²) * ( + Matrix<Scalar>::Identity(ddH[i].rows(), ddH[i].cols()) - z[i].conjugate() * z[i].transpose() / z² + ); + + Vector<Scalar> 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<class Gen> + std::vector<Vector<Scalar>> generateRandomδz(Gen& r) const { + std::vector<Vector<Scalar>> δz(z.size()); + + complex_normal_distribution<> d(0, 1, 0); + for (unsigned i = 1; i < z.size() - 1; i++) { + δz[i] = randomVector<Scalar>(z[0].size(), d, r); + } + + return δz; + } + + template<int p> + Real perturb(const Tensor<Scalar, p>& J, Real δ₀, const std::vector<Vector<Scalar>>& δ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<Vector<Scalar>> 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 <int p> - double step(const Tensor<Scalar, p>& J, double δ0) { - double δ = perturb(J, δ0); - spread(); - - return δ; + void relaxGradient(const Tensor<Scalar, p>& J, unsigned N, Real δ0) { + Real δ = δ0; + try { + for (unsigned i = 0; i < N; i++) { + std::vector<Vector<Scalar>> δz = generateGradientδz(J); + δ = 1.1 * perturb(J, δ, δz); + } + } catch (std::exception& e) { + } } template <int p> - void relax(const Tensor<Scalar, p>& J, unsigned N, double δ0) { - double δ = δ0; + void relaxDiscreteGradient(const Tensor<Scalar, p>& J, unsigned N, Real δ0, Real γ) { + Real δ = δ0; try { for (unsigned i = 0; i < N; i++) { - δ = 1.1 * step(J, δ); + std::vector<Vector<Scalar>> δz = generateDiscreteGradientδz(J, γ); + δ = 1.1 * perturb(J, δ, δz, γ); } } catch (std::exception& e) { } } + template <int p, class Gen> + void relaxRandom(const Tensor<Scalar, p>& J, unsigned N, Real δ0, Gen& r) { + Real δ = δ0; + for (unsigned i = 0; i < N; i++) { + try { + std::vector<Vector<Scalar>> δz = generateRandomδz(r); + δ = 1.1 * perturb(J, δ, δz); + } catch (std::exception& e) { + } + } + } + template <int p> - double cost(const Tensor<Scalar, p>& J) const { - double c = 0; + Real cost(const Tensor<Scalar, p>& J, Real γ = 0) const { + Real c = 0; for (unsigned i = 1; i < z.size() - 1; i++) { Vector<Scalar> 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<std::complex<double>> z(2), dz(2), ddz(2), dH(2); - Matrix<std::complex<double>> ddH(2, 2); + Vector<std::complex<Real>> z(2), dz(2), ddz(2), dH(2); + Matrix<std::complex<Real>> 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; |