From a074620c08d83eda9d7cb0b35f5b5706e922284f Mon Sep 17 00:00:00 2001 From: monty Date: Thu, 2 Nov 2023 06:51:02 -0600 Subject: [PATCH] ProgressMeter --- src-old/MadsEmcee.jl | 17 +++++++++++------ src/MadsMonteCarlo.jl | 7 ++++++- src/MadsSensitivityAnalysis.jl | 2 +- 3 files changed, 18 insertions(+), 8 deletions(-) diff --git a/src-old/MadsEmcee.jl b/src-old/MadsEmcee.jl index 3503c895..26818b1d 100644 --- a/src-old/MadsEmcee.jl +++ b/src-old/MadsEmcee.jl @@ -23,8 +23,10 @@ Example: Mads.emcee(llhood, numwalkers=10, numsamples_perwalker=100, thinning=1) ``` """ -function emcee(llhood::Function, numwalkers::Integer, x0::AbstractArray, numsamples_perwalker::Integer, thinning::Integer, a::Number=2.) - @assert length(size(x0)) == 2 +function sample(llhood::Function, numwalkers::Integer, x0::AbstractMatrix, numsamples_perwalker::Integer, thinning::Integer, a::Number=2.; filename::AbstractString="", load::Bool=true, save::Bool=true, rng::Random.AbstractRNG=Random.GLOBAL_RNG) + if filename != "" && isfile(filename) && load + chain, llhoodvals = JLD2.load(filename, "chain", "llhoods") + end x = copy(x0) chain = Array{Float64}(undef, size(x0, 1), numwalkers, div(numsamples_perwalker, thinning)) lastllhoodvals = RobustPmap.rpmap(llhood, map(i->x[:, i], 1:size(x, 2))) @@ -34,18 +36,18 @@ function emcee(llhood::Function, numwalkers::Integer, x0::AbstractArray, numsamp batch1 = 1:div(numwalkers, 2) batch2 = div(numwalkers, 2) + 1:numwalkers divisions = [(batch1, batch2), (batch2, batch1)] - @ProgressMeter.showprogress 4 for i = 1:numsamples_perwalker + @ProgressMeter.showprogress 1 for i = 1:numsamples_perwalker for ensembles in divisions active, inactive = ensembles - zs = map(u->((a - 1) * u + 1)^2 / a, rand(Mads.rng, length(active))) - proposals = map(i->zs[i] * x[:, active[i]] + (1 - zs[i]) * x[:, rand(Mads.rng, inactive)], 1:length(active)) + zs = map(u->((a - 1) * u + 1)^2 / a, rand(rng, length(active))) + proposals = map(i->zs[i] * x[:, active[i]] + (1 - zs[i]) * x[:, rand(rng, inactive)], 1:length(active)) newllhoods = RobustPmap.rpmap(llhood, proposals) for (j, walkernum) in enumerate(active) z = zs[j] newllhood = newllhoods[j] proposal = proposals[j] logratio = (size(x, 1) - 1) * log(z) + newllhood - lastllhoodvals[walkernum] - if log(rand(Mads.rng)) < logratio + if log(rand(rng)) < logratio lastllhoodvals[walkernum] = newllhood x[:, walkernum] = proposal end @@ -56,6 +58,9 @@ function emcee(llhood::Function, numwalkers::Integer, x0::AbstractArray, numsamp end end end + if save && filename != "" + JLD2.save(filename, "chain", chain, "llhoods", llhoodvals) + end return chain, llhoodvals end diff --git a/src/MadsMonteCarlo.jl b/src/MadsMonteCarlo.jl index f07c2e3b..ef2981c3 100644 --- a/src/MadsMonteCarlo.jl +++ b/src/MadsMonteCarlo.jl @@ -4,6 +4,7 @@ import JSON import AffineInvariantMCMC import DocumentFunction import BlackBoxOptim +import ProgressMeter import Random function p10_p50_p90(x::AbstractMatrix) @@ -113,11 +114,13 @@ function emceesampling(madsdata::AbstractDict, p0::AbstractMatrix; filename::Abs chain, llhoods, observations = loadecmeeresults(madsdata, filename) if isnothing(chain) if execute - @info("AffineInvariantMCMC will be executed!") + @info("AffineInvariantMCMC will be executed! No preexisting data to load ...") else + @warn("No preexisting data to load! AffineInvariantMCMC will be not executed!") return nothing, nothing, nothing end else + @info("AffineInvariantMCMC preexisting data loaded!") return chain, llhoods, observations end end @@ -132,11 +135,13 @@ function emceesampling(madsdata::AbstractDict, p0::AbstractMatrix; filename::Abs if newnsteps < 2 newnsteps = 2 end + @info("AffineInvariantMCMC burning stage ...") burninchain, _ = AffineInvariantMCMC.sample(arrayloglikelihood, numwalkers, p0, newnsteps, 1; filename="", rng=Mads.rng, save=false) newnsteps = div(nsteps, numwalkers) if newnsteps < 2 newnsteps = 2 end + @info("AffineInvariantMCMC exploration stage ...") chain, llhoods = AffineInvariantMCMC.sample(arrayloglikelihood, numwalkers, burninchain[:, :, end], newnsteps, thinning; filename="", rng=Mads.rng, save=false) chain, llhoods = AffineInvariantMCMC.flattenmcmcarray(chain, llhoods) observations = Mads.forward(madsdata, permutedims(chain)) diff --git a/src/MadsSensitivityAnalysis.jl b/src/MadsSensitivityAnalysis.jl index 6b1950e9..313e5708 100644 --- a/src/MadsSensitivityAnalysis.jl +++ b/src/MadsSensitivityAnalysis.jl @@ -505,7 +505,7 @@ function saltellibrute(madsdata::AbstractDict; N::Integer=1000, seed::Integer=-1 for i = eachindex(paramkeys) madsinfo("Parameter : $(paramkeys[i])") cond_means = Array{OrderedCollections.OrderedDict}(undef, numoneparamsamples) - @ProgressMeter.showprogress 1 "Computing ... " for j = 1:numoneparamsamples + @ProgressMeter.showprogress 1 "Computing ... " for j = 1:numoneparamsamples cond_means[j] = OrderedCollections.OrderedDict() for k = eachindex(obskeys) cond_means[j][obskeys[k]] = 0.