Work on Stokes lines, and new method involving parametric fits.

1 files changed, 211 insertions, 161 deletions

diff --git a/stokes.hpp b/stokes.hpp
index cbe5437..b2694a9 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -4,6 +4,8 @@
 #include "complex_normal.hpp"
 #include "dynamics.hpp"
 
+#include <iostream>
+
 class ropeRelaxationStallException: public std::exception {
   virtual const char* what() const throw() {
     return "Gradient descent stalled.";
@@ -11,49 +13,6 @@ class ropeRelaxationStallException: public std::exception {
 };
 
 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) {
-  Real z² = z.squaredNorm();
-  Real z´² = z´.squaredNorm();
-
-  Vector<Scalar> ż = zDot(z, dH);
-  Real ż² = ż.squaredNorm();
-
-  Real Reż·z´ = real(ż.dot(z´));
-
-  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 = - (
-        dżc * z´ + dż * z´.conjugate() - (dż * ż.conjugate() + dżc * ż) * Reż·z´ / ż²
-      ) / sqrt(ż² * z´²) / 2;
-
-  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²;
-
-  Real dReż·z´ = real(ż.dot(z´´) + ż´.dot(z´));
-
-  Vector<Scalar> ddtdLdz´ = - (
-      (
-        ż´.conjugate() - (
-          Reż·z´ * z´´.conjugate() + dReż·z´ * z´.conjugate()
-          - (Reż·z´ / z´²) * (z´´.dot(z´) + z´.dot(z´´)) * z´.conjugate() 
-          ) / z´²
-        )
-      - 0.5 * (
-        (ż.dot(ż´) + ż´.dot(ż)) / ż² + (z´´.dot(z´) + z´.dot(z´´)) / z´²
-        ) * (ż.conjugate() - Reż·z´ / z´² * z´.conjugate())
-      ) / sqrt(ż² * z´²) / 2;
-
-  return dLdz - ddtdLdz´;
-}
-
-template <class Scalar>
 class Rope {
   public:
     std::vector<Vector<Scalar>> z;
@@ -93,65 +52,10 @@ class Rope {
       return l;
     }
 
-    template <int p>
-    Real error(const Tensor<Scalar, p>& J) const {
-      Scalar H0, HN;
-      std::tie(H0, std::ignore, std::ignore) = hamGradHess(J, z.front());
-      std::tie(HN, std::ignore, std::ignore) = hamGradHess(J, z.back());
-
-      Real ImH = imag((H0 + HN) / 2.0);
-
-      Real err = 0;
-
-      for (unsigned i = 1; i < z.size() - 1; i++) {
-        Scalar Hi;
-        std::tie(Hi, std::ignore, std::ignore) = hamGradHess(J, z[i]);
-
-        err += pow(imag(Hi) - ImH, 2);
-      }
-
-      return sqrt(err);
-    }
-
     Vector<Scalar> dz(unsigned i) const {
       return z[i + 1] - z[i - 1];
     }
 
-    Vector<Scalar> ddz(unsigned i) const {
-      return 4.0 * (z[i + 1] + z[i - 1] - 2.0 * z[i]);
-    }
-
-    template <int p>
-    std::vector<Vector<Scalar>> generateGradientδz(const Tensor<Scalar, p>& J) const {
-      std::vector<Vector<Scalar>> δz(z.size());
-
-#pragma omp parallel for
-      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]);
-
-        δ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.
-      Real mag = 0;
-      for (unsigned i = 1; i < z.size() - 1; i++) {
-        mag += δz[i].norm();
-      }
-
-      for (unsigned i = 1; i < z.size() - 1; i++) {
-        δz[i] /= mag / n();
-      }
-
-      return δz;
-    }
-
     template <int p>
     std::vector<Vector<Scalar>> generateDiscreteGradientδz(const Tensor<Scalar, p>& J, Real γ) const {
       std::vector<Vector<Scalar>> δz(z.size());
@@ -185,29 +89,38 @@ class Rope {
           dC +=  0.5 * (ż[i + 1].conjugate() - dz(i + 1).conjugate() * real(ż[i + 1].dot(dz(i + 1))) / dz(i + 1).squaredNorm()) / (ż[i + 1].norm() * dz(i + 1).norm());
         }
 
-        dC += - γ * (z[i - 1] + z[i + 1]).conjugate();
+        dC += γ * (2 * z[i] - z[i - 1] - z[i + 1]).conjugate();
 
-        δz[i] = dC;
-      }
+        δz[i] = dC.conjugate();
 
-      Real size = 0;
-      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].conjugate() *  z[i].conjugate().dot(δz[i]) /  z²;
       }
 
       return δz;
     }
 
