Skip to content

Commit

Permalink
Klara removed
Browse files Browse the repository at this point in the history
  • Loading branch information
@montyvesselinov
montyvesselinov committed Aug 4, 2024
1 parent 97ad97e commit babed2e
Show file tree
Hide file tree
Showing 2 changed files with 23 additions and 108 deletions.
34 changes: 0 additions & 34 deletions examples/bayesian_sampling/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,42 +25,8 @@ if Mads.create_tests
JLD2.save(joinpath(d, "mcmcchains_emcee.jld2"), "mcmcchains_emcee", mcmcchains_emcee)
end

if !haskey(ENV, "MADS_NO_KLARA") && isdefined(Mads, :Klara) && isdefined(Klara, :BasicContMuvParameter)
Mads.bayessampling(md, 1; nsteps=1, burnin=1, thinning=1)
mcmcchains_bayes = Mads.bayessampling(md, 2; nsteps=10, burnin=1, thinning=1, seed=2016)
mcmcchain = Mads.bayessampling(md; nsteps=10, burnin=1, thinning=1, seed=2016)
Mads.savemcmcresults(permutedims(mcmcchain.value), rootname * "-test-mcmcchain1.json")
rm(rootname * "-test-mcmcchain1.json")
Mads.forward(md, mcmcchain.value; all=true)
end

md = Mads.loadmadsfile(joinpath(workdir, "w01.mads"))
rootname = Mads.getmadsrootname(md)
if !haskey(ENV, "MADS_NO_KLARA") && isdefined(Mads, :Klara) && isdefined(Klara, :BasicContMuvParameter)
mcmcchain = Mads.bayessampling(md; nsteps=10, burnin=1, thinning=1, seed=2016)
mcmcvalues = Mads.paramarray2dict(md, permutedims(mcmcchain.value)) # convert the parameters in the chain to a parameter dictionary of arrays
mcmcvalues_array = hcat(vcat(map(i->collect(mcmcvalues[i]), keys(mcmcvalues)))...)
Mads.forward(md, mcmcchain.value)
if !haskey(ENV, "MADS_NO_GADFLY")
Mads.scatterplotsamples(md, permutedims(mcmcchain.value), rootname * "-test-bayes-results.svg")
Mads.rmfile(rootname * "-test-bayes-results.svg")
Mads.spaghettiplots(md, mcmcvalues; keyword="", obs_plot_dots=false)
Mads.spaghettiplot(md, mcmcvalues; keyword="test")
Mads.spaghettiplots(md, 3; keyword="test", grayscale=true)
Mads.spaghettiplot(md, 3; keyword="test", grayscale=true)
s = splitdir(rootname)
for filesinadir in Mads.searchdir(Regex(string(s[2], "-", "[.]*", "spaghetti.svg")), path=s[1])
Mads.rmfile(filesinadir, path=s[1])
end
end
if Mads.create_tests
d = joinpath(workdir, "test_results")
Mads.mkdir(d)
JLD2.save(joinpath(d, "mcmcvalues.jld2"), "mcmcvalues", mcmcvalues)
end
good_mcmcvalues = JLD2.load(joinpath(workdir, "test_results", "mcmcvalues.jld2"), "mcmcvalues")
good_mcmcvalues_array = hcat(vcat(map(i->collect(good_mcmcvalues[i]), keys(good_mcmcvalues)))...)
end

good_mcmcchains_emcee = JLD2.load(joinpath(workdir, "test_results", "mcmcchains_emcee.jld2"), "mcmcchains_emcee")

Expand Down
97 changes: 23 additions & 74 deletions src/MadsMonteCarlo.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,16 @@ function loadecmeeresults(madsdata::AbstractDict, filename::AbstractString)
end

