A little bit of refactoring in the Rope class.

1 files changed, 26 insertions, 26 deletions

diff --git a/stokes.hpp b/stokes.hpp
index 117f4de..abf0c50 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -1,11 +1,5 @@
 #include "p-spin.hpp"
 
-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();
-}
-
 class ropeRelaxationStallException: public std::exception {
   virtual const char* what() const throw() {
     return "Gradient descent stalled.";
@@ -61,6 +55,10 @@ class Rope {
   public:
     std::vector<Vector<Scalar>> z;
 
+    unsigned n() const {
+      return z.size() - 2;
+    }
+
     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)));
@@ -85,7 +83,7 @@ class Rope {
         Matrix<Scalar> ddH;
         std::tie(std::ignore, dH, ddH) = hamGradHess(J, z[i]);
 
-        δz[i] = variation(z[i], this->dz(i), this->ddz(i), dH, ddH);
+        δz[i] = variation(z[i], dz(i), ddz(i), dH, ddH);
       }
 
       return δz;
@@ -93,12 +91,12 @@ class Rope {
 
     template <int p>
     double perturb(const Tensor<Scalar, p>& J, double δ0) {
-      std::vector<Vector<Scalar>> Δz = this->δz(J);
+      std::vector<Vector<Scalar>> Δz = δz(J);
 
       Rope rNew = *this;
       double δ = δ0;
 
-      while (rNew.cost(J) >= this->cost(J)) {
+      while (rNew.cost(J) >= cost(J)) {
         for (unsigned i = 1; i < z.size() - 1; i++) {
           rNew.z[i] = normalize(z[i] - δ * Δz[i].conjugate());
         }
@@ -144,11 +142,22 @@ class Rope {
     }
 
     template <int p>
-    double step(const Tensor<Scalar, p>& J, double δ) {
-      double δNew = this->perturb(J, δ);
-      this->spread();
+    double step(const Tensor<Scalar, p>& J, double δ0) {
+      double δ = perturb(J, δ0);
+      spread();
 
-      return δNew;
+      return δ;
+    }
+
+    template <int p>
+    void relax(const Tensor<Scalar, p>& J, unsigned N, double δ0) {
+      double δ = δ0;
+      try {
+        for (unsigned i = 0; i < N; i++) {
+          δ = 2 * step(J, δ);
+        }
+      } catch (std::exception& e) {
+      }
     }
 
     template <int p>
@@ -159,25 +168,16 @@ class Rope {
         Vector<Scalar> dH;
         std::tie(std::ignore, dH, std::ignore) = hamGradHess(J, z[i]);
 
-        c += segmentCost(z[i], this->dz(i), dH);
+        Vector<Scalar> zD = zDot(z[i], dH);
+
+        c += 1.0 - pow(real(zD.dot(dz(i))), 2) / zD.squaredNorm() / dz(i).squaredNorm();
       }
 
       return c;
     }
 
-    template <int p>
-    void relax(const Tensor<Scalar, p>& J, unsigned N, double δ0) {
-      double δ = δ0;
-      try {
-        for (unsigned i = 0; i < N; i++) {
-          δ = 2 * this->step(J, δ);
-        }
-      } catch (std::exception& e) {
-      }
-    }
-
     Rope<Scalar> interpolate() const {
-      Rope<Scalar> r(2 * (z.size() - 2) + 1, z[0], z[z.size() - 1]);
+      Rope<Scalar> r(2 * n() + 1, z[0], z[z.size() - 1]);
 
       for (unsigned i = 0; i < z.size(); i++) {
         r.z[2 * i] = z[i];