Revert "Work to debug log-fourier"

This reverts commit 77e8c94e5b41287e97e36e34a97478400637d102.

5 files changed, 40 insertions, 101 deletions

diff --git a/Makefile b/Makefile
index 0497df5..c5cac94 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator log-fourier_integrator_long log_get_energy log_integrator log-fourier_test
+all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator log-fourier_integrator_long log_get_energy
 
 CC := clang++ -std=c++17 -flto -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -DFFTW_THREADS=1
 
@@ -32,15 +32,9 @@ fourier_integrator: fourier_integrator.cpp fourier.hpp fourier.o p-spin.o
 log-fourier_integrator: log-fourier_integrator.cpp log-fourier.hpp log-fourier.o p-spin.o
 	$(CC) log-fourier_integrator.cpp log-fourier.o p-spin.o -lfftw3 -lgsl -o log-fourier_integrator
 
-log-fourier_test: log-fourier_test.cpp log-fourier.hpp log-fourier.o
-	$(CC) log-fourier_test.cpp log-fourier.o -lfftw3 -lgsl -o log-fourier_test
-
 log-fourier_integrator_long: log-fourier_integrator.cpp log-fourier.hpp log-fourier_long.o p-spin_long.o
 	$(CC) log-fourier_integrator.cpp log-fourier_long.o p-spin_long.o -D QUAD_PRECISION -lfftw3l -lgsl -o log-fourier_integrator_long
 
-log_integrator: log_integrator.cpp log-fourier.hpp log-fourier.o p-spin.o
-	$(CC) log_integrator.cpp log-fourier.o p-spin.o -lfftw3 -lgsl -o log_integrator
-
 get_energy: get_energy.cpp fourier.hpp fourier.o p-spin.o
 	$(CC) get_energy.cpp fourier.o p-spin.o -lfftw3 -lgsl -o get_energy
 
