Skip to content

Commit

Permalink
reformat
Browse files Browse the repository at this point in the history
  • Loading branch information
@montyvesselinov
montyvesselinov committed Nov 5, 2024
1 parent 168d422 commit 1a591a4
Show file tree
Hide file tree
Showing 4 changed files with 44 additions and 44 deletions.
22 changes: 11 additions & 11 deletions examples/anasol/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ if !haskey(ENV, "MADS_NO_GADFLY")
Mads.spaghettiplots(md, 2)
Mads.rmfile(joinpath(workdir, "w01short-ax-2-spaghetti.svg"))
Mads.rmfile(joinpath(workdir, "w01short-vx-2-spaghetti.svg"))
Mads.plotmass([1.,2.],[10.,20.],[5.,10.], joinpath(workdir, "plotmass-test.svg"))
Mads.plotmass([1.0, 2.0], [10.0, 20.0], [5.0, 10.0], joinpath(workdir, "plotmass-test.svg"))
Mads.rmfile(joinpath(workdir, "plotmass-test.svg"))
Mads.graphon()
end
Expand Down Expand Up @@ -99,7 +99,7 @@ Mads.wellon!(md, "w1a")
Mads.savemadsfile(md, joinpath(workdir, "test.mads"))
Mads.savemadsfile(md, Mads.getparamdict(md))
Mads.savemadsfile(md, Mads.getparamdict(md), joinpath(workdir, "test.mads"))
Mads.savemadsfile(md, Mads.getparamdict(md), joinpath(workdir, "test.mads"), explicit=true)
Mads.savemadsfile(md, Mads.getparamdict(md), joinpath(workdir, "test.mads"); explicit=true)
Mads.rmfile(joinpath(workdir, "w01shortexp-rerun.mads"))
Mads.rmfile(joinpath(workdir, "test.mads"))
Mads.setmadsinputfile("test.mads")
Expand Down Expand Up @@ -128,7 +128,7 @@ mdbad = deepcopy(md)
Mads.setparamsdistnormal!(mdbad, fill(1, length(m1)), fill(2, length(m1)))
Mads.setparamsdistuniform!(mdbad, fill(1, length(m1)), fill(2, length(m1)))

Mads.computemass("w01lambda", time=50, path=workdir)
Mads.computemass("w01lambda"; time=50, path=workdir)
Mads.rmdir(joinpath(workdir, "mass_reduced.svg"))

d = joinpath(workdir, "test_results")
Expand Down Expand Up @@ -163,35 +163,35 @@ good_targetkeys = JLD2.load(joinpath(d, "targetkeys.jld2"), "tk")

Test.@testset "Anasol" begin
# Test Mads.forward(md; all=true)
ssr = 0.
ssr = 0.0
for obskey in union(Set(keys(forward_results)), Set(keys(good_forward_results)))
ssr += (forward_results[obskey] - good_forward_results[obskey])^2
end
Test.@test isapprox(ssr, 0., atol=1e-8)
Test.@test isapprox(ssr, 0.0, atol=1e-8)

# Test Mads.getparamrandom(md, 5, init_dist=true)
# Test Mads.getparamrandom(md, 5; init_dist=true)
for obskey in union(Set(keys(paramrandom_true)), Set(keys(good_paramrandom_true)))
Test.@test isapprox(paramrandom_true[obskey], good_paramrandom_true[obskey], atol=1e-8)
end

# Test Mads.getparamrandom(md, 5, init_dist=false)
# Test Mads.getparamrandom(md, 5; init_dist=false)
for obskey in union(Set(keys(paramrandom_false)), Set(keys(good_paramrandom_false)))
Test.@test isapprox(paramrandom_false[obskey], good_paramrandom_false[obskey], atol=1e-8)
end

# Test Mads.getparamdict(md)
ssr = 0.
ssr = 0.0
for obskey in union(Set(keys(pd)), Set(keys(good_pd)))
ssr += (pd[obskey] - good_pd[obskey])^2
end
Test.@test isapprox(ssr, 0., atol=1e-8)
Test.@test isapprox(ssr, 0.0, atol=1e-8)

