Work on Stokes lines.

2 files changed, 74 insertions, 73 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;
 }
diff --git a/stokes.hpp b/stokes.hpp
index 3814cf4..6cc5c4a 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -1,4 +1,5 @@
 #include "p-spin.hpp"
+#include <iostream>
 
 class ropeRelaxationStallException: public std::exception {
   virtual const char* what() const throw() {
@@ -130,6 +131,11 @@ class Rope {
         δz[i] = variation(z[i], dz(i), ddz(i), dH, ddH);
       }
 
+      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();
+//        δz[i] = δz[i] - ((δz[i].conjugate().dot(dz(i))) / dz(i).squaredNorm()) * dz(i).conjugate();
+      }
+
       // We return a δz with average norm of one.
       double mag = 0;
       for (unsigned i = 1; i < z.size() - 1; i++) {
@@ -152,14 +158,18 @@ class Rope {
       // We rescale the step size by the distance between points.
       double Δl = length() / (n() + 1);
 
-      while (rNew.cost(J) >= cost(J)) {
+      while (true) {
         for (unsigned i = 1; i < z.size() - 1; i++) {
-          rNew.z[i] = normalize(z[i] - (δ * Δl) * Δz[i].conjugate());
+          rNew.z[i] = normalize(z[i] - (δ * Δl) * Δz[i]);
         }
 
-        δ /= 2;
+        if (rNew.cost(J) < cost(J)) {
+          break;
+        } else {
+          δ /= 2;
+        }
 
-        if (δ < 1e-30) {
+        if (δ < 1e-15) {
           throw ropeRelaxationStallException();
         }
       }
@@ -206,7 +216,7 @@ class Rope {
       double δ = δ0;
       try {
         for (unsigned i = 0; i < N; i++) {
-          δ = 2 * step(J, δ);
+          δ = 1.1 * step(J, δ);
         }
       } catch (std::exception& e) {
       }