1 files changed, 36 insertions, 13 deletions

diff --git a/langevin.cpp b/langevin.cpp
index cf61b85..acac6c7 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -1,6 +1,7 @@
 #include <getopt.h>
 
 #include <eigen3/Eigen/LU>
+#include <utility>
 
 #include "complex_normal.hpp"
 #include "p-spin.hpp"
@@ -8,6 +9,7 @@
 
 #include "pcg-cpp/include/pcg_random.hpp"
 #include "randutils/randutils.hpp"
+#include "unsupported/Eigen/CXX11/src/Tensor/TensorFFT.h"
 
 using Rng = randutils::random_generator<pcg32>;
 
@@ -129,7 +131,9 @@ int main(int argc, char* argv[]) {
 
   // simulation parameters
   double ε = 1e-4;
+  double εJ = 5e-2;
   double δ = 1e-2;   // threshold for determining saddle
+  double Δ = 1e-3;
   double γ = 1e-2;   // step size
   unsigned t = 1000; // number of Langevin steps
   unsigned M = 100;
@@ -182,24 +186,43 @@ int main(int argc, char* argv[]) {
   Vector z0 = normalize(randomVector(N, complex_normal_distribution<>(0, 1, 0), r.engine()));
 
   Vector zSaddle = findSaddle(J, z0, ε);
-
-  double W;
+  Vector zSaddlePrev = Vector::Zero(N);
   Vector z = zSaddle;
-  for (unsigned i = 0; i < n; i++) {
-    std::tie(W, z) = langevin(J, z, T, γ, M, r);
+
+  while (δ < (zSaddle - zSaddlePrev).norm()) { // Until we find two saddles sufficiently close...
+    std::tie(std::ignore, z) = langevin(J, z, T, γ, M, r);
     try {
-      Vector zNewSaddle = findSaddle(J, z, ε);
-      Scalar H;
-      Matrix ddH;
-      std::tie(H, std::ignore, ddH) = hamGradHess(J, zNewSaddle);
-      Eigen::SelfAdjointEigenSolver<Matrix> es(ddH.adjoint() * ddH);
-      std::cout << N << "\t" << M * (i+1) << "\t" << H << "\t"
-        << zNewSaddle.transpose() << "\t" << es.eigenvalues().transpose()
-        << std::endl;
-    std::cerr << M * (i+1) << " steps taken to move " << (zNewSaddle - zSaddle).norm() << ", saddle information saved." << std::endl;
+      Vector zSaddleNext = findSaddle(J, z, ε);
+      if (Δ < (zSaddleNext - zSaddle).norm()) { // Ensure we are finding distinct saddles.
+        zSaddlePrev = zSaddle;
+        zSaddle = zSaddleNext;
+      }
     } catch (std::exception& e) {
       std::cerr << "Could not find a saddle: " << e.what() << std::endl;
     }
+
+    std::cerr << "Current saddles are " << (zSaddle - zSaddlePrev).norm() << " apart." << std::endl;
+  }
+
+  std::cerr << "Found sufficiently nearby saddles, perturbing J." << std::endl;
+
+  complex_normal_distribution<> dJ(0, εJ * σ, 0);
+
+  std::function<void(Tensor&, std::array<unsigned, p>)> perturbJ =
+    [&dJ, &r] (Tensor& JJ, std::array<unsigned, p> ind) {
+      Scalar Ji = getJ<Scalar, p>(JJ, ind);
+      setJ<Scalar, p>(JJ, ind, Ji + dJ(r.engine()));
+    };
+
+  for (unsigned i = 0; i < n; i++) {
+    Tensor Jp = J;
+
+    iterateOver<Scalar, p>(Jp, perturbJ);
+
+    Vector zSaddleNew = findSaddle(Jp, zSaddle, ε);
+    Vector zSaddlePrevNew = findSaddle(Jp, zSaddlePrev, ε);
+
+    std::cout << (zSaddleNew - zSaddlePrevNew).norm() << std::endl;
   }
 
   return 0;