# Test computeconcentrations(paramdict)
ssr = 0.
ssr = 0.0
for obskey in union(Set(keys(forward_preds)), Set(keys(good_forward_preds)))
ssr += (forward_preds[obskey] - good_forward_preds[obskey])^2
end
Test.@test isapprox(ssr, 0., atol=1e-8)
Test.@test isapprox(ssr, 0.0, atol=1e-8)

# Test Mads.of(md)
Test.@test isapprox(madsOf, 628963.6972820368, atol=1e-8) #test
Expand Down
28 changes: 14 additions & 14 deletions examples/model_analysis/bayes_weight_analsis.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,29 +7,29 @@ Mads.mkdir("bayes_results")

@info("Bayesian analysis with different initial parameter sets and observation weights (standard deviation errors)")
@info("Bayesian analysis for the initial parameter guesses:")
for w = (1000000, 1000, 1)
for w in (1000000, 1000, 1)
Mads.setobsweights!(md, w)
mcmcchain = Mads.bayessampling(md; nsteps=1000000, burnin=100000, thinning=100, seed=2016)
Mads.scatterplotsamples(md, permutedims(mcmcchain.value), "bayes_results/bayes_init_w$w.png")
o = Mads.forward(md, mcmcchain.value)
Mads.spaghettiplot(md, o, filename="bayes_results/bayes_init_w$(w)_spaghetti.png")
Printf.@printf "Init: Observation Weight %d StdDev %f ->`o5` prediction: min = %f max = %f\n" w 1/w min(o[:,5]...) max(o[:,5]...)
f = Gadfly.plot(x=o[:,5], Gadfly.Guide.xlabel("o5"), Gadfly.Geom.histogram())
Mads.spaghettiplot(md, o; filename="bayes_results/bayes_init_w$(w)_spaghetti.png")
Printf.@printf "Init: Observation Weight %d StdDev %f ->`o5` prediction: min = %f max = %f\n" w 1 / w min(o[:, 5]...) max(o[:, 5]...)
f = Gadfly.plot(; x=o[:, 5], Gadfly.Guide.xlabel("o5"), Gadfly.Geom.histogram())
Gadfly.draw(Gadfly.PNG("bayes_results/bayes_init_w$(w)_o5.png", 6Gadfly.inch, 4Gadfly.inch), f)
end
pinit = Dict(zip(Mads.getparamkeys(md), Mads.getparamsinit(md)))

n = 100
@info("Calibration using $n random initial guesses for model parameters")
r = Mads.calibraterandom(md, n, all=true, seed=2016, save_results=false)
pnames = collect(keys(r[1,3]))
p = hcat(map(i->collect(values(r[i,3])), 1:n)...)'
r = Mads.calibraterandom(md, n; all=true, seed=2016, save_results=false)
pnames = collect(keys(r[1, 3]))
p = hcat(map(i -> collect(values(r[i, 3])), 1:n)...)'
np = length(pnames)
@info("Identify the 3 different global optima with different model parameter estimates")

ind_n0 = abs(p[:,4]) .< 0.1
ind_n0 = abs(p[:, 4]) .< 0.1
in0 = findall(ind_n0 .== true)[1]
ind_n1 = abs(p[:,4]-1) .< 0.1
ind_n1 = abs(p[:, 4] - 1) .< 0.1
in1 = findall(ind_n1 .== true)[1]
ind_n01 = !(ind_n0 | ind_n1)
in01 = findall(ind_n01 .== true)[1]
Expand All @@ -39,15 +39,15 @@ v = [in0, in1, in01]

