Skip to content

Commit

Permalink
format
Browse files Browse the repository at this point in the history
  • Loading branch information
@odunbar
odunbar committed Jul 11, 2024
1 parent 184141e commit 29bfe59
Show file tree
Hide file tree
Showing 11 changed files with 213 additions and 199 deletions.
38 changes: 18 additions & 20 deletions examples/EDMF_data/plot_posterior.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ using CalibrateEmulateSample.ParameterDistributions

# 5-parameter calibration exp
exp_name = "ent-det-tked-tkee-stab-calibration"
date_of_run = Date(2024,06,14)
date_of_run = Date(2024, 06, 14)

# Output figure read/write directory
figure_save_directory = joinpath(@__DIR__, "output", exp_name, string(date_of_run))
Expand Down Expand Up @@ -56,15 +56,12 @@ labels = get_name(posterior)

burnin = 50_000

data_rf= (; [(Symbol(labels[i]), posterior_samples[i, burnin:end]) for i in 1:length(labels)]...)
transformed_data_rf= (; [(Symbol(labels[i]), transformed_posterior_samples[i, burnin:end]) for i in 1:length(labels)]...)
data_rf = (; [(Symbol(labels[i]), posterior_samples[i, burnin:end]) for i in 1:length(labels)]...)
transformed_data_rf =
(; [(Symbol(labels[i]), transformed_posterior_samples[i, burnin:end]) for i in 1:length(labels)]...)

p = pairplot(
data_rf=> (PairPlots.Contourf(sigmas=1:1:3),),
)
trans_p = pairplot(
transformed_data_rf=> (PairPlots.Contourf(sigmas=1:1:3),),
)
p = pairplot(data_rf => (PairPlots.Contourf(sigmas = 1:1:3),))
trans_p = pairplot(transformed_data_rf => (PairPlots.Contourf(sigmas = 1:1:3),))

