From fb6b828d755dffd73b9168158b1d5c4d7d0a9923 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 1 Mar 2021 16:52:07 +0100 Subject: Work to find beyond-threshold points. --- langevin.cpp | 137 +++++++++++++++++++++++++++++++++++++++++------------------ p-spin.hpp | 39 +++++++++++++++++ 2 files changed, 134 insertions(+), 42 deletions(-) diff --git a/langevin.cpp b/langevin.cpp index e5b370c..dc7c5cf 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -41,12 +41,9 @@ std::tuple matchImaginaryEnergies(c 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); - Real prevdist = abs(imag(H1-H2)); - Real γ = prevdist / 10; + Real ΔH = abs(imag(getHamiltonian(J, z10) - getHamiltonian(J, z20))) / z10.size(); + Real γ = ΔH; std::function)> perturbJ = [&γ, &dJ, &r] (ComplexTensor& JJ, std::array ind) { @@ -54,7 +51,7 @@ std::tuple matchImaginaryEnergies(c setJ(JJ, ind, Ji + γ * dJ(r.engine())); }; - while (true) { + while (ΔH > ε) { ComplexTensor Jp = J; iterateOver(Jp, perturbJ); @@ -62,18 +59,17 @@ std::tuple matchImaginaryEnergies(c z1 = findSaddle(Jp, z10, Δ); z2 = findSaddle(Jp, z20, Δ); - Real dist = (z1 - z2).norm(); + Real Δz = (z1 - z2).norm(); - std::tie(H1, std::ignore, std::ignore) = hamGradHess(Jp, z1); - std::tie(H2, std::ignore, std::ignore) = hamGradHess(Jp, z2); + if (Δz > 1e-2) { + Real ΔHNew = abs(imag(getHamiltonian(Jp, z1) - getHamiltonian(Jp, z2))) / z1.size(); - if (abs(imag(H1 - H2)) < prevdist && dist > Real(1) / 100) { - J = Jp; - prevdist = abs(imag(H1 - H2)); - γ = prevdist / 10; + if (ΔHNew < ΔH) { + J = Jp; + ΔH = ΔHNew; + γ = ΔH; - if (abs(imag(H1 - H2)) < ε && dist > Real(1) / 100) { - break; + std::cerr << "Match imaginary energies: Found couplings with ΔH = " << ΔH << std::endl; } } } catch (std::exception& e) {} @@ -82,20 +78,19 @@ std::tuple matchImaginaryEnergies(c return {J, z1, z2}; } - int main(int argc, char* argv[]) { // model parameters unsigned N = 10; // number of spins Real T = 1; // temperature - Real Rκ = 1; // real part of distribution parameter + Real Rκ = 0; // real part of distribution parameter Real Iκ = 0; // imaginary part of distribution parameter // simulation parameters - Real ε = 1e-12; + Real ε = 1e-15; Real εJ = 1e-5; Real δ = 1e-2; // threshold for determining saddle Real Δ = 1e-3; - Real γ = 1e-2; // step size + Real γ = 1e-1; // step size unsigned t = 1000; // number of Langevin steps unsigned M = 100; unsigned n = 100; @@ -146,7 +141,6 @@ int main(int argc, char* argv[]) { complex_normal_distribution d(0, 1, 0); ComplexTensor J = generateCouplings(N, complex_normal_distribution(0, σ, κ), r.engine()); - ComplexVector z = normalize(randomVector(N, d, r.engine())); std::function energyNormGrad = [] (const ComplexTensor& J, const ComplexVector& z) { @@ -155,45 +149,104 @@ int main(int argc, char* argv[]) { return W; }; - std::unordered_map saddles; - std::list> nearbySaddles; - while (true) { // Until we find two saddles sufficiently close... - std::tie(std::ignore, z) = metropolis(J, z, energyNormGrad, T, γ, M, d, r.engine()); - try { - ComplexVector zSaddleNext = findSaddle(J, z, ε); - uint64_t saddleKey = round(1e2 * real(zSaddleNext(0))); - auto got = saddles.find(saddleKey); - - if (got == saddles.end()) { - saddles[saddleKey] = zSaddleNext; - nearbySaddles = saddlesCloserThan(saddles, δ); - if (nearbySaddles.size() > 0) { - break; - } - std::cerr << "Found " << saddles.size() << " distinct saddles." << std::endl; + std::function energyThres = [N, κ] + (const ComplexTensor& J, const ComplexVector& z) { + Real a = z.squaredNorm() / (Real)N; + /* + Complex ε = getHamiltonian(J, z) / (Real)N; + Real t = norm(ε) / getThresholdEnergyDensity(p, κ, ε, a); + */ + + Real amin = 1.1547; + Real amax = 1.46806; + + /* + Real tmin = 1; + Real tmax = 1.17514; + */ + + Real E = 0; + + if (a > amax) { + E += pow(a - amax, 2); + } else if (a < amin) { + E += pow(a - amin, 2); } - } catch (std::exception& e) { -// std::cerr << "Could not find a saddle: " << e.what() << std::endl; + + /* + if (t > tmax) { + E += pow(t - tmax, 2); + } else if (t < tmin) { + E += pow(t - tmin, 2); + } + */ + + return E; + }; + + Real largestThreshold = 0; + ComplexVector saddlePastThreshold; + bool foundSaddle = false; + +#pragma omp parallel default(none) shared(largestThreshold, saddlePastThreshold, foundSaddle, M, J, energyThres, d, N, γ, ε, κ) + { + Rng rPrivate; + ComplexVector z = normalize(randomVector(N, d, rPrivate.engine())); + + while (largestThreshold < 1) { // Until we find a saddle past the threshold... + std::tie(std::ignore, z) = metropolis(J, z, energyThres, (Real)0.1, γ, M, d, rPrivate.engine()); + try { + ComplexVector latestSaddle = findSaddle(J, z, ε); + Real threshold = getProportionOfThreshold(κ, J, latestSaddle); + +#pragma omp critical + { + if (threshold > largestThreshold + 1e-6) { + largestThreshold = threshold; + saddlePastThreshold = latestSaddle; + std::cerr << "Found saddle with threshold porportion " << largestThreshold << std::endl;; + } + } + } catch (std::exception& e) {} } + } + + std::cerr << "Found saddle with energy beyond threshold." << std::endl; + + ComplexVector zSaddleNext; + Real closestSaddle = std::numeric_limits::infinity(); + ComplexVector z = saddlePastThreshold; + while (closestSaddle > δ) { // Until we find two saddles sufficiently close... + std::tie(std::ignore, z) = metropolis(J, saddlePastThreshold, energyNormGrad, T, γ, 1e4, d, r.engine()); + try { + zSaddleNext = findSaddle(J, z, ε); + Real saddleDistance = (zSaddleNext - saddlePastThreshold).norm(); + if (saddleDistance < closestSaddle && saddleDistance > 1e-2) { + closestSaddle = saddleDistance; + std::cerr << "Nearby saddle: found saddle at distance " << saddleDistance << std::endl; + } + } catch (std::exception& e) {} } std::cerr << "Found sufficiently nearby saddles, perturbing J to equalize Im H." << std::endl; - ComplexVector z1 = nearbySaddles.front()[0]; - ComplexVector z2 = nearbySaddles.front()[1]; + ComplexVector z1 = saddlePastThreshold; + ComplexVector z2 = zSaddleNext; - std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, 1e-15, ε, r); + std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, 1e-14, ε, r); std::cerr << "Im H is now sufficently close, starting to relax rope." << std::endl; + std::cerr << "Threshold proportions of saddles are " << getProportionOfThreshold(κ, J, z1) << " and " << getProportionOfThreshold(κ, J, z2) << std::endl; + Rope stokes(10, z1, z2, J); for (unsigned i = 0; i < 9; i++) { stokes.relaxDiscreteGradient(J, 1e6, 0.1, 0); - std::cout << stokes.n() << " " << stokes.cost(J) << std::endl; + std::cout << stokes.n() << " " << stokes.cost(J) << " " << stokes.length() << std::endl; stokes = stokes.interpolate(); } diff --git a/p-spin.hpp b/p-spin.hpp index f1dc07f..b2bdf07 100644 --- a/p-spin.hpp +++ b/p-spin.hpp @@ -26,6 +26,45 @@ std::tuple, Matrix> hamGradHess(const Tensor +Scalar getHamiltonian(const Tensor& J, const Vector& z) { + Scalar H; + std::tie(H, std::ignore, std::ignore) = hamGradHess(J, z); + return H; +} + +template +Vector getGradient(const Tensor& J, const Vector& z) { + Vector dH; + std::tie(std::ignore, dH, std::ignore) = hamGradHess(J, z); + return dH; +} + +template +Matrix getHessian(const Tensor& J, const Vector& z) { + Matrix ddH; + std::tie(std::ignore, std::ignore, ddH) = hamGradHess(J, z); + return ddH; +} + +template +Real getThresholdEnergyDensity(unsigned p, Scalar κ, Scalar ε, Real a) { + Real apm2 = pow(a, p - 2); + Scalar δ = κ / apm2; + Real θ = arg(κ) + 2 * arg(ε); + + return (p - 1) * apm2 / (2 * p) * pow(1 - norm(δ), 2) / (1 + norm(δ) - 2 * abs(δ) * cos(θ)); +} + +template +Real getProportionOfThreshold(Scalar κ, const Tensor& J, const Vector& z) { + Real N = z.size(); + Scalar ε = getHamiltonian(J, z) / N; + Real a = z.squaredNorm() / N; + + return norm(ε) / getThresholdEnergyDensity(p, κ, ε, a); +} + template Vector zDot(const Vector& z, const Vector& dH) { return -dH.conjugate() + (dH.dot(z) / z.squaredNorm()) * z.conjugate(); -- cgit v1.2.3-70-g09d2