@info("Bayesian analysis for the 3 different global optima")
for i = 1:3
Mads.setparamsinit!(md, r[v[i],3])
for w = (1000000, 1000, 1)
Mads.setparamsinit!(md, r[v[i], 3])
for w in (1000000, 1000, 1)
Mads.setobsweights!(md, w)
mcmcchain = Mads.bayessampling(md; nsteps=1000000, burnin=100000, thinning=100, seed=2016)
Mads.scatterplotsamples(md, permutedims(mcmcchain.value), "bayes_results/bayes_opt_$(optnames[i])_w$w.png")
o = Mads.forward(md, mcmcchain.value)
Mads.spaghettiplot(md, o, filename="bayes_results/bayes_opt_$(optnames[i])_w$(w)_spaghetti.png")
Printf.@printf "O%-3s: Observation Weight %d StdDev %f -> `o5` prediction: min = %f max = %f\n" optnames[i] w 1/w min(o[:,5]...) max(o[:,5]...)
f = Gadfly.plot(x=o[:,5], Gadfly.Guide.xlabel("o5"), Gadfly.Geom.histogram())
Mads.spaghettiplot(md, o; filename="bayes_results/bayes_opt_$(optnames[i])_w$(w)_spaghetti.png")
Printf.@printf "O%-3s: Observation Weight %d StdDev %f -> `o5` prediction: min = %f max = %f\n" optnames[i] w 1 / w min(o[:, 5]...) max(o[:, 5]...)
f = Gadfly.plot(; x=o[:, 5], Gadfly.Guide.xlabel("o5"), Gadfly.Geom.histogram())
Gadfly.draw(Gadfly.PNG("bayes_results/bayes_opt_$(optnames[i])_w$(w)_o5.png", 6Gadfly.inch, 4Gadfly.inch), f)
end
end
Expand Down
10 changes: 5 additions & 5 deletions examples/model_analysis/emcee_walkers_analsis.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,13 +6,13 @@ md = Mads.loadmadsfile(joinpath("models", "internal-polynomial.mads"))
Mads.mkdir("emcee_results")

