From 5ee6815f0734b2089c5b4c068cc21f2983bdba24 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 5 Jan 2021 18:11:21 +0100
Subject: A lot of work, and fixed a huge bug regarding the meaning of .dot in
 the Eigen library for complex vectors.

---
 stereographic.hpp | 67 +++++++++++++++++++++++++++++++++++++------------------
 1 file changed, 45 insertions(+), 22 deletions(-)

(limited to 'stereographic.hpp')

diff --git a/stereographic.hpp b/stereographic.hpp
index 4109bc7..59432b7 100644
--- a/stereographic.hpp
+++ b/stereographic.hpp
@@ -1,12 +1,14 @@
 #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 = ζ.dot(ζ);
-  std::complex<double> b = 2.0 * sqrt(N) / (1.0 + a);
+  Scalar a = ζ.transpose() * ζ;
+  Scalar b = 2 * sqrt(N) / (1.0 + a);
 
   for (unsigned i = 0; i < N - 1; i++) {
     z(i) = b * ζ(i);
@@ -32,15 +34,14 @@ Matrix stereographicJacobian(const Vector& ζ) {
   unsigned N = ζ.size() + 1;
   Matrix J(N - 1, N);
 
-  Scalar b = ζ.dot(ζ);
+  Scalar b = ζ.transpose() * ζ;
 
   for (unsigned i = 0; i < N - 1; i++) {
     for (unsigned j = 0; j < N - 1; j++) {
-      J(i, j) = -4.0 * ζ(i) * ζ(j);
+      J(i, j) = - 4 * sqrt(N) * ζ(i) * ζ(j) / pow(1.0 + b, 2);
       if (i == j) {
-        J(i, j) += 2.0 * (1.0 + b);
+        J(i, j) += 2 * sqrt(N) * (1.0 + b) / pow(1.0 + b, 2);
       }
-      J(i, j) *= sqrt(N) / pow(1.0 + b, 2);
     }
 
     J(i, N - 1) = 4.0 * sqrt(N) * ζ(i) / pow(1.0 + b, 2);
@@ -49,25 +50,40 @@ Matrix stereographicJacobian(const Vector& ζ) {
   return J;
 }
 
-Matrix stereographicMetric(const Vector& ζ) {
-  unsigned N = ζ.size();
-  Matrix g(N, N);
+Matrix stereographicMetric(const Matrix& J) {
+  return J * J.adjoint();
+}
+
+Eigen::Tensor<Scalar, 3> stereographicChristoffel(const Vector& ζ, const Matrix& gInvJac) {
+  unsigned N = ζ.size() + 1;
+  Eigen::Tensor<Scalar, 3> dJ(N - 1, N - 1, N);
 
-  double a = ζ.cwiseAbs2().sum();
-  Scalar b = ζ.dot(ζ);
+  Scalar b = ζ.transpose() * ζ;
 
-  for (unsigned i = 0; i < N; i++) {
-    for (unsigned j = 0; j < N; j++) {
-      g(i, j) = 16.0 * (std::conj(ζ(i)) * ζ(j) * (1.0 + a) -
-                        std::real(std::conj(ζ(i) * ζ(j)) * (1.0 + b)));
+  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);
+        if (i == j) {
+          dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(k) / pow(1.0 + b, 3);
+        }
+        if (i == k) {
+          dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(j) / pow(1.0 + b, 3);
+        }
+        if (j == k) {
+          dJ(i, j, k) -= 4 * sqrt(N) * (1.0 + b) * ζ(i) / pow(1.0 + b, 3);
+        }
+      }
+      dJ(i, j, N - 1) = - 16 * sqrt(N) * ζ(i) * ζ(j) / pow(1.0 + b, 3);
       if (i == j) {
-        g(i, j) += 4.0 * std::abs(1.0 + b);
+        dJ(i, j, N - 1) += 4 * sqrt(N) * (1.0 + b) * ζ(i) / pow(1.0 + b, 3);
       }
-      g(i, j) *= (N + 1) / std::norm(pow(b, 2));
     }
   }
 
-  return g;
+  std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(2, 1)};
+
+  return dJ.contract(Eigen::TensorMap<Eigen::Tensor<const Scalar, 2>>(gInvJac.data(), N - 1, N), ip);
 }
 
 std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, const Vector& ζ) {
@@ -82,11 +98,18 @@ std::tuple<Scalar, Vector, Matrix> stereographicHamGradHess(const Tensor& J, con
   Matrix hessZ;
   std::tie(hamiltonian, gradZ, hessZ) = hamGradHess(J, z);
 
-  Matrix g = stereographicMetric(ζ);
-  Matrix gj = g.llt().solve(jacobian);
+  Matrix g = stereographicMetric(jacobian);
+  Matrix gInvJac = g.ldlt().solve(jacobian);
+
+  grad = gInvJac * gradZ;
+
+  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>(0, 0)};
+  Eigen::Tensor<Scalar, 2> H2 = Γ.contract(dfT, ip);
 
-  grad = gj * gradZ;
-  hess = gj * hessZ * gj.transpose();
+  hess = jacobian * hessZ * jacobian.transpose() + (Eigen::Map<Matrix>(H2.data(), ζ.size(), ζ.size())).transpose();
 
   return {hamiltonian, grad, hess};
 }
-- 
cgit v1.2.3-54-g00ecf