save(density_filepath, p)
save(transformed_density_filepath, trans_p)
Expand Down Expand Up @@ -93,19 +90,20 @@ density_filepath = joinpath(figure_save_directory, "$(case_rf)_$(case_gp)_poster
transformed_density_filepath = joinpath(figure_save_directory, "$(case_rf)_$(case_gp)_posterior_dist_phys.png")
labels = get_name(posterior)
data_gp = (; [(Symbol(labels[i]), posterior_samples[i, burnin:end]) for i in 1:length(labels)]...)
transformed_data_gp = (; [(Symbol(labels[i]), transformed_posterior_samples[i, burnin:end]) for i in 1:length(labels)]...)
transformed_data_gp =
(; [(Symbol(labels[i]), transformed_posterior_samples[i, burnin:end]) for i in 1:length(labels)]...)
#
#
#
gp_smoothing = 1 # >1 = smoothing KDE in plotting

p = pairplot(
data_rf=> (PairPlots.Contourf(sigmas=1:1:3),),
data_gp => (PairPlots.Contourf(sigmas=1:1:3,bandwidth=gp_smoothing),),
data_rf => (PairPlots.Contourf(sigmas = 1:1:3),),
data_gp => (PairPlots.Contourf(sigmas = 1:1:3, bandwidth = gp_smoothing),),
)
trans_p = pairplot(
transformed_data_rf=> (PairPlots.Contourf(sigmas=1:1:3),),
transformed_data_gp => (PairPlots.Contourf(sigmas=1:1:3,bandwidth=gp_smoothing),),
transformed_data_rf => (PairPlots.Contourf(sigmas = 1:1:3),),
transformed_data_gp => (PairPlots.Contourf(sigmas = 1:1:3, bandwidth = gp_smoothing),),
)

save(density_filepath, p)
Expand Down Expand Up @@ -135,22 +133,22 @@ transformed_density_filepath = joinpath(figure_save_directory, "all_posterior_di
labels = get_name(posterior)

data_prior = (; [(Symbol(labels[i]), posterior_samples[i, burnin:end]) for i in 1:length(labels)]...)
transformed_data_prior = (; [(Symbol(labels[i]), transformed_posterior_samples[i, burnin:end]) for i in 1:length(labels)]...)
transformed_data_prior =
(; [(Symbol(labels[i]), transformed_posterior_samples[i, burnin:end]) for i in 1:length(labels)]...)
#
#
#

p = pairplot(
data_rf=> (PairPlots.Contourf(sigmas=1:1:3),),
data_gp => (PairPlots.Contourf(sigmas=1:1:3,bandwidth=gp_smoothing),),
data_rf => (PairPlots.Contourf(sigmas = 1:1:3),),
data_gp => (PairPlots.Contourf(sigmas = 1:1:3, bandwidth = gp_smoothing),),
data_prior => (PairPlots.Scatter(),),
)
trans_p = pairplot(
transformed_data_rf=> (PairPlots.Contourf(sigmas=1:1:3),),
transformed_data_gp => (PairPlots.Contourf(sigmas=1:1:3,bandwidth=gp_smoothing),),
transformed_data_rf => (PairPlots.Contourf(sigmas = 1:1:3),),
transformed_data_gp => (PairPlots.Contourf(sigmas = 1:1:3, bandwidth = gp_smoothing),),
transformed_data_prior => (PairPlots.Scatter(),),
)

save(density_filepath, p)
save(transformed_density_filepath, trans_p)

13 changes: 5 additions & 8 deletions examples/EDMF_data/uq_for_edmf.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,14 +213,11 @@ function main()
"n_features_opt" => 200,
"localization" => SEC(0.05),
)
if case == "RF-prior"
overrides = Dict("verbose" => true,
"cov_sample_multiplier" => 0.01,
"n_iteration" => 0,
)
end
nugget = 1e-6
rng_seed = 99330
if case == "RF-prior"
overrides = Dict("verbose" => true, "cov_sample_multiplier" => 0.01, "n_iteration" => 0)
end
nugget = 1e-6
rng_seed = 99330
rng = Random.MersenneTwister(rng_seed)
input_dim = size(get_inputs(input_output_pairs), 1)
output_dim = size(get_outputs(input_output_pairs), 1)
Expand Down
138 changes: 81 additions & 57 deletions examples/Emulator/G-function/emulate.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ function main()
n_repeats = 3 # repeat exp with same data.
n_dimensions = 3
# To create the sampling
n_data_gen = 800
n_data_gen = 800

data =
SobolData(params = OrderedDict([Pair(Symbol("x", i), Uniform(0, 1)) for i in 1:n_dimensions]), N = n_data_gen)
Expand All @@ -68,7 +68,7 @@ function main()
CairoMakie.save(joinpath(output_directory, "GFunction_slices_truth_$(n_dimensions).png"), f1, px_per_unit = 3)
CairoMakie.save(joinpath(output_directory, "GFunction_slices_truth_$(n_dimensions).pdf"), f1, px_per_unit = 3)

n_train_pts = n_dimensions * 250
n_train_pts = n_dimensions * 250
ind = shuffle!(rng, Vector(1:n_data))[1:n_train_pts]
# now subsample the samples data
n_tp = length(ind)
Expand All @@ -91,7 +91,7 @@ function main()
TV = [(1 / (3 * (1 + ai)^2)) * prod_tmp2[i] / prod_tmp for (i, ai) in enumerate(a)]



cases = ["Prior", "GP", "RF-scalar"]
case = cases[3]
decorrelate = true
Expand All @@ -104,10 +104,10 @@ function main()
"train_fraction" => 0.8,#0.7
"n_iteration" => 20, # (=multiple of recompute_cov_at - 1 is most efficient)
"cov_sample_multiplier" => 1.0,
# "localization" => SEC(0.1),#,Doesnt help much tbh
# "accelerator" => NesterovAccelerator(),
# "localization" => SEC(0.1),#,Doesnt help much tbh
# "accelerator" => NesterovAccelerator(),
"n_ensemble" => 200, #40*n_dimensions,
# THIS currently needs to be set in src if used: "recompute_cov_at" => 10,
# THIS currently needs to be set in src if used: "recompute_cov_at" => 10,
)
if case == "Prior"
# don't do anything
Expand All @@ -130,7 +130,7 @@ function main()
rank = n_dimensions #<= 10 ? n_dimensions : 10
kernel_structure = SeparableKernel(LowRankFactor(rank, nugget), OneDimFactor())
n_features = n_dimensions <= 10 ? n_dimensions * 100 : 1000
if (n_features/n_train_pts > 0.9) && (n_features/n_train_pts < 1.1)
if (n_features / n_train_pts > 0.9) && (n_features / n_train_pts < 1.1)
@warn "The number of features similar to the number of training points, poor performance expected, change one or other of these"
end
mlt = ScalarRandomFeatureInterface(
Expand Down Expand Up @@ -159,55 +159,63 @@ function main()
GC.gc() #collect garbage

# PLotting:
if rep_idx == 1
if rep_idx == 1
f3, ax3, plt3 = scatter(
1:n_dimensions,
result_preds[1][:firstorder];
result_preds[1][:firstorder];
color = :red,
markersize =8,
markersize = 8,
marker = :cross,
label = "V-emulate",
title = "input dimension: $(n_dimensions)",
)
scatter!(ax3, result[:firstorder], color = :red, markersize = 8, label="V-approx")
scatter!(ax3, V, color = :red, markersize = 12, marker = :xcross, label="V-true")
scatter!(ax3, result[:firstorder], color = :red, markersize = 8, label = "V-approx")
scatter!(ax3, V, color = :red, markersize = 12, marker = :xcross, label = "V-true")
scatter!(
ax3,
1:n_dimensions,
result_preds[1][:totalorder];
color = :blue,
label = "TV-emulate",
markersize =8,
markersize = 8,
marker = :cross,
)
scatter!(ax3, result[:totalorder], color = :blue, markersize = 8,label="TV-approx")
scatter!(ax3, TV, color = :blue, markersize = 12, marker = :xcross, label="TV-true")
scatter!(ax3, result[:totalorder], color = :blue, markersize = 8, label = "TV-approx")
scatter!(ax3, TV, color = :blue, markersize = 12, marker = :xcross, label = "TV-true")
axislegend(ax3)

CairoMakie.save(joinpath(output_directory, "GFunction_sens_$(case)_$(n_dimensions).png"), f3, px_per_unit = 3)
CairoMakie.save(joinpath(output_directory, "GFunction_sens_$(case)_$(n_dimensions).pdf"), f3, px_per_unit = 3)

CairoMakie.save(
joinpath(output_directory, "GFunction_sens_$(case)_$(n_dimensions).png"),
f3,
px_per_unit = 3,
)
CairoMakie.save(
joinpath(output_directory, "GFunction_sens_$(case)_$(n_dimensions).pdf"),
f3,
px_per_unit = 3,
)
else
# get percentiles:
fo_mat = zeros(n_dimensions,rep_idx)
to_mat = zeros(n_dimensions,rep_idx)
for (idx,rp) in enumerate(result_preds)
fo_mat[:,idx] = rp[:firstorder]
to_mat[:,idx] = rp[:totalorder]
fo_mat = zeros(n_dimensions, rep_idx)
to_mat = zeros(n_dimensions, rep_idx)

for (idx, rp) in enumerate(result_preds)
fo_mat[:, idx] = rp[:firstorder]
to_mat[:, idx] = rp[:totalorder]
end
firstorder_med = percentile.(eachrow(fo_mat),50)
firstorder_low = percentile.(eachrow(fo_mat),5)
firstorder_up = percentile.(eachrow(fo_mat),95)
totalorder_med = percentile.(eachrow(to_mat),50)
totalorder_low = percentile.(eachrow(to_mat),5)
totalorder_up = percentile.(eachrow(to_mat),95)

firstorder_med = percentile.(eachrow(fo_mat), 50)
firstorder_low = percentile.(eachrow(fo_mat), 5)
firstorder_up = percentile.(eachrow(fo_mat), 95)

totalorder_med = percentile.(eachrow(to_mat), 50)
totalorder_low = percentile.(eachrow(to_mat), 5)
totalorder_up = percentile.(eachrow(to_mat), 95)

println("(50%) firstorder: ", firstorder_med)
println("(5%) firstorder: ", firstorder_low)
println("(95%) firstorder: ", firstorder_up)

println("(50%) totalorder: ", totalorder_med)
println("(5%) totalorder: ", totalorder_low)
println("(95%) totalorder: ", totalorder_up)
Expand All @@ -222,8 +230,8 @@ function main()
label = "V-emulate",
title = "input dimension: $(n_dimensions)",
)
scatter!(ax3, result[:firstorder], color = :red, markersize = 8, label="V-approx")
scatter!(ax3, V, color = :red, markersize = 12, marker = :xcross, label="V-true")
scatter!(ax3, result[:firstorder], color = :red, markersize = 8, label = "V-approx")
scatter!(ax3, V, color = :red, markersize = 12, marker = :xcross, label = "V-true")
errorbars!(
ax3,
1:n_dimensions,
Expand All @@ -234,13 +242,21 @@ function main()
color = :blue,
label = "TV-emulate",
)
scatter!(ax3, result[:totalorder], color = :blue, markersize = 8,label="TV-approx")
scatter!(ax3, TV, color = :blue, markersize = 12, marker = :xcross, label="TV-true")
scatter!(ax3, result[:totalorder], color = :blue, markersize = 8, label = "TV-approx")
scatter!(ax3, TV, color = :blue, markersize = 12, marker = :xcross, label = "TV-true")
axislegend(ax3)