diff --git a/log-fourier.cpp b/log-fourier.cpp
index bc2dd87..9d1f2cb 100644
--- a/log-fourier.cpp
+++ b/log-fourier.cpp
@@ -64,13 +64,15 @@ std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real
     for (unsigned n = 0; n < pad*N; n++) {
       if (n < N) {
         a[n] = c[n] * exp((1 - k) * τ(n));
+      } else if (n >= (pad - 1) * N) {
+        a[n] = c[pad*N-n-1] * exp((1 - k) * τ(pad*N-n-1));
       } else {
         a[n] = 0;
       }
     }
     FFTW_EXECUTE(a_to_â);
     for (unsigned n = 0; n < pad*N; n++) {
-      â[(pad*N / 2 + n) % (pad*N)] *= std::exp(II * s(n)) * std::pow(II * σ, II * s(n) - k) * Γ(k - II * s(n));
+      â[(pad*N / 2 + n) % (pad*N)] *= std::pow(II * σ, II * s(n) - k) * Γ(k - II * s(n));
     }
     FFTW_EXECUTE(â_to_a);
     for (unsigned n = 0; n < N; n++) {
@@ -85,16 +87,30 @@ std::vector<Real> LogarithmicFourierTransform::inverse(const std::vector<Complex
   std::vector<Real> c(N);
   std::vector<Real> σs = {1, -1};
   for (Real σ : σs) {
-    for (unsigned n = 0; n < N; n++) {
-      a[n] = (ĉ[n].real() + II * σ * ĉ[n].imag()) * std::exp((1 - k) * ω(n));
+    for (unsigned n = 0; n < pad * N; n++) {
+      if (n < N) {
+        if (σ < 0) {
+          a[n] = std::conj(ĉ[n]) * std::exp((1 - k) * ω(n));
+        } else {
+          a[n] = ĉ[n] * std::exp((1 - k) * ω(n));
+        }
+      } else if (n >= (pad - 1) * N) {
+        if (σ < 0) {
+          a[n] = ĉ[pad*N-n-1] * std::exp((1 - k) * ω(pad*N-n-1));
+        } else {
+          a[n] = std::conj(ĉ[pad*N-n-1]) * std::exp((1 - k) * ω(pad*N-n-1));
+        }
+      } else {
+        a[n] = 0;
+      }
     }
     FFTW_EXECUTE(a_to_â);
-    for (unsigned n = 0; n < N; n++) {
-      â[(N / 2 + n) % N] *= std::pow(-II * σ, II * s(n) - k) * Γ(k - II * s(n));
+    for (unsigned n = 0; n < pad*N; n++) {
+      â[(pad*N / 2 + n) % (pad*N)] *= std::pow(-II * σ, II * s(n) - k) * Γ(k - II * s(n));
     }
     FFTW_EXECUTE(â_to_a);
     for (unsigned n = 0; n < N; n++) {
-      c[n] += std::exp(-k * τ(n)) * a[n].real() / (2 * M_PI * N);
+      c[n] += std::exp(-k * τ(n)) * a[(pad - 1)*N+n].real() / (Real)(pad*N) / (2 * M_PI);
     }
   }
 
diff --git a/log-fourier.hpp b/log-fourier.hpp
index 675b1b0..f0b53ef 100644
--- a/log-fourier.hpp
+++ b/log-fourier.hpp
@@ -17,7 +17,6 @@ private:
   unsigned pad;
   Real k;
   Real Δτ;
-  Real Δs;
   Real τₛ;
   Real ωₛ;
   Real sₛ;
diff --git a/log-fourier_test.cpp b/log-fourier_test.cpp
deleted file mode 100644
index 678e4b1..0000000
--- a/log-fourier_test.cpp
+++ /dev/null
@@ -1,50 +0,0 @@
-#include "log-fourier.hpp"
-#include <getopt.h>
-#include <iostream>
-
-int main(int argc, char* argv[]) {
-  /* Log-Fourier parameters */
-  unsigned log2n = 8;
-  Real Δτ = 0.1;
-  Real k = 0.1;
-  unsigned pad = 4;
-
-  int opt;
-
-  while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:g:k:D:e:0:lS:x:P:")) != -1) {
-    switch (opt) {
-    case '2':
-      log2n = atoi(optarg);
-      break;
-    case 'k':
-      k = atof(optarg);
-      break;
-    case 'D':
-      Δτ = atof(optarg);
-      break;
-    case 'P':
-      pad = atoi(optarg);
-      break;
-    default:
-      exit(1);
-    }
-  }
-
-  unsigned N =log2n;
-
-  LogarithmicFourierTransform fft(N, k, Δτ, pad);
-
-  std::vector<Complex> a(N);
-
-  for (unsigned i = 0; i < N; i++) {
-    a[i] = 1 / (1 + pow(fft.ν(i), 2));
-  }
-
-  std::vector<Real> â = fft.inverse(a);
-
-  for (unsigned i = 0; i < N; i++) {
-    std::cout << fft.t(i) << " " << â[i] << std::endl;
-  }
-
-  return 0;
-}
diff --git a/log_integrator.cpp b/log_integrator.cpp
index 426af40..a8d9778 100644
--- a/log_integrator.cpp
+++ b/log_integrator.cpp
@@ -92,7 +92,6 @@ int main(int argc, char* argv[]) {
 
   std::vector<Real> Cₜ₋₁(N);
   std::vector<Real> Rₜ₋₁(N);
-  std::vector<Real> Γ(N);
   std::vector<Complex> Ĉₜ₋₁(N);
   std::vector<Complex> Ȓₜ₋₁(N);
 
@@ -109,7 +108,6 @@ int main(int argc, char* argv[]) {
         Cₜ₋₁[n] = Γ₀ * exp(-μₜ₋₁ * fft.t(n)) / μₜ₋₁;
       }
       Rₜ₋₁[n] = exp(-μₜ₋₁ * fft.t(n));
-      Γ[n] = (Γ₀ / τ₀) * exp(-fft.t(n) / τ₀);
 
       Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μₜ₋₁, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));
       Ȓₜ₋₁[n] = (Real)1.0 / (μₜ₋₁ + II * fft.ν(n));
@@ -144,47 +142,23 @@ int main(int argc, char* argv[]) {
 
       for (unsigned i = 0; i < N; i++) {
         dC[i] += -μₜ * Cₜ[i];
-        dR[i] += -μₜ * Rₜ[i];
-        for (unsigned j = 1; j < N - i; j++) {
-          Real Γᵢ₊ⱼ = (Γ[i + j] + Γ[i + j - 1]) / 2;
-          signed jR = round(log(fft.t(i+j) - fft.t(i)) / Δτ) + N / 2;
-          Real Rⱼ = jR >= N ? 0 : jR < 0 ? 1 : Rₜ[jR];
-          Real Cⱼ = jR >= N ? 0 : jR < 0 ? 1 : Cₜ[jR];
-          Real ΓR2 = Γᵢ₊ⱼ * Rⱼ;
-          Real dfCR = (dfC[i + j] * Rⱼ + dfC[i + j - 1] * Rⱼ) / 2;
-          Real RddfCC = (RddfC[i + j] * Cⱼ + RddfC[i + j - 1] * Cⱼ) / 2;
-          dC[i] += (fft.t(i+j) - fft.t(i+j - 1)) * (ΓR2 + pow(β, 2) * (dfCR + RddfCC));
+        Real ΓR;
+        for (unsigned j = 0; j < N; j++) {
+          
         }
-        for (unsigned j = 1; j < i; j++) {
-          signed jR = round(log(fft.t(i) - fft.t(i-j)) / Δτ) + N / 2;
-          Real Rⱼ = jR >= N ? 0 : jR < 0 ? 1 : Rₜ[jR];
-          Real Cⱼ = jR >= N ? 0 : jR < 0 ? 1 : Cₜ[jR];
-          dC[i] += (fft.t(i-j) - fft.t(i-j - 1)) * pow(β, 2) * (RddfC[i - 1 - j] * Cⱼ + RddfC[i - 1 - j+1] * Cⱼ) / 2;
-          dR[i] += (fft.t(i-j) - fft.t(i-j - 1)) * pow(β, 2) * (RddfC[i - 1 - j] * Rⱼ + RddfC[i - 1 - j+1] * Rⱼ) / 2;
-        }
-
-        if (dC[i] > 0) dC[i] = 0;
-        if (dR[i] > 0) dR[i] = 0;
       }
 
-      std::vector<Real> Rₜ₊₁(N);
-      std::vector<Real> Cₜ₊₁(N);
 
-      Rₜ₊₁[0] = 1;
-      Cₜ₊₁[N - 1] = 0;
-
-      for (unsigned i = 1; i < N; i++) {
-        Cₜ₊₁[N - 1 - i] = Cₜ₊₁[N - i] - (fft.t(N - i) - fft.t(N - i - 1)) * (dC[N - i] + dC[N - 1 - i]) / 2;
-        Rₜ₊₁[i] = Rₜ₊₁[i - 1] + (fft.t(i) - fft.t(i - 1)) * (dR[i - 1] + dR[i]) / 2;
+      std::vector<Complex> Ĉₜ₊₁(N);
+      std::vector<Complex> Ȓₜ₊₁(N);
+      for (unsigned n = 0; n < N; n++) {
+        Ȓₜ₊₁[n] = ((Real)1.0 + std::pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μₜ + II * fft.ν(n));
+        Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + II * fft.ν(n));
       }
+      std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
+      std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
 
-      μₜ /= pow(tanh(Cₜ₊₁[0]-1)+1, x);
-
-      Real C₀ = Cₜ₊₁[0];
-
-      for (unsigned i = 0; i < N; i++) {
-        Rₜ₊₁[i] = Rₜ₊₁[i] < 0 ? 0 : Rₜ₊₁[i];
-      }
+      μₜ *= pow(tanh(Cₜ₊₁[0]-1)+1, x);
 
       ΔCₜ = 0;
       for (unsigned i = 0; i < N; i++) {
@@ -209,6 +183,8 @@ int main(int argc, char* argv[]) {
       for (unsigned i = 0; i < N; i++) {
         Cₜ[i] += γ * (Cₜ₊₁[i] - Cₜ[i]);
         Rₜ[i] += γ * (Rₜ₊₁[i] - Rₜ[i]);
+        Ĉₜ[i] += γ * (Ĉₜ₊₁[i] - Ĉₜ[i]);
+        Ȓₜ[i] += γ * (Ȓₜ₊₁[i] - Ȓₜ[i]);
       }
 
       std::cerr << "\x1b[2K" << "\r";
@@ -219,6 +195,8 @@ int main(int argc, char* argv[]) {
       γ₀ /= 2;
       Cₜ = Cₜ₋₁;
       Rₜ = Rₜ₋₁;
+      Ĉₜ = Ĉₜ₋₁;
+      Ȓₜ = Ȓₜ₋₁;
       μₜ = μₜ₋₁;
     } else {
       Real E = energy(fft, Cₜ, Rₜ, p, s, λ, β);
@@ -231,6 +209,8 @@ int main(int argc, char* argv[]) {
       β += Δβ;
       Cₜ₋₁ = Cₜ;
       Rₜ₋₁ = Rₜ;
+      Ĉₜ₋₁ = Ĉₜ;
+      Ȓₜ₋₁ = Ȓₜ;
       μₜ₋₁ = μₜ;
     }
   }