Progress towards a working Stokes line finder.

2 files changed, 154 insertions, 67 deletions

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;
 }
diff --git a/stokes.hpp b/stokes.hpp
index b60ae0e..43a69e3 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -1,47 +1,96 @@
-#include "stereographic.hpp"
+#include "p-spin.hpp"
+#include <iostream>
+#include "dynamics.hpp"
 
 template <class Scalar>
-class Rope {
-  private:
-    std::vector<Vector<Scalar>> z;
+Vector<Scalar> zDot(const Vector<Scalar>& z, const Vector<Scalar>& dH) {
+  return -dH.conjugate() + (dH.dot(z) / z.squaredNorm()) * z.conjugate();
+}
 
-    template <int p>
-    Vector<Scalar> δz(const Tensor<Scalar, p>& J) const {
-      Vector<Scalar> dz(z.size());
-
-      for (unsigned i = 1; i < z.size() - 1; i++) {
-        Vector<Scalar> z12 = (z[i] + z[i - 1]) / 2.0;
-        Vector<Scalar> g;
-        std::tie(std::ignore, g, std::ignore) = stereographicHamGradHessHess(J, z12);
+template <class Scalar>
+double segmentCost(const Vector<Scalar>& z, const Vector<Scalar>& dz, const Vector<Scalar>& dH) {
+  Vector<Scalar> zD = zDot(z, dH);
+  return 1.0 - pow(real(zD.dot(dz)), 2) / zD.squaredNorm() / dz.squaredNorm();
+}
 
-        Vector<Scalar> Δ = z[i] - z[i - 1];
+template <class Scalar>
+Vector<Scalar> variation(const Vector<Scalar>& z, const Vector<Scalar>& z´, const Vector<Scalar>& z´´, const Vector<Scalar>& dH, const Matrix<Scalar>& ddH) {
+  double z² = z.squaredNorm();
+  double z´² = z´.squaredNorm();
+
+  Vector<Scalar> ż = zDot(z, dH);
+  double ż² = ż.squaredNorm();
+
+  double Reż·z´² = real(pow(ż.dot(z´), 2));
+
+  Matrix<Scalar> dż = (dH.conjugate() - (dH.dot(z) / z²) * z.conjugate()) * z.adjoint() / z²;
+  Matrix<Scalar> dżc = -ddH + (ddH * z.conjugate()) * z.transpose() / z²
+    + (z.dot(dH) / z²) * (Matrix<Scalar>::Identity(ddH.rows(), ddH.cols()) - z.conjugate() * z.transpose() / z²);
+
+  Vector<Scalar> dLdz = - 0.5 * (
+        ż.dot(z´) * (dżc * z´) + z´.dot(ż) * (dż * z´.conjugate())
+        - (dż * ż.conjugate() + dżc * ż) * Reż·z´² / ż²
+      ) / (ż² * z´²);
+
+  Vector<Scalar> ż´ = -(ddH * z´).conjugate() + ((ddH * z´).dot(z) / z²) * z.conjugate() + (
+        dH.dot(z) * z´.conjugate() + dH.dot(z´) * z.conjugate() - (dH.dot(z) * (z´.dot(z) + z.dot(z´)) / z²) * z.conjugate()
+      ) / z²;
+
+  Vector<Scalar> ddtdLdz´ = -0.5 * (
+      ż´.dot(z´) * ż.conjugate() + ż.dot(z´´) * ż.conjugate() + ż.dot(z´) * ż´.conjugate()
+      - (
+        Reż·z´² * z´´.conjugate() + real(2.0 * ż.dot(z´) * (ż.dot(z´´) + ż´.dot(z´))) * z´.conjugate()
+        - Reż·z´² * (z´´.dot(z´) + z´.dot(z´´)) * z´.conjugate() / z´²
+        ) / z´²
+      - (
+        (ż.dot(ż´) + ż´.dot(ż)) / ż² + (z´´.dot(z´) + z´.dot(z´´)) / z´²
+        ) * (
+        ż.dot(z´) * ż.conjugate() - Reż·z´² / z´² * z´.conjugate()
+        )
+      ) / (ż² * z´²);
+
+  return dLdz - ddtdLdz´;
+}
 
-        g.normalize();
-        Δ.normalize();
+template <class Scalar>
+class Rope {
+  public:
+    std::vector<Vector<Scalar>> z;
 
-        if (abs(arg(g.dot(Δ))) < M_PI / 2) {
-          dz[i] = g - Δ;
-        } else {
-          dz[i] = -g - Δ;
-        }
+    Rope(unsigned N, const Vector<Scalar>& z1, const Vector<Scalar>& z2) : z(N + 2) {
+      for (unsigned i = 0; i < N + 2; i++) {
+        z[i] = z1 + (z2 - z1) * ((double)i / (N + 1.0));
+        z[i] = normalize(z[i]);
       }
+    }
+
+    Vector<Scalar> dz(unsigned i) const {
+      return z[i + 1] - z[i - 1];
+    }
 
-      return dz;
+    Vector<Scalar> ddz(unsigned i) const {
+      return 4.0 * (z[i + 1] + z[i - 1] - 2.0 * z[i]);
     }
 
     template <int p>
     void perturb(const Tensor<Scalar, p>& J, double δ) {
-      z += δ * this->δz(J);
+      for (unsigned i = 1; i < z.size() - 1; i++) {
+        Vector<Scalar> dH;
+        Matrix<Scalar> ddH;
+        std::tie(std::ignore, dH, ddH) = hamGradHess(J, z[i]);
+
+        Vector<Scalar> δz = variation(z[i], this->dz(i), this->ddz(i), dH, ddH);
+
+        z[i] -= δ * δz.conjugate();
+        z[i] = normalize(z[i]);
+      }
     }
 
     void spread() {
-      std::vector<double> d(z.size() - 1);
       double l = 0;
 
       for (unsigned i= 0; i < z.size() - 1; i++) {
-        double Δ = (z[i + 1] - z[i]).norm();
-        d[i] = Δ;
-        l += Δ;
+        l += (z[i + 1] - z[i]).norm();
       }
 
       double a = 0;
@@ -49,19 +98,16 @@ class Rope {
 
       for (unsigned i = 1; i < z.size() - 1; i++) {
         double b = i * l / (z.size() - 1);
-        while (b >= a) {
-          a += d[pos];
+
+        while (b > a) {
           pos++;
+          a += (z[pos] - z[pos - 1]).norm();
         } 
 
-        z[i] = z[pos] - (a - b) * (z[pos] - z[pos - 1]).normalized();
-      }
-    }
+        Vector<Scalar> δz = z[pos] - z[pos - 1];
 
-  public:
-    Rope(unsigned N, const Vector<Scalar>& z1, const Vector<Scalar>& z2) : z(N + 2) {
-      for (unsigned i = 0; i < N + 2; i++) {
-        z[i] = z1 + (z2 - z1) * ((double)i / (N + 1.0));
+        z[i] = z[pos] - (a - b) / δz.norm() * δz;
+        z[i] = normalize(z[i]);
       }
     }
 
@@ -73,7 +119,16 @@ class Rope {
 
     template <int p>
     double cost(const Tensor<Scalar, p>& J) const {
-      return this->δz(J).norm();
+      double c = 0;
+
+      for (unsigned i = 1; i < z.size() - 1; i++) {
+        Vector<Scalar> dH;
+        std::tie(std::ignore, dH, std::ignore) = hamGradHess(J, z[i]);
+
+        c += segmentCost(z[i], this->dz(i), dH);
+      }
+
+      return c;
     }
 
     template <int p>
@@ -84,7 +139,7 @@ class Rope {
     }
 
     Rope<Scalar> interpolate() const {
-      Rope<Scalar> r(2 * z.size() - 1);
+      Rope<Scalar> r(2 * z.size() - 1, z[0], z[z.size() - 1]);
 
       for (unsigned i = 0; i < z.size(); i++) {
         r.z[2 * i] = z[i];