Work to find beyond-threshold points.

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();
   }
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<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();