Output full correlation files each step, and decrease γ on failure

1 files changed, 56 insertions, 47 deletions

diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 0805add..81e7634 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -2,6 +2,11 @@
 #include "p-spin.hpp"
 #include <getopt.h>
 #include <iostream>
+#include <fstream>
+
+std::string logFourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k) {
+  return prefix + "_" + std::to_string(p) + "_" + std::to_string(s) + "_" + std::to_string(λ) + "_" + std::to_string(τ₀) + "_" + std::to_string(β) + "_" + std::to_string(log2n) + "_" + std::to_string(Δτ)  + "_" + std::to_string(k) + ".dat";
+}
 
 int main(int argc, char* argv[]) {
   /* Model parameters */
@@ -65,7 +70,7 @@ int main(int argc, char* argv[]) {
   LogarithmicFourierTransform fft(N, k, Δτ, 4);
 
   Real Γ₀ = 1.0 + τ₀;
-  Real μ = 1.0;
+  Real μₜ₋₁ = 1.0;
 
   std::vector<Real> Cₜ₋₁(N);
   std::vector<Real> Rₜ₋₁(N);
@@ -75,20 +80,21 @@ int main(int argc, char* argv[]) {
   /* Start from the exact solution for β = 0 */
   for (unsigned n = 0; n < N; n++) {
     if (τ₀ != 1) {
-      Cₜ₋₁[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
+      Cₜ₋₁[n] = Γ₀ * (exp(-μₜ₋₁ * fft.t(n)) - μₜ₋₁ * τ₀ * exp(-fft.t(n) / τ₀)) / (μₜ₋₁ - pow(μₜ₋₁, 3) * pow(τ₀, 2));
     } else {
       Cₜ₋₁[n] = Γ₀ * exp(-fft.t(n)) * (1 + fft.t(n));
     }
-    Rₜ₋₁[n] = exp(-μ * fft.t(n));
+    Rₜ₋₁[n] = exp(-μₜ₋₁ * fft.t(n));
 
-    Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
-    Ȓₜ₋₁[n] = 1.0 / (μ + 1i * fft.ν(n));
+    Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μₜ₋₁, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
+    Ȓₜ₋₁[n] = 1.0 / (μₜ₋₁ + 1i * fft.ν(n));
   }
 
   std::vector<Real> Cₜ = Cₜ₋₁;
   std::vector<Real> Rₜ = Rₜ₋₁;
   std::vector<Complex> Ĉₜ = Ĉₜ₋₁;
   std::vector<Complex> Ȓₜ = Ȓₜ₋₁;
+  Real μₜ = μₜ₋₁;
 
   Real β = 0;
   while (β < βₘₐₓ) {
@@ -106,25 +112,13 @@ int main(int argc, char* argv[]) {
       std::vector<Complex> Ĉₜ₊₁(N);
       std::vector<Complex> Ȓₜ₊₁(N);
       for (unsigned n = 0; n < N; n++) {
-        Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μ + 1i * fft.ν(n));
-        Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μ + 1i * fft.ν(n));
+        Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μₜ + 1i * fft.ν(n));
+        Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + 1i * fft.ν(n));
       }
-
       std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
-
-      /*
-      for (unsigned n = 0; n < N; n++) {
-        RddfC[n] = Rₜ₊₁[n] * ddf(λ, p, s, Cₜ[n]);
-      }
-      RddfCt = fft.fourier(RddfC, false);
-
-      for (unsigned n = 0; n < N; n++) {
-        Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n));
-      }
-      */
       std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
 
-      μ *= pow(tanh(Cₜ₊₁[0]-1)+1, 0.05);
+      μₜ *= pow(tanh(Cₜ₊₁[0]-1)+1, 0.05);
 
       ΔC = 0;
       for (unsigned i = 0; i < N; i++) {
@@ -141,36 +135,51 @@ int main(int argc, char* argv[]) {
       }
 
       std::cerr << "\x1b[2K" << "\r";
-      std::cerr << β << " " << μ << " " << ΔC << " " << γ << " " << Cₜ[0];
+      std::cerr << β << " " << μₜ << " " << ΔC << " " << γ << " " << Cₜ[0];
     }
 
-    /* Integrate the energy using Simpson's rule */
-    Real E = 0;
-    for (unsigned n = 0; n < N/2-1; n++) {
-      Real h₂ₙ   = fft.t(2*n+1) - fft.t(2*n);
-      Real h₂ₙ₊₁ = fft.t(2*n+2) - fft.t(2*n+1);
-      Real f₂ₙ   = Rₜ[2*n]   * df(λ, p, s, Cₜ[2*n]);
-      Real f₂ₙ₊₁ = Rₜ[2*n+1] * df(λ, p, s, Cₜ[2*n+1]);
-      Real f₂ₙ₊₂ = Rₜ[2*n+2] * df(λ, p, s, Cₜ[2*n+2]);
-      E += (h₂ₙ + h₂ₙ₊₁) / 6 * (
-            (2 - h₂ₙ₊₁ / h₂ₙ) * f₂ₙ
-            + pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁
-            + (2 - h₂ₙ / h₂ₙ₊₁) * f₂ₙ₊₂
-          );
-    }
-    E *= β;
-
-    std::cerr << "\x1b[2K" << "\r";
-    std::cerr << β << " " << μ << " " << Ĉₜ[0].real() << " " << E << " " << γ << std::endl;
-    β += Δβ;
-    Cₜ₋₁ = Cₜ;
-    Rₜ₋₁ = Rₜ;
-    Ĉₜ₋₁ = Ĉₜ;
-    Ȓₜ₋₁ = Ȓₜ;
-  }
+    if (std::isnan(Cₜ[0])) {
+      Cₜ = Cₜ₋₁;
+      Rₜ = Rₜ₋₁;
+      Ĉₜ = Ĉₜ₋₁;
+      Ȓₜ = Ȓₜ₋₁;
+      μₜ = μₜ₋₁;
+      γ /= 2;
+    } else {
+      /* Integrate the energy using Simpson's rule */
+      Real E = 0;
+      for (unsigned n = 0; n < N/2-1; n++) {
+        Real h₂ₙ   = fft.t(2*n+1) - fft.t(2*n);
+        Real h₂ₙ₊₁ = fft.t(2*n+2) - fft.t(2*n+1);
+        Real f₂ₙ   = Rₜ[2*n]   * df(λ, p, s, Cₜ[2*n]);
+        Real f₂ₙ₊₁ = Rₜ[2*n+1] * df(λ, p, s, Cₜ[2*n+1]);
+        Real f₂ₙ₊₂ = Rₜ[2*n+2] * df(λ, p, s, Cₜ[2*n+2]);
+        E += (h₂ₙ + h₂ₙ₊₁) / 6 * (
+              (2 - h₂ₙ₊₁ / h₂ₙ) * f₂ₙ
+              + pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁
+              + (2 - h₂ₙ / h₂ₙ₊₁) * f₂ₙ₊₂
+            );
+      }
+      E *= β;
 
-  for (unsigned i = 0; i < N; i++) {
-    std::cout << fft.t(i) << " " << Cₜ[i] << std::endl;
+      std::cerr << "\x1b[2K" << "\r";
+      std::cerr << β << " " << μₜ << " " << Ĉₜ[0].real() << " " << E << " " << γ << std::endl;
+
+      std::ofstream outfile(logFourierFile("C", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::out | std::ios::binary);
+      outfile.write((const char*)(Cₜ.data()), N * sizeof(Real));
+      outfile.close();
+
+      std::ofstream outfileR(logFourierFile("R", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::out | std::ios::binary);
+      outfileR.write((const char*)(Rₜ.data()), N * sizeof(Real));
+      outfileR.close();
+
+      β += Δβ;
+      Cₜ₋₁ = Cₜ;
+      Rₜ₋₁ = Rₜ;
+      Ĉₜ₋₁ = Ĉₜ;
+      Ȓₜ₋₁ = Ȓₜ;
+      μₜ₋₁ = μₜ;
+    }
   }
 
   return 0;