2 files changed, 164 insertions, 0 deletions
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;
+}