Progress towards a working Stokes line finder.

1 files changed, 92 insertions, 37 deletions

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];