CairoMakie.save(joinpath(output_directory, "GFunction_sens_$(case)_$(n_dimensions).png"), f3, px_per_unit = 3)
CairoMakie.save(joinpath(output_directory, "GFunction_sens_$(case)_$(n_dimensions).pdf"), f3, px_per_unit = 3)


CairoMakie.save(
joinpath(output_directory, "GFunction_sens_$(case)_$(n_dimensions).png"),
f3,
px_per_unit = 3,
)
CairoMakie.save(
joinpath(output_directory, "GFunction_sens_$(case)_$(n_dimensions).pdf"),
f3,
px_per_unit = 3,
)

end
# plots - first 3 dimensions
if rep_idx == 1
Expand All @@ -250,8 +266,16 @@ function main()
scatter!(ax2, samples[:, i], y_preds[1][:], color = :blue)
scatter!(ax2, samples[ind, i], y[ind] + noise, color = :red, markersize = 8)
end
CairoMakie.save(joinpath(output_directory, "GFunction_slices_$(case)_$(n_dimensions).png"), f2, px_per_unit = 3)
CairoMakie.save(joinpath(output_directory, "GFunction_slices_$(case)_$(n_dimensions).pdf"), f2, px_per_unit = 3)
CairoMakie.save(
joinpath(output_directory, "GFunction_slices_$(case)_$(n_dimensions).png"),
f2,
px_per_unit = 3,
)
CairoMakie.save(
joinpath(output_directory, "GFunction_slices_$(case)_$(n_dimensions).pdf"),
f2,
px_per_unit = 3,
)
end
end

Expand All @@ -271,19 +295,19 @@ function main()
println("***************************************************************")


jldsave(
joinpath(output_directory, "Gfunction_$(case)_$(n_dimensions).jld2");
sobol_pts = samples,
train_idx = ind,
analytic_V=V,
analytic_TV=TV,
estimated_sobol = result,
mlt_sobol=result_preds,
mlt_pred_y=y_preds,
true_y=y,
noise_y=Γ,
observed_y = output,
)
jldsave(
joinpath(output_directory, "Gfunction_$(case)_$(n_dimensions).jld2");
sobol_pts = samples,
train_idx = ind,
analytic_V = V,
analytic_TV = TV,
estimated_sobol = result,
mlt_sobol = result_preds,
mlt_pred_y = y_preds,
true_y = y,
noise_y = Γ,
observed_y = output,
)


return y_preds, result_preds
Expand Down
Loading

0 comments on commit 29bfe59

Please sign in to comment.