diff options
| -rw-r--r-- | .gitignore | 11 | ||||
| -rw-r--r-- | Makefile | 37 | ||||
| -rw-r--r-- | fourier.cpp | 239 | ||||
| -rw-r--r-- | fourier.hpp | 70 | ||||
| -rw-r--r-- | fourier_integrator.cpp | 184 | ||||
| -rw-r--r-- | get_energy.cpp | 88 | ||||
| -rw-r--r-- | integrator.cpp | 139 | ||||
| -rw-r--r-- | log-fourier.cpp | 278 | ||||
| -rw-r--r-- | log-fourier.hpp | 58 | ||||
| -rw-r--r-- | log-fourier_integrator.cpp | 230 | ||||
| -rw-r--r-- | log_get_energy.cpp | 100 | ||||
| -rw-r--r-- | p-spin.cpp | 25 | ||||
| -rw-r--r-- | p-spin.hpp | 7 | ||||
| -rw-r--r-- | print_correlations.cpp | 86 | ||||
| -rw-r--r-- | types.hpp | 43 |
15 files changed, 795 insertions, 800 deletions
@@ -1,9 +1,10 @@ *.dat *.o -correlations -fourier_integrator vgcore.* -walk -get_energy -integrator fftw.wisdom +walk +correlations +log-fourier_integrator +log-fourier_integrator_long +log_get_energy +log_get_energy_long @@ -1,4 +1,4 @@ -all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator +all: walk correlations log-fourier_integrator log-fourier_integrator_long log_get_energy log_get_energy_long print_correlations print_correlations_long CC := clang++ -std=c++17 -flto -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -DFFTW_THREADS=1 @@ -8,17 +8,32 @@ walk: walk.cpp correlations: correlations.cpp $(CC) correlations.cpp -o correlations -integrator: integrator.cpp fourier.o fourier.hpp - $(CC) integrator.cpp fourier.o -lfftw3 -lgsl -o integrator +p-spin.o: p-spin.cpp p-spin.hpp types.hpp + $(CC) p-spin.cpp -c -o p-spin.o -fourier.o: fourier.cpp fourier.hpp - $(CC) fourier.cpp -c -o fourier.o +p-spin_long.o: p-spin.cpp p-spin.hpp types.hpp + $(CC) p-spin.cpp -D QUAD_PRECISION -c -o p-spin_long.o -fourier_integrator: fourier_integrator.cpp fourier.hpp fourier.o - $(CC) fourier_integrator.cpp fourier.o -lfftw3 -lgsl -o fourier_integrator +log-fourier.o: log-fourier.cpp log-fourier.hpp types.hpp + $(CC) log-fourier.cpp -c -o log-fourier.o -log-fourier_integrator: log-fourier_integrator.cpp fourier.hpp fourier.o - $(CC) log-fourier_integrator.cpp fourier.o -lfftw3 -lgsl -o log-fourier_integrator +log-fourier_long.o: log-fourier.cpp log-fourier.hpp types.hpp + $(CC) log-fourier.cpp -D QUAD_PRECISION -c -o log-fourier_long.o -get_energy: get_energy.cpp fourier.hpp fourier.o - $(CC) get_energy.cpp fourier.o -lfftw3 -lgsl -o get_energy +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_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_get_energy: log_get_energy.cpp log-fourier.hpp log-fourier.o p-spin.o + $(CC) log_get_energy.cpp log-fourier.o p-spin.o -lfftw3 -lgsl -o log_get_energy + +log_get_energy_long: log_get_energy.cpp log-fourier.hpp log-fourier_long.o p-spin_long.o + $(CC) log_get_energy.cpp log-fourier_long.o p-spin_long.o -D QUAD_PRECISION -lfftw3l -lgsl -o log_get_energy_long + +print_correlations: print_correlations.cpp log-fourier.hpp log-fourier.o p-spin.o + $(CC) print_correlations.cpp log-fourier.o p-spin.o -lfftw3 -lgsl -o print_correlations + +print_correlations_long: print_correlations.cpp log-fourier.hpp log-fourier_long.o p-spin_long.o + $(CC) print_correlations.cpp log-fourier_long.o p-spin_long.o -D QUAD_PRECISION -lfftw3l -lgsl -o print_correlations_long diff --git a/fourier.cpp b/fourier.cpp deleted file mode 100644 index bb1c92f..0000000 --- a/fourier.cpp +++ /dev/null @@ -1,239 +0,0 @@ -#include "fourier.hpp" -#include <fftw3.h> - -inline Real fP(unsigned p, Real q) { - return 0.5 * pow(q, p); -} - -inline Real dfP(unsigned p, Real q) { - return 0.5 * p * pow(q, p - 1); -} - -inline Real ddfP(unsigned p, Real q) { - return 0.5 * p * (p - 1) * pow(q, p - 2); -} - -Real f(Real λ, unsigned p, unsigned s, Real q) { - return (1 - λ) * fP(p, q) + λ * fP(s, q); -} - -Real df(Real λ, unsigned p, unsigned s, Real q) { - return (1 - λ) * dfP(p, q) + λ * dfP(s, q); -} - -Real ddf(Real λ, unsigned p, unsigned s, Real q) { - return (1 - λ) * ddfP(p, q) + λ * ddfP(s, 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_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); - fftw_export_wisdom_to_filename("fftw.wisdom"); -} - -FourierTransform::~FourierTransform() { - fftw_destroy_plan(plan_r2c); - fftw_destroy_plan(plan_c2r); - fftw_free(a); - fftw_free(â); - fftw_cleanup(); -} - -std::vector<Complex> FourierTransform::fourier(const std::vector<Real>& c) { - for (unsigned i = 0; i < 2 * n; i++) { - a[i] = c[i]; - } - fftw_execute(plan_r2c); - std::vector<Complex> ĉ(n + 1); - for (unsigned i = 0; i < n + 1; i++) { - ĉ[i] = â[i] * (Δτ * M_PI); - } - return ĉ; -} - -std::vector<Complex> FourierTransform::fourier() { - fftw_execute(plan_r2c); - std::vector<Complex> ĉ(n+1); - for (unsigned i = 0; i < n+1; i++) { - ĉ[i] = â[i] * (Δτ * M_PI); - } - return ĉ; -} - -std::vector<Real> FourierTransform::convolve(const std::vector<Real>& Γh, const std::vector<Real>& R) { - a[0] = 0; - for (unsigned i = 1; i < n; i++) { - a[i] = R[i]; - a[2 * n - i] = -R[i]; - } - fftw_execute(plan_r2c); - for (unsigned i = 1; i < n + 1; i++) { - â[i] *= Γh[i] * (Δτ * M_PI); - } - fftw_execute(plan_c2r); - std::vector<Real> dC(n); - for (unsigned i = 0; i < n; i++) { - dC[i] = a[i] * (Δω / (2 * M_PI)); - } - - return dC; -} - -std::vector<Real> FourierTransform::inverse(const std::vector<Complex>& ĉ) { - for (unsigned i = 0; i < n + 1; i++) { - â[i] = ĉ[i]; - } - fftw_execute(plan_c2r); - std::vector<Real> c(2*n); - for (unsigned i = 0; i < 2*n; i++) { - c[i] = a[i] * (Δω / (2 * M_PI)); - } - return c; -} - -void FourierTransform::writeToA(unsigned i, Real ai) { - a[i] = ai; -} - -LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad) : N(N), pad(pad), k(k), Δτ(Δτ) { - τₛ = -0.5 * N; - ωₛ = -0.5 * N; - sₛ = -0.5 * pad * N; - a = reinterpret_cast<Complex*>(fftw_alloc_complex(pad*N)); - â = reinterpret_cast<Complex*>(fftw_alloc_complex(pad*N)); - fftw_import_wisdom_from_filename("fftw.wisdom"); - a_to_â = fftw_plan_dft_1d(pad*N, reinterpret_cast<fftw_complex*>(a), reinterpret_cast<fftw_complex*>(â), FFTW_BACKWARD, 0); - â_to_a = fftw_plan_dft_1d(pad*N, reinterpret_cast<fftw_complex*>(â), reinterpret_cast<fftw_complex*>(a), FFTW_BACKWARD, 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 (n + sₛ) * 2*M_PI / (pad * N * Δτ); -} - -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 < pad*N; n++) { - if (n < N) { - a[n] = c[n] * exp((1 - k) * τ(n)); - } else { - a[n] = 0; - } - } - fftw_execute(a_to_â); - for (unsigned n = 0; n < pad*N; n++) { - â[(pad*N / 2 + n) % (pad*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)) * a[(pad - 1)*N+n].real() / (Real)(pad*N); - } - } - - return ĉ; -} - -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 < pad * N; n++) { - if (n < N) { - a[n] = ĉ[n] * exp((1 - k) * ω(n)); - } else { - a[n] = 0; - } - } - fftw_execute(a_to_â); - for (unsigned n = 0; n < pad*N; n++) { - â[(pad*N / 2 + n) % (pad*N)] *= pow(-1i * σ, 1i * s(n) - k) * gamma(k - 1i * s(n)); - } - fftw_execute(â_to_a); - for (unsigned n = 0; n < N; n++) { - c[n] += exp(-k * τ(n)) * a[(pad - 1)*N+n].real() / (Real)(pad*N) / (2 * M_PI); - } - } - - return c; -} - -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"; -} - -Real energy(const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real y, Real Δτ) { - Real e = 0; - - for (unsigned i = 0; i < C.size() / 2; i++) { - e += y * R[i] * df(λ, p, s, C[i]) * M_PI * Δτ; - } - - return e; -} - -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 λ) { - for (unsigned i = 0; i < C.size() / 2; i++) { - fft.writeToA(i, R[i] * ddf(λ, p, s, C[i])); - } - for (unsigned i = C.size() / 2; i < C.size(); i++) { - fft.writeToA(i, 0); - } - std::vector<Complex> RddfCt = fft.fourier(); - - for (unsigned i = 0; i < C.size(); i++) { - fft.writeToA(i, df(λ, p, s, C[i])); - } - std::vector<Complex> dfCt = fft.fourier(); - - return {RddfCt, dfCt}; -} - -Real estimateZ(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Complex>& Ct, const std::vector<Real>& R, const std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real y) { - auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ); - Real Γ₀ = 1 + τ₀ / 2; - - return ((Γ₀ * std::conj(Rt[0]) + pow(y, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real(); -} diff --git a/fourier.hpp b/fourier.hpp deleted file mode 100644 index 1708667..0000000 --- a/fourier.hpp +++ /dev/null @@ -1,70 +0,0 @@ -#pragma once -#include <cmath> -#include <fftw3.h> -#include <complex> -#include <vector> -#include <tuple> -#include <gsl/gsl_sf_gamma.h> - -#ifndef FFTW_THREADS -#define FFTW_THREADS 1 -#endif - -using Real = double; -using Complex = std::complex<Real>; -using namespace std::complex_literals; - -Real f(Real λ, unsigned p, unsigned s, Real q); -Real df(Real λ, unsigned p, unsigned s, Real q); -Real ddf(Real λ, unsigned p, unsigned s, Real q); - -class FourierTransform { -private: - Real* a; - Complex* â; - fftw_plan plan_r2c; - fftw_plan plan_c2r; - unsigned n; - Real Δω; - Real Δτ; -public: - FourierTransform(unsigned n, Real Δω, Real Δτ, unsigned flags = 0); - ~FourierTransform(); - std::vector<Complex> fourier(const std::vector<Real>& c); - 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); -}; - -class LogarithmicFourierTransform { -private: - Complex* a; - Complex* â; - fftw_plan a_to_â; - fftw_plan â_to_a; - unsigned N; - unsigned pad; - Real k; - Real Δτ; - Real τₛ; - Real ωₛ; - Real sₛ; -public: - LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad = 4); - ~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<Real> inverse(const std::vector<Complex>& ĉ); -}; - -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 λ); -Real estimateZ(FourierTransform& fft, const std::vector<Real>& C, const std::vector<Complex>& Ct, const std::vector<Real>& R, const std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real y); diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp deleted file mode 100644 index 4f4eff5..0000000 --- a/fourier_integrator.cpp +++ /dev/null @@ -1,184 +0,0 @@ -#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 y₀ = 0; - Real yₘₐₓ = 0.5; - Real Δy = 0.05; - - unsigned log2n = 8; - Real τₘₐₓ = 20; - - unsigned maxIterations = 1000; - Real ε = 1e-14; - Real γ = 1; - - bool loadData = false; - - 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': - y₀ = 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; - case 'l': - loadData = true; - break; - default: - exit(1); - } - } - - unsigned n = pow(2, log2n); - - Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / n; - Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ); - - Real z = 0.5; - Real Γ₀ = 1 + τ₀ / 2; - - std::vector<Real> C(2 * n); - std::vector<Real> R(2 * n); - - FourierTransform fft(n, Δω, Δτ); - std::vector<Complex> Ct; - std::vector<Complex> Rt; - - Real y = y₀; - - if (!loadData) { - // start from the exact solution for τ₀ = 0 - for (unsigned i = 0; i < n; i++) { - Real τ = i * Δτ * M_PI; - if (τ₀ > 0) { - C[i] = Γ₀ / 2 * (exp(-z * τ) - z * τ₀ * exp(-τ / τ₀)) / (z - pow(z, 3) * pow(τ₀, 2)); - } else { - C[i] = Γ₀ / 2 * exp(-z * τ) / z; - } - if (i > 0) { - C[2 * n - i] = C[i]; - } - R[i] = exp(-z * τ); - } - Ct = fft.fourier(C); - Rt = fft.fourier(R); - } else { - std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary); - cfile.read((char*)(C.data()), (C.size() / 2) * sizeof(Real)); - cfile.close(); - for (unsigned i = 1; i < n; i++) { - C[2 * n - i] = C[i]; - } - std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary); - rfile.read((char*)(R.data()), (R.size() / 2) * sizeof(Real)); - rfile.close(); - - Ct = fft.fourier(C); - Rt = fft.fourier(R); - - z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, y); - } - - while (y += Δy, y <= yₘₐₓ) { - Real ΔC = 1; - Real ΔCprev = 1000; - unsigned it = 0; - while (sqrt(2 * ΔC / C.size()) > ε) { - it++; - auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ); - - for (unsigned i = 0; i < Rt.size(); i++) { - Real ω = i * Δω; - Rt[i] = (1.0 + pow(y, 2) * RddfCt[i] * Rt[i]) / (z + 1i * ω); - } - - for (unsigned i = 0; i < Ct.size(); i++) { - Real ω = i * Δω; - Ct[i] = (Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(y, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (z + 1i * ω); - } - - std::vector<Real> Cnew = fft.inverse(Ct); - std::vector<Real> Rnew = fft.inverse(Rt); - for (unsigned i = n; i < 2 * n; i++) { - Rnew[i] = 0; - } - - ΔC = 0; - for (unsigned i = 0; i < Cnew.size() / 2; i++) { - ΔC += pow(Cnew[i] - C[i], 2); - ΔC += pow(Rnew[i] - R[i], 2); - } - - for (unsigned i = 0; i < Cnew.size(); i++) { - C[i] += γ * (Cnew[i] - C[i]); - } - - for (unsigned i = 0; i < Rnew.size() / 2; i++) { - R[i] += γ * (Rnew[i] - R[i]); - } - - z *= Cnew[0]; - - if (it % maxIterations == 0) { - if (ΔC < ΔCprev) { - γ = std::min(1.1 * γ, 1.0); - } else { - γ /= 2; - } - - ΔCprev = ΔC; - } - - std::cerr << it << " " << p << " " << s << " " << τ₀ << " " << y << " " << sqrt(2 * ΔC / C.size()) << " " << γ << std::endl; - - } - - Real e = energy(C, R, p, s, λ, y, Δτ); - - std::cerr << "y " << y << " " << e << " " << z << std::endl; - - std::ofstream outfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::out | std::ios::binary); - outfile.write((const char*)(C.data()), (C.size() / 2) * sizeof(Real)); - outfile.close(); - - std::ofstream outfileR(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::out | std::ios::binary); - outfileR.write((const char*)(R.data()), (R.size() / 2) * sizeof(Real)); - outfileR.close(); - } - - return 0; -} diff --git a/get_energy.cpp b/get_energy.cpp deleted file mode 100644 index c33c04e..0000000 --- a/get_energy.cpp +++ /dev/null @@ -1,88 +0,0 @@ -#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 y₀ = 0; - Real yₘₐₓ = 0.5; - Real Δy = 0.05; - - unsigned log2n = 8; - Real τₘₐₓ = 20; - - int opt; - - while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:")) != -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': - y₀ = atof(optarg); - break; - case 'y': - yₘₐₓ = atof(optarg); - break; - case 'd': - Δy = atof(optarg); - break; - default: - exit(1); - } - } - - unsigned n = pow(2, log2n); - - Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / n; - Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ); - - Real y = y₀; - - FourierTransform fft(n, Δω, Δτ); - - std::vector<Real> C(2 * n); - std::vector<Real> R(2 * n); - - while (y += Δy, y <= yₘₐₓ) { - std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary); - if (cfile.is_open()) { - cfile.read((char*)(C.data()), (C.size() / 2) * sizeof(Real)); - cfile.close(); - for (unsigned i = 1; i < n; i++) { - C[2 * n - i] = C[i]; - } - - std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, y, log2n, τₘₐₓ), std::ios::binary); - rfile.read((char*)(R.data()), (R.size() / 2) * sizeof(Real)); - rfile.close(); - - Real e = energy(C, R, p, s, λ, y, Δτ); - - std::vector<Complex> Ct = fft.fourier(C); - std::vector<Complex> Rt = fft.fourier(R); - - Real z = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, y); - - std::cout << y << " " << e << " " << Ct[0].real() << " " << z << std::endl; - } - } - - return 0; -} diff --git a/integrator.cpp b/integrator.cpp deleted file mode 100644 index 7e5c512..0000000 --- a/integrator.cpp +++ /dev/null @@ -1,139 +0,0 @@ -#include "fourier.hpp" -#include <fstream> -#include <fftw3.h> -#include <getopt.h> -#include <iostream> - -Real energy(const std::vector<Real>& C, std::vector<Real>& R, Real λ, unsigned p, unsigned s, Real Δτ) { - Real I = 0; - for (unsigned σ = 0; σ < C.size(); σ++) { - I += Δτ * df(λ, p, s, C[σ]) * R[σ]; - } - return I; -} - -int main(int argc, char* argv[]) { - unsigned p = 3; - unsigned s = 4; - Real λ = 0.5; - Real τₘₐₓ = 1e3; - Real τ₀ = 0; - Real β₀ = 0; - Real βₘₐₓ = 1; - Real Δβ = 1e-2; - unsigned iterations = 10; - unsigned log2n = 8; - Real ε = 1e-14; - - int opt; - - while ((opt = getopt(argc, argv, "T:2:t:0:b:d:I:")) != -1) { - switch (opt) { - case 'T': - τₘₐₓ = atof(optarg); - break; - case '2': - log2n = atof(optarg); - break; - case 't': - τ₀ = atof(optarg); - break; - case '0': - β₀ = atof(optarg); - break; - case 'b': - βₘₐₓ = atof(optarg); - break; - case 'd': - Δβ = atof(optarg); - break; - case 'I': - iterations = (unsigned)atof(optarg); - break; - default: - exit(1); - } - } - - unsigned N = pow(2, log2n); - - Real Δτ = (1 + τ₀ / 2) * τₘₐₓ / M_PI / N; - Real Δω = M_PI / ((1 + τ₀ / 2) * τₘₐₓ); - - FourierTransform fft(N, Δω, Δτ, FFTW_ESTIMATE); - - Real Γ₀ = 1; - Real μ = 1; - if (τ₀ > 0) { - μ = (sqrt(1+4*Γ₀*τ₀) - 1) / (2 * τ₀); - } - - Real τ = 0; - std::vector<Real> C(N); - std::vector<Real> R(N); - std::vector<Real> Γ(N); - std::vector<Real> Γh(N+1); - - Γh[0] = Γ₀; - - for (unsigned i = 0; i < N; i++) { - Real τ = i * Δτ; - Real ω = (i + 1) * Δω * M_PI; - if (τ₀ > 0) { - C[i] = (Γ₀ / μ) * (exp(-μ * τ) - μ * τ₀ * exp(-τ / τ₀)) / (1 - pow(μ * τ₀, 2)); - Γ[i] = (Γ₀ / τ₀) * exp(-τ / τ₀); - } else { - C[i] = (Γ₀ / μ) * exp(-μ * τ); - } - Γh[i+1] = Γ₀ / (1 + pow(ω * τ₀, 2)); - R[i] = exp(-μ * τ); - } - - for (Real β = β₀; β < βₘₐₓ; β += Δβ) { - Real Rerr = 100; - while (sqrt(Rerr / N) > ε) { - /* First step: integrate R from C */ - std::vector<Real> R₊(N); - R₊[0] = 1; - for (unsigned i = 1; i < N; i++) { - Real I = 0; - for (unsigned j = 0; j <= i; j++) { - I += R[i - j] * ddf(λ, p, s, C[i - j]) * R[j] * Δτ; - } - Real dR = -μ * R₊[i - 1] + pow(β, 2) * I; - R₊[i] = R₊[i - 1] + dR * Δτ; - } - - Rerr = 0; - for (unsigned i = 0; i < N; i++) { - Rerr += pow(R[i] - R₊[i], 2); - } - - R = R₊; - - /* Second step: integrate C from R */ - std::vector<Real> dC = fft.convolve(Γh, R); - Real Cₜ₊₁ = 0; - for (unsigned i = 0; i < N; i++) { - Real Cₜ = Cₜ₊₁ + dC[N - i - 1] * Δτ; - C[N - i - 1] = Cₜ; - Cₜ₊₁ = Cₜ; - } - - /* Third step: adjust μ */ - μ *= C[0]; - - std::cerr << β << " " << sqrt(Rerr / N) << std::endl; - } - - std::ofstream outfile(fourierFile("Ci", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary); - outfile.write((const char*)(C.data()), N * sizeof(Real)); - outfile.close(); - - std::ofstream outfileR(fourierFile("Ri", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary); - outfileR.write((const char*)(R.data()), N * sizeof(Real)); - outfileR.close(); - } - - return 0; -} diff --git a/log-fourier.cpp b/log-fourier.cpp new file mode 100644 index 0000000..16a7f3f --- /dev/null +++ b/log-fourier.cpp @@ -0,0 +1,278 @@ +#include "log-fourier.hpp" +#include "p-spin.hpp" + +#include <fstream> + +Complex Γ(Complex z) { + gsl_sf_result logΓ; + gsl_sf_result argΓ; + + gsl_sf_lngamma_complex_e((double)z.real(), (double)z.imag(), &logΓ, &argΓ); + + return std::exp((Real)logΓ.val + II * (Real)argΓ.val); +} + +LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad, Real shift) : N(N), pad(pad), k(k), Δτ(Δτ), ts(N), νs(N), Γs(pad * N), exp1kω(N), exp1kτ(N), expkω(N), expkτ(N), shift(shift) { + τₛ = -0.5 * N; + ωₛ = -0.5 * N; + sₛ = -0.5 * pad * N; + a = reinterpret_cast<Complex*>(FFTW_ALLOC_COMPLEX(pad*N)); + â = reinterpret_cast<Complex*>(FFTW_ALLOC_COMPLEX(pad*N)); + FFTW_IMPORT_WISDOM("fftw.wisdom"); + a_to_â = FFTW_PLAN_DFT_1D(pad*N, reinterpret_cast<FFTW_COMPLEX*>(a), reinterpret_cast<FFTW_COMPLEX*>(â), FFTW_BACKWARD, 0); + â_to_a = FFTW_PLAN_DFT_1D(pad*N, reinterpret_cast<FFTW_COMPLEX*>(â), reinterpret_cast<FFTW_COMPLEX*>(a), FFTW_BACKWARD, 0); + FFTW_EXPORT_WISDOM("fftw.wisdom"); + for (unsigned n = 0; n < N; n++) { + ts[n] = std::exp(τ(n)) / shift; + νs[n] = std::exp(ω(n)) * shift; + exp1kτ[n] = std::exp((1 - k) * τ(n)); + exp1kω[n] = std::exp((1 - k) * ω(n)); + expkτ[n] = std::exp(-k * τ(n)); + expkω[n] = std::exp(-k * ω(n)); + } + for (unsigned n = 0; n < pad * N; n++) { + Γs[n] = Γ(k - II * s(n)); + } +} + +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 (n + sₛ) * 2*M_PI / (pad * N * Δτ); +} + +Real LogarithmicFourierTransform::t(unsigned n) const { + return ts[n]; +} + +Real LogarithmicFourierTransform::ν(unsigned n) const { + return νs[n]; +} + +std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real>& c, bool symmetric) { + std::vector<Complex> ĉ(N); + std::vector<Real> σs = {1}; + /* c is either even or zero for negative arguments */ + if (symmetric){ + σs.push_back(-1); + } + for (Real σ : σs) { + for (unsigned n = 0; n < pad*N; n++) { + if (n < N) { + a[n] = c[n] * exp1kτ[n]; + } else if (n >= (pad - 1) * N) { + a[n] = c[pad*N-n-1] * exp1kτ[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::pow(II * σ, II * s(n) - k) * Γs[n]; + } + FFTW_EXECUTE(â_to_a); + for (unsigned n = 0; n < N; n++) { + ĉ[n] += expkω[n] * a[(pad - 1)*N+n] / (Real)(pad*N); + } + } + + for (unsigned n = 0; n < N; n++) { + ĉ[n] -= ĉ[N - 1]; + if (symmetric) ĉ[n] = ĉ[n].real(); + ĉ[n] /= shift; + } + + return ĉ; +} + +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 < pad * N; n++) { + if (n < N) { + a[n] = (ĉ[n].real() + II * σ * ĉ[n].imag()) * exp1kω[n]; + } else if (n >= (pad - 1) * N) { + a[n] = (ĉ[pad*N-n-1].real() - II * σ * ĉ[pad*N-n-1].imag()) * exp1kω[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::pow(-II * σ, II * s(n) - k) * Γs[n]; + } + FFTW_EXECUTE(â_to_a); + for (unsigned n = 0; n < N; n++) { + c[n] += expkτ[n] * a[(pad - 1)*N+n].real() / (Real)(pad*N) / (2 * M_PI); + } + } + + for (unsigned n = 0; n < N; n++) { + c[n] -= c[N - 1]; + c[n] *= shift; + } + + return c; +} + +std::string logFourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real shift) { + 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(shift) + ".dat"; +} + +void logFourierSave(const std::vector<Real>& C, const std::vector<Real>& R, const std::vector<Complex>& Ĉ, const std::vector<Complex>& Ȓ, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k) { + unsigned N = std::pow(2, log2n); + 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 outfileCt(logFourierFile("Ct", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::out | std::ios::binary); + outfileCt.write((const char*)(Ĉ.data()), N * sizeof(Complex)); + outfileCt.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(); + + std::ofstream outfileRt(logFourierFile("Rt", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::out | std::ios::binary); + outfileRt.write((const char*)(Ȓ.data()), N * sizeof(Complex)); + outfileRt.close(); +} + +bool logFourierLoad(std::vector<Real>& C, std::vector<Real>& R, std::vector<Complex>& Ĉ, std::vector<Complex>& Ȓ, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k) { + std::ifstream cfile(logFourierFile("C", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::binary); + std::ifstream rfile(logFourierFile("R", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::binary); + std::ifstream ctfile(logFourierFile("Ct", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::binary); + std::ifstream rtfile(logFourierFile("Rt", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::binary); + + if ((!cfile.is_open() || !rfile.is_open()) || (!ctfile.is_open() || !rtfile.is_open())) { + return false; + } + + unsigned N = std::pow(2, log2n); + + cfile.read((char*)(C.data()), N * sizeof(Real)); + cfile.close(); + + rfile.read((char*)(R.data()), N * sizeof(Real)); + rfile.close(); + + ctfile.read((char*)(Ĉ.data()), N * sizeof(Complex)); + ctfile.close(); + + rtfile.read((char*)(Ȓ.data()), N * sizeof(Complex)); + rtfile.close(); + + return true; +} + +std::tuple<std::vector<Complex>, std::vector<Complex>> ΣD(LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, Real β, unsigned p, unsigned s, Real λ) { + std::vector<Real> D(C.size()); + std::vector<Real> Σ(C.size()); + Real β² = std::pow(β, 2); + for (unsigned n = 0; n < C.size(); n++) { + D[n] = β² * ∂f(λ, p, s, C[n]); + Σ[n] = β² * R[n] * ∂∂f(λ, p, s, C[n]); + } + std::vector<Complex> Σhat = fft.fourier(Σ, false); + std::vector<Complex> Dhat = fft.fourier(D, true); + + return {Σhat, Dhat}; +} + +Real estimateZ(LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Complex>& Ĉ, const std::vector<Real>& R, const std::vector<Complex>& Ȓ, unsigned p, unsigned s, Real λ, Real τ₀, Real β) { + auto [Σhat, Dhat] = ΣD(fft, C, R, β, p, s, λ); + Real Γ₀ = 1.0; + + return (((2 * Γ₀ + Dhat[0]) * std::conj(Ȓ[0]) + Σhat[0] * Ĉ[0]) / Ĉ[0]).real(); +} + +Real energy(const LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real β) { + unsigned n₀ = 0; + for (unsigned n = 0; n < C.size(); n++) { + if (C[n] > 1 || R[n] > 1) n₀ = n % 2 == 0 ? n / 2 : (n + 1) / 2; + } + Real E = fft.t(2*n₀) * ∂f(λ, p, s, 1); + for (unsigned n = n₀; n < C.size()/2-1; n++) { + Real R₂ₙ = R[2*n]; + Real R₂ₙ₊₁ = R[2*n+1]; + Real R₂ₙ₊₂ = R[2*n+2]; + Real C₂ₙ = C[2*n]; + Real C₂ₙ₊₁ = C[2*n+1]; + Real C₂ₙ₊₂ = C[2*n+2]; + + if (C₂ₙ₊₂ < 0 || R₂ₙ₊₂ < 0) break; + + 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₂ₙ * ∂f(λ, p, s, C₂ₙ); + Real f₂ₙ₊₁ = R₂ₙ₊₁ * ∂f(λ, p, s, C₂ₙ₊₁); + Real f₂ₙ₊₂ = R₂ₙ₊₂ * ∂f(λ, p, s, C₂ₙ₊₂); + + E += (h₂ₙ + h₂ₙ₊₁) / 6 * ( + (2 - h₂ₙ₊₁ / h₂ₙ) * f₂ₙ + + std::pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁ + + (2 - h₂ₙ / h₂ₙ₊₁) * f₂ₙ₊₂ + ); + } + return β * E; +} + +Real C0(const LogarithmicFourierTransform& fft, const std::vector<Complex>& Ĉ) { + Real C = 0; + for (unsigned n = 0; n < Ĉ.size()/2-1; n++) { + Real Ĉ₂ₙ = Ĉ[2*n].real(); + Real Ĉ₂ₙ₊₁ = Ĉ[2*n+1].real(); + Real Ĉ₂ₙ₊₂ = Ĉ[2*n+2].real(); + + 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₂ₙ = Ĉ₂ₙ; + Real f₂ₙ₊₁ = Ĉ₂ₙ₊₁; + Real f₂ₙ₊₂ = Ĉ₂ₙ₊₂; + + C += (h₂ₙ + h₂ₙ₊₁) / 6 * ( + (2 - h₂ₙ₊₁ / h₂ₙ) * f₂ₙ + + std::pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁ + + (2 - h₂ₙ / h₂ₙ₊₁) * f₂ₙ₊₂ + ); + } + return C * pow(fft.shift, 2) / M_PI; +} + +void smooth(std::vector<Real>& C, Real ε) { + unsigned N = C.size(); + bool trigger0 = false; + bool trigger1 = false; + for (unsigned i = 0; i < N; i++) { + if (C[i] < ε || trigger0) { + C[i] = 0; + trigger0 = true; + } + if (C[N-1-i] > 1 - ε || trigger1) { + C[N-1-i] = 1; + trigger1 = true; + } + } + + Real Cmax = 0; + for (unsigned i = 0; i < N; i++) { + if (C[N-1-i] > Cmax) Cmax = C[N-1-i]; + C[N-1-i] = Cmax; + } +} + diff --git a/log-fourier.hpp b/log-fourier.hpp new file mode 100644 index 0000000..755f7e9 --- /dev/null +++ b/log-fourier.hpp @@ -0,0 +1,58 @@ +#pragma once + +#include "types.hpp" + +#include <vector> +#include <tuple> + +#include <fftw3.h> +#include <gsl/gsl_sf_gamma.h> + +class LogarithmicFourierTransform { +private: + Complex* a; + Complex* â; + FFTW_PLAN a_to_â; + FFTW_PLAN â_to_a; + unsigned N; + unsigned pad; + Real k; + Real Δτ; + Real τₛ; + Real ωₛ; + Real sₛ; + std::vector<Real> ts; + std::vector<Real> νs; + std::vector<Complex> Γs; + std::vector<Real> exp1kω; + std::vector<Real> exp1kτ; + std::vector<Real> expkω; + std::vector<Real> expkτ; +public: + Real shift; + LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad = 4, Real shift = 0.5); + ~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<Real> inverse(const std::vector<Complex>& ĉ); +}; + +std::string logFourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real shift); + +void logFourierSave(const std::vector<Real>& C, const std::vector<Real>& R, const std::vector<Complex>& Ĉ, const std::vector<Complex>& Ȓ, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real shift); + +bool logFourierLoad(std::vector<Real>& C, std::vector<Real>& R, std::vector<Complex>& Ĉ, std::vector<Complex>& Ȓ, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real shift); + +std::tuple<std::vector<Complex>, std::vector<Complex>> ΣD(LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, Real β, unsigned p, unsigned s, Real λ); + +Real estimateZ(LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Complex>& Ĉ, const std::vector<Real>& Ȓ, const std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β); + +Real energy(const LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real β); + +Real C0(const LogarithmicFourierTransform& fft, const std::vector<Complex>& Ĉ); + +void smooth(std::vector<Real>& C, Real ε); diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp index 177db46..e13cf1b 100644 --- a/log-fourier_integrator.cpp +++ b/log-fourier_integrator.cpp @@ -1,25 +1,36 @@ -#include "fourier.hpp" +#include "log-fourier.hpp" + #include <getopt.h> #include <iostream> +#include <iomanip> int main(int argc, char* argv[]) { + /* Model parameters */ unsigned p = 2; unsigned s = 2; Real λ = 0.5; Real τ₀ = 0; + /* Log-Fourier parameters */ unsigned log2n = 8; Real Δτ = 0.1; - Real k = 0.1; + Real k = -0.01; + Real logShift = 0; + bool shiftSquare = false; - Real ε = 1e-16; - Real γ = 1; - Real βₘₐₓ = 0.7; + /* Iteration parameters */ + Real ε = 1e-15; + Real γ₀ = 1; + Real x = 1; + Real β₀ = 0; + Real βₘₐₓ = 20; Real Δβ = 0.01; + bool loadData = false; + unsigned pad = 2; int opt; - while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:g:k:D:")) != -1) { + while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:k:D:e:0:lg:x:P:h:S")) != -1) { switch (opt) { case 'p': p = atoi(optarg); @@ -40,14 +51,35 @@ int main(int argc, char* argv[]) { Δβ = atof(optarg); break; case 'g': - γ = atof(optarg); + γ₀ = atof(optarg); break; case 'k': k = atof(optarg); break; + case 'h': + logShift = atof(optarg); + break; case 'D': Δτ = atof(optarg); break; + case 'e': + ε = atof(optarg); + break; + case '0': + β₀ = atof(optarg); + break; + case 'x': + x = atof(optarg); + break; + case 'P': + pad = atoi(optarg); + break; + case 'l': + loadData = true; + break; + case 'S': + shiftSquare = true; + break; default: exit(1); } @@ -55,81 +87,151 @@ int main(int argc, char* argv[]) { unsigned N = pow(2, log2n); - LogarithmicFourierTransform fft(N, k, Δτ, 4); + Real Γ₀ = 1; + Real μ₀ = τ₀ > 0 ? (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀) : Γ₀; - Real Γ₀ = 1.0; - Real μ = Γ₀; - if (τ₀ > 0) { - μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀); - } + Real shift = μ₀ * pow(10, logShift); + if (shiftSquare) shift *= μ₀; + LogarithmicFourierTransform fft(N, k, Δτ, pad, shift); - std::vector<Real> C(N); - std::vector<Real> R(N); - std::vector<Complex> Ct(N); - std::vector<Complex> Rt(N); + std::cerr << "Starting, μ₀ = " << μ₀ << ", range " << fft.t(0) << " " << fft.t(N-1) << std::endl; - // start from the exact solution for β = 0 - 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)); + Real μₜ₋₁ = μ₀; - Ct[n] = 2 * Γ₀ / (pow(μ, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2)); - Rt[n] = 1.0 / (μ + 1i * fft.ν(n)); + std::vector<Real> Cₜ₋₁(N); + std::vector<Real> Rₜ₋₁(N); + std::vector<Complex> Ĉₜ₋₁(N); + std::vector<Complex> Ȓₜ₋₁(N); + std::vector<Real> Γ(N); + + for (unsigned n = 0; n < N; n++) { + Γ[n] = Γ₀ / (1 + std::pow(τ₀ * fft.ν(n), 2)); } - Real β = 0; - while (β < βₘₐₓ) { - Real ΔC = 100; - while (ΔC > ε) { - std::vector<Real> RddfC(N); - std::vector<Real> dfC(N); + if (!loadData) { + /* Start from the exact solution for β = 0 */ for (unsigned n = 0; n < N; n++) { - RddfC[n] = R[n] * ddf(λ, p, s, C[n]); - dfC[n] = df(λ, p, s, C[n]); + if (τ₀ > 0) { + if (τ₀ == 2) { + Cₜ₋₁[n] = Γ₀ * std::exp(-fft.t(n) / 2) * (1 + fft.t(n) / 2); + } else { + Cₜ₋₁[n] = Γ₀ * (std::exp(-μ₀ * fft.t(n)) / μ₀ - τ₀ * std::exp(-fft.t(n) / τ₀)) / (1 - pow(μ₀ * τ₀, 2)); + } + } else { + Cₜ₋₁[n] = Γ₀ * std::exp(-μ₀ * fft.t(n)) / μ₀; + } + Rₜ₋₁[n] = std::exp(-μ₀ * fft.t(n)); + + Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μ₀, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2)); + Ȓₜ₋₁[n] = (Real)1 / (μ₀ + II * fft.ν(n)); } - std::vector<Complex> RddfCt = fft.fourier(RddfC, false); - std::vector<Complex> dfCt = fft.fourier(dfC, true); + } else { + if (!logFourierLoad(Cₜ₋₁, Rₜ₋₁, Ĉₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀, log2n, Δτ, logShift)) { + return 1; + } + μₜ₋₁ = estimateZ(fft, Cₜ₋₁, Ĉₜ₋₁, Rₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀); + } - std::vector<Complex> Rtnew(N); - std::vector<Complex> Ctnew(N); + std::vector<Real> Cₜ = Cₜ₋₁; + std::vector<Real> Rₜ = Rₜ₋₁; + std::vector<Complex> Ĉₜ = Ĉₜ₋₁; + std::vector<Complex> Ȓₜ = Ȓₜ₋₁; + Real μₜ = μₜ₋₁; - for (unsigned n = 0; n < N; n++) { - Rtnew[n] = (1.0 + pow(β, 2) * RddfCt[n] * Rt[n]) / (μ + 1i * fft.ν(n)); - Ctnew[n] = (2 * Γ₀ * std::conj(Rt[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ct[n] + dfCt[n] * std::conj(Rt[n]))) / (μ + 1i * fft.ν(n)); - } + Real β = β₀ + Δβ; + while (β < βₘₐₓ) { + Real γ = γ₀; + Real ΔCₜ = 100; + Real ΔCₜ₋₁ = 101; + while (ΔCₜ > ε) { + auto [Σ, D] = ΣD(fft, Cₜ, Rₜ, β, p, s, λ); - std::vector<Real> Cnew = fft.inverse(Ctnew); - std::vector<Real> Rnew = fft.inverse(Rtnew); + std::vector<Complex> Ĉₜ₊₁(N); + std::vector<Complex> Ȓₜ₊₁(N); - ΔC = 0; - for (unsigned i = 0; i < N; i++) { - ΔC += std::norm(Ct[i] - Ctnew[i]); - ΔC += std::norm(Rt[i] - Rtnew[i]); - } - ΔC = sqrt(ΔC) / (2*N); + Real C₀ = 0; + Real μ₊ = 0; + Real μ₋ = 0; + + while (std::abs(C₀ - 1) > ε) { + for (unsigned n = 0; n < N; n++) { + Ĉₜ₊₁[n] = (((2 * Γ[n] + D[n]) * std::conj(Ȓₜ[n]) + Σ[n] * Ĉₜ[n]) / (μₜ + II * fft.ν(n))).real(); + } + C₀ = C0(fft, Ĉₜ₊₁); + if (C₀ > 1) { + μ₋ = μₜ; + } else { + μ₊ = μₜ; + } + if (μ₋ > 0 && μ₊ > 0) { + μₜ = (μ₊ + μ₋) / 2; + } else { + μₜ *= pow(tanh(C₀-1)+1, x); + } + } + + ΔCₜ = 0; + for (unsigned n = 0; n < N; n++) { + ΔCₜ += std::norm(((2 * Γ[n] + D[n]) * std::conj(Ȓₜ[n]) + Σ[n] * Ĉₜ[n]) - Ĉₜ[n] * (μₜ + II * fft.ν(n))); + ΔCₜ += std::norm(((Real)1 + Σ[n] * Ȓₜ[n]) - Ȓₜ[n] * (μₜ + II * fft.ν(n))); + } + ΔCₜ = sqrt(ΔCₜ) / (2*N); + + for (unsigned n = 0; n < N; n++) { + Ȓₜ₊₁[n] = ((Real)1 + Σ[n] * Ȓₜ[n]) / (μₜ + II * fft.ν(n)); + } + + std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁); + std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁); + + smooth(Cₜ₊₁, ε); + smooth(Rₜ₊₁, ε); - for (unsigned i = 0; i < N; i++) { - C[i] += γ * (Cnew[i] - C[i]); - R[i] += γ * (Rnew[i] - R[i]); - Ct[i] += γ * (Ctnew[i] - Ct[i]); - Rt[i] += γ * (Rtnew[i] - Rt[i]); + 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]); + } + + if (ΔCₜ > ΔCₜ₋₁ * 2 && ΔCₜ < 1e-10) { + γ = std::max(γ / 2, (Real)1e-2); + } + + ΔCₜ₋₁ = ΔCₜ; + + std::cerr << "\x1b[2K" << "\r"; + std::cerr << std::setprecision(7) << β << " " << Δβ << " " << μₜ << " " << ΔCₜ << " " << γ; } - μ *= C[0]; + if (std::isnan(ΔCₜ)) { + β -= Δβ; + Δβ *= 0.1; + β += Δβ; + Cₜ = Cₜ₋₁; + Rₜ = Rₜ₋₁; + Ĉₜ = Ĉₜ₋₁; + Ȓₜ = Ȓₜ₋₁; + μₜ = μₜ₋₁; + } else { + Real E = energy(fft, Cₜ, Rₜ, p, s, λ, β); + + std::cerr << "\x1b[2K" << "\r"; + std::cerr << std::setprecision(7) << β << " " << Δβ << " " << μₜ << " " << Ĉₜ[0].real() << " " << E << std::endl; -// std::cerr << ΔC << std::endl; - } + logFourierSave(Cₜ, Rₜ, Ĉₜ, Ȓₜ, p, s, λ, τ₀, β, log2n, Δτ, logShift); - std::cerr << β << " " << μ << " " << Ct[0].real() << std::endl; - β += Δβ; - } + if (Ĉₜ[0].real() / Ĉₜ₋₁[0].real() > 1.25) { + Δβ = std::round(1e6 * Δβ / 2) / 1e6; + } - for (unsigned i = 0; i < N; i++) { - std::cout << fft.t(i) << " " << Rt[i].imag() << std::endl; + β = std::round(1e6 * (β + Δβ)) / 1e6; + Cₜ₋₁ = Cₜ; + Rₜ₋₁ = Rₜ; + Ĉₜ₋₁ = Ĉₜ; + Ȓₜ₋₁ = Ȓₜ; + μₜ₋₁ = μₜ; + } } return 0; diff --git a/log_get_energy.cpp b/log_get_energy.cpp new file mode 100644 index 0000000..d156fd4 --- /dev/null +++ b/log_get_energy.cpp @@ -0,0 +1,100 @@ +#include "log-fourier.hpp" + +#include <getopt.h> +#include <iostream> +#include <iomanip> +#include <filesystem> + +int main(int argc, char* argv[]) { + /* Model parameters */ + unsigned p = 2; + unsigned s = 2; + Real λ = 0.5; + Real τ₀ = 0; + + /* Log-Fourier parameters */ + unsigned log2n = 8; + Real Δτ = 0.1; + Real k = 0.1; + unsigned pad = 2; + Real logShift = 0; + bool shiftSquare = false; + + /* Iteration parameters */ + Real β₀ = 0; + Real βₘₐₓ = 0.7; + Real Δβ = 0.01; + + int opt; + + while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:k:h:D:0:S")) != -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 'b': + βₘₐₓ = atof(optarg); + break; + case 'd': + Δβ = atof(optarg); + break; + case 'k': + k = atof(optarg); + break; + case 'h': + logShift = atof(optarg); + break; + case 'D': + Δτ = atof(optarg); + break; + case '0': + β₀ = atof(optarg); + break; + case 'S': + shiftSquare = true; + break; + default: + exit(1); + } + } + + unsigned N = pow(2, log2n); + Real Γ₀ = 1; + Real μ₀ = τ₀ > 0 ? (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀) : Γ₀; + + Real shift = μ₀ * pow(10, logShift); + if (shiftSquare) shift *= μ₀; + LogarithmicFourierTransform fft(N, k, Δτ, pad, shift); + + std::vector<Real> C(N); + std::vector<Real> R(N); + std::vector<Complex> Ĉ(N); + std::vector<Complex> Ȓ(N); + + Real β = β₀; + + std::cout << std::setprecision(16); + + while (β = std::round(1e6 * (β + Δβ)) / 1e6, β <= βₘₐₓ) { + if (std::filesystem::exists(logFourierFile("C", p, s, λ, τ₀, β, log2n, Δτ, logShift))) { + logFourierLoad(C, R, Ĉ, Ȓ, p, s, λ, τ₀, β, log2n, Δτ, logShift); + + Real e = energy(fft, C, R, p, s, λ, β); + + Real μ = estimateZ(fft, C, Ĉ, R, Ȓ, p, s, λ, τ₀, β); + + std::cout << p << " " << s << " " << λ << " " << τ₀ << " " << β << " " << μ << " " << Ĉ[0].real() << " " << e << " " << Ȓ[0].real() << std::endl; + } + } + + return 0; +} diff --git a/p-spin.cpp b/p-spin.cpp new file mode 100644 index 0000000..ba18c6a --- /dev/null +++ b/p-spin.cpp @@ -0,0 +1,25 @@ +#include "p-spin.hpp" + +inline Real fₚ(unsigned p, Real q) { + return 0.5 * pow(q, p); +} + +inline Real ∂fₚ(unsigned p, Real q) { + return 0.5 * p * pow(q, p - 1); +} + +inline Real ∂∂fₚ(unsigned p, Real q) { + return 0.5 * p * (p - 1) * pow(q, p - 2); +} + +Real f(Real λ, unsigned p, unsigned s, Real q) { + return (1 - λ) * fₚ(p, q) + λ * fₚ(s, q); +} + +Real ∂f(Real λ, unsigned p, unsigned s, Real q) { + return (1 - λ) * ∂fₚ(p, q) + λ * ∂fₚ(s, q); +} + +Real ∂∂f(Real λ, unsigned p, unsigned s, Real q) { + return (1 - λ) * ∂∂fₚ(p, q) + λ * ∂∂fₚ(s, q); +} diff --git a/p-spin.hpp b/p-spin.hpp new file mode 100644 index 0000000..484bc17 --- /dev/null +++ b/p-spin.hpp @@ -0,0 +1,7 @@ +#pragma once + +#include "types.hpp" + +Real f(Real λ, unsigned p, unsigned s, Real q); +Real ∂f(Real λ, unsigned p, unsigned s, Real q); +Real ∂∂f(Real λ, unsigned p, unsigned s, Real q); diff --git a/print_correlations.cpp b/print_correlations.cpp new file mode 100644 index 0000000..4627119 --- /dev/null +++ b/print_correlations.cpp @@ -0,0 +1,86 @@ +#include "log-fourier.hpp" + +#include <getopt.h> +#include <iostream> +#include <iomanip> +#include <filesystem> + +int main(int argc, char* argv[]) { + /* Model parameters */ + unsigned p = 2; + unsigned s = 2; + Real λ = 0.5; + Real τ₀ = 0; + + /* Log-Fourier parameters */ + unsigned log2n = 8; + Real Δτ = 0.1; + Real k = 0.1; + unsigned pad = 2; + Real logShift = 0; + bool shiftSquare = false; + + /* Iteration parameters */ + Real β = 0; + + int opt; + + while ((opt = getopt(argc, argv, "p:s:2:T:t:k:h:D:0:S")) != -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 'k': + k = atof(optarg); + break; + case 'h': + logShift = atof(optarg); + break; + case 'D': + Δτ = atof(optarg); + break; + case '0': + β = atof(optarg); + break; + case 'S': + shiftSquare = true; + break; + default: + exit(1); + } + } + + unsigned N = pow(2, log2n); + Real Γ₀ = 1; + Real μ₀ = τ₀ > 0 ? (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀) : Γ₀; + + Real shift = μ₀ * pow(10, logShift); + if (shiftSquare) shift *= μ₀; + LogarithmicFourierTransform fft(N, k, Δτ, pad, shift); + + std::vector<Real> C(N); + std::vector<Real> R(N); + std::vector<Complex> Ĉ(N); + std::vector<Complex> Ȓ(N); + + std::cout << std::setprecision(16); + + if (std::filesystem::exists(logFourierFile("C", p, s, λ, τ₀, β, log2n, Δτ, logShift))) { + logFourierLoad(C, R, Ĉ, Ȓ, p, s, λ, τ₀, β, log2n, Δτ, logShift); + + for (unsigned n = 0; n < N; n++) { + std::cout << fft.t(n) << " " << C[n] << " " << R[n] << " " << fft.ν(n) << " " << Ĉ[n].real() << " " << Ȓ[n].real() << " " << Ȓ[n].imag() << std::endl; + } + } + + return 0; +} diff --git a/types.hpp b/types.hpp new file mode 100644 index 0000000..71556fd --- /dev/null +++ b/types.hpp @@ -0,0 +1,43 @@ +#pragma once +#include <complex> +#include <fftw3.h> + +using namespace std::complex_literals; + +#ifndef QUAD_PRECISION + +#define FFTW_PLAN fftw_plan +#define FFTW_ALLOC_COMPLEX fftw_alloc_complex +#define FFTW_IMPORT_WISDOM fftw_import_wisdom_from_filename +#define FFTW_EXPORT_WISDOM fftw_export_wisdom_to_filename +#define FFTW_PLAN_DFT_1D fftw_plan_dft_1d +#define FFTW_DESTROY_PLAN fftw_destroy_plan +#define FFTW_FREE fftw_free +#define FFTW_CLEANUP fftw_cleanup +#define FFTW_EXECUTE fftw_execute +#define FFTW_COMPLEX fftw_complex + +using Real = double; + +#define II 1i + +#else + +#define FFTW_PLAN fftwl_plan +#define FFTW_ALLOC_COMPLEX fftwl_alloc_complex +#define FFTW_IMPORT_WISDOM fftwl_import_wisdom_from_filename +#define FFTW_EXPORT_WISDOM fftwl_export_wisdom_to_filename +#define FFTW_PLAN_DFT_1D fftwl_plan_dft_1d +#define FFTW_DESTROY_PLAN fftwl_destroy_plan +#define FFTW_FREE fftwl_free +#define FFTW_CLEANUP fftwl_cleanup +#define FFTW_EXECUTE fftwl_execute +#define FFTW_COMPLEX fftwl_complex + +using Real = long double; + +#define II 1il + +#endif + +using Complex = std::complex<Real>; |