+log-fourier_integrator

+log_get_energy

+all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator log_get_energy

+integrator: integrator.cpp fourier.o p-spin.o fourier.hpp

+ $(CC) integrator.cpp fourier.o p-spin.o -lfftw3 -lgsl -o integrator

+fourier.o: fourier.cpp fourier.hpp types.hpp p-spin.o

+p-spin.o: p-spin.cpp p-spin.hpp types.hpp

+ $(CC) p-spin.cpp -c -o p-spin.o

+log-fourier.o: log-fourier.cpp log-fourier.hpp types.hpp

+ $(CC) log-fourier.cpp -c -o log-fourier.o

+fourier_integrator: fourier_integrator.cpp fourier.hpp fourier.o p-spin.o

+ $(CC) fourier_integrator.cpp fourier.o p-spin.o -lfftw3 -lgsl -o fourier_integrator

+

+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

+

+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

+

+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

+#include "p-spin.hpp"

+ for (unsigned n = 0; n < C.size() / 4 -1; n++) {

+ Real h₂ₙ = M_PI * Δτ;

+ Real h₂ₙ₊₁ = M_PI * Δτ;

+ 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₂ₙ₊₂

+ );

+ return y * e;

+#include "types.hpp"

+

+ Real β₀ = 0;

+ Real βₘₐₓ = 0.5;

+ Real Δβ = 0.05;

+ β₀ = atof(optarg);

+ βₘₐₓ = atof(optarg);

+ Δβ = atof(optarg);

+ Real Γ₀ = 1.0;

+ Real μ = Γ₀;

+ if (τ₀ > 0) {

+ μ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);

+ }

+ Real β = β₀;

+ C[i] = Γ₀ * (exp(-μ * τ) - μ * τ₀ * exp(-τ / τ₀)) / (μ - pow(μ, 3) * pow(τ₀, 2));

+ C[i] = Γ₀ * exp(-μ * τ) / μ;

+ R[i] = exp(-μ * τ);

