From f47dd594893995aa669de416f189676e7e97caba Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 23 Nov 2021 11:04:45 +0100
Subject: Trying new method for distinguishing Stoke's lines.

---
 collectStokesData.hpp | 52 ++++++++++++++++++++++++++++++---------------------
 1 file changed, 31 insertions(+), 21 deletions(-)

diff --git a/collectStokesData.hpp b/collectStokesData.hpp
index eae81fd..9492b25 100644
--- a/collectStokesData.hpp
+++ b/collectStokesData.hpp
@@ -6,8 +6,9 @@ using Complex = std::complex<Real>;
 
 template<int ...ps, class Generator, typename... T>
 void collectStokesData(std::string tag, unsigned N, Generator& r, double ε, double δz, bool minimum, T... μs) {
-  unsigned nGs = 8;
-  unsigned nTs = 32;
+  unsigned nIts = 5;
+  unsigned nGs = 4;
+  unsigned coDim = 16;
   Real newSaddleThres = 1e-4;
 
   pSpinModel<Real, ps...> ReM(N, r, μs...);
@@ -38,7 +39,7 @@ void collectStokesData(std::string tag, unsigned N, Generator& r, double ε, dou
   Vector<Complex> zSaddle = zMin;
 
   while ((zSaddle - zMin).norm() < newSaddleThres * N || abs(imag(M.getHamiltonian(zSaddle))) < 1e-10) {
-    Vector<Complex> z0 = normalize(zSaddle + δz * randomVector<Complex>(N, d, r));
+    Vector<Complex> z0 = normalize(zSaddle + N * δz * randomVector<Complex>(N, d, r));
     zSaddle = findSaddle(M, z0, ε);
   }
 
@@ -54,9 +55,6 @@ void collectStokesData(std::string tag, unsigned N, Generator& r, double ε, dou
 
   M *= exp(Complex(0, φ));
 
-  Cord c(M, zMin, zSaddle, nGs);
-  c.relaxNewton(nTs, 1, 1e-10, 1e3);
-
   std::ofstream file("stokes_info_" + tag + ".dat");
 
   file.precision(15);
@@ -78,23 +76,35 @@ void collectStokesData(std::string tag, unsigned N, Generator& r, double ε, dou
   file << eigenSz.eigenvalues().transpose() << std::endl;
   file << φ << " " << (xMin - zSaddle).norm() << std::endl;
 
-  Real reConstraintError = 0;
-  Real imConstraintError = 0;
-  Real imEnergyError = 0;
+  unsigned nTest = 1000;
 
-  unsigned nTest = nTs * 10;
+  file << nIts << " " << nGs << " " << coDim << " " << nTest << std::endl;
 
-  for (unsigned i = 0; i < nTest; i++) {
-    Vector<Complex> zi = c.f((i + 1.0) / (nTest + 1.0));
-    imEnergyError += pow(std::imag(M.getHamiltonian(zi) - M.getHamiltonian(zMin)), 2);
-    Complex constraintError = (Complex)(zi.transpose() * zi) - (Complex)N;
-    reConstraintError += pow(real(constraintError), 2);
-    imConstraintError += pow(imag(constraintError), 2);
-  }
+  std::vector<Vector<Complex>> oldGs;
+
+  for (int j = 0; j < nIts; j++) {
+    Cord c(M, zMin, zSaddle, nGs + j * 2);
+    for (unsigned i = 0; i < oldGs.size(); i++) {
+      c.gs[i] = oldGs[i];
+    }
+    c.relaxNewton(c.gs.size() * 4, 1, 1e-10, 1e4);
+    oldGs = c.gs;
+
+    Real reConstraintError = 0;
+    Real imConstraintError = 0;
+    Real imEnergyError = 0;
+
+    for (unsigned i = 0; i < nTest; i++) {
+      Vector<Complex> zi = c.f((i + 1.0) / (nTest + 1.0));
+      imEnergyError += pow(std::imag(M.getHamiltonian(zi) - M.getHamiltonian(zMin)), 2);
+      Complex constraintError = (Complex)(zi.transpose() * zi) - (Complex)N;
+      reConstraintError += pow(real(constraintError), 2);
+      imConstraintError += pow(imag(constraintError), 2);
+    }
 
-  file << nGs << " " << nTest << std::endl;;
-  file << c.totalCost(nTest) / (nTest) << " " << sqrt(imEnergyError / (Real)(nTest)) << " " << sqrt(reConstraintError / (Real)(nTest)) << " " << sqrt(imConstraintError / (Real)(nTest)) << std::endl;
-  for (const Vector<Complex>& gi : c.gs) {
-    file << gi.transpose() << std::endl;
+    file << oldGs.size() << " " <<  c.totalCost(nTest) / (nTest) << " " << sqrt(imEnergyError / (Real)(nTest)) << " " << sqrt(reConstraintError / (Real)(nTest)) << " " << sqrt(imConstraintError / (Real)(nTest)) << std::endl;
+    for (const Vector<Complex>& gi : c.gs) {
+      file << gi.transpose() << std::endl;
+    }
   }
 }
-- 
cgit v1.2.3-54-g00ecf