Started implementing log-Fourier method, not working right now.

4 files changed, 228 insertions, 6 deletions

diff --git a/Makefile b/Makefile
index 688a99c..ad67e7e 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-all: walk correlations integrator fourier_integrator get_energy
+all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator
 
 CC := clang++ -std=c++17 -flto -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -DFFTW_THREADS=1
 
@@ -9,13 +9,16 @@ correlations: correlations.cpp
 	$(CC) correlations.cpp -o correlations
 
 integrator: integrator.cpp fourier.o fourier.hpp
-	$(CC) integrator.cpp fourier.o -lfftw3 -lfftw3_omp -o integrator
+	$(CC) integrator.cpp fourier.o -lfftw3 -lgsl -o integrator
 
 fourier.o: fourier.cpp fourier.hpp
 	$(CC) fourier.cpp -c -o fourier.o
 
 fourier_integrator: fourier_integrator.cpp fourier.hpp fourier.o
-	$(CC) fourier_integrator.cpp fourier.o -lfftw3 -lfftw3_omp -o fourier_integrator
+	$(CC) fourier_integrator.cpp fourier.o -lfftw3 -lgsl -o fourier_integrator
+
+log-fourier_integrator: log-fourier_integrator.cpp fourier.hpp fourier.o
+	$(CC) log-fourier_integrator.cpp fourier.o -lfftw3 -lgsl -o log-fourier_integrator
 
 get_energy: get_energy.cpp fourier.hpp fourier.o
-	$(CC) get_energy.cpp fourier.o -lfftw3 -lfftw3_omp -o get_energy
+	$(CC) get_energy.cpp fourier.o -lfftw3 -lgsl -o get_energy
diff --git a/fourier.cpp b/fourier.cpp
index bc4b633..cf45ad1 100644
--- a/fourier.cpp
+++ b/fourier.cpp
@@ -1,4 +1,5 @@
 #include "fourier.hpp"
