From 7cc214d81264722e5766885fefa477efd4311e3c Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 4 Apr 2023 18:16:43 +0200
Subject: Collect only data on stationary points with adjacent indices.

---
 collectMorseData.hpp | 32 +++++++++++++++++++++++++-------
 1 file changed, 25 insertions(+), 7 deletions(-)

diff --git a/collectMorseData.hpp b/collectMorseData.hpp
index d63050d..eddb825 100644
--- a/collectMorseData.hpp
+++ b/collectMorseData.hpp
@@ -31,6 +31,12 @@ void collectMorseData(std::string tag, unsigned N, Generator& r, Real ε, Real 
   xHess -= (xHess * xMin) * xMin.transpose() / (Real)N;
   Eigen::EigenSolver<Matrix<Real>> eigenS(xHess);
 
+  unsigned indx = 0;
+
+  for (Real λ : eigenS.eigenvalues().real()) {
+    if (λ < 0) indx++;
+  }
+
   pSpinModel M = ReM;
   Vector<Real> zMin = xMin;
 
@@ -38,18 +44,30 @@ void collectMorseData(std::string tag, unsigned N, Generator& r, Real ε, Real 
 
   Real δz = δz₀;
 
-  while ((zSaddle - zMin).norm() < newSaddleThres * N) {
+  Real Hz = -std::numeric_limits<Real>::infinity();
+  Vector<Real> dHz;
+  Matrix<Real> ddHz;
+  Matrix<Real> zHess;
+
+  unsigned indz = std::numeric_limits<unsigned>::infinity();
+
+  while ((zSaddle - zMin).norm() < newSaddleThres * N || !(indz == indx + 1 || indz == indx - 1) ) {
     Vector<Real> z0 = normalize(zSaddle + δz * randomVector<Real>(N, Red, r));
     zSaddle = findSaddle(M, z0, ε);
     δz *= 1.01;
+
+    std::tie(Hz, dHz, ddHz, std::ignore) = M.hamGradHess(zSaddle);
+    zHess = ddHz - (zSaddle.transpose() * dHz) * Matrix<Real>::Identity(N, N) / (Real)N;
+    zHess -= (zHess * zSaddle) * zSaddle.adjoint() / zSaddle.squaredNorm();
+    Eigen::EigenSolver<Matrix<Real>> eigenSz(zHess);
+
+    indz = 0;
+
+    for (Real λ : eigenSz.eigenvalues().real()) {
+      if (λ < 0) indz++;
+    }
   }
 
-  Real Hz;
-  Vector<Real> dHz;
-  Matrix<Real> ddHz;
-  std::tie(Hz, dHz, ddHz, std::ignore) = M.hamGradHess(zSaddle);
-  Matrix<Real> zHess = ddHz - (zSaddle.transpose() * dHz) * Matrix<Real>::Identity(N, N) / (Real)N;
-  zHess -= (zHess * zSaddle) * zSaddle.adjoint() / zSaddle.squaredNorm();
   Eigen::EigenSolver<Matrix<Real>> eigenSz(zHess);
 
   std::ofstream file("stokes_info_" + tag + ".dat");
-- 
cgit v1.2.3-54-g00ecf