Moved some repeated code into the pSpinModel class.

3 files changed, 77 insertions, 79 deletions

diff --git a/langevin.cpp b/langevin.cpp
index 8ecd69c..0e5a616 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -81,17 +81,10 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  Complex κ(Rκ, Iκ);
-  Real σp = sqrt(factorial(p) / ((Real)2 * pow(N, p - 1)));
-  Real σ2 = sqrt(factorial(2) / ((Real)2 * pow(N, 2 - 1)));
-  Real σ4 = sqrt(factorial(4) / ((Real)2 * pow(N, 4 - 1)));
-
   Rng r;
 
-  pSpinModel<Real, 2, 4> pSpin;
-
-  std::get<0>(pSpin.Js) = generateCouplings<Real, 2>(N, std::normal_distribution<Real>(0, σ2), r.engine());
-  std::get<1>(pSpin.Js) = generateCouplings<Real, 4>(N, std::normal_distribution<Real>(0, σ4 / 10), r.engine());
+//  pSpinModel<Real, 2, 4> pSpin({J2, J4});
+  pSpinModel<Real, 2, 4>pSpin(N, r.engine(), 1, 0.1);
 
   std::normal_distribution<Real> Red(0, 1);
 
@@ -100,18 +93,17 @@ int main(int argc, char* argv[]) {
   RealVector dHr;
   RealMatrix ddHr;
   std::tie(Hr, dHr, ddHr, std::ignore) = pSpin.hamGradHess(zMin);
-  Eigen::EigenSolver<Matrix<Real>> eigenS(ddHr - (dHr * zMin.transpose() + zMin.dot(dHr) * Matrix<Real>::Identity(zMin.size(), zMin.size()) + (ddHr * zMin) * zMin.transpose()) / (Real)zMin.size() + 2.0 * zMin * zMin.transpose());
-  std::cout << eigenS.eigenvalues().transpose() << std::endl;
+  Eigen::EigenSolver<Matrix<Real>> eigenS(ddHr - (dHr * zMin.transpose() + zMin.dot(dHr) * Matrix<Real>::Identity(zMin.size(), zMin.size()) + (ddHr * zMin) * zMin.transpose()) / (Real)zMin.size() + zMin * zMin.transpose() / (Real)N);
   for (unsigned i = 0; i < N; i++) {
-    RealVector zNew = normalize(zMin + 0.01 * eigenS.eigenvectors().col(i).real());
-    std::cout << pSpin.getHamiltonian(zNew) - Hr << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << std::endl;
+    RealVector zNew = normalize(zMin + 1e-3 * eigenS.eigenvectors().col(i).real());
+    std::cout << eigenS.eigenvalues()(i) << " " << 2.0 * (pSpin.getHamiltonian(zNew) - Hr) / 1e-6 << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << " " << eigenS.eigenvectors().col(0).dot(eigenS.eigenvectors().col(i)) << std::endl;
   }
   std::cout << std::endl;
   getchar();
 
   complex_normal_distribution<Real> d(0, 1, 0);
 
-  pSpinModel<Complex, 2, 4> complexPSpin = pSpin.cast<Complex>();;
+  pSpinModel complexPSpin = pSpin.cast<Complex>();;
   ComplexVector zSaddle = zMin.cast<Complex>();
 
   ComplexVector zSaddleNext;
@@ -134,8 +126,7 @@ int main(int argc, char* argv[]) {
   Real φ = atan2( H2.imag() - H1.imag(), H1.real() - H2.real());
   std::cerr << (zSaddle - zSaddleNext).norm() / (Real)N << " " << φ << " " << H1 * exp(Complex(0, φ)) << " " << H2 * exp(Complex(0, φ)) << std::endl;
 
-  std::get<0>(complexPSpin.Js) = exp(Complex(0, φ)) * std::get<0>(complexPSpin.Js);
-  std::get<1>(complexPSpin.Js) = exp(Complex(0, φ)) * std::get<1>(complexPSpin.Js);
+  complexPSpin *= exp(Complex(0, φ));
 
   Cord test(complexPSpin, zSaddle, zSaddleNext, 3);
   test.relaxNewton(10, 1, 1e4);
diff --git a/p-spin.hpp b/p-spin.hpp
index 7a49bca..544abf1 100644
--- a/p-spin.hpp
+++ b/p-spin.hpp
@@ -12,73 +12,92 @@ Vector<typename Derived::Scalar> normalize(const Eigen::MatrixBase<Derived>& z)
   return z * sqrt((Real)z.size() / (typename Derived::Scalar)(z.transpose() * z));
 }
 
-template <class Scalar>
-std::tuple<Matrix<Scalar>, Tensor<Scalar, 3>> hamGradTensorHelper(const Vector<Scalar>& z, const Tensor<Scalar, 2>& J) {
-  Tensor<Scalar, 3> J3(z.size(), z.size(), z.size());;
-  J3.setZero();
-  Matrix<Scalar> Jz = Eigen::Map<const Matrix<Scalar>>(J.data(), z.size(), z.size());
-
-  return {Jz, J3};
-}
+template <typename Scalar, int... ps>
+class pSpinModel {
+private:
+  std::tuple<Matrix<Scalar>, Tensor<Scalar, 3>> hamGradTensorHelper(const Vector<Scalar>& z, const Tensor<Scalar, 2>& J) const {
+    Tensor<Scalar, 3> J3(z.size(), z.size(), z.size());;
+    J3.setZero();
+    Matrix<Scalar> Jz = Eigen::Map<const Matrix<Scalar>>(J.data(), z.size(), z.size());
 
-template <class Scalar, int p>
-std::tuple<Matrix<Scalar>, Tensor<Scalar, 3>> hamGradTensorHelper(const Vector<Scalar>& z, const Tensor<Scalar, p>& J) {
-  Tensor<Scalar, 3> J3 = contractDown(J, z);
-  Tensor<Scalar, 1> zT = Eigen::TensorMap<Tensor<const Scalar, 1>>(z.data(), {z.size()});
-  Tensor<Scalar, 2> J3zT = J3.contract(zT, ip00);
-  Matrix<Scalar> Jz = Eigen::Map<const Matrix<Scalar>>(J3zT.data(), z.size(), z.size());
+    return {Jz, J3};
+  }
 
-  return {Jz, J3};
-}
+  template <int p>
+  std::tuple<Matrix<Scalar>, Tensor<Scalar, 3>> hamGradTensorHelper(const Vector<Scalar>& z, const Tensor<Scalar, p>& J) const {
+    Tensor<Scalar, 3> J3 = contractDown(J, z);
+    Tensor<Scalar, 1> zT = Eigen::TensorMap<Tensor<const Scalar, 1>>(z.data(), {z.size()});
+    Tensor<Scalar, 2> J3zT = J3.contract(zT, ip00);
+    Matrix<Scalar> Jz = Eigen::Map<const Matrix<Scalar>>(J3zT.data(), z.size(), z.size());
 
-template <class Scalar, int p, int... ps>
-std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>, Tensor<Scalar, 3>> hamGradHessHelper(const Vector<Scalar>& z, const Tensor<Scalar, p>& J, const Tensor<Scalar, ps>& ...Js) {
-  auto [Jz, J3] = hamGradTensorHelper(z, J);
-  
-  Vector<Scalar> Jzz = Jz * z;
-  Scalar Jzzz = Jzz.transpose() * z;
-
-  Real pBang = factorial(p);
-
-  Tensor<Scalar, 3> dddH = ((p - 2) * (p - 1) * p / pBang) * J3;
-  Matrix<Scalar> ddH = ((p - 1) * p / pBang) * Jz;
-  Vector<Scalar> dH = (p / pBang) * Jzz;
-  Scalar H = Jzzz / pBang;
-
-  if constexpr (sizeof...(Js) > 0) {
-    auto [Hs, dHs, ddHs, dddHs] = hamGradHessHelper(z, Js...);
-    H += Hs;
-    dH += dHs;
-    ddH += ddHs;
-    dddH += dddHs;
+    return {Jz, J3};
   }
 
-  return {H, dH, ddH, dddH};
-}
+  template <int p, int... qs>
+  std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>, Tensor<Scalar, 3>> hamGradHessHelper(const Vector<Scalar>& z, const Tensor<Scalar, p>& J, const Tensor<Scalar, qs>& ...Js) const {
+    auto [Jz, J3] = hamGradTensorHelper(z, J);
+    
+    Vector<Scalar> Jzz = Jz * z;
+    Scalar Jzzz = Jzz.transpose() * z;
+
+    Real pBang = factorial(p);
+
+    Tensor<Scalar, 3> dddH = ((p - 2) * (p - 1) * p / pBang) * J3;
+    Matrix<Scalar> ddH = ((p - 1) * p / pBang) * Jz;
+    Vector<Scalar> dH = (p / pBang) * Jzz;
+    Scalar H = Jzzz / pBang;
+
+    if constexpr (sizeof...(Js) > 0) {
+      auto [Hs, dHs, ddHs, dddHs] = hamGradHessHelper(z, Js...);
+      H += Hs;
+      dH += dHs;
+      ddH += ddHs;
+      dddH += dddHs;
+    }
+
+    return {H, dH, ddH, dddH};
+  }
 
-template <typename Scalar, int... ps>
-class pSpinModel {
 public:
   std::tuple<Tensor<Scalar, ps>...> Js;
 
+  pSpinModel() {}
+
+  pSpinModel(const std::tuple<Tensor<Scalar, ps>...>& Js) : Js(Js) {}
+
+  template <class Generator, typename... T>
+  pSpinModel<Real>(unsigned N, Generator& r, T... μs) {
+    Js = std::make_tuple(μs * generateRealPSpinCouplings<Real, ps>(N, r)...);
+  }
+
   unsigned dimension() const {
     return std::get<0>(Js).dimension(0);
   }
 
   template <typename NewScalar>
   pSpinModel<NewScalar, ps...> cast() const {
-    pSpinModel<NewScalar, ps...> M;
-    M.Js = std::apply(
+    return std::apply(
         [] (const Tensor<Scalar, ps>& ...Ks) -> std::tuple<Tensor<NewScalar, ps>...> {
           return std::make_tuple(Ks.template cast<NewScalar>()...);
         }, Js
-        );
+      );
+  }
+
+  template <typename T>
+  pSpinModel<Scalar, ps...>& operator*=(T x) {
+    std::tuple<Tensor<Scalar, ps>...> newJs = std::apply(
+        [x] (const Tensor<Scalar, ps>& ...Ks) -> std::tuple<Tensor<Scalar, ps>...> {
+          return std::make_tuple((x * Ks)...);
+        }, Js
+      );
+
+    Js = newJs;
 
-    return M;
+    return *this;
   }
 
   std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>, Tensor<Scalar, 3>> hamGradHess(const Vector<Scalar>& z) const {
-    return std::apply([&z](const Tensor<Scalar, ps>& ...Ks) -> std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>, Tensor<Scalar, 3>> { return hamGradHessHelper(z, Ks...); }, Js);
+    return std::apply([&z, this](const Tensor<Scalar, ps>& ...Ks) -> std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>, Tensor<Scalar, 3>> { return hamGradHessHelper(z, Ks...); }, Js);
   }
 
   Scalar getHamiltonian(const Vector<Scalar>& z) const {
diff --git a/tensor.hpp b/tensor.hpp
index fc99042..1212bda 100644
--- a/tensor.hpp
+++ b/tensor.hpp
@@ -5,6 +5,8 @@
 
 #include <eigen3/unsupported/Eigen/CXX11/Tensor>
 
+#include "factorial.hpp"
+
 template <class Scalar>
 using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
 
@@ -89,25 +91,11 @@ Tensor<Scalar, p> generateCouplings(unsigned N, Distribution d, Generator& r) {
   return J;
 }
 
-template <class Scalar, int p, class Distribution, class Generator>
-Tensor<Scalar, p> plantState(const Tensor<Scalar, p>& J, const Vector<Scalar>& z, double β) {
-  Tensor<Scalar, p> JPlanted = J;
-
-  std::function<void(Tensor<Scalar, p>&, std::array<unsigned, p>)> plant =
-    [&z, β] (Tensor<Scalar, p>& JJ, std::array<unsigned, p> ind) {
-      Scalar Ji = getJ<Scalar, p>(JJ, ind);
-      Scalar zzz = 1;
-
-      for (unsigned i : ind) {
-        zzz *= z(i);
-      }
-
-      setJ<Scalar, p>(JJ, ind, Ji - β * zzz / pow(zzz.size(), p));
-    };
-
-  iterateOver<Scalar, p>(JPlanted, plant);
+template <class Real, int p, class Generator>
+Tensor<Real, p> generateRealPSpinCouplings(unsigned N, Generator& r) {
+  Real σp = sqrt(factorial(p) / ((Real)2 * pow(N, p - 1)));
 
-  return JPlanted;
+  return generateCouplings<Real, p>(N, std::normal_distribution<Real>(0, σp), r);
 }
 
 template <class Scalar>