1 files changed, 59 insertions, 68 deletions
diff --git a/langevin.cpp b/langevin.cpp
index c3f8436..51cb2a0 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -35,6 +35,55 @@ std::list<std::array<ComplexVector, 2>> saddlesCloserThan(const std::unordered_m
   return pairs;
 }
 
+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);
+
+  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;
+
+  std::function<void(ComplexTensor&, std::array<unsigned, p>)> perturbJ =
+    [&γ, &dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) {
+      Complex Ji = getJ<Complex, p>(JJ, ind);
+      setJ<Complex, p>(JJ, ind, Ji + γ * dJ(r.engine()));
+    };
+
+  while (true) {
+    ComplexTensor Jp = J;
+    iterateOver<Complex, p>(Jp, perturbJ);
+
+    try {
+      z1 = findSaddle(Jp, z10, Δ);
+      z2 = findSaddle(Jp, z20, Δ);
+
+      double dist = (z1 - z2).norm();
+
+      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) {
+        J = Jp;
+        prevdist = abs(imag(H1 - H2));
+        γ = 0.1 * prevdist;
+        std::cerr << prevdist << std::endl;
+
+        if (abs(imag(H1 - H2)) < ε && dist > 1e-2) {
+          break;
+        }
+      }
+    } catch (std::exception& e) {}
+  }
+
+  return {J, z1, z2};
+}
+
+
 int main(int argc, char* argv[]) {
   // model parameters
   unsigned N = 10; // number of spins
@@ -43,7 +92,7 @@ int main(int argc, char* argv[]) {
   double Iκ = 0;   // imaginary part of distribution parameter
 
   // simulation parameters
-  double ε = 1e-4;
+  double ε = 1e-12;
   double εJ = 1e-5;
   double δ = 1e-2;   // threshold for determining saddle
   double Δ = 1e-3;
@@ -133,80 +182,22 @@ int main(int argc, char* argv[]) {
 
   std::cerr << "Found sufficiently nearby saddles, perturbing J to equalize Im H." << std::endl;
 
-  complex_normal_distribution<> dJ(0, εJ * σ, 0);
-
-  std::function<void(ComplexTensor&, std::array<unsigned, p>)> perturbJ =
-    [&εJ, &dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) {
-      Complex Ji = getJ<Complex, p>(JJ, ind);
-      setJ<Complex, p>(JJ, ind, Ji + εJ * dJ(r.engine()));
-    };
-
-  ComplexTensor Jp = J;
-  Complex H1, H2;
-  ComplexVector z1, z2;
-  std::tie(H1, std::ignore, std::ignore) = hamGradHess(Jp, nearbySaddles.front()[0]);
-  std::tie(H2, std::ignore, std::ignore) = hamGradHess(Jp, nearbySaddles.front()[1]);
-  double prevdist = abs(imag(H1-H2));
-  εJ = 1e4 * prevdist;
-
-  while (true) {
-    ComplexTensor Jpp = Jp;
-    iterateOver<Complex, p>(Jpp, perturbJ);
+  ComplexVector z1 = nearbySaddles.front()[0];
+  ComplexVector z2 = nearbySaddles.front()[1];
 
-    try {
-      z1 = findSaddle(Jpp, nearbySaddles.front()[0], ε);
-      z2 = findSaddle(Jpp, 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);
 
-      double dist = (z1 - z2).norm();
-
-      std::tie(H1, std::ignore, std::ignore) = hamGradHess(Jpp, z1);
-      std::tie(H2, std::ignore, std::ignore) = hamGradHess(Jpp, z2);
-
-      if (abs(imag(H1 - H2)) < prevdist && dist > 1e-2) {
-        Jp = Jpp;
-        prevdist = abs(imag(H1 - H2));
-        εJ = 1e4 * prevdist;
-
-        if (abs(imag(H1 - H2)) < 1e-13 && dist > 1e-2) {
-          std::cerr << "Found distinct critical points with sufficiently similar Im H." << std::endl;
-          break;
-        }
-      }
+    Rope<Complex> stokes(10, z1, z2, J);
 
+    for (unsigned i = 0; i < 8; i++) {
+      stokes.relax(J, 1e5, 0.1);
 
+      std::cout << stokes.n() << " " << stokes.cost(J) << " " << stokes.length() << " " << stokes.error(J) << std::endl;
 
-    } catch (std::exception& e) {
-      std::cerr << "Couldn't find a saddle with new couplings, skipping." << std::endl;
+      stokes = stokes.interpolate();
     }
   }
 
-  Rope<Complex> stokes(10, z1, z2, Jp);
-
-  std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl;
-
-  stokes.relax(Jp, 1e5, 1e-1);
-
-  std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl;
-
-  stokes = stokes.interpolate();
-  stokes.relax(Jp, 1e5, 1e-1);
-
-  std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl;
-
-  stokes = stokes.interpolate();
-  stokes.relax(Jp, 1e5, 1e-1);
-
-  std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl;
-
-  stokes = stokes.interpolate();
-  stokes.relax(Jp, 1e5, 1e-1);
-
-  std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl;
-
-  stokes = stokes.interpolate();
-  stokes.relax(Jp, 1e5, 1e-1);
-
-  std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl;
-
   return 0;
 }