From c9cb8dd499c2b2648a78e63f2a0810451a69f8fd Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 12 Nov 2021 16:41:14 +0100
Subject: Some tweaking of the Stokes finder, how the lines are recorded, and
 fewer executables.

---
 collectStokesData.hpp | 37 ++++++++++++++++++++++---------------
 1 file changed, 22 insertions(+), 15 deletions(-)

(limited to 'collectStokesData.hpp')

diff --git a/collectStokesData.hpp b/collectStokesData.hpp
index 433f5ad..b32ccf0 100644
--- a/collectStokesData.hpp
+++ b/collectStokesData.hpp
@@ -5,10 +5,11 @@
 using Complex = std::complex<Real>;
 
 template<int ...ps, class Generator, typename... T>
-void collectStokesData(std::ofstream& file, unsigned N, Generator& r, double ε, double δz, bool minimum, T... μs) {
+void collectStokesData(std::string tag, unsigned N, Generator& r, double ε, double δz, bool minimum, T... μs) {
+  std::ofstream file("stokes_info_" + tag + ".dat");
   unsigned nGs = 8;
-  unsigned nTs = 20;
-  Real newSaddleThres = 1e-3;
+  unsigned nTs = 32;
+  Real newSaddleThres = 1e-4;
 
   pSpinModel<Real, ps...> ReM(N, r, μs...);
 
@@ -57,35 +58,41 @@ void collectStokesData(std::ofstream& file, unsigned N, Generator& r, double ε,
   file.precision(15);
 
   file << N << std::endl;
-  ((file << ps), ...) << std::endl;;
-  ((file << μs), ...) << std::endl;;
+  ((file << ps << " "), ...) << std::endl;;
+  ((file << μs << " "), ...) << std::endl;;
 
-  std::apply([&file](const Tensor<Real, ps>&... Js) -> void {
-        ((file << Js << 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 {
+      std::make_tuple(tensorFile.write(Js.data(), Js.size() * sizeof(Complex))...);
       } , 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 << std::endl;
+  file << Hz << " " << zSaddle.squaredNorm() << std::endl;
   file << eigenSz.eigenvalues().transpose() << std::endl;
-  file << φ << std::endl;
+  file << φ << " " << (xMin - zSaddle).norm() << std::endl;
 
   Cord c(M, zMin, zSaddle, nGs);
   c.relaxNewton(nTs, 1, 1e-10, 1e3);
 
-  Complex constraintError = 0;
+  Real reConstraintError = 0;
+  Real imConstraintError = 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));
+  unsigned nTest = nTs * 10;
+
+  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);
-    constraintError += pow(((Complex)(zi.transpose() * zi) - (Complex)N), 2);
+    Complex constraintError = (Complex)(zi.transpose() * zi) - (Complex)N;
+    reConstraintError += pow(real(constraintError), 2);
+    imConstraintError += pow(imag(constraintError), 2);
   }
 
-  file << nGs << " " << nTs << std::endl;;
-  file << c.totalCost(100 * nTs) / (100 * nTs) << " " << sqrt(imEnergyError) / (100 * nTs) << " " << sqrt(constraintError) / (100.0 * nTs) << std::endl;
+  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;
   }
-- 
cgit v1.2.3-54-g00ecf