Cleaned up code and algorithmic improvements for data collection.

5 files changed, 141 insertions, 189 deletions

diff --git a/collectStokesData.hpp b/collectStokesData.hpp
index 50aad24..433f5ad 100644
--- a/collectStokesData.hpp
+++ b/collectStokesData.hpp
@@ -1,21 +1,34 @@
 #include "stokes.hpp"
-#include <iostream>
+#include "complex_normal.hpp"
+#include <fstream>
 
 using Complex = std::complex<Real>;
 
 template<int ...ps, class Generator, typename... T>
-void collectStokesData(unsigned N, Generator& r, double ε, T... μs) {
+void collectStokesData(std::ofstream& file, unsigned N, Generator& r, double ε, double δz, bool minimum, T... μs) {
+  unsigned nGs = 8;
+  unsigned nTs = 20;
+  Real newSaddleThres = 1e-3;
+
   pSpinModel<Real, ps...> ReM(N, r, μs...);
 
   std::normal_distribution<Real> Red(0, 1);
 
-  Vector<Real> xMin = randomMinimum(ReM, Red, r, ε);
+  Vector<Real> xMin(N);
+
+  if (minimum) {
+    xMin = randomMinimum<Real>(ReM, Red, r, ε);
+  } else {
+    xMin = randomSaddle<Real, Real>(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));
+  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);
 
   complex_normal_distribution<Real> d(0, 1, 0);
 
@@ -24,29 +37,56 @@ void collectStokesData(unsigned N, Generator& r, double ε, T... μs) {
 
   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));
+  while ((zSaddle - zMin).norm() < newSaddleThres * N || abs(imag(M.getHamiltonian(zSaddle))) < 1e-10) {
+    Vector<Complex> z0 = normalize(zSaddle + δz * randomVector<Complex>(N, d, r));
     zSaddle = findSaddle(M, z0, ε);
   }
 
-  Complex H1 = M.getHamiltonian(zMin);
-  Complex H2 = M.getHamiltonian(zSaddle);
+  Complex Hz;
+  Vector<Complex> dHz;
+  Matrix<Complex> ddHz;
+  std::tie(Hz, dHz, ddHz, std::ignore) = M.hamGradHess(zSaddle);
+  Matrix<Complex> zHess = ddHz - (zSaddle.transpose() * dHz) * Matrix<Complex>::Identity(N, N) / (Real)N;
+  zHess -= (zHess * zSaddle) * zSaddle.adjoint() / zSaddle.squaredNorm();
+  Eigen::ComplexEigenSolver<Matrix<Complex>> eigenSz(zHess);
 
-  Real φ = atan2( H2.imag() - H1.imag(), H1.real() - H2.real());
+  Real φ = atan2(Hz.imag(), Hx - Hz.real());
 
   M *= exp(Complex(0, φ));
 
-  Cord c(M, zMin, zSaddle, 3);
-  c.relaxNewton(10, 1, 1e4);
+  file.precision(15);
 
-  ((std::cout << ps), ...) << std::endl;;
-  ((std::cout << μs), ...) << std::endl;;
+  file << N << std::endl;
+  ((file << ps), ...) << std::endl;;
+  ((file << μs), ...) << std::endl;;
 