-    template<class Gen>
-    std::vector<Vector<Scalar>> generateRandomδz(Gen& r) const {
-      std::vector<Vector<Scalar>> δz(z.size());
+    void spread() {
+      Real l = length();
+
+      Real a = 0;
+      unsigned pos = 0;
+
+      std::vector<Vector<Scalar>> zNew = z;
 
-      complex_normal_distribution<> d(0, 1, 0);
       for (unsigned i = 1; i < z.size() - 1; i++) {
-        δz[i] = randomVector<Scalar>(z[0].size(), d, r);
+        Real b = i * l / (z.size() - 1);
+
+        while (b > a) {
+          pos++;
+          a += (z[pos] - z[pos - 1]).norm();
+        }
+
+        Vector<Scalar> δz = z[pos] - z[pos - 1];
+
+        zNew[i] = normalize(z[pos] - (a - b) / δz.norm() * δz);
       }
 
-      return δz;
+      z = zNew;
     }
 
     template<int p>
@@ -222,8 +135,6 @@ class Rope {
           rNew.z[i] = normalize(z[i] - (δ * Δl) * δz[i]);
         }
 
-        rNew.spread();
-
         if (rNew.cost(J, γ) < cost(J, γ)) {
           break;
         } else {
@@ -235,71 +146,33 @@ class Rope {
         }
       }
 
+//      rNew.spread();
+
       z = rNew.z;
 
       return δ;
     }
 
-    void spread() {
-      Real l = length();
-
-      Real a = 0;
-      unsigned pos = 0;
-
-      std::vector<Vector<Scalar>> zNew = z;
-
-      for (unsigned i = 1; i < z.size() - 1; i++) {
-        Real b = i * l / (z.size() - 1);
-
-        while (b > a) {
-          pos++;
-          a += (z[pos] - z[pos - 1]).norm();
-        } 
-
-        Vector<Scalar> δz = z[pos] - z[pos - 1];
-
-        zNew[i] = normalize(z[pos] - (a - b) / δz.norm() * δz);
-      }
-
-      z = zNew;
-    }
-
-    template <int p>
-    void relaxGradient(const Tensor<Scalar, p>& J, unsigned N, Real δ0) {
-      Real δ = δ0;
-      try {
-        for (unsigned i = 0; i < N; i++) {
-          std::vector<Vector<Scalar>> δz = generateGradientδz(J);
-          δ = 1.1 * perturb(J, δ, δz);
-        }
-      } catch (std::exception& e) {
-      }
-    }
-
     template <int p>
     void relaxDiscreteGradient(const Tensor<Scalar, p>& J, unsigned N, Real δ0, Real γ) {
       Real δ = δ0;
       try {
         for (unsigned i = 0; i < N; i++) {
           std::vector<Vector<Scalar>> δz = generateDiscreteGradientδz(J, γ);
-          δ = 1.1 * perturb(J, δ, δz, γ);
+          double stepSize = 0;
+          for (const Vector<Scalar>& v : δz) {
+            stepSize += v.norm();
+          }
+          if (stepSize / δz.size() < 1e-6) {
+            break;
+          }
+          std::cout << cost(J) << " " << stepSize / δz.size() << std::endl;
+          δ = 2 * perturb(J, δ, δz, γ);
         }
       } catch (std::exception& e) {
       }
     }
 
-    template <int p, class Gen>
-    void relaxRandom(const Tensor<Scalar, p>& J, unsigned N, Real δ0, Gen& r) {
-      Real δ = δ0;
-        for (unsigned i = 0; i < N; i++) {
-          try {
-            std::vector<Vector<Scalar>> δz = generateRandomδz(r);
-            δ = 1.1 * perturb(J, δ, δz);
-          } catch (std::exception& e) {
-        }
-      }
-    }
-
     template <int p>
     Real cost(const Tensor<Scalar, p>& J, Real γ = 0) const {
       Real c = 0;
@@ -334,3 +207,180 @@ class Rope {
       return r;
     }
 };
