From 1d30476d67b1dcfaa4240dc8569a05c6caa84fed Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 4 Apr 2023 12:31:06 +0200
Subject: Added utilities for the study of real Morse lines.

Made real analogues of the library collectStokesData.hpp and
pureStokes.cpp (now titled collectMorseData.hpp and pureMorse.cpp) for
the study of Morse lines on ordinary (real) p-spin models.
---
 collectMorseData.hpp | 108 +++++++++++++++++++++++++++++++++++++++++++++++++++
 pureMorse.cpp        |  56 ++++++++++++++++++++++++++
 2 files changed, 164 insertions(+)
 create mode 100644 collectMorseData.hpp
 create mode 100644 pureMorse.cpp

diff --git a/collectMorseData.hpp b/collectMorseData.hpp
new file mode 100644
index 0000000..d63050d
--- /dev/null
+++ b/collectMorseData.hpp
@@ -0,0 +1,108 @@
+#include "stokes.hpp"
+#include "complex_normal.hpp"
+#include <fstream>
+
+template<int ...ps, class Generator, typename... T>
+void collectMorseData(std::string tag, unsigned N, Generator& r, Real ε, Real δz₀, bool minimum, bool spike, T... μs) {
+  unsigned nIts = 5;
+  unsigned nGs = 4;
+  unsigned coDim = 16;
+  Real newSaddleThres = 1e-4;
+
+  pSpinModel<Real, ps...> ReM(N, r, μs...);
+
+  std::normal_distribution<Real> Red(0, 1);
+
+  Vector<Real> xMin(N);
+
+  if (minimum) {
+    xMin = randomMinimum<Real>(ReM, Red, r, ε);
+  } else {
+    xMin = randomSaddle<Real, Real>(ReM, Red, r, ε);
+  }
+
+  Real Hx;
+  std::tie(Hx, std::ignore, std::ignore, std::ignore) = ReM.hamGradHess(xMin);
+
+  Vector<Real> dHx;
+  Matrix<Real> ddHx;
+  std::tie(Hx, dHx, ddHx, std::ignore) = ReM.hamGradHess(xMin);
+  Matrix<Real> xHess = ddHx - xMin.dot(dHx) * Matrix<Real>::Identity(xMin.size(), xMin.size()) / (Real)N;
+  xHess -= (xHess * xMin) * xMin.transpose() / (Real)N;
+  Eigen::EigenSolver<Matrix<Real>> eigenS(xHess);
+
+  pSpinModel M = ReM;
+  Vector<Real> zMin = xMin;
+
+  Vector<Real> zSaddle = zMin;
+
+  Real δz = δz₀;
+
+  while ((zSaddle - zMin).norm() < newSaddleThres * N) {
+    Vector<Real> z0 = normalize(zSaddle + δz * randomVector<Real>(N, Red, r));
+    zSaddle = findSaddle(M, z0, ε);
+    δz *= 1.01;
+  }
+
+  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");
+
+  file.precision(15);
+
+  file << N << std::endl;
+  ((file << ps << " "), ...) << std::endl;;
+  ((file << μs << " "), ...) << std::endl;;
+
+  std::ofstream tensorFile("stokes_tensor_" + tag + ".dat",  std::ios::out | std::ios::binary | std::ios::trunc);
+  std::apply([&tensorFile](const Tensor<Real, ps>&... Js) -> void {
+      (tensorFile.write((char*)Js.data(), Js.size() * sizeof(Real)), ...);
+      } , ReM.Js);
+
+  file << xMin.transpose() << std::endl;
+  file << Hx << std::endl;
+  file << eigenS.eigenvalues().real().transpose() << std::endl;
+  file << zSaddle.transpose() << std::endl;
+  file << Hz << " " << zSaddle.squaredNorm() << std::endl;
+  file << eigenSz.eigenvalues().transpose() << std::endl;
+  file << 0 << " " << (xMin - zSaddle).norm() << std::endl;
+
+  unsigned nTest = 1000;
+
+  file << nIts << " " << nGs << " " << coDim << " " << nTest << std::endl;
+
+  std::vector<Vector<Real>> oldGs;
+
+  for (int j = 0; j < nIts; j++) {
+    Cord c(M, zMin, zSaddle, nGs + pow(2, j));
+    for (unsigned i = 0; i < oldGs.size(); i++) {
+      c.gs[i] = oldGs[i];
+    }
+    c.relaxNewton(c.gs.size() * 2, 1, 1e-10, 1e4);
+    oldGs = c.gs;
+
+    Real reConstraintError = 0;
+    Real imConstraintError = 0;
+    Real imEnergyError = 0;
+    Real ell = 0;
+
+    for (unsigned i = 0; i < nTest; i++) {
+      Real t = (i + 1.0) / (nTest + 1.0);
+      Vector<Real> zi = c.f(t);
+      Real constraintError = (Real)(zi.transpose() * zi) - (Real)N;
+      reConstraintError += pow(constraintError, 2);
+      ell += c.df(t).norm();
+    }
+
+    file << oldGs.size() << " " << ell / nTest << " " <<  c.totalCost(nTest) / (nTest) << " " << sqrt(imEnergyError / (Real)(nTest)) << " " << sqrt(reConstraintError / (Real)(nTest)) << " " << sqrt(imConstraintError / (Real)(nTest)) << std::endl;
+    for (const Vector<Real>& gi : c.gs) {
+      file << gi.transpose() << std::endl;
+    }
+  }
+}
diff --git a/pureMorse.cpp b/pureMorse.cpp
new file mode 100644
index 0000000..f157c35
--- /dev/null
+++ b/pureMorse.cpp
@@ -0,0 +1,56 @@
+#include <getopt.h>
+#include <chrono>
+
+#include "collectMorseData.hpp"
+
+#include "pcg-cpp/include/pcg_random.hpp"
+#include "randutils/randutils.hpp"
+
+using Rng = randutils::random_generator<pcg32>;
+
+int main(int argc, char* argv[]) {
+  // model parameters
+  unsigned N = 10; // number of spins
+  // simulation parameters
+  Real ε = 1e-15;
+  Real δz₀ = 1e-3;
+  unsigned n = 10;
+  bool useMinima = false;
+  bool useSpike = false;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "N:e:d:n:ms")) != -1) {
+    switch (opt) {
+    case 'N':
+      N = (unsigned)atof(optarg);
+      break;
+    case 'e':
+      ε = atof(optarg);
+      break;
+    case 'd':
+      δz₀ = atof(optarg);
+      break;
+    case 'n':
+      n = atof(optarg);
+      break;
+    case 'm':
+      useMinima = true;
+      break;
+    case 's':
+      useSpike = true;
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  Rng r;
+
+  for (unsigned i = 0; i < n; i++) {
+    auto tag = std::chrono::high_resolution_clock::now();
+    collectMorseData<3>(std::to_string(tag.time_since_epoch().count()), N, r.engine(), ε, δz₀, useMinima, useSpike, 1.0);
+  }
+
+  return 0;
+}
-- 
cgit v1.2.3-54-g00ecf