More changes.

4 files changed, 69 insertions, 41 deletions

diff --git a/dynamics.hpp b/dynamics.hpp
index 85eef71..b373659 100644
--- a/dynamics.hpp
+++ b/dynamics.hpp
@@ -92,26 +92,26 @@ Vector<Scalar> randomVector(unsigned N, Distribution d, Generator& r) {
 }
 
 template <class Scalar, int p, class Distribution, class Generator>
-std::tuple<double, Vector<Scalar>> langevin(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, double T, double γ, unsigned N, Distribution d, Generator& r) {
+std::tuple<double, Vector<Scalar>> metropolis(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0,
+    std::function<double(const Tensor<Scalar, p>&, const Vector<Scalar>&)>& energy,
+    double T, double γ, unsigned N, Distribution d, Generator& r) {
   Vector<Scalar> z = z0;
 
-  double W;
-  std::tie(W, std::ignore) = WdW(J, z);
+  double E = energy(J, z);
 
   std::uniform_real_distribution<double> D(0, 1);
 
   for (unsigned i = 0; i < N; i++) {
-    Vector<Scalar> zNewTmp = z + randomVector<Scalar>(z.size(), d, r);
+    Vector<Scalar> zNewTmp = z + γ * randomVector<Scalar>(z.size(), d, r);
     Vector<Scalar> zNew = normalize(zNewTmp);
 
-    double WNew;
-    std::tie(WNew, std::ignore) = WdW(J, zNew);
+    double ENew = energy(J, zNew);
 
-    if (exp((W - WNew) / T) > D(r)) {
+    if (E - ENew > T * log(D(r))) {
       z = zNew;
-      W = WNew;
+      E = ENew;
     }
   }
 
-  return {W, z};
+  return {E, z};
 }
diff --git a/langevin.cpp b/langevin.cpp
index 870879b..5ca2d72 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -83,11 +83,41 @@ int main(int argc, char* argv[]) {
   ComplexVector z0 = normalize(randomVector<Complex>(N, d, r.engine()));
 
   ComplexVector zSaddle = findSaddle(J, z0, ε);
-  ComplexVector zSaddlePrev = ComplexVector::Zero(N);
   ComplexVector z = zSaddle;
 
+  std::function<double(const ComplexTensor&, const ComplexVector&)> energyNormGrad = []
+    (const ComplexTensor& J, const ComplexVector& z) {
+      double W;
+      std::tie(W, std::ignore) = WdW(J, z);
+      return W;
+    };
+
+  double aGoal = 1e3;
+
+  std::function<double(const ComplexTensor&, const ComplexVector&)> energyInvA = [aGoal]
+    (const ComplexTensor& J, const ComplexVector& z) {
+      double a = z.squaredNorm();
+      if (a > aGoal) {
+        return -aGoal;
+      } else {
+        return -z.squaredNorm();
+      }
+    };
+
+  while (zSaddle.squaredNorm() < aGoal) {
+    std::tie(std::ignore, z) = metropolis(J, z, energyInvA, T, γ, 100, d, r.engine());
+    try {
+      std::cerr << "Starting descent from " << z.squaredNorm() << "." << std::endl;
+      zSaddle = findSaddle(J, z, ε);
+    } catch (std::exception& e) {
+    }
+    std::cerr << "Current saddle is of size " << zSaddle.squaredNorm() << "." << std::endl;
+  }
+
+  ComplexVector zSaddlePrev = ComplexVector::Zero(N);
+
   while (δ < (zSaddle - zSaddlePrev).norm()) { // Until we find two saddles sufficiently close...
-    std::tie(std::ignore, z) = langevin(J, z, T, γ, M, d, r.engine());
+    std::tie(std::ignore, z) = metropolis(J, z, energyNormGrad, T, γ, M, d, r.engine());
     try {
       ComplexVector zSaddleNext = findSaddle(J, z, ε);
       if (Δ < (zSaddleNext - zSaddle).norm()) { // Ensure we are finding distinct saddles.
diff --git a/p-spin.hpp b/p-spin.hpp
index 4621db6..e5d4f94 100644
--- a/p-spin.hpp
+++ b/p-spin.hpp
@@ -5,15 +5,6 @@
 #include "tensor.hpp"
 #include "factorial.hpp"
 
-template <class Scalar>
-using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
-
-template <class Scalar>
-using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
-
-template <class Scalar, int p>
-using Tensor = Eigen::Tensor<Scalar, p>;
-
 template <class Scalar, int p>
 std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>> hamGradHess(const Tensor<Scalar, p>& J, const Vector<Scalar>& z) {
   Matrix<Scalar> Jz = contractDown(J, z); // Contracts J into p - 2 copies of z.
diff --git a/tensor.hpp b/tensor.hpp
index 21f2a89..6a5b2c2 100644
--- a/tensor.hpp
+++ b/tensor.hpp
@@ -5,43 +5,52 @@
 
 #include <eigen3/unsupported/Eigen/CXX11/Tensor>
 
+template <class Scalar>
+using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
+
+template <class Scalar>
+using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
+
+template <class Scalar, int p>
+using Tensor = Eigen::Tensor<Scalar, p>;
+
 template <class Scalar, int p, std::size_t... Indices>
-Eigen::Tensor<Scalar, p> initializeJHelper(unsigned N, std::index_sequence<Indices...>) {
+Tensor<Scalar, p> initializeJHelper(unsigned N, std::index_sequence<Indices...>) {
   std::array<unsigned, p> Ns;
   std::fill_n(Ns.begin(), p, N);
-  return Eigen::Tensor<Scalar, p>(std::get<Indices>(Ns)...);
+  return Tensor<Scalar, p>(std::get<Indices>(Ns)...);
 }
 
 template <class Scalar, int p>
-Eigen::Tensor<Scalar, p> initializeJ(unsigned N) {
+Tensor<Scalar, p> initializeJ(unsigned N) {
   return initializeJHelper<Scalar, p>(N, std::make_index_sequence<p>());
 }
 
 template <class Scalar, int p, std::size_t... Indices>
-void setJHelper(Eigen::Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind, Scalar val, std::index_sequence<Indices...>) {
+void setJHelper(Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind, Scalar val, std::index_sequence<Indices...>) {
   J(std::get<Indices>(ind)...) = val;
 }
 
 template <class Scalar, int p>
-void setJ(Eigen::Tensor<Scalar, p>& J, std::array<unsigned, p> ind, Scalar val) {
+void setJ(Tensor<Scalar, p>& J, std::array<unsigned, p> ind, Scalar val) {
   do {
     setJHelper<Scalar, p>(J, ind, val, std::make_index_sequence<p>());
   } while (std::next_permutation(ind.begin(), ind.end()));
 }
 
 template <class Scalar, int p, std::size_t... Indices>
-Scalar getJHelper(const Eigen::Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind, std::index_sequence<Indices...>) {
+Scalar getJHelper(const Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind, std::index_sequence<Indices...>) {
   return J(std::get<Indices>(ind)...);
 }
 
 template <class Scalar, int p>
-Scalar getJ(const Eigen::Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind) {
+Scalar getJ(const Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind) {
   return getJHelper<Scalar, p>(J, ind, std::make_index_sequence<p>());
 }
 
 template <class Scalar, int p>
-void iterateOverHelper(Eigen::Tensor<Scalar, p>& J,
-    std::function<void(Eigen::Tensor<Scalar, p>&, std::array<unsigned, p>)>& f,
+void iterateOverHelper(Tensor<Scalar, p>& J,
+    std::function<void(Tensor<Scalar, p>&, std::array<unsigned, p>)>& f,
     unsigned l, std::array<unsigned, p> is) {
   if (l == 0) {
     f(J, is);
@@ -61,17 +70,17 @@ void iterateOverHelper(Eigen::Tensor<Scalar, p>& J,
 }
 
 template <class Scalar, int p>
-void iterateOver(Eigen::Tensor<Scalar, p>& J, std::function<void(Eigen::Tensor<Scalar, p>&, std::array<unsigned, p>)>& f) {
+void iterateOver(Tensor<Scalar, p>& J, std::function<void(Tensor<Scalar, p>&, std::array<unsigned, p>)>& f) {
   std::array<unsigned, p> is;
   iterateOverHelper<Scalar, p>(J, f, p, is);
 }
 
 template <class Scalar, int p, class Distribution, class Generator>
-Eigen::Tensor<Scalar, p> generateCouplings(unsigned N, Distribution d, Generator& r) {
-  Eigen::Tensor<Scalar, p> J = initializeJ<Scalar, p>(N);
+Tensor<Scalar, p> generateCouplings(unsigned N, Distribution d, Generator& r) {
+  Tensor<Scalar, p> J = initializeJ<Scalar, p>(N);
 
-  std::function<void(Eigen::Tensor<Scalar, p>&, std::array<unsigned, p>)> setRandom =
-    [&d, &r] (Eigen::Tensor<Scalar, p>& JJ, std::array<unsigned, p> ind) {
+  std::function<void(Tensor<Scalar, p>&, std::array<unsigned, p>)> setRandom =
+    [&d, &r] (Tensor<Scalar, p>& JJ, std::array<unsigned, p> ind) {
       setJ<Scalar, p>(JJ, ind, d(r));
     };
 
@@ -81,17 +90,15 @@ Eigen::Tensor<Scalar, p> generateCouplings(unsigned N, Distribution d, Generator
 }
 
 template <class Scalar>
-Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>
-contractDown(const Eigen::Tensor<Scalar, 2>& J, const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>& z) {
-  return Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>>(J.data(), z.size(), z.size());
+Matrix<Scalar> contractDown(const Tensor<Scalar, 2>& J, const Vector<Scalar>& z) {
+  return Eigen::Map<const Matrix<Scalar>>(J.data(), z.size(), z.size());
 }
 
 const std::array<Eigen::IndexPair<int>, 1> ip00 = {Eigen::IndexPair<int>(0, 0)};
 
 template <class Scalar, int r>
-Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>
-contractDown(const Eigen::Tensor<Scalar, r>& J, const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>& z) {
-  Eigen::Tensor<Scalar, 1> zT = Eigen::TensorMap<Eigen::Tensor<const Scalar, 1>>(z.data(), {z.size()});
-  Eigen::Tensor<Scalar, r - 1> Jz = J.contract(zT, ip00);
+Matrix<Scalar> contractDown(const Tensor<Scalar, r>& J, const Vector<Scalar>& z) {
+  Tensor<Scalar, 1> zT = Eigen::TensorMap<Tensor<const Scalar, 1>>(z.data(), {z.size()});
+  Tensor<Scalar, r - 1> Jz = J.contract(zT, ip00);
   return contractDown(Jz, z);
 }