Cleaned up code and removed unused portions.

2 files changed, 37 insertions, 96 deletions

diff --git a/langevin.cpp b/langevin.cpp
index 49efab4..e5ed0b8 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -69,41 +69,36 @@ Vector findSaddle(const Tensor& J, const Vector& z0, double ε, double δW, doub
   double W;
   std::tie(W, std::ignore) = WdW(J, z0);
 
-  Vector ζ = euclideanToStereographic(z0);
+  Vector z = z0;
+  Vector ζ = euclideanToStereographic(z);
+
   Vector dH;
   Matrix ddH;
-  std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ);
-
-  unsigned steps = 0;
-  unsigned gradSteps = 0;
+  std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ, z);
 
   while (W > ε) {
-    // ddH is complex symmetric, which is (almost always) invertible
+    // ddH is complex symmetric, which is (almost always) invertible, so a
+    // partial pivot LU decomposition can be used.
     Vector dζ = ddH.partialPivLu().solve(dH);
     Vector ζNew = ζ - dζ;
+    Vector zNew = stereographicToEuclidean(ζNew);
 
     double WNew;
-    std::tie(WNew, std::ignore) = WdW(J, stereographicToEuclidean(ζNew));
+    std::tie(WNew, std::ignore) = WdW(J, zNew);
 
     if (WNew < W) { // If Newton's step lowered the objective, accept it!
       ζ = ζNew;
+      z = zNew;
       W = WNew;
     } else { // Otherwise, do gradient descent until W is a factor δW smaller.
-      Vector z;
-      std::tie(W, z) = gradientDescent(J, stereographicToEuclidean(ζ), γ0, W / δW);
+      std::tie(W, z) = gradientDescent(J, z, γ0, W / δW);
       ζ = euclideanToStereographic(z);
-      gradSteps++;
     }
 
-    std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ);
-    steps++;
-
-    if (steps % 100 == 0) {
-      std::cerr << steps << " minimization steps, W is " << W << " " << gradSteps << "." << std::endl;
-    }
+    std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ, z);
   }
 
