diff options
| -rw-r--r-- | langevin.cpp | 137 | ||||
| -rw-r--r-- | 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<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(c 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); - 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<void(ComplexTensor&, std::array<unsigned, p>)> perturbJ = [&γ, &dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) { @@ -54,7 +51,7 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(c setJ<Complex, p>(JJ, ind, Ji + γ * dJ(r.engine())); }; - while (true) { + while (ΔH > ε) { ComplexTensor Jp = J; iterateOver<Complex, p>(Jp, perturbJ); @@ -62,18 +59,17 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> 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<ComplexTensor, ComplexVector, ComplexVector> 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<Real> d(0, 1, 0); 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<Real(const ComplexTensor&, const ComplexVector&)> energyNormGrad = [] (const ComplexTensor& J, const ComplexVector& z) { @@ -155,45 +149,104 @@ int main(int argc, char* argv[]) { return W; }; - std::unordered_map<uint64_t, ComplexVector> saddles; - std::list<std::array<ComplexVector, 2>> 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<Real(const ComplexTensor&, const ComplexVector&)> 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<Complex>(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<Real>::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(); } @@ -26,6 +26,45 @@ std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>> hamGradHess(const Tensor<Scal return {hamiltonian, gradient, hessian}; } +template <class Scalar, int p> +Scalar getHamiltonian(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) { + Scalar H; + std::tie(H, std::ignore, std::ignore) = hamGradHess(J, z); + return H; +} + +template <class Scalar, int p> +Vector<Scalar> getGradient(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) { + Vector<Scalar> dH; + std::tie(std::ignore, dH, std::ignore) = hamGradHess(J, z); + return dH; +} + +template <class Scalar, int p> +Matrix<Scalar> getHessian(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) { + Matrix<Scalar> ddH; + std::tie(std::ignore, std::ignore, ddH) = hamGradHess(J, z); + return ddH; +} + +template <class Scalar> +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 <class Scalar, int p> +Real getProportionOfThreshold(Scalar κ, const Tensor<Scalar, p>& J, const Vector<Scalar>& z) { + Real N = z.size(); + Scalar ε = getHamiltonian(J, z) / N; + Real a = z.squaredNorm() / N; + + return norm(ε) / getThresholdEnergyDensity(p, κ, ε, a); +} + template <class Scalar> Vector<Scalar> zDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) { return -dH.conjugate() + (dH.dot(z) / z.squaredNorm()) * z.conjugate(); |