+ std::ifstream cfile(fourierFile("C", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::binary);

+ std::ifstream rfile(fourierFile("R", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::binary);

+ μ = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, β);

+ std::vector<Real> Cb = C;

+ std::vector<Real> Rb = R;

+ std::vector<Complex> Ctb = Ct;

+ std::vector<Complex> Rtb = Rt;

+ Real μb = μ;

+

+ Real fac = 1;

+ while (β <= βₘₐₓ) {

+ Rt[i] = (1.0 + pow(β, 2) * RddfCt[i] * Rt[i]) / (μ + 1i * ω);

+ Ct[i] = (2 * Γ₀ * std::conj(Rt[i]) / (1 + pow(τ₀ * ω, 2)) + pow(β, 2) * (RddfCt[i] * Ct[i] + dfCt[i] * std::conj(Rt[i]))) / (μ + 1i * ω);

+ μ *= pow(tanh(C[0]-1)+1, 0.05);

+

+ std::cerr << "\x1b[2K" << "\r";

+ std::cerr << μ << " " << p << " " << s << " " << τ₀ << " " << β << " " << sqrt(2 * ΔC / C.size()) << " " << γ << " " << C[0];

+ if (std::isnan(C[0])) {

+ Δβ /= 2;

+ β -= Δβ;

+ C = Cb;

+ R = Rb;

+ Ct = Ctb;

+ Rt = Rtb;

+ μ = μb;

+ } else {

+

+ std::cerr << "\x1b[2K" << "\r";

+ Real e = energy(C, R, p, s, λ, β, Δτ);

+ std::cerr << "y " << β << " " << e << " " << μ << std::endl;

+

+ std::ofstream outfile(fourierFile("C", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);

+ std::ofstream outfileR(fourierFile("R", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);

+ β += Δβ;

+ Cb = C;

+ Rb = R;

+ Rtb = Rt;

+ Ctb = Ct;

+ μb = μ;

+ }

+#include <iomanip>

+ std::cout << std::setprecision(15);

+

+#include "p-spin.hpp"

+#include "log-fourier.hpp"

+#include "p-spin.hpp"

+#include <complex>

+#include <fstream>

+

+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 Γ(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};

+ /* 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] * 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) * Γ(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)(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) {

+ if (σ < 0) {

+ a[n] = std::conj(ĉ[n]) * exp((1 - k) * ω(n));

+ } else {

+ 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) * Γ(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 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";

+}

+

+void logFourierSave(const std::vector<Real>& C, const std::vector<Real>& R, const std::vector<Complex>& Ct, const std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k) {

+ unsigned N = 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*)(Ct.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*)(Rt.data()), N * sizeof(Complex));

+ outfileRt.close();

+}

+

+bool logFourierLoad(std::vector<Real>& C, std::vector<Real>& R, std::vector<Complex>& Ct, std::vector<Complex>& Rt, 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 = 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*)(Ct.data()), N * sizeof(Complex));

+ ctfile.close();

+

+ rtfile.read((char*)(Rt.data()), N * sizeof(Complex));

+ rtfile.close();

+

+ return true;

+}

+

+std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ) {

+ std::vector<Real> dfC(C.size());

+ std::vector<Real> RddfC(C.size());

+ for (unsigned n = 0; n < C.size(); n++) {

+ RddfC[n] = R[n] * ddf(λ, p, s, C[n]);

+ dfC[n] = df(λ, p, s, C[n]);

+ }

+ std::vector<Complex> RddfCt = fft.fourier(RddfC, false);

+ std::vector<Complex> dfCt = fft.fourier(dfC, true);

+

+ return {RddfCt, dfCt};

+}

+

+Real estimateZ(LogarithmicFourierTransform& 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 β) {

+ auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ);

+ Real Γ₀ = 1.0;

+

+ return ((2 * Γ₀ * std::conj(Rt[0]) + pow(β, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real();

+}

+

+Real energy(const LogarithmicFourierTransform& fft, std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real β) {

+ Real E = 0;

+ for (unsigned n = 0; n < C.size()/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₂ₙ₊₂

+ );

+ }

+ return β * E;

+}

+

+#pragma once

+#include "types.hpp"

+

+#include <cmath>

+#include <fftw3.h>

+#include <vector>

+#include <tuple>

+#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ₛ;

+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>& ĉ);

+};

+

+std::string logFourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k);

+

+void logFourierSave(const std::vector<Real>& C, const std::vector<Real>& R, const std::vector<Complex>& Ct, const std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k);

+

+bool logFourierLoad(std::vector<Real>& C, std::vector<Real>& R, std::vector<Complex>& Ct, std::vector<Complex>& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k);

+

+std::tuple<std::vector<Complex>, std::vector<Complex>> RddfCtdfCt(LogarithmicFourierTransform& fft, const std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ);

+

+Real estimateZ(LogarithmicFourierTransform& 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 β);

+

+Real energy(const LogarithmicFourierTransform& fft, std::vector<Real>& C, const std::vector<Real>& R, unsigned p, unsigned s, Real λ, Real β);

+#include "log-fourier.hpp"

+ /* Model parameters */

+ /* Log-Fourier parameters */

+ /* Iteration parameters */

+ Real ε = 1e-14;

+ Real γ₀ = 1;

+ Real β₀ = 0;

+ bool loadData = false;

+ unsigned stepsToRespond = 1000;

+ while ((opt = getopt(argc, argv, "p:s:2:T:t:b:d:g:k:D:e:0:lS:")) != -1) {

+ γ₀ = atof(optarg);

+ case 'e':

+ ε = atof(optarg);

+ break;

+ case '0':

+ β₀ = atof(optarg);

+ break;

+ case 'l':

+ loadData = true;

+ break;

+ case 'S':

+ stepsToRespond = atoi(optarg);

+ break;

+ Real μₜ₋₁ = Γ₀;

+ μₜ₋₁ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀);

+ std::vector<Real> Cₜ₋₁(N);

+ std::vector<Real> Rₜ₋₁(N);

+ std::vector<Complex> Ĉₜ₋₁(N);

+ std::vector<Complex> Ȓₜ₋₁(N);

+ if (!loadData) {

+ /* Start from the exact solution for β = 0 */

+ for (unsigned n = 0; n < N; n++) {

+ if (τ₀ > 0) {

+ if (τ₀ == 2) {

+ Cₜ₋₁[n] = Γ₀ * exp(-fft.t(n) / 2) * (1 + fft.t(n) / 2);

+ } else {

+ 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));

+

+ Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μₜ₋₁, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2));

+ Ȓₜ₋₁[n] = 1.0 / (μₜ₋₁ + 1i * fft.ν(n));

+ } else {

+ logFourierLoad(Cₜ₋₁, Rₜ₋₁, Ĉₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀, log2n, Δτ, k);

+ μₜ₋₁ = estimateZ(fft, Cₜ₋₁, Ĉₜ₋₁, Rₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀);

+ std::vector<Real> Cₜ = Cₜ₋₁;

+ std::vector<Real> Rₜ = Rₜ₋₁;

+ std::vector<Complex> Ĉₜ = Ĉₜ₋₁;

+ std::vector<Complex> Ȓₜ = Ȓₜ₋₁;

+ Real μₜ = μₜ₋₁;

+

+ Real β = β₀ + Δβ;

+ Real γ = γ₀;

+ Real ΔCmin = 1000;

+ Real ΔCₜ = 100;

+ unsigned stepsUp = 0;

+ while (ΔCₜ > ε) {

+ auto [RddfCt, dfCt] = RddfCtdfCt(fft, Cₜ, Rₜ, p, s, λ);

+

+ 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));

+ }

+ std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);

+ std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);