-  return stereographicToEuclidean(ζ);
+  return z;
 }
 
 std::tuple<double, Vector> langevin(const Tensor& J, const Vector& z0, double T, double γ, unsigned N, Rng& r) {
@@ -189,13 +184,7 @@ int main(int argc, char* argv[]) {
   Rng r;
 
   Tensor J = generateCouplings<Scalar, PSPIN_P>(N, complex_normal_distribution<>(0, σ, κ), r.engine());
-  Vector z0 = randomVector(N, complex_normal_distribution<>(0, 1, 0), r.engine());
-  z0 *= sqrt(N) / sqrt((Scalar)(z0.transpose() * z0)); // Normalize.
-
-  std::function<bool(double, unsigned)> f = [δ](double W, unsigned) {
-    std::cout << W << std::endl;
-    return W < δ;
-  };
+  Vector z0 = normalize(randomVector(N, complex_normal_distribution<>(0, 1, 0), r.engine()));
 
   Vector zSaddle = findSaddle(J, z0, ε, δ, γ);
 
@@ -209,6 +198,7 @@ int main(int argc, char* argv[]) {
     std::tie(H, std::ignore, ddH) = hamGradHess(J, zNewSaddle);
     Eigen::SelfAdjointEigenSolver<Matrix> es(ddH.adjoint() * ddH);
     std::cout << (zNewSaddle - zSaddle).norm() << " " << real(H) << " " << imag(H) << " " << es.eigenvalues().transpose() << std::endl;
+    std::cerr << M * (i+1) << " steps taken to move " << (zNewSaddle - zSaddle).norm() << ", saddle information saved." << std::endl;
   }
 
   return 0;
diff --git a/stereographic.hpp b/stereographic.hpp
index a3f2cc6..8313f25 100644
--- a/stereographic.hpp
+++ b/stereographic.hpp
@@ -1,22 +1,21 @@
 #include <eigen3/Eigen/Cholesky>
 
-#include "Eigen/src/Core/util/Meta.h"
 #include "p-spin.hpp"
-#include "unsupported/Eigen/CXX11/src/Tensor/TensorMeta.h"
 
 Vector stereographicToEuclidean(const Vector& ζ) {
   unsigned N = ζ.size() + 1;
   Vector z(N);
+
   Scalar a = ζ.transpose() * ζ;
   Scalar b = 2 * sqrt(N) / (1.0 + a);
 
   for (unsigned i = 0; i < N - 1; i++) {
-    z(i) = b * ζ(i);
+    z(i) = ζ(i);
   }
 
-  z(N - 1) = b * (a - 1.0) / 2.0;
+  z(N - 1) = (a - 1.0) / 2.0;
 
-  return z;
+  return b * z;
 }
 
 Vector euclideanToStereographic(const Vector& z) {
@@ -24,94 +23,46 @@ Vector euclideanToStereographic(const Vector& z) {
   Vector ζ(N - 1);
 
   for (unsigned i = 0; i < N - 1; i++) {
-    ζ(i) = z(i) / (sqrt(N) - z(N - 1));
+    ζ(i) = z(i);
   }
 
-  return ζ;
+  return ζ / (sqrt(N) - z(N - 1));
 }
 
 Matrix stereographicJacobian(const Vector& ζ) {
-  unsigned N = ζ.size() + 1;
-  Matrix J(N - 1, N);
-
-  Scalar b = ζ.transpose() * ζ;
-
-  for (unsigned i = 0; i < N - 1; i++) {
-    for (unsigned j = 0; j < N - 1; j++) {
-      J(i, j) = - 4 * sqrt(N) * ζ(i) * ζ(j) / pow(1.0 + b, 2);
-      if (i == j) {
-        J(i, j) += 2 * sqrt(N) * (1.0 + b) / pow(1.0 + b, 2);
-      }
-    }
-
-    J(i, N - 1) = 4.0 * sqrt(N) * ζ(i) / pow(1.0 + b, 2);
-  }
-
-  return J;
-}
-
-Matrix stereographicMetric(const Matrix& J) {
-  return J * J.adjoint();
-}
-
-// Gives the Christoffel symbol, with Γ_(i1, i2)^(i3).
-Eigen::Tensor<Scalar, 3> stereographicChristoffel(const Vector& ζ, const Matrix& gInvJacConj) {
-  unsigned N = ζ.size() + 1;
-  Eigen::Tensor<Scalar, 3> dJ(N - 1, N - 1, N);
+  unsigned N = ζ.size();
+  Matrix J(N, N + 1);
 
   Scalar b = 1.0 + (Scalar)(ζ.transpose() * ζ);
 
-  for (unsigned i = 0; i < N - 1; i++) {
-    for (unsigned j = 0; j < N - 1; j++) {
-      for (unsigned k = 0; k < N - 1; k++) {
-        dJ(i, j, k) = 16 * sqrt(N) * ζ(i) * ζ(j) * ζ(k) / pow(b, 3);
-        if (i == j) {
-          dJ(i, j, k) -= 4 * sqrt(N) * ζ(k) / pow(b, 2);
-        }
-        if (i == k) {
-          dJ(i, j, k) -= 4 * sqrt(N) * ζ(j) / pow(b, 2);
-        }
-        if (j == k) {
-          dJ(i, j, k) -= 4 * sqrt(N) * ζ(i) / pow(b, 2);
-        }
-      }
-      dJ(i, j, N - 1) = - 16 * sqrt(N) * ζ(i) * ζ(j) / pow(b, 3);
+  for (unsigned i = 0; i < N; i++) {
+    for (unsigned j = 0; j < N; j++) {
+      J(i, j) = - ζ(i) * ζ(j);
+
       if (i == j) {
-        dJ(i, j, N - 1) += 4 * sqrt(N) * ζ(i) / pow(b, 2);
+        J(i, j) += 0.5 * b;
       }
     }
-  }
 
-  std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(2, 1)};
+    J(i, N) = ζ(i);
+  }
 
-  return dJ.contract(Eigen::TensorMap<Eigen::Tensor<const Scalar, 2>>(gInvJacConj.data(), N - 1, N), ip);
+  return 4 * sqrt(N + 1) * J / pow(b, 2);
 }
 
-std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, const Vector& ζ) {
-  Vector grad;
-  Matrix hess;
-
-  Matrix jacobian = stereographicJacobian(ζ);
-  Vector z = stereographicToEuclidean(ζ);
-
-  std::complex<double> hamiltonian;
+std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, const Vector& ζ, const Vector& z) {
+  Scalar hamiltonian;
   Vector gradZ;
   Matrix hessZ;
   std::tie(hamiltonian, gradZ, hessZ) = hamGradHess(J, z);
 
-  Matrix metric = stereographicMetric(jacobian);
-
-  grad = metric.llt().solve(jacobian) * gradZ;
+  Matrix jacobian = stereographicJacobian(ζ);
 
-  /* This is much slower to calculate than the marginal speedup it offers...
-  Eigen::Tensor<Scalar, 3> Γ = stereographicChristoffel(ζ, gInvJac.conjugate());
-  Vector df = jacobian * gradZ;
-  Eigen::Tensor<Scalar, 1> dfT = Eigen::TensorMap<Eigen::Tensor<Scalar, 1>>(df.data(), {df.size()});
-  std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(2, 0)};
-  Eigen::Tensor<Scalar, 2> H2 = Γ.contract(dfT, ip);
-  */
+  Matrix metric = jacobian * jacobian.adjoint();
 
-  hess = jacobian * hessZ * jacobian.transpose(); // - Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size());
+  // The metric is Hermitian and positive definite, so a Cholesky decomposition can be used.
+  Vector grad = metric.llt().solve(jacobian) * gradZ;
+  Matrix hess = jacobian * hessZ * jacobian.transpose();
 
   return {hamiltonian, grad, hess};
 }