From 99b8f462ef45f46234d3eeb58d7ee041ad63d643 Mon Sep 17 00:00:00 2001 From: GianlucaFuwa Date: Wed, 11 Dec 2024 00:56:38 +0100 Subject: [PATCH] Fix some flow stuff and add adaptive sigma0 for opes --- MetaAnalysis/src/bias.jl | 68 ++++++++++++------------- src/bias/Bias.jl | 12 ++++- src/bias/metadynamics.jl | 2 + src/bias/opes.jl | 10 +++- src/bias/parametric.jl | 2 + src/bias/weights.jl | 2 + src/main/Main.jl | 3 +- src/main/runbuild.jl | 30 +++++++---- src/measurements/measurement_methods.jl | 9 ++-- test/test_reversibility.jl | 65 +++++++++++++---------- 10 files changed, 123 insertions(+), 80 deletions(-) diff --git a/MetaAnalysis/src/bias.jl b/MetaAnalysis/src/bias.jl index 159e0ae..5f12822 100644 --- a/MetaAnalysis/src/bias.jl +++ b/MetaAnalysis/src/bias.jl @@ -178,6 +178,40 @@ struct Kernel σ::Float64 end +function (k::Kernel)(s, cutoff², penalty) + return evaluate_kernel(s, k.height, k.center, k.σ, cutoff², penalty) +end + +@inline function evaluate_kernel(s, height, center, σ, cutoff², penalty) + diff = (center - s) / σ + diff² = diff^2 + out = ifelse(diff² >= cutoff², 0.0, height * (exp(-0.5diff²) - penalty)) + return out +end + +function derivative(k::Kernel, s, cutoff², penalty) + return kernel_derivative(s, k.height, k.center, k.σ, cutoff², penalty) +end + +@inline function kernel_derivative(s, height, center, σ, cutoff², penalty) + diff = (center - s) / σ + diff² = diff^2 + val = ifelse(diff² >= cutoff², 0.0, height * (exp(-0.5diff²) - penalty)) + out = -diff / σ * val + return out +end + +Base.:*(c::Real, k::Kernel) = Kernel(c * k.height, k.center, k.σ) + +function merge(k::Kernel, other::Kernel) # Kernel merger + h = k.height + other.height + c = (k.height * k.center + other.height * other.center) / h + s_my_part = k.height * (k.σ^2 + k.center^2) + s_other_part = other.height * (other.σ^2 + other.center^2) + s² = (s_my_part + s_other_part) / h - c^2 + return Kernel(h, c, sqrt(s²)) +end + mutable struct OPES is_first_step::Bool @@ -362,37 +396,3 @@ function opes_from_file!(dict, usebias) return kernels, length(kernels) end end - -function (k::Kernel)(s, cutoff², penalty) - return evaluate_kernel(s, k.height, k.center, k.σ, cutoff², penalty) -end - -@inline function evaluate_kernel(s, height, center, σ, cutoff², penalty) - diff = (center - s) / σ - diff² = diff^2 - out = ifelse(diff² >= cutoff², 0.0, height * (exp(-0.5diff²) - penalty)) - return out -end - -function derivative(k::Kernel, s, cutoff², penalty) - return kernel_derivative(s, k.height, k.center, k.σ, cutoff², penalty) -end - -@inline function kernel_derivative(s, height, center, σ, cutoff², penalty) - diff = (center - s) / σ - diff² = diff^2 - val = ifelse(diff² >= cutoff², 0.0, height * (exp(-0.5diff²) - penalty)) - out = -diff / σ * val - return out -end - -Base.:*(c::Real, k::Kernel) = Kernel(c * k.height, k.center, k.σ) - -function merge(k::Kernel, other::Kernel) # Kernel merger - h = k.height + other.height - c = (k.height * k.center + other.height * other.center) / h - s_my_part = k.height * (k.σ^2 + k.center^2) - s_other_part = other.height * (other.σ^2 + other.center^2) - s² = (s_my_part + s_other_part) / h - c^2 - return Kernel(h, c, sqrt(s²)) -end diff --git a/src/bias/Bias.jl b/src/bias/Bias.jl index d6e0b3c..90422b3 100644 --- a/src/bias/Bias.jl +++ b/src/bias/Bias.jl @@ -55,7 +55,9 @@ function Bias(p::ParameterSet, U; mpi_multi_sim=false, instance=mpi_myrank(), du @level1("- Constructing Bias instance $(inum)...") kind_of_bias = Unicode.normalize(p.kind_of_bias; casefold=true) TCV = get_cvtype_from_parameters(p) - smearing = StoutSmearing(U; numlayers=p.numsmears_for_cv, rho=p.rhostout_for_cv) + numsmears = p.numsmears_for_cv + rho = p.rhostout_for_cv + smearing = StoutSmearing(U; numlayers=maximum(numsmears), rho=rho) is_static = dummy ? true : (inum==0 ? false : p.is_static[inum]) sstr = (is_static || kind_of_bias == "parametric") ? "static" : "dynamic" @level1("| Type: $(sstr) $(kind_of_bias)") @@ -70,7 +72,11 @@ function Bias(p::ParameterSet, U; mpi_multi_sim=false, instance=mpi_myrank(), du error("kind_of_bias $(kind_of_bias) not supported. Try metad, opes or parametric") end - @level1("| CV: $(string(TCV)) with $(string(smearing))") + for i in eachindex(p.numsmears_for_cv) + @level1 """ + | CV$(i): $(string(TCV)) with StoutSmearing(numlayers=$(numsmears[i]), rho=$(rho)) + """ + end if !(bias isa Parametric) is_opes = bias isa OPES @@ -157,6 +163,8 @@ kind_of_cv(b::Bias) = b.kind_of_cv update_bias!(::NoBias, args...; kwargs...) = nothing update_bias!(::Nothing, args...; kwargs...) = nothing write_to_file(::AbstractBias, args...) = nothing +is_adaptive(b::Bias) = is_adaptive(b.bias) +set_σ₀!(b::Bias, val) = set_σ₀!(b.bias, val) include("metadynamics.jl") include("opes.jl") diff --git a/src/bias/metadynamics.jl b/src/bias/metadynamics.jl index 77fae09..03f6356 100644 --- a/src/bias/metadynamics.jl +++ b/src/bias/metadynamics.jl @@ -115,6 +115,8 @@ end Base.length(m::Metadynamics) = length(m.values) Base.eachindex(m::Metadynamics) = eachindex(m.values) Base.lastindex(m::Metadynamics) = lastindex(m.values) +is_adaptive(::Metadynamics) = false +set_σ₀!(::Metadynamics, ::Any) = nothing function Base.setindex!(m::Metadynamics, v, i) return m.values[i] = v diff --git a/src/bias/opes.jl b/src/bias/opes.jl index 1c23b75..3b37bab 100644 --- a/src/bias/opes.jl +++ b/src/bias/opes.jl @@ -131,7 +131,8 @@ function OPES(; @level1("| SIGMA0: $(σ₀)") @assert σ₀ >= 0 "SIGMA0 must be >= 0" @level1("| SIGMA_MIN: $(σ_min)") - @assert σ_min >= 0 "SIGMA_MIN must be > 0" + @assert σ_min > 0 "SIGMA_MIN must be > 0" + σ₀ > 0 && (@assert σ₀ >= σ_min "If SIGMA0 is > 0 then it must be bigger than SIGMA_MIN") @level1("| FIXED_SIGMA: $(fixed_σ)") @level1("| EPSILON: $(ϵ)") @assert ϵ > 0 "EPSILON must be > 0, maybe your BARRIER is to high?" @@ -312,6 +313,13 @@ end get_kernels(o::OPES) = o.kernels get_δkernels(o::OPES) = o.δkernels +is_adaptive(o::OPES) = (o.σ₀ == 0) + +function set_σ₀!(o::OPES, val) + o.σ₀ = val + @level1("| sigma0 is adaptively set to $(o.σ₀)!") + return nothing +end update!(o::OPES, cv, itrj) = update_opes!(o, cv, itrj) diff --git a/src/bias/parametric.jl b/src/bias/parametric.jl index d31c0b7..cd40bf4 100644 --- a/src/bias/parametric.jl +++ b/src/bias/parametric.jl @@ -37,6 +37,8 @@ function Parametric(p::ParameterSet; dummy=false) return Parametric(cvlims, penalty_weight, Q, A, Z) end +is_adaptive(::Parametric) = false +set_σ₀!(::Parametric, ::Any) = nothing update!(::Parametric, cv, args...) = nothing clear!(::Parametric) = nothing diff --git a/src/bias/weights.jl b/src/bias/weights.jl index a3a7bf3..e27b9d4 100644 --- a/src/bias/weights.jl +++ b/src/bias/weights.jl @@ -66,6 +66,8 @@ function calc_weight(o::OPES, cv, args...) end function calc_weight(m::Metadynamics, cv, weight_method) + w = 0.0 + if weight_method == "tiwari" # average over exp(V) in denom w = calc_weight_tiwari(m, cv) elseif weight_method == "balanced_exp" # average over V in denom diff --git a/src/main/Main.jl b/src/main/Main.jl index 60ab216..00c087f 100644 --- a/src/main/Main.jl +++ b/src/main/Main.jl @@ -4,10 +4,11 @@ using Dates using DelimitedFiles using InteractiveUtils using Random +using Statistics using ..MetaIO using ..Utils -import ..BiasModule: NoBias, calc_weights, recalc_CV!, update_bias! +import ..BiasModule: NoBias, calc_weights, is_adaptive, recalc_CV!, update_bias!, set_σ₀! import ..DiracOperators: QuenchedFermionAction, fermaction_from_str import ..Fields: calc_gauge_action, is_distributed, normalize! import ..Measurements: MYEXT_str diff --git a/src/main/runbuild.jl b/src/main/runbuild.jl index 6d3e9fa..39700b1 100644 --- a/src/main/runbuild.jl +++ b/src/main/runbuild.jl @@ -84,13 +84,15 @@ function build_bias!(univ, parameters, updatemethod; mpi_multi_sim=false) additional_string=additional_string, ) - measurements_with_flow = MeasurementMethods( - U, - parameters.measure_dir, - parameters.measurements_with_flow; - additional_string=additional_string, - flow=true, - ) + measurements_with_flow = ntuple(length(gflow)) do i + MeasurementMethods( + U, + parameters.measure_dir, + parameters.measurements_with_flow; + additional_string=additional_string, + flow=gflow[i], + ) + end checkpointer = Checkpointer( parameters.ensemble_dir, parameters.save_checkpoint_every @@ -124,7 +126,8 @@ function metabuild!( bias = univ.bias comm = mpi_comm() starting_Q = parameters.starting_Q - mpi_barrier() + therm_cv = Vector{Float64}(undef, parameters.numtherm) + adaptive_σ = is_adaptive(bias) @level1("- Thermalization:") _, runtime_therm = @timed begin @@ -140,10 +143,13 @@ function metabuild!( bias=NoBias(), metro_test=itrj>10, # So we dont get stuck at the beginning therm=true, - mpi_multi_sim=mpi_multi_sim, ) end + if adaptive_σ + recalc_CV!(U, bias) + therm_cv[itrj] = U.CV + end @level1("| Elapsed time:\t$(updatetime) [s] @ $(string(current_time()))") end end @@ -153,6 +159,11 @@ function metabuild!( mpi_barrier() + if adaptive_σ + std_cv = mpi_allgather(std(therm_cv)::Float64, comm) + set_σ₀!(bias, mean(std_cv)) + end + @level1("- Production:") _, runtime_prod = @timed begin numaccepts = 0.0 @@ -166,7 +177,6 @@ function metabuild!( fermion_action=fermion_action, bias=bias, metro_test=true, - mpi_multi_sim=mpi_multi_sim, ) numaccepts += accepted end diff --git a/src/measurements/measurement_methods.jl b/src/measurements/measurement_methods.jl index 327352a..eabf53b 100644 --- a/src/measurements/measurement_methods.jl +++ b/src/measurements/measurement_methods.jl @@ -88,11 +88,12 @@ end function calc_measurements_flowed( - m::MeasurementMethods, flow::Tuple, U, itrj, myinstance=mpi_myrank(); - mpi_multi_sim=false + m::Tuple, flow::Tuple, U, itrj, myinstance=mpi_myrank(); mpi_multi_sim=false ) - for f in flow - calc_measurements_flowed(m, f, U, itrj, myinstance; mpi_multi_sim=mpi_multi_sim) + for i in eachindex(flow) + calc_measurements_flowed( + m[i], flow[i], U, itrj, myinstance; mpi_multi_sim=mpi_multi_sim + ) end return nothing diff --git a/test/test_reversibility.jl b/test/test_reversibility.jl index 7cb552a..91c61bb 100644 --- a/test/test_reversibility.jl +++ b/test/test_reversibility.jl @@ -3,35 +3,40 @@ using Random using Printf using TimerOutputs -function test_reversibility(; - integrator=Leapfrog(), gaction=WilsonGaugeAction, faction=nothing, Nf=2, with_bias=false +function test_reversibility( + ; + integrator=Leapfrog(), + gaction=WilsonGaugeAction, + faction=QuenchedFermionAction, + Nf=2, + with_bias=false, + hmc_trajectory=1, + hmc_steps=100, ) - println("┌ Testing reversibility of $integrator") - println("| gauge action: $gaction") - println("| fermion action: $faction") - println("| bias enabled: $(with_bias)") - println("└\n") - MetaQCD.MetaIO.set_global_logger!(1, devnull; tc=true) + MetaQCD.MetaIO.set_global_logger!(1, nothing; tc=true) Random.seed!(123) - N = 4 + N = 12 U = Gaugefield{CPU,Float64,gaction}(N, N, N, N, 6.0) - fermion = if faction === nothing - nothing + fermion = if isnothing(faction) || faction === QuenchedFermionAction + QuenchedFermionAction() else (faction(U, 0.01; Nf=Nf, cg_tol_md=1e-7, cg_tol_action=1e-9),) end + println("┌ Testing reversibility of $integrator") + println("| gauge action: $gaction") + println("| fermion action: $faction") + println("| bias enabled: $(with_bias)") + println("└\n") to = TimerOutput() - @timeit to "load_config!" load_config!(BridgeFormat(), U, "./test/testconf.txt") - - hmc_trajectory = 1 - hmc_steps = 100 + identity_gauges!(U) + # @timeit to "load_config!" load_config!(BridgeFormat(), U, "./test/testconf.txt") bias = if with_bias Bias( Clover(), - StoutSmearing(U, 4, 0.12), + StoutSmearing(U; numlayers=4, rho=0.12), true, Parametric((-5, 5), 10, 0, 100, 1.4), nothing, @@ -40,7 +45,7 @@ function test_reversibility(; 0, ) else - nothing + NoBias() end hmc = HMC( @@ -52,23 +57,23 @@ function test_reversibility(; nothing end - dH_n = reversibility_test(hmc, U, fermion, nothing, "without bias", to) + dH_n = reversibility_test(hmc, U, fermion, NoBias(), "without bias", to) dH_b = with_bias ? reversibility_test(hmcB, U, fermion, bias, "with bias", to) : nothing # show(to) return dH_n, dH_b end function reversibility_test(hmc::HMC{TI}, U, fermion, bias, str, to) where {TI} - @timeit to "sample_pseudofermions!" if isnothing(fermion) + @timeit to "sample_pseudofermions!" if isnothing(hmc.ϕ) nothing else - sample_pseudofermions!(hmc.ϕ[1], fermion[1], U) + sample_pseudofermions!(hmc.ϕ, fermion, U) end @timeit to "gaussian_TA!" gaussian_TA!(hmc.P, hmc.friction) @timeit to "gauge action" Sg_old = calc_gauge_action(U) @timeit to "kinetic energy" trP²_old = -calc_kinetic_energy(hmc.P) @timeit to "fermion action" Sf_old = - isnothing(fermion) ? 0.0 : calc_fermion_action(fermion[1], U, hmc.ϕ[1]) + isnothing(hmc.ϕ) ? 0.0 : calc_fermion_action(fermion, U, hmc.ϕ) @timeit to "CV" CV_old = isnothing(bias) ? 0.0 : calc_CV(U, bias) V_old = isnothing(bias) ? 0.0 : bias(CV_old) H_old = Sg_old + trP²_old + Sf_old + V_old @@ -76,7 +81,10 @@ function reversibility_test(hmc::HMC{TI}, U, fermion, bias, str, to) where {TI} @timeit to "evolve!" evolve!(hmc.integrator, U, hmc, fermion, bias) @timeit to "normalize!" MetaQCD.normalize!(U) - @show (calc_gauge_action(U) - calc_kinetic_energy(hmc.P)) - H_old + CV_mid = isnothing(bias) ? 0.0 : calc_CV(U, bias) + V_mid = isnothing(bias) ? 0.0 : bias(CV_mid) + H_mid = calc_gauge_action(U) - calc_kinetic_energy(hmc.P) + V_mid + ΔH_mid = H_mid - H_old @timeit to "invert momenta" MetaQCD.Fields.mul!(hmc.P, -1) @timeit to "evolve!" evolve!(hmc.integrator, U, hmc, fermion, bias) @@ -85,18 +93,19 @@ function reversibility_test(hmc::HMC{TI}, U, fermion, bias, str, to) where {TI} @timeit to "gauge action" Sg_new = calc_gauge_action(U) @timeit to "kinetic energy" trP²_new = -calc_kinetic_energy(hmc.P) @timeit to "fermion action" Sf_new = - isnothing(fermion) ? 0.0 : calc_fermion_action(fermion[1], U, hmc.ϕ[1]) + isnothing(hmc.ϕ) ? 0.0 : calc_fermion_action(fermion, U, hmc.ϕ) @timeit to "CV" CV_new = isnothing(bias) ? 0.0 : calc_CV(U, bias) V_new = isnothing(bias) ? 0.0 : bias(CV_new) H_new = Sg_new + trP²_new + Sf_new + V_new ΔH = H_new - H_old println("┌ $str") - @printf("| ΔSg = %+.5e (%+.10e)\n", Sg_new - Sg_old, Sg_old) - @printf("| ΔtrP² = %+.5e (%+.10e)\n", trP²_new - trP²_old, trP²_old) - @printf("| ΔSf = %+.5e (%+.10e)\n", Sf_new - Sf_old, Sf_old) - @printf("| ΔV = %+.5e (%+.10e)\n", V_new - V_old, V_old) - @printf("| ΔH = %+.5e (%+.10e)\n", ΔH, H_old) + @printf("| ΔSg = %+.5e (%+.10e)\n", Sg_new - Sg_old, Sg_old) + @printf("| ΔtrP² = %+.5e (%+.10e)\n", trP²_new - trP²_old, trP²_old) + @printf("| ΔSf = %+.5e (%+.10e)\n", Sf_new - Sf_old, Sf_old) + @printf("| ΔV = %+.5e (%+.10e)\n", V_new - V_old, V_old) + @printf("| ΔH_mid = %+.5e (%+.10e)\n", ΔH_mid, H_old) + @printf("| ΔH = %+.5e (%+.10e)\n", ΔH, H_old) @printf("└\n") return ΔH end