diff options
| -rw-r--r-- | complex_normal.hpp | 2 | ||||
| -rw-r--r-- | langevin.cpp | 17 | ||||
| -rw-r--r-- | stereographic.hpp | 10 | ||||
| -rw-r--r-- | stokes.hpp | 14 |
4 files changed, 18 insertions, 25 deletions
diff --git a/complex_normal.hpp b/complex_normal.hpp index 2ffcea1..65b054f 100644 --- a/complex_normal.hpp +++ b/complex_normal.hpp @@ -15,7 +15,7 @@ public: complex_normal_distribution(std::complex<T> μ, T σ, std::complex<T> κ) : μ(μ), d(0, σ / sqrt(2)) { δx = sqrt(1 + std::abs(κ)); δy = sqrt(1 - std::abs(κ)); - eθ = std::polar(1.0, std::arg(κ)); + eθ = std::polar((T)1, std::arg(κ)); } template <class Generator> std::complex<T> operator()(Generator& g) { diff --git a/langevin.cpp b/langevin.cpp index 8c191f3..e5b370c 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -37,7 +37,7 @@ 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, Real ε, Real Δ, Generator& r) { - Real σ = sqrt(factorial(p) / (2.0 * pow(z10.size(), p - 1))); + Real σ = sqrt(factorial(p) / (Real(2) * pow(z10.size(), p - 1))); complex_normal_distribution<Real> dJ(0, σ, 0); ComplexTensor J = J0; @@ -46,7 +46,7 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(c std::tie(H1, std::ignore, std::ignore) = hamGradHess(J, z10); std::tie(H2, std::ignore, std::ignore) = hamGradHess(J, z20); Real prevdist = abs(imag(H1-H2)); - Real γ = 0.1 * prevdist; + Real γ = prevdist / 10; std::function<void(ComplexTensor&, std::array<unsigned, p>)> perturbJ = [&γ, &dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) { @@ -67,13 +67,12 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(c std::tie(H1, std::ignore, std::ignore) = hamGradHess(Jp, z1); std::tie(H2, std::ignore, std::ignore) = hamGradHess(Jp, z2); - if (abs(imag(H1 - H2)) < prevdist && dist > 1e-2) { + if (abs(imag(H1 - H2)) < prevdist && dist > Real(1) / 100) { J = Jp; prevdist = abs(imag(H1 - H2)); - γ = 0.1 * prevdist; - std::cerr << prevdist << std::endl; + γ = prevdist / 10; - if (abs(imag(H1 - H2)) < ε && dist > 1e-2) { + if (abs(imag(H1 - H2)) < ε && dist > Real(1) / 100) { break; } } @@ -140,7 +139,7 @@ int main(int argc, char* argv[]) { } Complex κ(Rκ, Iκ); - Real σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1))); + Real σ = sqrt(factorial(p) / ((Real)2 * pow(N, p - 1))); Rng r; @@ -185,7 +184,9 @@ int main(int argc, char* argv[]) { ComplexVector z1 = nearbySaddles.front()[0]; ComplexVector z2 = nearbySaddles.front()[1]; - std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, 1e-14, ε, r); + std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, 1e-15, ε, r); + + std::cerr << "Im H is now sufficently close, starting to relax rope." << std::endl; Rope stokes(10, z1, z2, J); diff --git a/stereographic.hpp b/stereographic.hpp index 2dfa735..58c6b31 100644 --- a/stereographic.hpp +++ b/stereographic.hpp @@ -10,13 +10,13 @@ Vector<Scalar> stereographicToEuclidean(const Vector<Scalar>& ζ) { Vector<Scalar> z(N); Scalar a = ζ.transpose() * ζ; - Scalar b = 2 * sqrt(N) / (1.0 + a); + Scalar b = 2 * sqrt((Real)N) / ((Real)1 + a); for (unsigned i = 0; i < N - 1; i++) { z(i) = ζ(i); } - z(N - 1) = (a - 1.0) / 2.0; + z(N - 1) = (a - (Real)1) / (Real)2; return b * z; } @@ -30,7 +30,7 @@ Vector<Scalar> euclideanToStereographic(const Vector<Scalar>& z) { ζ(i) = z(i); } - return ζ / (sqrt(N) - z(N - 1)); + return ζ / (sqrt((Real)N) - z(N - 1)); } template <class Scalar> @@ -38,14 +38,14 @@ Matrix<Scalar> stereographicJacobian(const Vector<Scalar>& ζ) { unsigned N = ζ.size(); Matrix<Scalar> J(N, N + 1); - Scalar b = 1.0 + (Scalar)(ζ.transpose() * ζ); + Scalar b = (Real)1 + (Scalar)(ζ.transpose() * ζ); for (unsigned i = 0; i < N; i++) { for (unsigned j = 0; j < N; j++) { J(i, j) = - ζ(i) * ζ(j); if (i == j) { - J(i, j) += 0.5 * b; + J(i, j) += b / (Real)2; } } @@ -1,7 +1,8 @@ +#pragma once + #include "p-spin.hpp" #include "complex_normal.hpp" #include "dynamics.hpp" -#include <iostream> class ropeRelaxationStallException: public std::exception { virtual const char* what() const throw() { @@ -189,20 +190,11 @@ class Rope { δz[i] = dC; } + Real size = 0; 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; } |