From 2ede6db86243e59223cd89e6debd7107d02eabd5 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Sun, 20 Apr 2025 12:44:38 -0300 Subject: Standardized saving and loading of files --- .gitignore | 1 + Makefile | 5 ++- log-fourier.cpp | 89 ++++++++++++++++++++++++++++++++++++++++++++++ log-fourier.hpp | 11 ++++++ log-fourier_integrator.cpp | 61 ++----------------------------- log_get_energy.cpp | 82 ++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 190 insertions(+), 59 deletions(-) create mode 100644 log_get_energy.cpp diff --git a/.gitignore b/.gitignore index f6c3bea..fe5f9c3 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,4 @@ walk get_energy integrator fftw.wisdom +log_get_energy diff --git a/Makefile b/Makefile index a1d0907..6cf8a52 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,4 @@ -all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator +all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator log_get_energy CC := clang++ -std=c++17 -flto -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -DFFTW_THREADS=1 @@ -28,3 +28,6 @@ log-fourier_integrator: log-fourier_integrator.cpp log-fourier.hpp log-fourier.o 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 diff --git a/log-fourier.cpp b/log-fourier.cpp index 7461a70..f7d6f4b 100644 --- a/log-fourier.cpp +++ b/log-fourier.cpp @@ -1,5 +1,7 @@ #include "log-fourier.hpp" +#include "p-spin.hpp" #include +#include LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad) : N(N), pad(pad), k(k), Δτ(Δτ) { τₛ = -0.5 * N; @@ -106,3 +108,90 @@ std::vector LogarithmicFourierTransform::inverse(const std::vector& C, const std::vector& R, const std::vector& Ct, const std::vector& 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& C, std::vector& R, std::vector& Ct, std::vector& 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> RddfCtdfCt(LogarithmicFourierTransform& fft, const std::vector& C, const std::vector& R, unsigned p, unsigned s, Real λ) { + std::vector dfC(C.size()); + std::vector 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 RddfCt = fft.fourier(RddfC, false); + std::vector dfCt = fft.fourier(dfC, true); + + return {RddfCt, dfCt}; +} + +Real estimateZ(LogarithmicFourierTransform& fft, const std::vector& C, const std::vector& Ct, const std::vector& R, const std::vector& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β) { + auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ); + Real Γ₀ = 1 + τ₀; + + 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& C, const std::vector& 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; +} + diff --git a/log-fourier.hpp b/log-fourier.hpp index f57c6bc..b11db0f 100644 --- a/log-fourier.hpp +++ b/log-fourier.hpp @@ -31,3 +31,14 @@ public: std::vector inverse(const std::vector& ĉ); }; +std::string logFourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k); + +void logFourierSave(const std::vector& C, const std::vector& R, const std::vector& Ct, const std::vector& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k); + +bool logFourierLoad(std::vector& C, std::vector& R, std::vector& Ct, std::vector& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k); + +std::tuple, std::vector> RddfCtdfCt(LogarithmicFourierTransform& fft, const std::vector& C, const std::vector& R, unsigned p, unsigned s, Real λ); + +Real estimateZ(LogarithmicFourierTransform& fft, const std::vector& C, const std::vector& Ct, const std::vector& R, const std::vector& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β); + +Real energy(const LogarithmicFourierTransform& fft, std::vector& C, const std::vector& R, unsigned p, unsigned s, Real λ, Real β); diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp index a032401..cf9819a 100644 --- a/log-fourier_integrator.cpp +++ b/log-fourier_integrator.cpp @@ -1,32 +1,6 @@ #include "log-fourier.hpp" -#include "p-spin.hpp" #include #include -#include - -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"; -} - -std::tuple, std::vector> RddfCtdfCt(LogarithmicFourierTransform& fft, const std::vector& C, const std::vector& R, unsigned p, unsigned s, Real λ) { - std::vector dfC(C.size()); - std::vector 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 RddfCt = fft.fourier(RddfC, false); - std::vector dfCt = fft.fourier(dfC, true); - - return {RddfCt, dfCt}; -} - -Real estimateZ(LogarithmicFourierTransform& fft, const std::vector& C, const std::vector& Ct, const std::vector& R, const std::vector& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β) { - auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ); - Real Γ₀ = 1 + τ₀ / 2; - - return ((2 * Γ₀ * std::conj(Rt[0]) + pow(β, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real(); -} int main(int argc, char* argv[]) { /* Model parameters */ @@ -119,16 +93,7 @@ int main(int argc, char* argv[]) { Ȓₜ₋₁[n] = 1.0 / (μₜ₋₁ + 1i * fft.ν(n)); } } else { - std::ifstream cfile(logFourierFile("C", p, s, λ, τ₀, β₀, log2n, Δτ, k), std::ios::binary); - cfile.read((char*)(Cₜ₋₁.data()), N * sizeof(Real)); - cfile.close(); - std::ifstream rfile(logFourierFile("R", p, s, λ, τ₀, β₀, log2n, Δτ, k), std::ios::binary); - rfile.read((char*)(Rₜ₋₁.data()), N * sizeof(Real)); - rfile.close(); - - Ĉₜ₋₁ = fft.fourier(Cₜ₋₁, true); - Ȓₜ₋₁ = fft.fourier(Rₜ₋₁, false); - + logFourierLoad(Cₜ₋₁, Rₜ₋₁, Ĉₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀, log2n, Δτ, k); μₜ₋₁ = estimateZ(fft, Cₜ₋₁, Ĉₜ₋₁, Rₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀); } @@ -181,32 +146,12 @@ int main(int argc, char* argv[]) { μₜ = μₜ₋₁; γ /= 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 *= β; + Real E = energy(fft, Cₜ, Rₜ, p, s, λ, β); 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(); + logFourierSave(Cₜ, Rₜ, Ĉₜ, Ȓₜ, p, s, λ, τ₀, β, log2n, Δτ, k); β += Δβ; Cₜ₋₁ = Cₜ; diff --git a/log_get_energy.cpp b/log_get_energy.cpp new file mode 100644 index 0000000..7a531fd --- /dev/null +++ b/log_get_energy.cpp @@ -0,0 +1,82 @@ +#include "log-fourier.hpp" +#include +#include +#include + +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 C(N); + std::vector R(N); + std::vector Ct(N); + std::vector Rt(N); + + Real β = β₀; + + 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 << β << " " << μ << " " << Ct[0].real() << " " << e << std::endl; + } + } + + return 0; +} -- cgit v1.2.3-70-g09d2