From 58e62d50243bc45dbd4ae466642e801a1b46dd69 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 9 Nov 2021 17:59:53 +0100
Subject: Progress towards a systematized approach.

---
 collectStokesData.hpp | 52 +++++++++++++++++++++++++++++++++++++++++++++++++++
 dynamics.hpp          |  4 ++--
 langevin.cpp          | 17 +++++++++++------
 3 files changed, 65 insertions(+), 8 deletions(-)
 create mode 100644 collectStokesData.hpp

diff --git a/collectStokesData.hpp b/collectStokesData.hpp
new file mode 100644
index 0000000..50aad24
--- /dev/null
+++ b/collectStokesData.hpp
@@ -0,0 +1,52 @@
+#include "stokes.hpp"
+#include <iostream>
+
+using Complex = std::complex<Real>;
+
+template<int ...ps, class Generator, typename... T>
+void collectStokesData(unsigned N, Generator& r, double ε, T... μs) {
+  pSpinModel<Real, ps...> ReM(N, r, μs...);
+
+  std::normal_distribution<Real> Red(0, 1);
+
+  Vector<Real> xMin = randomMinimum(ReM, Red, r, ε);
+
+  Real Hx;
+  Vector<Real> dHx;
+  Matrix<Real> ddHx;
+  std::tie(Hx, dHx, ddHx, std::ignore) = ReM.hamGradHess(xMin);
+  Eigen::EigenSolver<Matrix<Real>> eigenS(ddHx - xMin.dot(dHx) * Matrix<Real>::Identity(xMin.size(), xMin.size()) / Real(N));
+
+  complex_normal_distribution<Real> d(0, 1, 0);
+
+  pSpinModel M = ReM.template cast<Complex>();;
+  Vector<Complex> zMin = xMin.cast<Complex>();
+
+  Vector<Complex> zSaddle = zMin;
+
+  while ((zSaddle - zMin).norm() < 1e-3 * N || abs(imag(M.getHamiltonian(zSaddle))) < 1e-10) {
+    Vector<Complex> z0 = normalize(zSaddle + 0.5 * randomVector<Complex>(N, d, r));
+    zSaddle = findSaddle(M, z0, ε);
+  }
+
+  Complex H1 = M.getHamiltonian(zMin);
+  Complex H2 = M.getHamiltonian(zSaddle);
+
+  Real φ = atan2( H2.imag() - H1.imag(), H1.real() - H2.real());
+
+  M *= exp(Complex(0, φ));
+
+  Cord c(M, zMin, zSaddle, 3);
+  c.relaxNewton(10, 1, 1e4);
+
+  ((std::cout << ps), ...) << std::endl;;
+  ((std::cout << μs), ...) << std::endl;;
+
+  std::apply([](const Tensor<Real, ps>&... Js) -> void {
+        ((std::cout << Js << std::endl), ...);
+      } , ReM.Js);
+
+  std::cout << xMin.transpose() << std::endl;
+  std::cout << Hx << std::endl;
+  std::cout << eigenS.eigenvalues().real().transpose() << std::endl;
+}
diff --git a/dynamics.hpp b/dynamics.hpp
index 36b2dfa..4f33c1a 100644
--- a/dynamics.hpp
+++ b/dynamics.hpp
@@ -9,7 +9,7 @@
 template <class Real, int ...ps>
 Vector<Real> findMinimum(const pSpinModel<Real, ps...>& M, const Vector<Real>& z0, Real ε) {
   Vector<Real> z = z0;
-  Real λ = 10;
+  Real λ = 100;
 
   Real H;
   Vector<Real> dH;
@@ -91,7 +91,7 @@ Vector<Real> randomMinimum(const pSpinModel<Real, ps...>& M, Distribution d, Gen
   bool foundSaddle = false;
 
   while (!foundSaddle) {
-    Vector<Real> z0 = normalize(randomVector<Real>(M.dimension(), d, r.engine()));
+    Vector<Real> z0 = normalize(randomVector<Real>(M.dimension(), d, r));
 
     try {
       zSaddle = findMinimum(M, z0, ε);
diff --git a/langevin.cpp b/langevin.cpp
index 0e5a616..ea28d65 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -11,6 +11,8 @@
 #include "dynamics.hpp"
 #include "stokes.hpp"
 
+#include "collectStokesData.hpp"
+
 #include "pcg-cpp/include/pcg_random.hpp"
 #include "randutils/randutils.hpp"
 #include "unsupported/Eigen/CXX11/src/Tensor/TensorFFT.h"
@@ -83,20 +85,23 @@ int main(int argc, char* argv[]) {
 
   Rng r;
 
-//  pSpinModel<Real, 2, 4> pSpin({J2, J4});
-  pSpinModel<Real, 2, 4>pSpin(N, r.engine(), 1, 0.1);
+  collectStokesData<2, 4>(N, r.engine(), 1e-15, 1.0, 0.01);
+
+  pSpinModel<Real, 2, 4>pSpin(N, r.engine(), 1, 0.01);
 
   std::normal_distribution<Real> Red(0, 1);
 
-  RealVector zMin = randomMinimum(pSpin, Red, r, ε);
+  RealVector zMin = randomMinimum(pSpin, Red, r.engine(), ε);
   Real Hr;
   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() + zMin * zMin.transpose() / (Real)N);
+  RealMatrix M1 = ddHr - (dHr * zMin.transpose() + zMin.dot(dHr) * Matrix<Real>::Identity(zMin.size(), zMin.size()) + (ddHr * zMin) * zMin.transpose()) / (Real)zMin.size();
+  RealMatrix M2 = ddHr - zMin.dot(dHr) * Matrix<Real>::Identity(zMin.size(), zMin.size()) / Real(N);
+  Eigen::EigenSolver<Matrix<Real>> eigenS(M2);
   for (unsigned i = 0; i < N; i++) {
-    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;
+    RealVector zNew = normalize(zMin + 1e-5 * eigenS.eigenvectors().col(i).real());
+    std::cout << eigenS.eigenvalues()(i) << " " << 2.0 * (pSpin.getHamiltonian(zNew) - Hr) / 1e-10 << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << " " << eigenS.eigenvectors().col(0).dot(eigenS.eigenvectors().col(i)) << std::endl;
   }
   std::cout << std::endl;
   getchar();
-- 
cgit v1.2.3-54-g00ecf