From 7bf6e952b53699977f5091a78f0f9f48f7b359c5 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 9 Nov 2021 10:13:59 +0100
Subject: Generalized code to easily allow mixed p-spins.

---
 stokes.hpp | 70 ++++++++++++++++----------------------------------------------
 1 file changed, 18 insertions(+), 52 deletions(-)

(limited to 'stokes.hpp')

diff --git a/stokes.hpp b/stokes.hpp
index b3ae90b..c2b0ec5 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -4,17 +4,18 @@
 #include "complex_normal.hpp"
 #include "dynamics.hpp"
 
-template <class Scalar>
+template <class Scalar, int ...ps>
 class Cord {
 public:
+  const pSpinModel<Scalar, ps...>& M;
   std::vector<Vector<Scalar>> gs;
   Vector<Scalar> z0;
   Vector<Scalar> z1;
+  double φ;
 
-  template <int p>
-  Cord(const Tensor<Scalar, p>& J, const Vector<Scalar>& z2, const Vector<Scalar>& z3, unsigned ng) : gs(ng, Vector<Scalar>::Zero(z2.size())) {
-    Scalar H2 = getHamiltonian(J, z2);
-    Scalar H3 = getHamiltonian(J, z3);
+  Cord(const pSpinModel<Scalar, ps...>& M, const Vector<Scalar>& z2, const Vector<Scalar>& z3, unsigned ng) : M(M), gs(ng, Vector<Scalar>::Zero(z2.size())) {
+    Scalar H2 = M.getHamiltonian(z2);
+    Scalar H3 = M.getHamiltonian(z3);
 
     if (real(H2) > real(H3)) {
       z0 = z2;
@@ -53,64 +54,31 @@ public:
     return z;
   }
 
-  template <int p>
-  Real cost(const Tensor<Scalar, p>& J, Real t) const {
+  Real cost(Real t) const {
     Vector<Scalar> z = f(t);
     Scalar H;
     Vector<Scalar> dH;
-    std::tie(H, dH, std::ignore) = hamGradHess(J, z);
+    std::tie(H, dH, std::ignore, std::ignore) = M.hamGradHess(z);
     Vector<Scalar> ż = zDot(z, dH);
     Vector<Scalar> dz = df(t);
 
     return 1 - real(ż.dot(dz)) / ż.norm() / dz.norm();
   }
 
-  template <int p>
-  Real totalCost(const Tensor<Scalar, p>& J, unsigned nt) const {
+  Real totalCost(unsigned nt) const {
     Real tc = 0;
 
     for (unsigned i = 0; i < nt; i++) {
       Real t = (i + 1.0) / (nt + 1.0);
-      tc += cost(J, t);
+      tc += cost(t);
     }
 
     return tc;
   }
 
-  template <int p>
-  std::vector<Vector<Scalar>> dgs(const Tensor<Scalar, p>& J, Real t) const {
+  std::pair<Vector<Scalar>, Matrix<Scalar>> gsGradHess(Real t) const {
     Vector<Scalar> z = f(t);
-    auto [H, dH, ddH] = hamGradHess(J, z);
-    Vector<Scalar> ż = zDot(z, dH);
-    Vector<Scalar> dz = df(t);
-    Matrix<Scalar> dż = dzDot(z, dH);
-    Matrix<Scalar> dżc = dzDotConjugate(z, dH, ddH);
-
-    std::vector<Vector<Scalar>> x;
-    x.reserve(gs.size());
-
-    for (unsigned i = 0; i < gs.size(); i++) {
-      Real fdg = gCoeff(i, t);
-      Real dfdg = dgCoeff(i, t);
-      Vector<Scalar> dC = - 0.5 / ż.norm() / dz.norm() * (
-          dfdg * ż.conjugate() + fdg * dżc * dz + fdg * dż * dz.conjugate()
-          - real(dz.dot(ż)) * (
-              dfdg * dz.conjugate() / dz.squaredNorm() +
-              fdg * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm()
-            )
-        );
-
-      x.push_back(dC.conjugate());
-    }
-
-    return x;
-  }
-
-  template <int p>
-  std::pair<Vector<Scalar>, Matrix<Scalar>> gsGradHess(const Tensor<Scalar, p>& J, Real t) const {
-    Vector<Scalar> z = f(t);
-    auto [H, dH, ddH] = hamGradHess(J, z);
-    Tensor<Scalar, 3> dddH = ((p - 2) * (p - 1) * p / (Real)factorial(p)) * J; // BUG IF P != 3
+    auto [H, dH, ddH, dddH] = M.hamGradHess(z);
     Vector<Scalar> ż = zDot(z, dH);
     Tensor<Scalar, 1> żT = Eigen::TensorMap<Tensor<const Scalar, 1>>(ż.data(), {z.size()});
     Vector<Scalar> dz = df(t);
@@ -241,14 +209,13 @@ public:
     return {x, M};
   }
 
-  template <int p>
-  std::pair<Real, Real> relaxStepNewton(const Tensor<Scalar, p>& J, unsigned nt, Real δ₀) {
+  std::pair<Real, Real> relaxStepNewton(unsigned nt, Real δ₀) {
     Vector<Scalar> dgsTot = Vector<Scalar>::Zero(2 * z0.size() * gs.size());
     Matrix<Scalar> HessTot = Matrix<Scalar>::Zero(2 * z0.size() * gs.size(), 2 * z0.size() * gs.size());
 
     for (unsigned i = 0; i < nt; i++) {
       Real t = (i + 1.0) / (nt + 1.0);
-      auto [dgsi, Hessi] = gsGradHess(J, t);
+      auto [dgsi, Hessi] = gsGradHess(t);
 
       dgsTot += dgsi / nt;
       HessTot += Hessi / nt;
@@ -259,7 +226,7 @@ public:
     Cord cNew(*this);
 
     Real δ = δ₀;
-    Real oldCost = totalCost(J, nt);
+    Real oldCost = totalCost(nt);
     Real newCost = std::numeric_limits<Real>::infinity();
 
     Matrix<Scalar> I = abs(HessTot.diagonal().array()).matrix().asDiagonal();
@@ -273,7 +240,7 @@ public:
         }
       }
 
-      newCost = cNew.totalCost(J, nt);
+      newCost = cNew.totalCost(nt);
 
       δ *= 2;
     }
@@ -283,13 +250,12 @@ public:
     return {δ / 2, stepSize};
   }
 
-  template <int p>
-  void relaxNewton(const Tensor<Scalar, p>& J, unsigned nt, Real δ₀, unsigned maxSteps) {
+  void relaxNewton(unsigned nt, Real δ₀, unsigned maxSteps) {
     Real δ = δ₀;
     Real stepSize = std::numeric_limits<Real>::infinity();
     unsigned steps = 0;
     while (stepSize > 1e-7 && steps < maxSteps) {
-      std::tie(δ, stepSize) = relaxStepNewton(J, nt, δ);
+      std::tie(δ, stepSize) = relaxStepNewton(nt, δ);
       δ /= 3;
       steps++;
     }
-- 
cgit v1.2.3-54-g00ecf