From 4069a0f073f0de40c2e8b41c2d2d3558bb0d8be1 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Thu, 18 Feb 2021 08:41:19 +0100
Subject: More Stokes tweaks.

---
 langevin.cpp | 28 ++++++++++++----------------
 stokes.hpp   | 50 ++++++++++++++++++++++++++++++++++++++++++++------
 2 files changed, 56 insertions(+), 22 deletions(-)

diff --git a/langevin.cpp b/langevin.cpp
index e1464ec..e466860 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -182,30 +182,26 @@ int main(int argc, char* argv[]) {
 
   Rope<Complex> stokes(10, z1, z2);
 
-  std::cout << stokes.cost(Jp) << std::endl;
+  std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl;
 
-  stokes.relax(Jp, 10000, 1e-4);
+  stokes.relax(Jp, 1e5, 1e-1);
 
-  std::cout << stokes.cost(Jp) << std::endl;
+  std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl;
 
-  Rope<Complex> stokes2 = stokes.interpolate();
+  stokes = stokes.interpolate();
+  stokes.relax(Jp, 1e5, 1e-1);
 
-  stokes2.relax(Jp, 10000, 1e-4);
+  std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl;
 
-  std::cout << stokes2.cost(Jp) << std::endl;
+  stokes = stokes.interpolate();
+  stokes.relax(Jp, 1e5, 1e-1);
 
-  stokes2 = stokes2.interpolate();
+  std::cout << stokes.cost(Jp) << " " << stokes.length() << " " << stokes.error(Jp) << std::endl;
 
-  stokes2.relax(Jp, 10000, 1e-4);
-
-  std::cout << stokes2.cost(Jp) << std::endl;
-
-  stokes2 = stokes2.interpolate();
-
-  stokes2.relax(Jp, 10000, 1e-4);
-
-  std::cout << stokes2.cost(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 abf0c50..0756dc3 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -59,6 +59,36 @@ class Rope {
       return z.size() - 2;
     }
 
+    double length() const {
+      double l = 0;
+
+      for (unsigned i = 0; i < z.size() - 1; i++) {
+        l += (z[i + 1] - z[i]).norm();
+      }
+
+      return l;
+    }
+
+    template <int p>
+    double 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());
+
+      double ImH = imag((H0 + HN) / 2.0);
+
+      double 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);
+    }
+
     Rope(unsigned N, const Vector<Scalar>& z1, const Vector<Scalar>& z2) : z(N + 2) {
       for (unsigned i = 0; i < N + 2; i++) {
         z[i] = normalize(z1 + (z2 - z1) * ((double)i / (N + 1.0)));
@@ -86,6 +116,16 @@ class Rope {
         δz[i] = variation(z[i], dz(i), ddz(i), dH, ddH);
       }
 
+      // We return a δz with average norm of one.
+      double 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;
     }
 
@@ -95,10 +135,12 @@ class Rope {
 
       Rope rNew = *this;
       double δ = δ0;
+      // We rescale the step size by the distance between points.
+      double Δl = length() / (n() + 1);
 
       while (rNew.cost(J) >= cost(J)) {
         for (unsigned i = 1; i < z.size() - 1; i++) {
-          rNew.z[i] = normalize(z[i] - δ * Δz[i].conjugate());
+          rNew.z[i] = normalize(z[i] - (δ * Δl) * Δz[i].conjugate());
         }
 
         δ /= 2;
@@ -114,11 +156,7 @@ class Rope {
     }
 
     void spread() {
-      double l = 0;
-
-      for (unsigned i = 0; i < z.size() - 1; i++) {
-        l += (z[i + 1] - z[i]).norm();
-      }
+      double l = length();
 
       double a = 0;
       unsigned pos = 0;
-- 
cgit v1.2.3-70-g09d2