From 9c3ac9c97abeca3ebcb11a1a2c8a6e3cb3791735 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 16 Feb 2021 16:42:53 +0100
Subject: Progress towards a working Stokes line finder.

---
 langevin.cpp | 92 ++++++++++++++++++++++++++++++++++++++++--------------------
 1 file changed, 62 insertions(+), 30 deletions(-)

(limited to 'langevin.cpp')

diff --git a/langevin.cpp b/langevin.cpp
index 5ca2d72..5d9acc9 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -1,8 +1,10 @@
 #include <getopt.h>
+#include <stereographic.hpp>
 
 #include "complex_normal.hpp"
 #include "p-spin.hpp"
 #include "dynamics.hpp"
+#include "stokes.hpp"
 
 #include "pcg-cpp/include/pcg_random.hpp"
 #include "randutils/randutils.hpp"
@@ -26,7 +28,7 @@ int main(int argc, char* argv[]) {
 
   // simulation parameters
   double ε = 1e-4;
-  double εJ = 1e-2;
+  double εJ = 1e-5;
   double δ = 1e-2;   // threshold for determining saddle
   double Δ = 1e-3;
   double γ = 1e-2;   // step size
@@ -92,28 +94,6 @@ int main(int argc, char* argv[]) {
       return W;
     };
 
-  double aGoal = 1e3;
-
-  std::function<double(const ComplexTensor&, const ComplexVector&)> energyInvA = [aGoal]
-    (const ComplexTensor& J, const ComplexVector& z) {
-      double a = z.squaredNorm();
-      if (a > aGoal) {
-        return -aGoal;
-      } else {
-        return -z.squaredNorm();
-      }
-    };
-
-  while (zSaddle.squaredNorm() < aGoal) {
-    std::tie(std::ignore, z) = metropolis(J, z, energyInvA, T, γ, 100, d, r.engine());
-    try {
-      std::cerr << "Starting descent from " << z.squaredNorm() << "." << std::endl;
-      zSaddle = findSaddle(J, z, ε);
-    } catch (std::exception& e) {
-    }
-    std::cerr << "Current saddle is of size " << zSaddle.squaredNorm() << "." << std::endl;
-  }
-
   ComplexVector zSaddlePrev = ComplexVector::Zero(N);
 
   while (δ < (zSaddle - zSaddlePrev).norm()) { // Until we find two saddles sufficiently close...
@@ -131,7 +111,11 @@ int main(int argc, char* argv[]) {
     std::cerr << "Current saddles are " << (zSaddle - zSaddlePrev).norm() << " apart." << std::endl;
   }
 
-  std::cerr << "Found sufficiently nearby saddles, perturbing J." << std::endl;
+  Rope<Complex> stokest(20, zSaddle, zSaddlePrev);
+
+  stokest.relax(J, 10000, 0.0001);
+
+  std::cerr << "Found sufficiently nearby saddles, perturbing J to equalize Im H." << std::endl;
 
   complex_normal_distribution<> dJ(0, εJ * σ, 0);
 
@@ -141,20 +125,68 @@ int main(int argc, char* argv[]) {
       setJ<Complex, p>(JJ, ind, Ji + dJ(r.engine()));
     };
 
-  for (unsigned i = 0; i < n; i++) {
-    ComplexTensor Jp = J;
+  ComplexTensor Jp = J;
+  ComplexVector z1, z2;
+  Complex H1, H2;
+  double prevdist = 100;
 
-    iterateOver<Complex, p>(Jp, perturbJ);
+  while (true) {
+    ComplexTensor Jpp = Jp;
+    iterateOver<Complex, p>(Jpp, perturbJ);
 
     try {
-      ComplexVector zSaddleNew = findSaddle(Jp, zSaddle, ε);
-      ComplexVector zSaddlePrevNew = findSaddle(Jp, zSaddlePrev, ε);
+      z1 = findSaddle(Jpp, zSaddle, ε);
+      z2 = findSaddle(Jpp, zSaddlePrev, ε);
+
+      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));
+
+        std::cout << abs(imag(H1 - H2)) << " " << dist << std::endl;
+        if (abs(imag(H1 - H2)) < 1e-8 && dist > 1e-2) {
+          std::cout << "Found distinct critical points with sufficiently similar Im H." << std::endl;
+          break;
+        }
+      }
+
+
 
-      std::cout << zSaddleNew.transpose() << " " << zSaddlePrevNew.transpose() << std::endl;
     } catch (std::exception& e) {
       std::cerr << "Couldn't find a saddle with new couplings, skipping." << std::endl;
     }
   }
 
+  Rope<Complex> stokes(20, z1, z2);
+
+  std::cout << stokes.cost(Jp) << std::endl;
+
+  stokes.relax(Jp, 10000, 0.0001);
+
+  std::cout << stokes.cost(Jp) << std::endl;
+
+  Rope<Complex> stokes2 = stokes.interpolate();
+
+  stokes2.relax(Jp, 10000, 0.001);
+
+  std::cout << stokes2.cost(Jp) << std::endl;
+
+  stokes2 = stokes2.interpolate();
+
+  stokes2.relax(Jp, 10000, 0.001);
+
+  std::cout << stokes2.cost(Jp) << std::endl;
+
+  stokes2 = stokes2.interpolate();
+
+  stokes2.relax(Jp, 10000, 0.001);
+
+  std::cout << stokes2.cost(Jp) << std::endl;
+
+
   return 0;
 }
-- 
cgit v1.2.3-70-g09d2