-  std::apply([](const Tensor<Real, ps>&... Js) -> void {
-        ((std::cout << Js << std::endl), ...);
+  std::apply([&file](const Tensor<Real, ps>&... Js) -> void {
+        ((file << Js << std::endl), ...);
       } , ReM.Js);
 
-  std::cout << xMin.transpose() << std::endl;
-  std::cout << Hx << std::endl;
-  std::cout << eigenS.eigenvalues().real().transpose() << std::endl;
+  file << xMin.transpose() << std::endl;
+  file << Hx << std::endl;
+  file << eigenS.eigenvalues().real().transpose() << std::endl;
+  file << zSaddle.transpose() << std::endl;
+  file << Hz << std::endl;
+  file << eigenSz.eigenvalues().transpose() << std::endl;
+  file << φ << std::endl;
+
+  Cord c(M, zMin, zSaddle, nGs);
+  c.relaxNewton(nTs, 1, 1e-10, 1e3);
+
+  Complex constraintError = 0;
+  Real imEnergyError = 0;
+
+  for (unsigned i = 0; i < 100 * nTs; i++) {
+    Vector<Complex> zi = c.f((i + 1.0) / (100.0 * nTs + 1.0));
+    imEnergyError += pow(std::imag(M.getHamiltonian(zi) - M.getHamiltonian(zMin)), 2);
+    constraintError += pow(((Complex)(zi.transpose() * zi) - (Complex)N), 2);
+  }
+
+  file << nGs << " " << nTs << std::endl;;
+  file << c.totalCost(100 * nTs) / (100 * nTs) << " " << sqrt(imEnergyError) / (100 * nTs) << " " << sqrt(constraintError) / (100.0 * nTs) << std::endl;
+  for (const Vector<Complex>& gi : c.gs) {
+    file << gi.transpose() << std::endl;
+  }
 }
diff --git a/dynamics.hpp b/dynamics.hpp
index 4f33c1a..e578cdb 100644
--- a/dynamics.hpp
+++ b/dynamics.hpp
@@ -108,7 +108,7 @@ Vector<Scalar> randomSaddle(const pSpinModel<Real, ps...>& M, Distribution d, Ge
   bool foundSaddle = false;
 
   while (!foundSaddle) {
-    Vector<Scalar> z0 = normalize(randomVector<Scalar>(M.dimension(), d, r.engine()));
+    Vector<Scalar> z0 = normalize(randomVector<Scalar>(M.dimension(), d, r));
 
     try {
       zSaddle = findSaddle(M, z0, ε);
diff --git a/langevin.cpp b/langevin.cpp
deleted file mode 100644
index ea28d65..0000000
--- a/langevin.cpp
+++ /dev/null
@@ -1,142 +0,0 @@
-#include <getopt.h>
-#include <limits>
-#include <unordered_map>
-#include <list>
-
-#include "Eigen/Dense"
-#include "Eigen/src/Eigenvalues/ComplexEigenSolver.h"
-#include "Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h"
-#include "complex_normal.hpp"
-#include "p-spin.hpp"
-#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"
-
-#define PSPIN_P 3
-const int p = PSPIN_P; // polynomial degree of Hamiltonian
-using Complex = std::complex<Real>;
-using RealVector = Vector<Real>;
-using RealMatrix = Matrix<Real>;
-using RealTensor = Tensor<Real, p>;
-using ComplexVector = Vector<Complex>;
-using ComplexMatrix = Matrix<Complex>;
-using ComplexTensor = Tensor<Complex, p>;
-
-using Rng = randutils::random_generator<pcg32>;
-
-int main(int argc, char* argv[]) {
-  // model parameters
-  unsigned N = 10; // number of spins
-  Real T = 1;    // temperature
-  Real Rκ = 0;   // real part of distribution parameter
-  Real Iκ = 0;   // imaginary part of distribution parameter
-  // simulation parameters
-  Real ε = 1e-15;
-  Real εJ = 1e-5;
-  Real δ = 1;   // threshold for determining saddle
-  Real Δ = 1e-3;
-  Real γ = 1e-1;   // step size
-  unsigned t = 1000; // number of Langevin steps
-  unsigned M = 100;
-  unsigned n = 100;
-
-  int opt;
-
-  while ((opt = getopt(argc, argv, "N:M:n:T:e:r:i:g:t:E:")) != -1) {
-    switch (opt) {
-    case 'N':
-      N = (unsigned)atof(optarg);
-      break;
-    case 't':
-      t = (unsigned)atof(optarg);
-      break;
-    case 'T':
-      T = atof(optarg);
-      break;
-    case 'e':
-      δ = atof(optarg);
-      break;
-    case 'E':
-      ε = atof(optarg);
-      break;
-    case 'g':
-      γ = atof(optarg);
-    case 'r':
-      Rκ = atof(optarg);
-      break;
-    case 'i':
-      Iκ = atof(optarg);
-      break;
-    case 'n':
-      n = (unsigned)atof(optarg);
-      break;
-    case 'M':
-      M = (unsigned)atof(optarg);
-      break;
-    default:
-      exit(1);
-    }
-  }
-
-  Rng r;
-
-  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.engine(), ε);
-  Real Hr;
-  RealVector dHr;
-  RealMatrix ddHr;
-  std::tie(Hr, dHr, ddHr, std::ignore) = pSpin.hamGradHess(zMin);
-  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-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();
-
-  complex_normal_distribution<Real> d(0, 1, 0);
-
-  pSpinModel complexPSpin = pSpin.cast<Complex>();;
-  ComplexVector zSaddle = zMin.cast<Complex>();
-
-  ComplexVector zSaddleNext;
-  bool foundSaddle = false;
-  while (!foundSaddle) {
-    ComplexVector z0 = normalize(zSaddle + δ * randomVector<Complex>(N, d, r.engine()));
-    try {
-      zSaddleNext = findSaddle(complexPSpin, z0, ε);
-      Real saddleDistance = (zSaddleNext - zSaddle).norm();
-      if (saddleDistance / N > 1e-2) {
-        std::cout << saddleDistance << std::endl;
-        foundSaddle = true;
-      }
-    } catch (std::exception& e) {}
-  }
-
-  Complex H1 = complexPSpin.getHamiltonian(zSaddle);
-  Complex H2 = complexPSpin.getHamiltonian(zSaddleNext);
-
-  Real φ = atan2( H2.imag() - H1.imag(), H1.real() - H2.real());
-  std::cerr << (zSaddle - zSaddleNext).norm() / (Real)N << " " << φ << " " << H1 * exp(Complex(0, φ)) << " " << H2 * exp(Complex(0, φ)) << std::endl;
-
-  complexPSpin *= exp(Complex(0, φ));
-
-  Cord test(complexPSpin, zSaddle, zSaddleNext, 3);
-  test.relaxNewton(10, 1, 1e4);
-
-  std::cout << test.totalCost(10) << " " << test.totalCost(20) << std::endl;
-
-  return 0;
-}
diff --git a/pureStokesFromMinima.cpp b/pureStokesFromMinima.cpp
new file mode 100644
index 0000000..8f815cf
--- /dev/null
+++ b/pureStokesFromMinima.cpp
@@ -0,0 +1,54 @@
+#include <getopt.h>
+#include <chrono>
+#include <fstream>
+
+#include "collectStokesData.hpp"
+
+#include "pcg-cpp/include/pcg_random.hpp"
+#include "randutils/randutils.hpp"
+#include "unsupported/Eigen/CXX11/src/Tensor/TensorFFT.h"
+
+#define PSPIN_P 3
+const int p = PSPIN_P; // polynomial degree of Hamiltonian
+
+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 δ = 1;
+  unsigned n = 10;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "N:e:d:n:")) != -1) {
+    switch (opt) {
+    case 'N':
+      N = (unsigned)atof(optarg);
+      break;
+    case 'e':
+      ε = atof(optarg);
+      break;
+    case 'd':
+      δ = atof(optarg);
+      break;
+    case 'n':
+      n = atof(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  Rng r;
+
+  for (unsigned i = 0; i < n; i++) {
+    auto tag = std::chrono::high_resolution_clock::now();
+    std::ofstream output("stokes_" + std::to_string(tag.time_since_epoch().count()) + ".dat");
+    collectStokesData<3>(output, N, r.engine(), ε, δ, true, 1.0);
+  }
+
+  return 0;
+}
diff --git a/stokes.hpp b/stokes.hpp
index 9fd6f01..8c0c1fb 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -1,9 +1,10 @@
 #pragma once
 
 #include "p-spin.hpp"
-#include "complex_normal.hpp"
 #include "dynamics.hpp"
 
+#include <iostream>
+
 template <class Scalar, int ...ps>
 class Cord {
 public:
@@ -17,7 +18,7 @@ public:
     Scalar H2 = M.getHamiltonian(z2);
     Scalar H3 = M.getHamiltonian(z3);
 
-    if (real(H2) > real(H3)) {
+    if (std::real(H2) > std::real(H3)) {
       z0 = z2;
       z1 = z3;
     } else {
@@ -62,7 +63,7 @@ public:
     Vector<Scalar> ż = zDot(z, dH);
     Vector<Scalar> dz = df(t);
 
-    return 1 - real(ż.dot(dz)) / ż.norm() / dz.norm();
+    return 1 - std::real(ż.dot(dz)) / ż.norm() / dz.norm();
   }
 
   Real totalCost(unsigned nt) const {
@@ -117,24 +118,22 @@ public:
       Real dfdgn = dgCoeff(i, t);
       Vector<Scalar> dC = - 0.5 / ż.norm() / dz.norm() * (
           dfdgn * ż.conjugate() + fdgn * dżc * dz + fdgn * dż * dz.conjugate()
-          - real(dz.dot(ż)) * (
+          - std::real(dz.dot(ż)) * (
               dfdgn * dz.conjugate() / dz.squaredNorm() +
               fdgn * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm()
             )
         );
 
-      for (unsigned j = 0; j < dC.size(); j++) {
-        x(z.size() * i + j) = dC(j);
-        x(z.size() * gs.size() + z.size() * i + j) = conj(dC(j));
-      }
+      x.segment(z.size() * i, z.size()) = dC;
+      x.segment(z.size() * gs.size() + z.size() * i, z.size()) = dC.conjugate();
 
       for (unsigned j = 0; j < gs.size(); j++) {
         Real fdgm = gCoeff(j, t);
         Real dfdgm = dgCoeff(j, t);
-        Scalar CC = - real(dz.dot(ż)) / ż.norm() / dz.norm();
+        Scalar CC = - std::real(dz.dot(ż)) / ż.norm() / dz.norm();
         Vector<Scalar> dCm = - 0.5 / ż.norm() / dz.norm() * (
             dfdgm * ż.conjugate() + fdgm * dżc * dz + fdgm * dż * dz.conjugate()
-            - real(dz.dot(ż)) * (
+            - std::real(dz.dot(ż)) * (
                 dfdgm * dz.conjugate() / dz.squaredNorm() +
                 fdgm * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm()
               )
@@ -195,21 +194,17 @@ public:
              )
             ;
 
-        for (unsigned k = 0; k < z.size(); k++) {
-          for (unsigned l = 0; l < z.size(); l++) {
-            M(z.size() * i + k, z.size() * j + l) = dcdC(k, l);
-            M(z.size() * gs.size() + z.size() * i + k, z.size() * j + l) = conj(ddC(k, l));
-            M(z.size() * i + k, z.size() * gs.size() + z.size() * j + l) = ddC(k, l);
-            M(z.size() * gs.size() + z.size() * i + k, z.size() * gs.size() + z.size() * j + l) = conj(dcdC(k, l));
-          }
-        }
+        M.block(z.size() * i, z.size() * j, z.size(), z.size()) = dcdC;
+        M.block(z.size() * gs.size() + z.size() * i, z.size() * j, z.size(), z.size()) = ddC.conjugate();
+        M.block(z.size() * i, z.size() * gs.size() + z.size() * j, z.size(), z.size()) = ddC;
+        M.block(z.size() * gs.size() + z.size() * i, z.size() * gs.size() + z.size() * j, z.size(), z.size()) = dcdC.conjugate();
       }
     }
 
     return {x, M};
   }
 
-  std::pair<Real, Real> relaxStepNewton(unsigned nt, Real δ₀) {
+  std::tuple<Real, Real> relaxStepNewton(unsigned nt, Real δ₀) {
     Vector<Scalar> dgsTot = Vector<Scalar>::Zero(2 * z0.size() * gs.size());
     Matrix<Scalar> HessTot = Matrix<Scalar>::Zero(2 * z0.size() * gs.size(), 2 * z0.size() * gs.size());
 
@@ -231,32 +226,37 @@ public:
 
     Matrix<Scalar> I = abs(HessTot.diagonal().array()).matrix().asDiagonal();
 
+    Vector<Scalar> δg(gs.size() * z0.size());
+
     while (newCost > oldCost) {
-      Vector<Scalar> δg = (HessTot + δ * I).partialPivLu().solve(dgsTot);
+      δg = (HessTot + δ * I).partialPivLu().solve(dgsTot).tail(gs.size() * z0.size());
 
       for (unsigned i = 0; i < gs.size(); i++) {
-        for (unsigned j = 0; j < z0.size(); j++) {
-          cNew.gs[i](j) = gs[i](j) - δg[z0.size() * gs.size() + z0.size() * i + j];
-        }
+        cNew.gs[i] = gs[i] - δg.segment(z0.size() * i, z0.size());
       }
 
       newCost = cNew.totalCost(nt);
 
-      δ *= 2;
+      δ *= 1.5;
     }
 
+    std::cerr << newCost << " " << stepSize << " " << δ << std::endl;
+
     gs = cNew.gs;
 
-    return {δ / 2, stepSize};
+    return {δ / 1.5, stepSize};
   }
 
-  void relaxNewton(unsigned nt, Real δ₀, unsigned maxSteps) {
+  void relaxNewton(unsigned nt, Real δ₀, Real ε, unsigned maxSteps) {
     Real δ = δ₀;
     Real stepSize = std::numeric_limits<Real>::infinity();
     unsigned steps = 0;
-    while (stepSize > 1e-7 && steps < maxSteps) {
+    while (stepSize > ε && steps < maxSteps) {
       std::tie(δ, stepSize) = relaxStepNewton(nt, δ);
-      δ /= 3;
+      if (δ > 100) {
+        nt++;
+      }
+      δ /= 2;
       steps++;
     }
   }