@info("AffineInvariantMCMC (EMCEE) Bayesian analysis with different number of walkers:")
for nw = (100, 50, 20, 10, 5, 1)
for nw in (100, 50, 20, 10, 5, 1)
Mads.setobsweights!(md, 1000000)
@time chain, llhoods = Mads.emceesampling(md; numwalkers=nw, nsteps=1000000, burnin=100000, thinning=100, seed=2016)
Mads.scatterplotsamples(md, chain', "emcee_results/emcee_init_nw$nw.png")
o = Mads.forward(md, chain')
Mads.spaghettiplot(md, o, filename="emcee_results/emcee_init_nw$(nw)_spaghetti.png")
Printf.@printf "Init: Number of walkers %d ->`o5` prediction: min = %f max = %f\n" nw min(o[:,5]...) max(o[:,5]...)
f = Gadfly.plot(x=o[:,5], Gadfly.Guide.xlabel("o5"), Gadfly.Geom.histogram())
Mads.spaghettiplot(md, o; filename="emcee_results/emcee_init_nw$(nw)_spaghetti.png")
Printf.@printf "Init: Number of walkers %d ->`o5` prediction: min = %f max = %f\n" nw min(o[:, 5]...) max(o[:, 5]...)
f = Gadfly.plot(; x=o[:, 5], Gadfly.Guide.xlabel("o5"), Gadfly.Geom.histogram())
Gadfly.draw(Gadfly.PNG("emcee_results/emcee_init_nw$(nw)_o5.png", 6Gadfly.inch, 4Gadfly.inch), f)
end
end
28 changes: 14 additions & 14 deletions examples/model_analysis/emcee_weight_analsis.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,29 +7,29 @@ Mads.mkdir("emcee_results")

@info("AffineInvariantMCMC (EMCEE) Bayesian analysis with different initial parameter sets and observation weights (standard deviation errors)")
@info("AffineInvariantMCMC (EMCEE) Bayesian analysis for the initial parameter guesses:")
for w = (1000000, 1000, 1)
for w in (1000000, 1000, 1)
Mads.setobsweights!(md, w)
chain, llhoods = Mads.emceesampling(md; numwalkers=100, nsteps=1000000, burnin=100000, thinning=100, seed=2016)
Mads.scatterplotsamples(md, chain', "emcee_results/emcee_init_w$w.png")
o = Mads.forward(md, chain')
Mads.spaghettiplot(md, o, filename="emcee_results/emcee_init_w$(w)_spaghetti.png")
Printf.@printf "Init: Observation Weight %d StdDev %f ->`o5` prediction: min = %f max = %f\n" w 1/w min(o[:,5]...) max(o[:,5]...)
f = Gadfly.plot(x=o[:,5], Gadfly.Guide.xlabel("o5"), Gadfly.Geom.histogram())
Mads.spaghettiplot(md, o; filename="emcee_results/emcee_init_w$(w)_spaghetti.png")
Printf.@printf "Init: Observation Weight %d StdDev %f ->`o5` prediction: min = %f max = %f\n" w 1 / w min(o[:, 5]...) max(o[:, 5]...)
f = Gadfly.plot(; x=o[:, 5], Gadfly.Guide.xlabel("o5"), Gadfly.Geom.histogram())
Gadfly.draw(Gadfly.PNG("emcee_results/emcee_init_w$(w)_o5.png", 6Gadfly.inch, 4Gadfly.inch), f)
end
pinit = Dict(zip(Mads.getparamkeys(md), Mads.getparamsinit(md)))

n = 100
@info("Calibration using $n random initial guesses for model parameters")
r = Mads.calibraterandom(md, n, all=true, seed=2016, save_results=false)
pnames = collect(keys(r[1,3]))
p = hcat(map(i->collect(values(r[i,3])), 1:n)...)'
r = Mads.calibraterandom(md, n; all=true, seed=2016, save_results=false)
pnames = collect(keys(r[1, 3]))
p = hcat(map(i -> collect(values(r[i, 3])), 1:n)...)'
np = length(pnames)
@info("Identify the 3 different global optima with different model parameter estimates")

ind_n0 = abs(p[:,4]) .< 0.1
ind_n0 = abs(p[:, 4]) .< 0.1
in0 = findall(ind_n0 .== true)[1]
ind_n1 = abs(p[:,4]-1) .< 0.1
ind_n1 = abs(p[:, 4] - 1) .< 0.1
in1 = findall(ind_n1 .== true)[1]
ind_n01 = !(ind_n0 | ind_n1)
in01 = findall(ind_n01 .== true)[1]
Expand All @@ -39,15 +39,15 @@ v = [in0, in1, in01]

@info("AffineInvariantMCMC (EMCEE) Bayesian analysis for the 3 different global optima")
for i = 1:3
Mads.setparamsinit!(md, r[v[i],3])
for w = (1000000, 1000, 1)
Mads.setparamsinit!(md, r[v[i], 3])
for w in (1000000, 1000, 1)
Mads.setobsweights!(md, w)
chain, llhoods = Mads.emceesampling(md; numwalkers=100, nsteps=1000000, burnin=100000, thinning=100, seed=2016)
Mads.scatterplotsamples(md, chain', "emcee_results/emcee_opt_$(optnames[i])_w$w.png")
o = Mads.forward(md, chain')
Mads.spaghettiplot(md, o, filename="emcee_results/emcee_opt_$(optnames[i])_w$(w)_spaghetti.png")
Printf.@printf "O%-3s: Observation Weight %d StdDev %f -> `o5` prediction: min = %f max = %f\n" optnames[i] w 1/w min(o[:,5]...) max(o[:,5]...)
f = Gadfly.plot(x=o[:,5], Gadfly.Guide.xlabel("o5"), Gadfly.Geom.histogram())
Mads.spaghettiplot(md, o; filename="emcee_results/emcee_opt_$(optnames[i])_w$(w)_spaghetti.png")
Printf.@printf "O%-3s: Observation Weight %d StdDev %f -> `o5` prediction: min = %f max = %f\n" optnames[i] w 1 / w min(o[:, 5]...) max(o[:, 5]...)
f = Gadfly.plot(; x=o[:, 5], Gadfly.Guide.xlabel("o5"), Gadfly.Geom.histogram())
Gadfly.draw(Gadfly.PNG("emcee_results/emcee_opt_$(optnames[i])_w$(w)_o5.png", 6Gadfly.inch, 4Gadfly.inch), f)
end
end
Expand Down

0 comments on commit 1a591a4

Please sign in to comment.