Generalized code to easily allow mixed p-spins.

1 files changed, 28 insertions, 22 deletions

diff --git a/dynamics.hpp b/dynamics.hpp
index 7d04b1a..36b2dfa 100644
--- a/dynamics.hpp
+++ b/dynamics.hpp
@@ -6,22 +6,28 @@
 
 #include "p-spin.hpp"
 
-template <class Real, int p>
-Vector<Real> findMinimum(const Tensor<Real, p>& J, const Vector<Real>& z0, Real ε) {
+template <class Real, int ...ps>
+Vector<Real> findMinimum(const pSpinModel<Real, ps...>& M, const Vector<Real>& z0, Real ε) {
   Vector<Real> z = z0;
   Real λ = 10;
 
-  auto [H, dH, ddH] = hamGradHess(J, z);
+  Real H;
+  Vector<Real> dH;
+  Matrix<Real> ddH;
+  std::tie(H, dH, ddH, std::ignore) = M.hamGradHess(z);
 
   Vector<Real> g = dH - z.dot(dH) * z / (Real)z.size();
-  Matrix<Real> M = ddH - (dH * z.transpose() + z.dot(dH) * Matrix<Real>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + 2.0 * z * z.transpose();
+  Matrix<Real> m = ddH - (dH * z.transpose() + z.dot(dH) * Matrix<Real>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + 2.0 * z * z.transpose();
 
   while (g.norm() / z.size() > ε && λ < 1e8) {
-    Vector<Real> dz = (M + λ * (Matrix<Real>)abs(M.diagonal().array()).matrix().asDiagonal()).partialPivLu().solve(g);
+    Vector<Real> dz = (m + λ * (Matrix<Real>)abs(m.diagonal().array()).matrix().asDiagonal()).partialPivLu().solve(g);
     dz -= z.dot(dz) * z / (Real)z.size();
     Vector<Real> zNew = normalize(z - dz);
 
-    auto [HNew, dHNew, ddHNew] = hamGradHess(J, zNew);
+    Real HNew;
+    Vector<Real> dHNew;
+    Matrix<Real> ddHNew;
+    std::tie(HNew, dHNew, ddHNew, std::ignore) = M.hamGradHess(zNew);
 
     if (HNew * 1.0001 <= H) {
       z = zNew;
@@ -30,7 +36,7 @@ Vector<Real> findMinimum(const Tensor<Real, p>& J, const Vector<Real>& z0, Real
       ddH = ddHNew;
 
       g = dH - z.dot(dH) * z / (Real)z.size();
-      M = ddH - (dH * z.transpose() + z.dot(dH) * Matrix<Real>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + 2.0 * z * z.transpose();
+      m = ddH - (dH * z.transpose() + z.dot(dH) * Matrix<Real>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + 2.0 * z * z.transpose();
 
       λ /= 2;
     } else {
@@ -41,28 +47,28 @@ Vector<Real> findMinimum(const Tensor<Real, p>& J, const Vector<Real>& z0, Real
   return z;
 }
 
-template <class Real, class Scalar, int p>
-Vector<Scalar> findSaddle(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, Real ε) {
+template <class Real, class Scalar, int ...ps>
+Vector<Scalar> findSaddle(const pSpinModel<Scalar, ps...>& M, const Vector<Scalar>& z0, Real ε) {
   Vector<Scalar> z = z0;
 
   Vector<Scalar> dH;
   Matrix<Scalar> ddH;
-  std::tie(std::ignore, dH, ddH) = hamGradHess(J, z);
+  std::tie(std::ignore, dH, ddH, std::ignore) = M.hamGradHess(z);
 
   Scalar zz = z.transpose() * z;
   Vector<Scalar> g = dH - (z.transpose() * dH) * z / (Real)z.size() + z * (zz - (Real)z.size());
-  Matrix<Scalar> M = ddH - (dH * z.transpose() + (z.transpose() * dH) * Matrix<Scalar>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size()  + Matrix<Scalar>::Identity(z.size(), z.size()) * (zz - (Real)z.size()) + 2.0 * z * z.transpose();
+  Matrix<Scalar> m = ddH - (dH * z.transpose() + (z.transpose() * dH) * Matrix<Scalar>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size()  + Matrix<Scalar>::Identity(z.size(), z.size()) * (zz - (Real)z.size()) + 2.0 * z * z.transpose();
 
   while (g.norm() / z.size() > ε) {
-    Vector<Scalar> dz = M.partialPivLu().solve(g);
+    Vector<Scalar> dz = m.partialPivLu().solve(g);
     dz -= z.conjugate().dot(dz) / z.squaredNorm() * z.conjugate();
     z = normalize(z - dz);
 
-    std::tie(std::ignore, dH, ddH) = hamGradHess(J, z);
+    std::tie(std::ignore, dH, ddH, std::ignore) = M.hamGradHess(z);
 
     zz = z.transpose() * z;
     g = dH - (z.transpose() * dH) * z / (Real)z.size() + z * (zz - (Real)z.size());
-    M = ddH - (dH * z.transpose() + (z.transpose() * dH) * Matrix<Scalar>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size()  + Matrix<Scalar>::Identity(z.size(), z.size()) * (zz - (Real)z.size()) + 2.0 * z * z.transpose();
+    m = ddH - (dH * z.transpose() + (z.transpose() * dH) * Matrix<Scalar>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size()  + Matrix<Scalar>::Identity(z.size(), z.size()) * (zz - (Real)z.size()) + 2.0 * z * z.transpose();
   }
 
   return z;
@@ -79,16 +85,16 @@ Vector<Scalar> randomVector(unsigned N, Distribution d, Generator& r) {
   return z;
 }
 
-template <class Real, int p, class Distribution, class Generator>
-Vector<Real> randomMinimum(const Tensor<Real, p>& J, Distribution d, Generator& r, Real ε) {
+template <class Real, class Distribution, class Generator, int ...ps>
+Vector<Real> randomMinimum(const pSpinModel<Real, ps...>& M, Distribution d, Generator& r, Real ε) {
   Vector<Real> zSaddle;
   bool foundSaddle = false;
 
   while (!foundSaddle) {
-    Vector<Real> z0 = normalize(randomVector<Real>(J.dimension(0), d, r.engine()));
+    Vector<Real> z0 = normalize(randomVector<Real>(M.dimension(), d, r.engine()));
 
     try {
-      zSaddle = findMinimum(J, z0, ε);
+      zSaddle = findMinimum(M, z0, ε);
       foundSaddle = true;
     } catch (std::exception& e) {}
   }
@@ -96,16 +102,16 @@ Vector<Real> randomMinimum(const Tensor<Real, p>& J, Distribution d, Generator&
   return zSaddle;
 }
 
-template <class Real, class Scalar, int p, class Distribution, class Generator>
-Vector<Scalar> randomSaddle(const Tensor<Scalar, p>& J, Distribution d, Generator& r, Real ε) {
+template <class Real, class Scalar, class Distribution, class Generator, int ... ps>
+Vector<Scalar> randomSaddle(const pSpinModel<Real, ps...>& M, Distribution d, Generator& r, Real ε) {
   Vector<Scalar> zSaddle;
   bool foundSaddle = false;
 
   while (!foundSaddle) {
-    Vector<Scalar> z0 = normalize(randomVector<Scalar>(J.dimension(0), d, r.engine()));
+    Vector<Scalar> z0 = normalize(randomVector<Scalar>(M.dimension(), d, r.engine()));
 
     try {
-      zSaddle = findSaddle(J, z0, ε);
+      zSaddle = findSaddle(M, z0, ε);
       foundSaddle = true;
     } catch (std::exception& e) {}
   }