function emceesampling(madsdata::AbstractDict; filename::AbstractString="", load::Bool=true, save::Bool=true, execute::Bool=true, numwalkers::Integer=10, sigma::Number=0.01, seed::Integer=-1, rng::Union{Nothing,Random.AbstractRNG,DataType}=nothing, kw...)
if load && filename != ""
if filename != ""
load = save = true
end
if (load || save) && filename == ""
np = length(Mads.getoptparams(madsdata))
no = length(Mads.getobskeys(madsdata))
filename = joinpath(Mads.getmadsproblemdir(madsdata), Mads.getmadsrootname(madsdata) * "_emcee_results_$(np)_$(no)_$(numwalkers)_$(nexecutions)_$(burnin)_$(thinning).jld2")
madsinfo("Filename not provided! AffineInvariantMCMC results will be saved in $(filename) ...")
end
if load
chain, llhoods, observations = loadecmeeresults(madsdata, filename)
if isnothing(chain)
if execute
Expand Down Expand Up @@ -113,13 +122,21 @@ function emceesampling(madsdata::AbstractDict; filename::AbstractString="", load
p0[i, j] = pmin[i] + rand(Mads.rng, d) * (pmax[i] - pmin[i])
end
end
chain, llhoods, observations = emceesampling(madsdata, p0; filename=filename, numwalkers=numwalkers, seed=seed, rng=rng, kw...)
chain, llhoods, observations = emceesampling(madsdata, p0; filename=filename, load=load, save=save, execute=execute, numwalkers=numwalkers, seed=seed, rng=rng, kw...)
return chain, llhoods, observations
end

function emceesampling(madsdata::AbstractDict, p0::AbstractMatrix; filename::AbstractString="", load::Bool=true, save::Bool=true, execute::Bool=true, numwalkers::Integer=10, nexecutions::Integer=100, burnin::Integer=numwalkers, thinning::Integer=10, seed::Integer=-1, weightfactor::Number=1.0, rng::Union{Nothing,Random.AbstractRNG,DataType}=nothing, distributed_function::Bool=false)
function emceesampling(madsdata::AbstractDict, p0::AbstractMatrix; filename::AbstractString="", load::Bool=false, save::Bool=false, execute::Bool=true, numwalkers::Integer=10, nexecutions::Integer=100, burnin::Integer=numwalkers, thinning::Integer=10, seed::Integer=-1, weightfactor::Number=1.0, rng::Union{Nothing,Random.AbstractRNG,DataType}=nothing, distributed_function::Bool=false)
if filename != ""
load = save = true
end
np = length(Mads.getoptparams(madsdata))
no = length(Mads.getobskeys(madsdata))
if (load || save) && filename == ""
filename = joinpath(Mads.getmadsproblemdir(madsdata), Mads.getmadsrootname(madsdata) * "_emcee_results_$(np)_$(no)_$(numwalkers)_$(nexecutions)_$(burnin)_$(thinning).jld2")
madsinfo("Filename not provided! AffineInvariantMCMC results will be saved in $(filename) ...")
end
numsamples_perwalker = div(nexecutions, numwalkers)
if load && filename != ""
if load
bad_data = false
@info("AffineInvariantMCMC preexisting data loading $(filename) ...")
chain, llhoods, observations = loadecmeeresults(madsdata, filename)
Expand All @@ -130,7 +147,7 @@ function emceesampling(madsdata::AbstractDict, p0::AbstractMatrix; filename::Abs
if size(chain, 2) != div(numsamples_perwalker, thinning) * numwalkers
@warn("Preexisting data does not match the number of walkers and steps!")
bad_data = true
elseif size(chain, 1) != length(Mads.getoptparamkeys(madsdata))
elseif size(chain, 1) != no
@warn("Preexisting data does not match the number of parameters!")
bad_data = true
end
Expand All @@ -147,10 +164,6 @@ function emceesampling(madsdata::AbstractDict, p0::AbstractMatrix; filename::Abs
return chain, llhoods, observations
end
end
if save && filename == ""
filename = joinpath(Mads.getmadsproblemdir(madsdata), Mads.getmadsrootname(madsdata) * "_emcee_results_$(numwalkers)_$(nexecutions)_$(burnin)_$(thinning).jld2")
madsinfo("Filename not provided! AffineInvariantMCMC results will be saved in $(filename) ...")
end
Mads.setseed(seed; rng=rng)
madsloglikelihood = makemadsloglikelihood(madsdata; weightfactor=weightfactor)
arrayloglikelihood = Mads.makearrayloglikelihood(madsdata, madsloglikelihood)
Expand Down Expand Up @@ -178,7 +191,6 @@ function emceesampling(madsdata::AbstractDict, p0::AbstractMatrix; filename::Abs
end
return chain, llhoods, observations
end

@doc """
Bayesian sampling with Goodman & Weare's Affine Invariant Markov chain Monte Carlo (MCMC) Ensemble sampler (aka Emcee)
Expand Down Expand Up @@ -206,62 +218,6 @@ Mads.emceesampling(madsdata, p0; numwalkers=10, nsteps=100, burnin=10, thinning=
```
""" emceesampling

if isdefined(Mads, :Klara)
# && isdefined(Klara, :BasicContMuvParameter)
function bayessampling(madsdata::AbstractDict; nsteps::Integer=1000, burnin::Integer=100, thinning::Integer=1, seed::Integer=-1, rng::Union{Nothing,Random.AbstractRNG,DataType}=nothing)
Mads.setseed(seed; rng=rng)
madsloglikelihood = makemadsloglikelihood(madsdata)
arrayloglikelihood = makearrayloglikelihood(madsdata, madsloglikelihood)
optparamkeys = getoptparamkeys(madsdata)
initvals = Array{Float64}(undef, length(optparamkeys))
for i = eachindex(optparamkeys)
initvals[i] = madsdata["Parameters"][optparamkeys[i]]["init"]
end
# mcparams = Klara.BasicContMuvParameter(:p, logtarget=arrayloglikelihood)
# model = Klara.likelihood_model(mcparams, false)
# # sampler = Klara.MH(fill(1e-1, length(initvals)))
# sampler = Klara.RAM(fill(1e-1, length(initvals)))
# mcrange = Klara.BasicMCRange(nsteps=nsteps + burnin, burnin=burnin, thinning=thinning)
# mcparams0 = Dict(:p=>initvals)
# outopts = Dict{Symbol, Any}(:monitor=>[:value, :logtarget, :loglikelihood], :diagnostics=>[:accept])
# job = Klara.BasicMCJob(model, sampler, mcrange, mcparams0, outopts=outopts, tuner=Klara.VanillaMCTuner())
# Klara.run(job)
# chain = Klara.output(job)
# return chain
end
function bayessampling(madsdata::AbstractDict, numsequences::Integer; nsteps::Integer=1000, burnin::Integer=100, thinning::Integer=1, seed::Integer=-1)
if seed != 0
mcmcchains = RobustPmap.rpmap(i->bayessampling(madsdata; nsteps=nsteps, burnin=burnin, thinning=thinning, seed=seed+i), 1:numsequences)
else
mcmcchains = RobustPmap.rpmap(i->bayessampling(madsdata; nsteps=nsteps, burnin=burnin, thinning=thinning), 1:numsequences)
end
return mcmcchains
end

@doc """
Bayesian Sampling
$(DocumentFunction.documentfunction(bayessampling;
argtext=Dict("madsdata"=>"MADS problem dictionary",
"numsequences"=>"number of sequences executed in parallel"),
keytext=Dict("nsteps"=>"number of final realizations in the chain [default=`1000`]",
"burnin"=>"number of initial realizations before the MCMC are recorded [default=`100`]",
"thinning"=>"removal of any `thinning` realization [default=`1`]",
"seed"=>"random seed [default=`0`]")))
Returns:
- MCMC chain
Examples:
```julia
Mads.bayessampling(madsdata; nsteps=1000, burnin=100, thinning=1, seed=2016)
Mads.bayessampling(madsdata, numsequences; nsteps=1000, burnin=100, thinning=1, seed=2016)
```
""" bayessampling
end

"""
Save MCMC chain in a file
Expand All @@ -286,13 +242,6 @@ function savemcmcresults(chain::AbstractArray, filename::AbstractString)
println(f, "")
print(f, "Variance: ")
JSON.print(f, Statistics.var(chain; dims=1)[:])
#println(f, "")
#print(f, "MCse: ")
#JSON.print(f, Klara.mcse(chain))
#println(f, "")
#print(f, "MCvar: ")
#JSON.print(f, Klara.mcvar(chain))
#println(f, "")
close(f)
end

Expand Down

0 comments on commit babed2e

Please sign in to comment.