Many changes.

1 files changed, 122 insertions, 32 deletions

diff --git a/stokes.hpp b/stokes.hpp
index 6cc5c4a..cbf5ebb 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -1,4 +1,6 @@
 #include "p-spin.hpp"
+#include "complex_normal.hpp"
+#include "dynamics.hpp"
 #include <iostream>
 
 class ropeRelaxationStallException: public std::exception {
@@ -9,13 +11,13 @@ 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) {
-  double z² = z.squaredNorm();
-  double z´² = z´.squaredNorm();
+  Real z² = z.squaredNorm();
+  Real z´² = z´.squaredNorm();
 
   Vector<Scalar> ż = zDot(z, dH);
-  double ż² = ż.squaredNorm();
+  Real ż² = ż.squaredNorm();
 
-  double Reż·z´ = real(ż.dot(z´));
+  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²
@@ -33,7 +35,7 @@ Vector<Scalar> variation(const Vector<Scalar>& z, const Vector<Scalar>& z´, con
         ) * z.conjugate()
       ) / z²;
 
-  double dReż·z´ = real(ż.dot(z´´) + ż´.dot(z´));
+  Real dReż·z´ = real(ż.dot(z´´) + ż´.dot(z´));
 
   Vector<Scalar> ddtdLdz´ = - (
       (
@@ -72,7 +74,7 @@ class Rope {
       }
 
       for (unsigned i = 0; i < N + 2; i++) {
-        z[i] = normalize(z.front() + (z.back() - z.front()) * ((double)i / (N + 1.0)));
+        z[i] = normalize(z.front() + (z.back() - z.front()) * ((Real)i / (N + 1.0)));
       }
     }
 
@@ -80,8 +82,8 @@ class Rope {
       return z.size() - 2;
     }
 
-    double length() const {
-      double l = 0;
+    Real length() const {
+      Real l = 0;
 
       for (unsigned i = 0; i < z.size() - 1; i++) {
         l += (z[i + 1] - z[i]).norm();
@@ -91,14 +93,14 @@ class Rope {
     }
 
     template <int p>
-    double error(const Tensor<Scalar, p>& J) const {
+    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());
 
-      double ImH = imag((H0 + HN) / 2.0);
+      Real ImH = imag((H0 + HN) / 2.0);
 
-      double err = 0;
+      Real err = 0;
 
       for (unsigned i = 1; i < z.size() - 1; i++) {
         Scalar Hi;
@@ -119,7 +121,7 @@ class Rope {
     }
 
     template <int p>
-    std::vector<Vector<Scalar>> δz(const Tensor<Scalar, p>& J) const {
+    std::vector<Vector<Scalar>> generateGradientδz(const Tensor<Scalar, p>& J) const {
       std::vector<Vector<Scalar>> δz(z.size());
 
 #pragma omp parallel for
@@ -137,7 +139,7 @@ class Rope {
       }
 
       // We return a δz with average norm of one.
-      double mag = 0;
+      Real mag = 0;
       for (unsigned i = 1; i < z.size() - 1; i++) {
         mag += δz[i].norm();
       }
@@ -150,20 +152,87 @@ class Rope {
     }
 
     template <int p>
-    double perturb(const Tensor<Scalar, p>& J, double δ0) {
-      std::vector<Vector<Scalar>> Δz = δz(J);
+    std::vector<Vector<Scalar>> generateDiscreteGradientδz(const Tensor<Scalar, p>& J, Real γ) const {
+      std::vector<Vector<Scalar>> δz(z.size());
+
+      std::vector<Vector<Scalar>> ż(z.size());
+      std::vector<Vector<Scalar>> dH(z.size());
+      std::vector<Matrix<Scalar>> ddH(z.size());
+
+#pragma omp parallel for
+      for (unsigned i = 1; i < z.size() - 1; i++) {
+        std::tie(std::ignore, dH[i], ddH[i]) = hamGradHess(J, z[i]);
+        ż[i] = zDot(z[i], dH[i]);
+      }
+
+#pragma omp parallel for
+      for (unsigned i = 1; i < z.size() - 1; i++) {
+        Real z² = z[i].squaredNorm();
+        Matrix<Scalar> dż = (dH[i].conjugate() - (dH[i].dot(z[i]) / z²) * z[i].conjugate()) * z[i].adjoint() / z²;
+        Matrix<Scalar> dżc = -ddH[i] + (ddH[i] * z[i].conjugate()) * z[i].transpose() / z²
+          + (z[i].dot(dH[i]) / z²) * (
+              Matrix<Scalar>::Identity(ddH[i].rows(), ddH[i].cols()) - z[i].conjugate() * z[i].transpose() / z²
+            );
+
+        Vector<Scalar> dC = -0.5 * (dżc * dz(i) + dż * dz(i).conjugate() - (dżc * ż[i] + dż * ż[i].conjugate()) * real(ż[i].dot(dz(i))) / ż[i].squaredNorm()) / (ż[i].norm() * dz(i).norm());
+
+        if (i > 1) {
+          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());
+        }
+
+        if (i < z.size() - 2) {
+          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();
 
+        δz[i] = dC;
+      }
+
+      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();
+      }
+
+      // 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<class Gen>
+    std::vector<Vector<Scalar>> generateRandomδz(Gen& r) const {
+      std::vector<Vector<Scalar>> δz(z.size());
+
+      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);
+      }
+
+      return δz;
+    }
+
+    template<int p>
+    Real perturb(const Tensor<Scalar, p>& J, Real δ₀, const std::vector<Vector<Scalar>>& δz, Real γ = 0) {
       Rope rNew = *this;
-      double δ = δ0;
+      Real δ = δ₀;
       // We rescale the step size by the distance between points.
-      double Δl = length() / (n() + 1);
+      Real Δl = length() / (n() + 1);
 
       while (true) {
         for (unsigned i = 1; i < z.size() - 1; i++) {
-          rNew.z[i] = normalize(z[i] - (δ * Δl) * Δz[i]);
+          rNew.z[i] = normalize(z[i] - (δ * Δl) * δz[i]);
         }
 
-        if (rNew.cost(J) < cost(J)) {
+        rNew.spread();
+
+        if (rNew.cost(J, γ) < cost(J, γ)) {
           break;
         } else {
           δ /= 2;
@@ -180,15 +249,15 @@ class Rope {
     }
 
     void spread() {
-      double l = length();
+      Real l = length();
 
-      double a = 0;
+      Real a = 0;
       unsigned pos = 0;
 
       std::vector<Vector<Scalar>> zNew = z;
 
       for (unsigned i = 1; i < z.size() - 1; i++) {
-        double b = i * l / (z.size() - 1);
+        Real b = i * l / (z.size() - 1);
 
         while (b > a) {
           pos++;
@@ -204,27 +273,44 @@ class Rope {
     }
 
     template <int p>
-    double step(const Tensor<Scalar, p>& J, double δ0) {
-      double δ = perturb(J, δ0);
-      spread();
-
-      return δ;
+    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 relax(const Tensor<Scalar, p>& J, unsigned N, double δ0) {
-      double δ = δ0;
+    void relaxDiscreteGradient(const Tensor<Scalar, p>& J, unsigned N, Real δ0, Real γ) {
+      Real δ = δ0;
       try {
         for (unsigned i = 0; i < N; i++) {
-          δ = 1.1 * step(J, δ);
+          std::vector<Vector<Scalar>> δz = generateDiscreteGradientδz(J, γ);
+          δ = 1.1 * 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>
-    double cost(const Tensor<Scalar, p>& J) const {
-      double c = 0;
+    Real cost(const Tensor<Scalar, p>& J, Real γ = 0) const {
+      Real c = 0;
 
       for (unsigned i = 1; i < z.size() - 1; i++) {
         Vector<Scalar> dH;
@@ -235,6 +321,10 @@ class Rope {
         c += 1.0 - real(zD.dot(dz(i))) / zD.norm() / dz(i).norm();
       }
 
+      for (unsigned i = 0; i < z.size() - 1; i++) {
+        c += γ * (z[i + 1] - z[i]).squaredNorm();
+      }
+
       return c;
     }