diff --git a/examples/EDMF_data/plot_posterior.jl b/examples/EDMF_data/plot_posterior.jl index 712fd241..30b71e50 100644 --- a/examples/EDMF_data/plot_posterior.jl +++ b/examples/EDMF_data/plot_posterior.jl @@ -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)) @@ -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) @@ -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) @@ -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) - diff --git a/examples/EDMF_data/uq_for_edmf.jl b/examples/EDMF_data/uq_for_edmf.jl index bcb525ed..1f2a1f9f 100644 --- a/examples/EDMF_data/uq_for_edmf.jl +++ b/examples/EDMF_data/uq_for_edmf.jl @@ -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) diff --git a/examples/Emulator/G-function/emulate.jl b/examples/Emulator/G-function/emulate.jl index b4358283..84706b2c 100644 --- a/examples/Emulator/G-function/emulate.jl +++ b/examples/Emulator/G-function/emulate.jl @@ -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) @@ -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) @@ -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 @@ -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 @@ -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( @@ -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) @@ -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, @@ -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 @@ -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 @@ -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 diff --git a/examples/Emulator/G-function/plot_result.jl b/examples/Emulator/G-function/plot_result.jl index c98f93c5..17cc0230 100644 --- a/examples/Emulator/G-function/plot_result.jl +++ b/examples/Emulator/G-function/plot_result.jl @@ -8,8 +8,8 @@ function main() cases = ["Prior", "GP", "RF-scalar"] case = cases[3] n_dimensions = 3 - filename = joinpath(output_directory,"Gfunction_$(case)_$(n_dimensions).jld2") - + filename = joinpath(output_directory, "Gfunction_$(case)_$(n_dimensions).jld2") + ( sobol_pts, train_idx, @@ -17,11 +17,11 @@ function main() mlt_sobol, analytic_V, analytic_TV, - true_y, + true_y, noise_y, observed_y, estimated_sobol, - )= load( + ) = load( filename, "sobol_pts", "train_idx", @@ -36,56 +36,56 @@ function main() ) n_repeats = length(mlt_sobol) - + if n_repeats == 1 f3, ax3, plt3 = scatter( 1:n_dimensions, - mlt_sobol[1][:firstorder]; + mlt_sobol[1][:firstorder]; color = :red, - markersize =8, + markersize = 8, marker = :cross, label = "V-emulate", title = "input dimension: $(n_dimensions)", ) - scatter!(ax3, estimated_sobol[:firstorder], color = :red, markersize = 8, label="V-approx") - scatter!(ax3, analytic_V, color = :red, markersize = 12, marker = :xcross, label="V-true") + scatter!(ax3, estimated_sobol[:firstorder], color = :red, markersize = 8, label = "V-approx") + scatter!(ax3, analytic_V, color = :red, markersize = 12, marker = :xcross, label = "V-true") scatter!( ax3, 1:n_dimensions, mlt_sobol[1][:totalorder]; color = :blue, label = "TV-emulate", - markersize =8, + markersize = 8, marker = :cross, ) - scatter!(ax3, estimated_sobol[:totalorder], color = :blue, markersize = 8,label="TV-approx") - scatter!(ax3, analytic_TV, color = :blue, markersize = 12, marker = :xcross, label="TV-true") + scatter!(ax3, estimated_sobol[:totalorder], color = :blue, markersize = 8, label = "TV-approx") + scatter!(ax3, analytic_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) else # get percentiles: - fo_mat = zeros(n_dimensions,n_repeats) - to_mat = zeros(n_dimensions,n_repeats) - - for (idx,rp) in enumerate(mlt_sobol) - fo_mat[:,idx] = rp[:firstorder] - to_mat[:,idx] = rp[:totalorder] + fo_mat = zeros(n_dimensions, n_repeats) + to_mat = zeros(n_dimensions, n_repeats) + + for (idx, rp) in enumerate(mlt_sobol) + 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) @@ -100,8 +100,8 @@ function main() label = "V-emulate", title = "input dimension: $(n_dimensions)", ) - scatter!(ax3, estimated_sobol[:firstorder], color = :red, markersize = 8, label="V-approx") - scatter!(ax3, analytic_V, color = :red, markersize = 12, marker = :xcross, label="V-true") + scatter!(ax3, estimated_sobol[:firstorder], color = :red, markersize = 8, label = "V-approx") + scatter!(ax3, analytic_V, color = :red, markersize = 12, marker = :xcross, label = "V-true") errorbars!( ax3, 1:n_dimensions, @@ -112,13 +112,13 @@ function main() color = :blue, label = "TV-emulate", ) - scatter!(ax3, estimated_sobol[:totalorder], color = :blue, markersize = 8,label="TV-approx") - scatter!(ax3, analytic_TV, color = :blue, markersize = 12, marker = :xcross, label="TV-true") + scatter!(ax3, estimated_sobol[:totalorder], color = :blue, markersize = 8, label = "TV-approx") + scatter!(ax3, analytic_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) - + end # plots - first 3 dimensions plot_dim = n_dimensions >= 3 ? 3 : n_dimensions diff --git a/examples/Emulator/Ishigami/emulate.jl b/examples/Emulator/Ishigami/emulate.jl index f64ad7c8..84d9b19a 100644 --- a/examples/Emulator/Ishigami/emulate.jl +++ b/examples/Emulator/Ishigami/emulate.jl @@ -1,4 +1,5 @@ + using GlobalSensitivityAnalysis const GSA = GlobalSensitivityAnalysis using Distributions @@ -127,8 +128,8 @@ function main() result_pred = analyze(data, y_pred') push!(y_preds, y_pred) push!(result_preds, result_pred) - - jldsave(joinpath(output_directory,"emulator_repeat_$(rep_idx)_$(case).jld2"); emulator) + + jldsave(joinpath(output_directory, "emulator_repeat_$(rep_idx)_$(case).jld2"); emulator) end @@ -144,15 +145,15 @@ function main() VT3 = 8 * b^2 * π^8 / 225 jldsave( - joinpath(output_directory,"results_$case.jld2"); + joinpath(output_directory, "results_$case.jld2"); sobol_pts = samples, train_idx = ind, mlt_pred_y = y_preds, - mlt_sobol= result_preds, - analytic_sobol = [V1,V2,V3,VT1,VT2,VT3], - true_y=y, + mlt_sobol = result_preds, + analytic_sobol = [V1, V2, V3, VT1, VT2, VT3], + true_y = y, noise_y = Γ, - estimated_sobol = result + estimated_sobol = result, ) println(" ") diff --git a/examples/Emulator/Ishigami/plot_result.jl b/examples/Emulator/Ishigami/plot_result.jl index 4bcf2e98..a86f1f96 100644 --- a/examples/Emulator/Ishigami/plot_result.jl +++ b/examples/Emulator/Ishigami/plot_result.jl @@ -7,18 +7,9 @@ function main() output_directory = "output" cases = ["Prior", "GP", "RF-scalar"] case = cases[3] - filename = joinpath(output_directory,"results_$case.jld2") - - ( - sobol_pts, - train_idx, - mlt_pred_y, - mlt_sobol, - analytic_sobol, - true_y, - noise_y, - estimated_sobol, - )= load( + filename = joinpath(output_directory, "results_$case.jld2") + + (sobol_pts, train_idx, mlt_pred_y, mlt_sobol, analytic_sobol, true_y, noise_y, estimated_sobol) = load( filename, "sobol_pts", "train_idx", @@ -31,7 +22,7 @@ function main() ) n_repeats = length(mlt_sobol) - (V,V1,V2,V3,VT1,VT2,VT3) = analytic_sobol + (V, V1, V2, V3, VT1, VT2, VT3) = analytic_sobol f1 = Figure(resolution = (1.618 * 900, 300), markersize = 4) axx = Axis(f1[1, 1], xlabel = "x1", ylabel = "f") @@ -45,7 +36,7 @@ function main() save(joinpath(output_directory, "ishigami_slices_truth.png"), f1, px_per_unit = 3) save(joinpath(output_directory, "ishigami_slices_truth.pdf"), f1, px_per_unit = 3) - + # display some info println(" ") println("True Sobol Indices") @@ -93,7 +84,7 @@ function main() save(joinpath(output_directory, "ishigami_slices_$(case).pdf"), f2, px_per_unit = 3) - + end diff --git a/examples/Emulator/L63/emulate.jl b/examples/Emulator/L63/emulate.jl index 699e5617..5c12bd1b 100644 --- a/examples/Emulator/L63/emulate.jl +++ b/examples/Emulator/L63/emulate.jl @@ -27,7 +27,7 @@ function main() # rng rng = MersenneTwister(1232434) - n_repeats = 20 # repeat exp with same data. + n_repeats = 20 # repeat exp with same data. println("run experiment $n_repeats times") @@ -108,12 +108,12 @@ function main() "n_features_opt" => 150, "n_iteration" => 20, "accelerator" => NesterovAccelerator(), -# "localization" => EnsembleKalmanProcesses.Localizers.SECNice(0.05,1.0), # localization / s + # "localization" => EnsembleKalmanProcesses.Localizers.SECNice(0.05,1.0), # localization / s "n_ensemble" => 200, "verbose" => true, ) - + # Build ML tools if case == "GP" gppackage = Emulators.GPJL() @@ -134,7 +134,7 @@ function main() kernel_structure = kernel_structure, optimizer_options = rf_optimizer_overrides, ) - + elseif case ∈ ["RF-scalar", "RF-scalar-diagin"] n_features = 10 * Int(floor(sqrt(3 * n_tp))) kernel_structure = @@ -187,12 +187,15 @@ function main() if !(case == "GP") && (rep_idx == 1) JLD2.save( joinpath(output_directory, case * "_l63_config.jld2"), - "rf_optimizer_overrides", rf_optimizer_overrides, - "n_features", n_features, - "kernel_structure", kernel_structure, + "rf_optimizer_overrides", + rf_optimizer_overrides, + "n_features", + n_features, + "kernel_structure", + kernel_structure, ) end - + # Emulate if case ∈ ["RF-nosvd-nonsep", "RF-nosvd-sep"] decorrelate = false diff --git a/examples/Emulator/L63/plot_results.jl b/examples/Emulator/L63/plot_results.jl index fdf5386a..1ae772d9 100644 --- a/examples/Emulator/L63/plot_results.jl +++ b/examples/Emulator/L63/plot_results.jl @@ -10,13 +10,13 @@ function main() case = cases[2] # Load from saved files -#= - config_file = JLD2.load(joinpath(output_directory, case * "_l63_config.jld2")) - rf_optimizer_overrides = config_file["rf_optimizer_overrides"] - n_features = config_file["n_features"] - kernel_structure = config_file["kernel_structure"] - =# - + #= + config_file = JLD2.load(joinpath(output_directory, case * "_l63_config.jld2")) + rf_optimizer_overrides = config_file["rf_optimizer_overrides"] + n_features = config_file["n_features"] + kernel_structure = config_file["kernel_structure"] + =# + histogram_data = JLD2.load(joinpath(output_directory, case * "_l63_histdata.jld2")) solhist = histogram_data["solhist"] uhist = histogram_data["uhist"] @@ -27,13 +27,13 @@ function main() # plotting trajectories for just first repeat uplot_tmp = uplot[1] uhist_tmp = uhist[1] - + # copied from emulate.jl dt = 0.01 tmax_test = 20 #100 xx = 0.0:dt:tmax_test - + f = Figure(size = (900, 450)) axx = Axis(f[1, 1], xlabel = "time", ylabel = "x") axy = Axis(f[2, 1], xlabel = "time", ylabel = "y") diff --git a/src/RandomFeature.jl b/src/RandomFeature.jl index 8fcff07a..441ad61d 100644 --- a/src/RandomFeature.jl +++ b/src/RandomFeature.jl @@ -489,13 +489,13 @@ function calculate_mean_cov_and_coeffs( # sizes (output_dim x n_test), (output_dim x output_dim x n_test) ## TODO - the theory states that the following should be set: -# scaled_coeffs = sqrt(1 / (n_features)) * RF.Methods.get_coeffs(fitted_features) - scaled_coeffs = 1e-3*rand(size(scaled_coeffs))#overwrite with noise... + # scaled_coeffs = sqrt(1 / (n_features)) * RF.Methods.get_coeffs(fitted_features) + scaled_coeffs = 1e-3 * rand(size(scaled_coeffs))#overwrite with noise... # However the convergence is much improved with setting this to zero: #scaled_coeffs = 0 if decomp_type == "cholesky" chol_fac = RF.Methods.get_decomposition(RF.Methods.get_feature_factors(fitted_features)).L - + complexity = 2 * sum(log(chol_fac[i, i]) for i in 1:size(chol_fac, 1)) else svd_singval = RF.Methods.get_decomposition(RF.Methods.get_feature_factors(fitted_features)).S @@ -545,7 +545,7 @@ $(DocStringExtensions.TYPEDSIGNATURES) Calculate the empirical covariance, additionally applying the Noise Informed Covariance Estimator (NICE) Vishnay et al. 2024. """ -function nice_cov(sample_mat::AA, n_samples=400, δ::FT=1.0) where {AA <: AbstractMatrix, FT <: Real} +function nice_cov(sample_mat::AA, n_samples = 400, δ::FT = 1.0) where {AA <: AbstractMatrix, FT <: Real} n_sample_cov = size(sample_mat, 2) Γ = cov(sample_mat, dims = 2) @@ -562,7 +562,7 @@ function nice_cov(sample_mat::AA, n_samples=400, δ::FT=1.0) where {AA <: Abstra interp_steps = 100 # use find the variability in the corr coeff matrix entries std_corrs = approximate_corr_std.(corr, n_sample_cov, n_samples) # Found in EKP.Localizers !! slowest part of code -> could speed up by precomputing an interpolation of [-1,1] - + std_tol = sqrt(sum(std_corrs .^ 2)) α_min_exceeded = [max_exponent] for α in 2:2:max_exponent # even exponents give a PSD @@ -575,14 +575,14 @@ function nice_cov(sample_mat::AA, n_samples=400, δ::FT=1.0) where {AA <: Abstra end corr_psd = corr .^ α_min_exceeded[1] corr_psd_prev = corr .^ (α_min_exceeded[1] - 2) # previous PSD correction - + for α in LinRange(1.0, 0.0, interp_steps) corr_interp = ((1 - α) * (corr_psd_prev) + α * corr_psd) .* corr if norm(corr_interp - corr) < δ * std_tol - corr[:,:] = corr_interp #update the correlation matrix block + corr[:, :] = corr_interp #update the correlation matrix block break end - end + end return posdef_correct(V * corr * V) # rebuild the cov matrix @@ -672,11 +672,11 @@ function estimate_mean_and_coeffnorm_covariance( approx_σ2[id, id] = mean_of_covs[i, :, :] # this ordering, so we can take a mean/cov in dims = 2. blockmeans[id, :] = permutedims(means[i, :, :], (2, 1)) end - + sample_mat = vcat(blockmeans, coeffl2norm, complexity) correction = "shrinkage" -# correction = "nice" + # correction = "nice" if correction == "shrinkage" Γ = shrinkage_cov(sample_mat) elseif correction == "nice" @@ -886,9 +886,9 @@ function estimate_mean_and_coeffnorm_covariance( sample_mat = vcat(blockmeans, coeffl2norm, complexity) - + correction = "shrinkage" -# correction = "nice" + # correction = "nice" if correction == "shrinkage" Γ = shrinkage_cov(sample_mat) elseif correction == "nice" diff --git a/src/ScalarRandomFeature.jl b/src/ScalarRandomFeature.jl index dceb1e10..3e0f44a3 100644 --- a/src/ScalarRandomFeature.jl +++ b/src/ScalarRandomFeature.jl @@ -156,7 +156,6 @@ function ScalarRandomFeatureInterface( "accelerator" => EKP.DefaultAccelerator(), # acceleration with momentum "recompute_cov_at" => Inf, # whether to recompute the output_covariance "localization" => EKP.Localizers.NoLocalization(), # localization / sample error correction for small ensembles - ) if !isnothing(optimizer_options) @@ -402,32 +401,34 @@ function build_models!( initial_params = construct_initial_ensemble(rng, prior, n_ensemble) data = vcat(get_outputs(io_pairs_opt)[(n_train + 1):end], 0.0, 0.0) - ekiobj = [EKP.EnsembleKalmanProcess( - initial_params, - data, - Γ, - Inversion(), - scheduler = scheduler, - rng = rng, - accelerator = accelerator, - verbose = opt_verbose_flag, - localization_method = localization, - )] + ekiobj = [ + EKP.EnsembleKalmanProcess( + initial_params, + data, + Γ, + Inversion(), + scheduler = scheduler, + rng = rng, + accelerator = accelerator, + verbose = opt_verbose_flag, + localization_method = localization, + ), + ] err = zeros(n_iteration) # [4.] optimize with EKP Γ_old = [Γ] for i in 1:n_iteration - + #get parameters: lvec = transform_unconstrained_to_constrained(prior, get_u_final(ekiobj[1])) if i % recompute_cov_at == 0 - + internal_Γ_new, approx_σ2_new = estimate_mean_and_coeffnorm_covariance( srfi, rng, - transform_unconstrained_to_constrained(prior, mean(get_u_final(ekiobj[1]),dims=2)), + transform_unconstrained_to_constrained(prior, mean(get_u_final(ekiobj[1]), dims = 2)), regularization, n_features_opt, n_train, @@ -438,10 +439,10 @@ function build_models!( decomp_type, multithread_type, ) - + Γ_new = internal_Γ_new Γ_new[1:n_test, 1:n_test] += regularization # + approx_σ2 - Γ_new[(n_test + 1):end, (n_test + 1):end] += I + Γ_new[(n_test + 1):end, (n_test + 1):end] += I if opt_verbose_flag @info "updated algorithm covariance difference: $(norm(Γ_new - Γ_old[1]))" end @@ -453,15 +454,15 @@ function build_models!( data, Γ_new, Inversion(), - scheduler = DataMisfitController(terminate_at=1e4), + scheduler = DataMisfitController(terminate_at = 1e4), rng = rng, accelerator = DefaultAccelerator(), localization_method = Localizers.SECNice(200), verbose = opt_verbose_flag, ) - + end - + g_ens, _ = calculate_ensemble_mean_and_coeffnorm( srfi, rng, @@ -478,8 +479,7 @@ function build_models!( inflation = optimizer_options["inflation"] if inflation > 0 - terminated = - EKP.update_ensemble!(ekiobj[1], g_ens, additive_inflation = true, s = inflation) # small regularizing inflation + terminated = EKP.update_ensemble!(ekiobj[1], g_ens, additive_inflation = true, s = inflation) # small regularizing inflation else terminated = EKP.update_ensemble!(ekiobj[1], g_ens) # small regularizing inflation end diff --git a/src/VectorRandomFeature.jl b/src/VectorRandomFeature.jl index 87473b16..b6fffcaf 100644 --- a/src/VectorRandomFeature.jl +++ b/src/VectorRandomFeature.jl @@ -408,13 +408,13 @@ function build_models!( else # think of the regularization_matrix as the observational noise covariance, or a related quantity if !isposdef(regularization_matrix) - regularization= posdef_correct(regularization_matrix) + regularization = posdef_correct(regularization_matrix) println("RF regularization matrix is not positive definite, correcting") else regularization = regularization_matrix end - + end idx_shuffle = randperm(rng, n_data) @@ -458,7 +458,8 @@ function build_models!( if tikhonov_opt_val == 0 # Build the covariance Γ = internal_Γ - Γ[1:(n_test * output_dim), 1:(n_test * output_dim)] += isa(regularization, UniformScaling) ? regularization : kron(I(n_test), regularization) # + approx_σ2 + Γ[1:(n_test * output_dim), 1:(n_test * output_dim)] += + isa(regularization, UniformScaling) ? regularization : kron(I(n_test), regularization) # + approx_σ2 Γ[(n_test * output_dim + 1):end, (n_test * output_dim + 1):end] += I @info Γ[(n_test * output_dim + 1):end, (n_test * output_dim + 1):end] data = vcat(reshape(get_outputs(io_pairs_opt)[:, (n_train + 1):end], :, 1), 0.0, 0.0) #flatten data @@ -468,7 +469,7 @@ function build_models!( outsize = size(internal_Γ, 1) Γ = zeros(outsize + n_hp, outsize + n_hp) Γ[1:outsize, 1:outsize] = internal_Γ - Γ[1:(n_test * output_dim), 1:(n_test * output_dim)] += kron(I(n_test),regularization) # block diag regularization + Γ[1:(n_test * output_dim), 1:(n_test * output_dim)] += kron(I(n_test), regularization) # block diag regularization Γ[(n_test * output_dim + 1):outsize, (n_test * output_dim + 1):outsize] += I Γ[(outsize + 1):end, (outsize + 1):end] = tikhonov_opt_val .* cov(prior) @@ -545,8 +546,7 @@ function build_models!( end inflation = optimizer_options["inflation"] if inflation > 0 - terminated = - EKP.update_ensemble!(ekiobj, g_ens, additive_inflation = true, s = inflation) # small regularizing inflation + terminated = EKP.update_ensemble!(ekiobj, g_ens, additive_inflation = true, s = inflation) # small regularizing inflation else terminated = EKP.update_ensemble!(ekiobj, g_ens) # small regularizing inflation end