More changes.

4 files changed, 18 insertions, 25 deletions

diff --git a/complex_normal.hpp b/complex_normal.hpp
index 2ffcea1..65b054f 100644
--- a/complex_normal.hpp
+++ b/complex_normal.hpp
@@ -15,7 +15,7 @@ public:
   complex_normal_distribution(std::complex<T> μ, T σ, std::complex<T> κ) : μ(μ), d(0, σ / sqrt(2)) {
     δx = sqrt(1 + std::abs(κ));
     δy = sqrt(1 - std::abs(κ));
-    eθ = std::polar(1.0, std::arg(κ));
+    eθ = std::polar((T)1, std::arg(κ));
   }
 
   template <class Generator> std::complex<T> operator()(Generator& g) {
diff --git a/langevin.cpp b/langevin.cpp
index 8c191f3..e5b370c 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -37,7 +37,7 @@ std::list<std::array<ComplexVector, 2>> saddlesCloserThan(const std::unordered_m
 
 template <class Generator>
 std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(const ComplexTensor& J0, const ComplexVector& z10, const ComplexVector& z20, Real ε, Real Δ, Generator& r) {
-  Real σ = sqrt(factorial(p) / (2.0 * pow(z10.size(), p - 1)));
+  Real σ = sqrt(factorial(p) / (Real(2) * pow(z10.size(), p - 1)));
   complex_normal_distribution<Real> dJ(0, σ, 0);
 
   ComplexTensor J = J0;
@@ -46,7 +46,7 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(c
   std::tie(H1, std::ignore, std::ignore) = hamGradHess(J, z10);
   std::tie(H2, std::ignore, std::ignore) = hamGradHess(J, z20);
   Real prevdist = abs(imag(H1-H2));
-  Real γ = 0.1 * prevdist;
+  Real γ = prevdist / 10;
 
   std::function<void(ComplexTensor&, std::array<unsigned, p>)> perturbJ =
     [&γ, &dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) {
@@ -67,13 +67,12 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(c
       std::tie(H1, std::ignore, std::ignore) = hamGradHess(Jp, z1);
       std::tie(H2, std::ignore, std::ignore) = hamGradHess(Jp, z2);
 
-      if (abs(imag(H1 - H2)) < prevdist && dist > 1e-2) {
+      if (abs(imag(H1 - H2)) < prevdist && dist > Real(1) / 100) {
         J = Jp;
         prevdist = abs(imag(H1 - H2));
-        γ = 0.1 * prevdist;
-        std::cerr << prevdist << std::endl;
+        γ = prevdist / 10;
 
-        if (abs(imag(H1 - H2)) < ε && dist > 1e-2) {
+        if (abs(imag(H1 - H2)) < ε && dist > Real(1) / 100) {
           break;
         }
       }
@@ -140,7 +139,7 @@ int main(int argc, char* argv[]) {
   }
 
   Complex κ(Rκ, Iκ);
-  Real σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1)));
+  Real σ = sqrt(factorial(p) / ((Real)2 * pow(N, p - 1)));
 
   Rng r;
 
@@ -185,7 +184,9 @@ int main(int argc, char* argv[]) {
   ComplexVector z1 = nearbySaddles.front()[0];
   ComplexVector z2 = nearbySaddles.front()[1];
 
-  std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, 1e-14, ε, r);
+  std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, 1e-15, ε, r);
+
+  std::cerr << "Im H is now sufficently close, starting to relax rope." << std::endl;
 
   Rope stokes(10, z1, z2, J);
 
diff --git a/stereographic.hpp b/stereographic.hpp
index 2dfa735..58c6b31 100644
--- a/stereographic.hpp
+++ b/stereographic.hpp
@@ -10,13 +10,13 @@ Vector<Scalar> stereographicToEuclidean(const Vector<Scalar>& ζ) {
   Vector<Scalar> z(N);
 
   Scalar a = ζ.transpose() * ζ;
-  Scalar b = 2 * sqrt(N) / (1.0 + a);
+  Scalar b = 2 * sqrt((Real)N) / ((Real)1 + a);
 
   for (unsigned i = 0; i < N - 1; i++) {
     z(i) = ζ(i);
   }
 
-  z(N - 1) = (a - 1.0) / 2.0;
+  z(N - 1) = (a - (Real)1) / (Real)2;
 
   return b * z;
 }
@@ -30,7 +30,7 @@ Vector<Scalar> euclideanToStereographic(const Vector<Scalar>& z) {
     ζ(i) = z(i);
   }
 
-  return ζ / (sqrt(N) - z(N - 1));
+  return ζ / (sqrt((Real)N) - z(N - 1));
 }
 
 template <class Scalar>
@@ -38,14 +38,14 @@ Matrix<Scalar> stereographicJacobian(const Vector<Scalar>& ζ) {
   unsigned N = ζ.size();
   Matrix<Scalar> J(N, N + 1);
 
-  Scalar b = 1.0 + (Scalar)(ζ.transpose() * ζ);
+  Scalar b = (Real)1 + (Scalar)(ζ.transpose() * ζ);
 
   for (unsigned i = 0; i < N; i++) {
     for (unsigned j = 0; j < N; j++) {
       J(i, j) = - ζ(i) * ζ(j);
 
       if (i == j) {
-        J(i, j) += 0.5 * b;
+        J(i, j) += b / (Real)2;
       }
     }
 
diff --git a/stokes.hpp b/stokes.hpp
index cbf5ebb..cbe5437 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -1,7 +1,8 @@
+#pragma once
+
 #include "p-spin.hpp"
 #include "complex_normal.hpp"
 #include "dynamics.hpp"
-#include <iostream>
 
 class ropeRelaxationStallException: public std::exception {
   virtual const char* what() const throw() {
@@ -189,20 +190,11 @@ class Rope {
         δz[i] = dC;
       }
 
+      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();
       }
 
-      // 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;
     }