+

+ μₜ *= pow(tanh(Cₜ₊₁[0]-1)+1, 0.5);

+

+ ΔCₜ = 0;

+ for (unsigned i = 0; i < N; i++) {

+ ΔCₜ += std::norm(Cₜ[i] - Cₜ₊₁[i]);

+ ΔCₜ += std::norm(Rₜ[i] - Rₜ₊₁[i]);

+ }

+ ΔCₜ = sqrt(ΔCₜ) / (2*N);

+

+ if (ΔCₜ < 0.9 * ΔCmin) {

+ ΔCmin = ΔCₜ;

+ stepsUp = 0;

+ } else {

+ stepsUp++;

+ }

+

+ if (stepsUp > stepsToRespond) {

+ γ = std::max(γ/2, 1e-4);

+ stepsUp = 0;

+ ΔCmin = ΔCₜ;

+ }

+

+ 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";

+ std::cerr << β << " " << μₜ << " " << ΔCₜ << " " << γ << " " << Cₜ[0];

+ if (std::isnan(Cₜ[0])) {

+ γ₀ /= 2;

+ Cₜ = Cₜ₋₁;

+ Rₜ = Rₜ₋₁;

+ Ĉₜ = Ĉₜ₋₁;

+ Ȓₜ = Ȓₜ₋₁;

+ μₜ = μₜ₋₁;

+ } else {

+ Real E = energy(fft, Cₜ, Rₜ, p, s, λ, β);

+ std::cerr << "\x1b[2K" << "\r";

+ std::cerr << β << " " << μₜ << " " << Ĉₜ[0].real() << " " << E << " " << γ << std::endl;

+ logFourierSave(Cₜ, Rₜ, Ĉₜ, Ȓₜ, p, s, λ, τ₀, β, log2n, Δτ, k);

+ β += Δβ;

+ Cₜ₋₁ = Cₜ;

+ Rₜ₋₁ = Rₜ;

+ Ĉₜ₋₁ = Ĉₜ;

+ Ȓₜ₋₁ = Ȓₜ;

+ μₜ₋₁ = μₜ;

+#include "log-fourier.hpp"

+#include <getopt.h>

+#include <iomanip>

+#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;

+

+ /* 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:D:0:")) != -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 'D':

+ Δτ = atof(optarg);

+ break;

+ case '0':

+ β₀ = atof(optarg);

+ break;

+ default:

+ exit(1);

+ }

+ }

+

+ unsigned N = pow(2, log2n);

+

+ LogarithmicFourierTransform fft(N, k, Δτ, 4);

+

+ std::vector<Real> C(N);

+ std::vector<Real> R(N);

+ std::vector<Complex> Ct(N);

+ std::vector<Complex> Rt(N);

+

+ Real β = β₀;

+

+ std::cout << std::setprecision(16);

+

+ while (β += Δβ, β <= βₘₐₓ) {

+ if (logFourierLoad(C, R, Ct, Rt, p, s, λ, τ₀, β, log2n, Δτ, k)) {

+

+ Real e = energy(fft, C, R, p, s, λ, β);

+

+ Real μ = estimateZ(fft, C, Ct, R, Rt, p, s, λ, τ₀, β);

+

+ std::cout << p << " " << s << " " << λ << " " << τ₀ << " " << β << " " << μ << " " << Ct[0].real() << " " << e << " " << Rt[0].real() << std::endl;

+ }

+ }

+

+ return 0;

+}

+#include "p-spin.hpp"

+

+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);

+}

+#pragma once

+#include "types.hpp"

+

+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);

+#pragma once

+#include <complex>

+

+using Real = double;

+using Complex = std::complex<Real>;

+using namespace std::complex_literals;