+#include <boost/multiprecision/fwd.hpp>
 
 inline Real fP(unsigned p, Real q) {
   return 0.5 * pow(q, p);
@@ -27,8 +28,8 @@ Real ddf(Real λ, unsigned p, unsigned s, Real q) {
 FourierTransform::FourierTransform(unsigned n, Real Δω, Real Δτ, unsigned flags) : n(n), Δω(Δω), Δτ(Δτ) {
   a = fftw_alloc_real(2 * n);
   â = reinterpret_cast<Complex*>(fftw_alloc_complex(n + 1));
-  fftw_init_threads();
-  fftw_plan_with_nthreads(FFTW_THREADS);
+//  fftw_init_threads();
+//  fftw_plan_with_nthreads(FFTW_THREADS);
   fftw_import_wisdom_from_filename("fftw.wisdom");
   plan_r2c = fftw_plan_dft_r2c_1d(2 * n, a, reinterpret_cast<fftw_complex*>(â), flags);
   plan_c2r = fftw_plan_dft_c2r_1d(2 * n, reinterpret_cast<fftw_complex*>(â), a, flags);
@@ -99,6 +100,78 @@ void FourierTransform::writeToA(unsigned i, Real ai) {
   a[i] = ai;
 }
 
+LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δω, Real Δτ, Real Δs) : N(N), k(k), Δτ(Δτ), Δω(Δω), Δs(Δs) {
+  τₛ = -0.5 * N;
+  ωₛ = -0.5 * N;
+  sₛ = -0.5 * N;
+  a = reinterpret_cast<Complex*>(fftw_alloc_complex(N));
+  â = reinterpret_cast<Complex*>(fftw_alloc_complex(N));
+  fftw_import_wisdom_from_filename("fftw.wisdom");
+  a_to_â = fftw_plan_dft_1d(N, reinterpret_cast<fftw_complex*>(a), reinterpret_cast<fftw_complex*>(â), 1, 0);
+  â_to_a = fftw_plan_dft_1d(N, reinterpret_cast<fftw_complex*>(â), reinterpret_cast<fftw_complex*>(a), 1, 0);
+  fftw_export_wisdom_to_filename("fftw.wisdom");
+}
+
+LogarithmicFourierTransform::~LogarithmicFourierTransform() {
+  fftw_destroy_plan(a_to_â);
+  fftw_destroy_plan(â_to_a);
+  fftw_free(a);
+  fftw_free(â);
+  fftw_cleanup();
+}
+
+Real LogarithmicFourierTransform::τ(unsigned n) const {
+  return Δτ * (n + τₛ);
+}
+
+Real LogarithmicFourierTransform::ω(unsigned n) const {
+  return Δω * (n + ωₛ);
+}
+
+Real LogarithmicFourierTransform::s(unsigned n) const {
+  return Δs * (n + sₛ);
+}
+
+Real LogarithmicFourierTransform::t(unsigned n) const {
+  return exp(τ(n));
+}
+
+Real LogarithmicFourierTransform::ν(unsigned n) const {
+  return exp(ω(n));
+}
+
+Complex gamma(Complex z) {
+  gsl_sf_result logΓ;
+  gsl_sf_result argΓ;
+
+  gsl_sf_lngamma_complex_e(z.real(), z.imag(), &logΓ, &argΓ);
+
+  return exp(logΓ.val + 1i * argΓ.val);
+}
+
+std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real>& c, bool symmetric){ 
+  std::vector<Complex> ĉ(N);
+  std::vector<Real> σs = {1};
+  if (symmetric){
+    σs.push_back(-1);
+  }
+  for (Real σ : σs) {
+    for (unsigned n = 0; n < N; n++) {
+      a[n] = c[n] * exp((1 - k) * τ(n));
+    }
+    fftw_execute(a_to_â);
+    for (unsigned n = 0; n < N; n++) {
+      â[(n + N / 2) % N] *= pow(-1i * σ, 1i * (s(n)) - k) * gamma(k - 1i * (s(n)));
+    }
+    fftw_execute(â_to_a);
+    for (unsigned n = 0; n < N; n++) {
+      ĉ[n] += exp(-k * ω(n) * 4 * M_PI) * a[n] / pow(2 * M_PI, 1);
+    }
+  }
+
+  return ĉ;
+}
+
 std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ) {
   return prefix + "_" + std::to_string(p) + "_" + std::to_string(s) + "_" + std::to_string(λ) + "_" + std::to_string(τ₀) + "_" + std::to_string(y) + "_" + std::to_string(log2n) + "_" + std::to_string(τₘₐₓ) + ".dat";
 }
diff --git a/fourier.hpp b/fourier.hpp
index 791953b..22a57e7 100644
--- a/fourier.hpp
+++ b/fourier.hpp
@@ -4,6 +4,7 @@
 #include <complex>
 #include <vector>
 #include <tuple>
+#include <gsl/gsl_sf_gamma.h>
 
 #ifndef FFTW_THREADS
 #define FFTW_THREADS 1
@@ -36,6 +37,37 @@ public:
   std::vector<Real> convolve(const std::vector<Real>& Γh, const std::vector<Real>& R);
 };
 