+
+template <class Real, class Scalar, int p>
+bool stokesLineTest(const Tensor<Scalar, p>& J, const Vector<Scalar>& z1, const Vector<Scalar>& z2, unsigned n0, unsigned steps) {
+  Rope stokes(n0, z1, z2, J);
+
+  Real oldCost = stokes.cost(J);
+
+  for (unsigned i = 0; i < steps; i++) {
+    stokes.relaxDiscreteGradient(J, 1e6, 1, pow(2, steps));
+
+    Real newCost = stokes.cost(J);
+
+    if (newCost > oldCost) {
+      return false;
+    }
+
+    oldCost = newCost;
+
+    stokes = stokes.interpolate();
+  }
+  return true;
+}
+
+template <class Scalar>
+class Cord {
+public:
+  std::vector<Vector<Scalar>> gs;
+  Vector<Scalar> z0;
+  Vector<Scalar> z1;
+
+  template <int p>
+  Cord(const Tensor<Scalar, p>& J, const Vector<Scalar>& z2, const Vector<Scalar>& z3, unsigned ng) : gs(ng, Vector<Scalar>::Zero(z2.size())) {
+    Scalar H2 = getHamiltonian(J, z2);
+    Scalar H3 = getHamiltonian(J, z3);
+
+    if (real(H2) > real(H3)) {
+      z0 = z2;
+      z1 = z3;
+    } else {
+      z0 = z3;
+      z1 = z2;
+    }
+  }
+
+  Real gCoeff(unsigned i, Real t) const {
+    return (1 - t) * t * pow(t, i);
+  }
+
+  Real dgCoeff(unsigned i, Real t) const {
+    return (i + 1) * (1 - t) * pow(t, i) - pow(t, i + 1);
+  }
+
+  Vector<Scalar> f(Real t) const {
+    Vector<Scalar> z = (1 - t) * z0 + t * z1;
+
+    for (unsigned i = 0; i < gs.size(); i++) {
+      z += gCoeff(i, t) * gs[i];
+    }
+
+    return z;
+  }
+
+  Vector<Scalar> df(Real t) const {
+    Vector<Scalar> z = z1 - z0;
+
+    for (unsigned i = 0; i < gs.size(); i++) {
+      z += dgCoeff(i, t) * gs[i];
+    }
+
+    return z;
+  }
+
+  template <int p>
+  Real cost(const Tensor<Scalar, p>& J, Real t) const {
+    Vector<Scalar> z = f(t);
+    Scalar H;
+    Vector<Scalar> dH;
+    std::tie(H, dH, std::ignore) = hamGradHess(J, z);
+    Vector<Scalar> ż = zDot(z, dH);
+    Vector<Scalar> dz = df(t);
+
+    return 1 - real(ż.dot(dz)) / ż.norm() / dz.norm();
+  }
+
+  template <int p>
+  Real totalCost(const Tensor<Scalar, p>& J, unsigned nt) const {
+    Real tc = 0;
+
+    for (unsigned i = 0; i < nt; i++) {
+      Real t = (i + 1.0) / (nt + 1.0);
+      tc += cost(J, t);
+    }
+
+    return tc;
+  }
+
+  template <int p>
+  std::vector<Vector<Scalar>> dgs(const Tensor<Scalar, p>& J, Real t) const {
+    Vector<Scalar> z = f(t);
+    auto [H, dH, ddH] = hamGradHess(J, z);
+    Vector<Scalar> ż = zDot(z, dH);
+    Vector<Scalar> dz = df(t);
+    Matrix<Scalar> dż = dzDot(z, dH);
+    Matrix<Scalar> dżc = dzDotConjugate(z, dH, ddH);
+
+    std::vector<Vector<Scalar>> x;
+    x.reserve(gs.size());
+
+    for (unsigned i = 0; i < gs.size(); i++) {
+      Real fdg = gCoeff(i, t);
+      Real dfdg = dgCoeff(i, t);
+      Vector<Scalar> dC = - 0.5 / ż.norm() / dz.norm() * (
+          dfdg * ż.conjugate() + fdg * dżc * dz + fdg * dż * dz.conjugate()
+          - real(dz.dot(ż)) * (
+              dfdg * dz.conjugate() / dz.squaredNorm() +
+              fdg * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm()
+            )
+        );
+
+      x.push_back(dC.conjugate());
+    }
+
+    return x;
+  }
+
+  template <int p>
+  std::pair<Real, Real> relaxStep(const Tensor<Scalar, p>& J, unsigned nt, Real δ₀) {
+    std::vector<Vector<Scalar>> dgsTot(gs.size(), Vector<Scalar>::Zero(z0.size()));
+
+    for (unsigned i = 0; i < nt; i++) {
+      Real t = (i + 1.0) / (nt + 1.0);
+      std::vector<Vector<Scalar>> dgsI = dgs(J, t);
+
+      for (unsigned j = 0; j < gs.size(); j++) {
+        dgsTot[j] += dgsI[j] / nt;
+      }
+    }
+
+    Real stepSize = 0;
+    for (const Vector<Scalar>& dgi : dgsTot) {
+      stepSize += dgi.squaredNorm();
+    }
+    stepSize = sqrt(stepSize);
+
+    Cord cNew(*this);
+
+    Real δ = δ₀;
+    Real oldCost = totalCost(J, nt);
+    Real newCost = std::numeric_limits<Real>::infinity();
+
+    while (newCost > oldCost) {
+      for (unsigned i = 0; i < gs.size(); i++) {
+        cNew.gs[i] = gs[i] - δ * dgsTot[i];
+      }
+
+      newCost = cNew.totalCost(J, nt);
+
+      δ /= 2;
+    }
+
+    gs = cNew.gs;
+
+    return {2 * δ, stepSize};
+  }
+
+  template <int p>
+  void relax(const Tensor<Scalar, p>& J, unsigned nt, Real δ₀, unsigned maxSteps) {
+    Real δ = δ₀;
+    Real stepSize = std::numeric_limits<Real>::infinity();
+    unsigned steps = 0;
+    while (stepSize > 1e-7 && steps < maxSteps) {
+      std::tie(δ, stepSize) = relaxStep(J, nt, δ);
+      std::cout << totalCost(J, nt) << " " << δ << " " << stepSize << std::endl;
+      steps++;
+    }
+  }
+};