From 1094eaa434d791d3a511f853df9d6d66d169b2ab Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 6 Jan 2021 11:16:20 +0100
Subject: Lots of work, and fixed bug stemming from choice of matrix algorithm.

---
 langevin.cpp      |  5 ++++-
 stereographic.hpp | 20 ++++++++++----------
 2 files changed, 14 insertions(+), 11 deletions(-)

diff --git a/langevin.cpp b/langevin.cpp
index 3bc12ed..05b436e 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -7,6 +7,8 @@
 #include "pcg-cpp/include/pcg_random.hpp"
 #include "randutils/randutils.hpp"
 
+#include <eigen3/Eigen/LU>
+
 #include "complex_normal.hpp"
 #include "p-spin.hpp"
 #include "stereographic.hpp"
@@ -57,7 +59,8 @@ Vector findSaddle(const Tensor& J, const Vector& z0, double γ0, double δ, doub
   while (W > ε) {
     double γ = pow(r.variate<double, std::normal_distribution>(0, γ0), 2);
 
-    Vector ζNew = ζ - ddH.ldlt().solve(dH);
+    Vector dζ = ddH.partialPivLu().solve(dH);
+    Vector ζNew = ζ - dζ;
 
     double WNew;
     Vector dWNew;
diff --git a/stereographic.hpp b/stereographic.hpp
index 59432b7..4c3c831 100644
--- a/stereographic.hpp
+++ b/stereographic.hpp
@@ -58,25 +58,25 @@ Eigen::Tensor<Scalar, 3> stereographicChristoffel(const Vector& ζ, const Matrix
   unsigned N = ζ.size() + 1;
   Eigen::Tensor<Scalar, 3> dJ(N - 1, N - 1, N);
 
-  Scalar b = ζ.transpose() * ζ;
+  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(1.0 + b, 3);
+        dJ(i, j, k) = 16 * sqrt(N) * ζ(i) * ζ(j) * ζ(k) / pow(b, 3);
         if (i == j) {
-          dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(k) / pow(1.0 + b, 3);
+          dJ(i, j, k) -= 4 * sqrt(N) * ζ(k) / pow(b, 2);
         }
         if (i == k) {
-          dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(j) / pow(1.0 + b, 3);
+          dJ(i, j, k) -= 4 * sqrt(N) * ζ(j) / pow(b, 2);
         }
         if (j == k) {
-          dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(i) / pow(1.0 + b, 3);
+          dJ(i, j, k) -= 4 * sqrt(N) * ζ(i) / pow(b, 2);
         }
       }
-      dJ(i, j, N - 1) = - 16 * sqrt(N) * ζ(i) * ζ(j) / pow(1.0 + b, 3);
+      dJ(i, j, N - 1) = - 16 * sqrt(N) * ζ(i) * ζ(j) / pow(b, 3);
       if (i == j) {
-        dJ(i, j, N - 1) += 4 * sqrt(N) * (1.0 + b) * ζ(i) / pow(1.0 + b, 3);
+        dJ(i, j, N - 1) += 4 * sqrt(N) * ζ(i) / pow(b, 2);
       }
     }
   }
@@ -103,13 +103,13 @@ std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, con
 
   grad = gInvJac * gradZ;
 
-  Eigen::Tensor<Scalar, 3> Γ = stereographicChristoffel(ζ, gInvJac.conjugate());
+  Eigen::Tensor<Scalar, 3> Γ = stereographicChristoffel(ζ, gInvJac);
   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>(0, 0)};
+  std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(1, 0)};
   Eigen::Tensor<Scalar, 2> H2 = Γ.contract(dfT, ip);
 
-  hess = jacobian * hessZ * jacobian.transpose() + (Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size())).transpose();
+  hess = jacobian * hessZ * jacobian.transpose() + Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size());
 
   return {hamiltonian, grad, hess};
 }
-- 
cgit v1.2.3-70-g09d2