+class LogarithmicFourierTransform {
+private:
+  Complex* a;
+  Complex* â;
+  fftw_plan a_to_â;
+  fftw_plan â_to_a;
+  unsigned N;
+  Real k;
+  Real Δτ;
+  Real Δω;
+  Real Δs;
+  Real τₛ;
+  Real ωₛ;
+  Real sₛ;
+public:
+  LogarithmicFourierTransform(unsigned N, Real k, Real Δω, Real Δτ, Real Δs);
+  ~LogarithmicFourierTransform();
+  Real τ(unsigned n) const;
+  Real ω(unsigned n) const;
+  Real t(unsigned n) const;
+  Real ν(unsigned n) const;
+  Real s(unsigned n) const;
+  std::vector<Complex> fourier(const std::vector<Real>& c, bool symmetric);
+  std::vector<Complex> fourier();
+  std::vector<Real> inverse(const std::vector<Complex>& ĉ);
+  void writeToA(unsigned i, Real ai);
+  std::vector<Real> convolve(const std::vector<Real>& Γh, const std::vector<Real>& R);
+};
+
+Complex gamma(Complex z);
+
 std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ);
 Real energy(const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real y, Real Δτ);
 std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ);
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
new file mode 100644
index 0000000..247270b
--- /dev/null
+++ b/log-fourier_integrator.cpp
@@ -0,0 +1,114 @@
+#include "fourier.hpp"
+#include <getopt.h>
+#include <iostream>
+#include <fstream>
+
+int main(int argc, char* argv[]) {
+  unsigned p = 2;
+  unsigned s = 2;
+  Real λ = 0.5;
+  Real τ₀ = 0;
+  Real β₀ = 0;
+  Real yₘₐₓ = 0.5;
+  Real Δy = 0.05;
+
+  unsigned log2n = 8;
+  Real τₘₐₓ = 20;
+
+  unsigned maxIterations = 1000;
+  Real ε = 1e-14;
+  Real γ = 1;
+
+  int opt;
+
+  while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:I:g:l")) != -1) {
+    switch (opt) {
+    case 'p':
+      p = atoi(optarg);
+      break;
+    case 's':
+      s = atoi(optarg);
+      break;
+    case '2':
+      log2n = atoi(optarg);
+      break;
+    case 'T':
+      τₘₐₓ = atof(optarg);
+      break;
+    case 't':
+      τ₀ = atof(optarg);
+      break;
+    case '0':
+      β₀ = atof(optarg);
+      break;
+    case 'y':
+      yₘₐₓ = atof(optarg);
+      break;
+    case 'd':
+      Δy = atof(optarg);
+      break;
+    case 'I':
+      maxIterations = (unsigned)atof(optarg);
+      break;
+    case 'g':
+      γ = atof(optarg);
+      break;
+    default:
+      exit(1);
+    }
+  }
+
+  unsigned N = pow(2, log2n);
+
+  Real Δτ = 1e-2;
+  Real Δω = 1e-2;
+  Real Δs = 1e-2;
+  Real k = 0.1;
+
+  Real Γ₀ = 1.0;
+  Real μ = Γ₀;
+  if (τ₀ > 0) {
+    μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);
+  }
+
+  std::vector<Real> C(N);
+  std::vector<Real> R(N);
+
+  LogarithmicFourierTransform fft(N, k, Δω, Δτ, Δs);
+
+  for (unsigned n = 0; n < N; n++) {
+    if (τ₀ > 0) {
+      C[n] = Γ₀ * (exp(-μ * fft.t(n)) - μ * τ₀ * exp(-fft.t(n) / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));
+    } else {
+      C[n] = Γ₀ * exp(-μ * fft.t(n)) / μ;
+    }
+    R[n] = exp(-μ * fft.t(n));
+  }
+
+  std::vector<Complex> Ct = fft.fourier(C, true);
+  std::vector<Complex> Rt = fft.fourier(R, false);
+
+  /*
+  for (unsigned n = 0; n < N; n++) {
+    std::cout << fft.t(n) << " " << C[n] << std::endl;
+  }
+  */
+
+  for (unsigned n = 0; n < N; n++) {
+    std::cout << fft.ν(n) << " " << Ct[n].real() << std::endl;
+  }
+
+  /*
+  // start from the exact solution for τ₀ = 0
+  for (unsigned i = 0; i < N + 1; i++) {
+    Real ω = i * Δω;
+    Ct[i] = 2 * Γ₀ / (pow(μ, 2) + pow(ω, 2)) / (1 + pow(τ₀ * ω, 2));
+    Rt[i] = 1.0 / (μ + 1i * ω);
+    Γt[i] = Γ₀ / (1 + pow(τ₀ * ω, 2));
+  }
+  C = fft.inverse(Ct);
+  R = fft.inverse(Rt);
+  */
+
+  return 0;
+}