diff --git a/Project.toml b/Project.toml index c431ee511..74ed7f67d 100644 --- a/Project.toml +++ b/Project.toml @@ -15,9 +15,12 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +RandomFeatures = "36c3bae2-c0c3-419d-b3b4-eebadd35c5e5" ScikitLearn = "3646fa90-6ef7-5e7e-9f22-8aca16db6324" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" @@ -31,9 +34,10 @@ EnsembleKalmanProcesses = "0.8, 0.9, 0.10, 0.11, 0.12, 0.13, 0.14" GaussianProcesses = "0.12" MCMCChains = "4.14, 5" PyCall = "1.93" +RandomFeatures = "0.2" ScikitLearn = "0.6" StatsBase = "0.33" -julia = "1.6" +julia = "1.6, 1.7, 1.8" [extras] Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" diff --git a/examples/Emulator/GaussianProcess/plot_GP.jl b/examples/Emulator/GaussianProcess/plot_GP.jl index 5a815969c..a0ee9751f 100644 --- a/examples/Emulator/GaussianProcess/plot_GP.jl +++ b/examples/Emulator/GaussianProcess/plot_GP.jl @@ -68,9 +68,9 @@ if !isdir(output_directory) end #create the machine learning tools: Gaussian Process -gppackage = GPJL() +gppackage = SKLJL() pred_type = YType() -gaussian_process = GaussianProcess(gppackage, noise_learn = true) +gaussian_process = GaussianProcess(gppackage, noise_learn = false) # Generate training data (x-y pairs, where x ∈ ℝ ᵖ, y ∈ ℝ ᵈ) # x = [x1, x2]: inputs/predictors/features/parameters @@ -92,7 +92,7 @@ gx[2, :] = g2x # Add noise η μ = zeros(d) -Σ = 0.1 * [[0.8, 0.0] [0.0, 0.5]] # d x d +Σ = 0.1 * [[0.8, 0.1] [0.1, 0.5]] # d x d noise_samples = rand(MvNormal(μ, Σ), n) # y = G(x) + η Y = gx .+ noise_samples @@ -182,9 +182,9 @@ println("GP trained") # Plot mean and variance of the predicted observables y1 and y2 # For this, we generate test points on a x1-x2 grid. -n_pts = 50 -x1 = range(0.0, stop = 2 * π, length = n_pts) -x2 = range(0.0, stop = 2 * π, length = n_pts) +n_pts = 200 +x1 = range(0.0, stop = (4.0 / 5.0) * 2 * π, length = n_pts) +x2 = range(0.0, stop = (4.0 / 5.0) * 2 * π, length = n_pts) X1, X2 = meshgrid(x1, x2) # Input for predict has to be of size N_samples x input_dim inputs = permutedims(hcat(X1[:], X2[:]), (2, 1)) diff --git a/examples/Emulator/RandomFeature/Project.toml b/examples/Emulator/RandomFeature/Project.toml new file mode 100644 index 000000000..9d41c3993 --- /dev/null +++ b/examples/Emulator/RandomFeature/Project.toml @@ -0,0 +1,15 @@ +[deps] +CalibrateEmulateSample = "95e48a1f-0bec-4818-9538-3db4340308e3" +Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[compat] +FiniteDiff = "~2.10" +julia = "~1.6" diff --git a/examples/Emulator/RandomFeature/scalar_optimize_and_plot_RF.jl b/examples/Emulator/RandomFeature/scalar_optimize_and_plot_RF.jl new file mode 100644 index 000000000..a1a0fce93 --- /dev/null +++ b/examples/Emulator/RandomFeature/scalar_optimize_and_plot_RF.jl @@ -0,0 +1,257 @@ +# Reference the in-tree version of CalibrateEmulateSample on Julias load path +include(joinpath(@__DIR__, "..", "..", "ci", "linkfig.jl")) + +# Import modules +using Random +using StableRNGs +using Distributions +using Statistics +using LinearAlgebra +using CalibrateEmulateSample.Emulators +using CalibrateEmulateSample.DataContainers +using CalibrateEmulateSample.ParameterDistributions + +case = "scalar" +println("running case $case") + +plot_flag = true +if plot_flag + using Plots + gr(size = (1500, 700)) + Plots.scalefontsizes(1.3) + font = Plots.font("Helvetica", 18) + fontdict = Dict(:guidefont => font, :xtickfont => font, :ytickfont => font, :legendfont => font) + +end + +function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where {T} + m, n = length(vy), length(vx) + gx = reshape(repeat(vx, inner = m, outer = 1), m, n) + gy = reshape(repeat(vy, inner = 1, outer = n), m, n) + + return gx, gy +end +rng_seed = 41 +Random.seed!(rng_seed) +output_directory = joinpath(@__DIR__, "output") +if !isdir(output_directory) + mkdir(output_directory) +end + +#problem +n = 100 # number of training points +p = 2 # input dim +d = 2 # output dim +X = 2.0 * π * rand(p, n) +# G(x1, x2) +g1x = sin.(X[1, :]) .+ cos.(X[2, :]) +g2x = sin.(X[1, :]) .- cos.(X[2, :]) +gx = zeros(2, n) +gx[1, :] = g1x +gx[2, :] = g2x + +# Add noise η +μ = zeros(d) +Σ = 0.1 * [[0.8, 0.1] [0.1, 0.5]] # d x d +noise_samples = rand(MvNormal(μ, Σ), n) +# y = G(x) + η +Y = gx .+ noise_samples + +iopairs = PairedDataContainer(X, Y, data_are_columns = true) +@assert get_inputs(iopairs) == X +@assert get_outputs(iopairs) == Y + +#plot training data with and without noise +if plot_flag + p1 = plot( + X[1, :], + X[2, :], + g1x, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) + + figpath = joinpath(output_directory, "RF_" * case * "_observed_y1nonoise.png") + savefig(figpath) + + p2 = plot( + X[1, :], + X[2, :], + g2x, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) + figpath = joinpath(output_directory, "RF_" * case * "_observed_y2nonoise.png") + savefig(figpath) + + p3 = plot( + X[1, :], + X[2, :], + Y[1, :], + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) + figpath = joinpath(output_directory, "RF_" * case * "_observed_y1.png") + savefig(figpath) + + p4 = plot( + X[1, :], + X[2, :], + Y[2, :], + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) + figpath = joinpath(output_directory, "RF_" * case * "_observed_y2.png") + savefig(figpath) + +end + +# setup random features +n_features = 200 + +srfi = ScalarRandomFeatureInterface(n_features, p) +emulator = Emulator(srfi, iopairs, obs_noise_cov = Σ, normalize_inputs = true) +println("build RF with $n training points and $(n_features) random features.") + +optimize_hyperparameters!(emulator) # although RF already optimized + +# Plot mean and variance of the predicted observables y1 and y2 +# For this, we generate test points on a x1-x2 grid. +n_pts = 200 +x1 = range(0.0, stop = 4.0 / 5.0 * 2 * π, length = n_pts) +x2 = range(0.0, stop = 4.0 / 5.0 * 2 * π, length = n_pts) +X1, X2 = meshgrid(x1, x2) +# Input for predict has to be of size N_samples x input_dim +inputs = permutedims(hcat(X1[:], X2[:]), (2, 1)) + +rf_mean, rf_cov = predict(emulator, inputs, transform_to_real = true) +println("end predictions at ", n_pts * n_pts, " points") + + +#plot predictions +for y_i in 1:d + rf_var_temp = [diag(rf_cov[j]) for j in 1:length(rf_cov)] # (40000,) + rf_var = permutedims(vcat([x' for x in rf_var_temp]...), (2, 1)) # 2 x 40000 + + mean_grid = reshape(rf_mean[y_i, :], n_pts, n_pts) # 2 x 40000 + if plot_flag + p5 = plot( + x1, + x2, + mean_grid, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zlabel = "mean of y" * string(y_i), + zguidefontrotation = 90, + ) + end + var_grid = reshape(rf_var[y_i, :], n_pts, n_pts) + if plot_flag + p6 = plot( + x1, + x2, + var_grid, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zlabel = "var of y" * string(y_i), + zguidefontrotation = 90, + ) + + plot(p5, p6, layout = (1, 2), legend = false) + + savefig(joinpath(output_directory, "RF_" * case * "_y" * string(y_i) * "_predictions.png")) + end +end + +# Plot the true components of G(x1, x2) +g1_true = sin.(inputs[1, :]) .+ cos.(inputs[2, :]) +g1_true_grid = reshape(g1_true, n_pts, n_pts) +if plot_flag + p7 = plot( + x1, + x2, + g1_true_grid, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zlabel = "sin(x1) + cos(x2)", + zguidefontrotation = 90, + ) + savefig(joinpath(output_directory, "RF_" * case * "_true_g1.png")) +end + +g2_true = sin.(inputs[1, :]) .- cos.(inputs[2, :]) +g2_true_grid = reshape(g2_true, n_pts, n_pts) +if plot_flag + p8 = plot( + x1, + x2, + g2_true_grid, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zlabel = "sin(x1) - cos(x2)", + zguidefontrotation = 90, + ) + g_true_grids = [g1_true_grid, g2_true_grid] + + savefig(joinpath(output_directory, "RF_" * case * "_true_g2.png")) + +end + +# Plot the difference between the truth and the mean of the predictions +for y_i in 1:d + + # Reshape rf_cov to size N_samples x output_dim + rf_var_temp = [diag(rf_cov[j]) for j in 1:length(rf_cov)] # (40000,) + rf_var = permutedims(vcat([x' for x in rf_var_temp]...), (2, 1)) # 40000 x 2 + + mean_grid = reshape(rf_mean[y_i, :], n_pts, n_pts) + var_grid = reshape(rf_var[y_i, :], n_pts, n_pts) + # Compute and plot 1/variance * (truth - prediction)^2 + + if plot_flag + zlabel = "1/var * (true_y" * string(y_i) * " - predicted_y" * string(y_i) * ")^2" + + p9 = plot( + x1, + x2, + sqrt.(1.0 ./ var_grid .* (g_true_grids[y_i] .- mean_grid) .^ 2), + st = :surface, + camera = (30, 60), + c = :magma, + zlabel = zlabel, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) + + savefig(joinpath(output_directory, "RF_" * case * "_y" * string(y_i) * "_difference_truth_prediction.png")) + end +end diff --git a/examples/Emulator/RandomFeature/vector_optimize_and_plot_RF.jl b/examples/Emulator/RandomFeature/vector_optimize_and_plot_RF.jl new file mode 100644 index 000000000..05e4f348f --- /dev/null +++ b/examples/Emulator/RandomFeature/vector_optimize_and_plot_RF.jl @@ -0,0 +1,286 @@ +# Reference the in-tree version of CalibrateEmulateSample on Julias load path +include(joinpath(@__DIR__, "..", "..", "ci", "linkfig.jl")) + +# Import modules +using Random +using StableRNGs +using Distributions +using Statistics +using LinearAlgebra +using CalibrateEmulateSample.Emulators +using CalibrateEmulateSample.DataContainers +using CalibrateEmulateSample.ParameterDistributions +plot_flag = true +if plot_flag + using Plots + gr(size = (1500, 700)) + Plots.scalefontsizes(1.3) + font = Plots.font("Helvetica", 18) + fontdict = Dict(:guidefont => font, :xtickfont => font, :ytickfont => font, :legendfont => font) + +end + + +function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where {T} + m, n = length(vy), length(vx) + gx = reshape(repeat(vx, inner = m, outer = 1), m, n) + gy = reshape(repeat(vy, inner = 1, outer = n), m, n) + + return gx, gy +end + +function main() + + rng_seed = 41 + Random.seed!(rng_seed) + output_directory = joinpath(@__DIR__, "output") + if !isdir(output_directory) + mkdir(output_directory) + end + + cases = ["svd-diag", "svd-nondiag", "nosvd-diag", "nosvd-nondiag"] + case_mask = 1:4 # which cases to do + + #problem + n = 150 # number of training points + p = 2 # input dim + d = 2 # output dim + X = 2.0 * π * rand(p, n) + # G(x1, x2) + g1x = sin.(X[1, :]) .+ cos.(X[2, :]) + g2x = sin.(X[1, :]) .- cos.(X[2, :]) + gx = zeros(2, n) + gx[1, :] = g1x + gx[2, :] = g2x + + # Add noise η + μ = zeros(d) + Σ = 0.1 * [[0.8, 0.1] [0.1, 0.5]] # d x d + noise_samples = rand(MvNormal(μ, Σ), n) + # y = G(x) + η + Y = gx .+ noise_samples + + iopairs = PairedDataContainer(X, Y, data_are_columns = true) + @assert get_inputs(iopairs) == X + @assert get_outputs(iopairs) == Y + + #plot training data with and without noise + if plot_flag + p1 = plot( + X[1, :], + X[2, :], + g1x, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) + + figpath = joinpath(output_directory, "RF_observed_y1nonoise.png") + savefig(figpath) + + p2 = plot( + X[1, :], + X[2, :], + g2x, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) + figpath = joinpath(output_directory, "RF_observed_y2nonoise.png") + savefig(figpath) + + p3 = plot( + X[1, :], + X[2, :], + Y[1, :], + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) + figpath = joinpath(output_directory, "RF_observed_y1.png") + savefig(figpath) + + p4 = plot( + X[1, :], + X[2, :], + Y[2, :], + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) + figpath = joinpath(output_directory, "RF_observed_y2.png") + savefig(figpath) + + end + + for case in cases[case_mask] + println("running case $case") + + + + + # setup random features + n_features = 200 + + if case == "svd-diag" + vrfi = VectorRandomFeatureInterface(n_features, p, d, diagonalize_output = true) + emulator = Emulator(vrfi, iopairs, obs_noise_cov = Σ, normalize_inputs = true) + elseif case == "svd-nondiag" + vrfi = VectorRandomFeatureInterface(n_features, p, d) + emulator = Emulator(vrfi, iopairs, obs_noise_cov = Σ, normalize_inputs = true) + elseif case == "nosvd-diag" + vrfi = VectorRandomFeatureInterface(n_features, p, d, diagonalize_output = true) + emulator = Emulator(vrfi, iopairs, obs_noise_cov = Σ, normalize_inputs = true, decorrelate = false) + elseif case == "nosvd-nondiag" + vrfi = VectorRandomFeatureInterface(n_features, p, d) + emulator = Emulator(vrfi, iopairs, obs_noise_cov = Σ, normalize_inputs = true, decorrelate = false) + end + println("build RF with $n training points and $(n_features) random features.") + + optimize_hyperparameters!(emulator) # although RF already optimized + + # Plot mean and variance of the predicted observables y1 and y2 + # For this, we generate test points on a x1-x2 grid. + n_pts = 200 + x1 = range(0.0, stop = 4.0 / 5.0 * 2 * π, length = n_pts) + x2 = range(0.0, stop = 4.0 / 5.0 * 2 * π, length = n_pts) + X1, X2 = meshgrid(x1, x2) + # Input for predict has to be of size input_dim x N_samples + inputs = permutedims(hcat(X1[:], X2[:]), (2, 1)) + + rf_mean, rf_cov = predict(emulator, inputs, transform_to_real = true) + println("end predictions at ", n_pts * n_pts, " points") + + + #plot predictions + for y_i in 1:d + + rf_var_temp = [diag(rf_cov[j]) for j in 1:length(rf_cov)] # (40000,) + rf_var = permutedims(vcat([x' for x in rf_var_temp]...), (2, 1)) # 2 x 40000 + + mean_grid = reshape(rf_mean[y_i, :], n_pts, n_pts) # 2 x 40000 + if plot_flag + p5 = plot( + x1, + x2, + mean_grid, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zlabel = "mean of y" * string(y_i), + zguidefontrotation = 90, + ) + end + var_grid = reshape(rf_var[y_i, :], n_pts, n_pts) + if plot_flag + p6 = plot( + x1, + x2, + var_grid, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zlabel = "var of y" * string(y_i), + zguidefontrotation = 90, + ) + + plot(p5, p6, layout = (1, 2), legend = false) + + savefig(joinpath(output_directory, "RF_" * case * "_y" * string(y_i) * "_predictions.png")) + end + end + + # Plot the true components of G(x1, x2) + g1_true = sin.(inputs[1, :]) .+ cos.(inputs[2, :]) + g1_true_grid = reshape(g1_true, n_pts, n_pts) + if plot_flag + p7 = plot( + x1, + x2, + g1_true_grid, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zlabel = "sin(x1) + cos(x2)", + zguidefontrotation = 90, + ) + savefig(joinpath(output_directory, "RF_" * case * "_true_g1.png")) + end + + g2_true = sin.(inputs[1, :]) .- cos.(inputs[2, :]) + g2_true_grid = reshape(g2_true, n_pts, n_pts) + if plot_flag + p8 = plot( + x1, + x2, + g2_true_grid, + st = :surface, + camera = (30, 60), + c = :cividis, + xlabel = "x1", + ylabel = "x2", + zlabel = "sin(x1) - cos(x2)", + zguidefontrotation = 90, + ) + g_true_grids = [g1_true_grid, g2_true_grid] + + savefig(joinpath(output_directory, "RF_" * case * "_true_g2.png")) + + end + MSE = 1 / size(rf_mean, 2) * sqrt(sum((rf_mean[1, :] - g1_true) .^ 2 + (rf_mean[2, :] - g2_true) .^ 2)) + println("L^2 error of mean and latent truth:", MSE) + + # Plot the difference between the truth and the mean of the predictions + for y_i in 1:d + + # Reshape rf_cov to size N_samples x output_dim + rf_var_temp = [diag(rf_cov[j]) for j in 1:length(rf_cov)] # (40000,) + rf_var = permutedims(vcat([x' for x in rf_var_temp]...), (2, 1)) # 40000 x 2 + + mean_grid = reshape(rf_mean[y_i, :], n_pts, n_pts) + var_grid = reshape(rf_var[y_i, :], n_pts, n_pts) + # Compute and plot 1/variance * (truth - prediction)^2 + + if plot_flag + zlabel = "1/var * (true_y" * string(y_i) * " - predicted_y" * string(y_i) * ")^2" + + p9 = plot( + x1, + x2, + sqrt.(1.0 ./ var_grid .* (g_true_grids[y_i] .- mean_grid) .^ 2), + st = :surface, + camera = (30, 60), + c = :magma, + zlabel = zlabel, + xlabel = "x1", + ylabel = "x2", + zguidefontrotation = 90, + ) + + savefig( + joinpath(output_directory, "RF_" * case * "_y" * string(y_i) * "_difference_truth_prediction.png"), + ) + end + end + end +end + +main() diff --git a/examples/GCM/Manifest.toml b/examples/GCM/Manifest.toml new file mode 100644 index 000000000..b5f8553ff --- /dev/null +++ b/examples/GCM/Manifest.toml @@ -0,0 +1,1568 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.8.0" +manifest_format = "2.0" +project_hash = "0cb8f6abcf6e1f799ac049a17e1b0cde938f0a63" + +[[deps.AMD]] +deps = ["Libdl", "LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "00163dc02b882ca5ec032400b919e5f5011dbd31" +uuid = "14f7f29c-3bd6-536c-9a0b-7339e30b5a3e" +version = "0.5.0" + +[[deps.AbstractFFTs]] +deps = ["ChainRulesCore", "LinearAlgebra"] +git-tree-sha1 = "69f7020bd72f069c219b5e8c236c1fa90d2cb409" +uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" +version = "1.2.1" + +[[deps.AbstractMCMC]] +deps = ["BangBang", "ConsoleProgressMonitor", "Distributed", "LogDensityProblems", "Logging", "LoggingExtras", "ProgressLogging", "Random", "StatsBase", "TerminalLoggers", "Transducers"] +git-tree-sha1 = "02b9f1388a7f7a7540ffc6f058c8397a0469add7" +uuid = "80f14c24-f653-4e6a-9b94-39d6b0f70001" +version = "4.2.1" + +[[deps.AbstractTrees]] +git-tree-sha1 = "52b3b436f8f73133d7bc3a6c71ee7ed6ab2ab754" +uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" +version = "0.4.3" + +[[deps.Adapt]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "195c5505521008abea5aee4f96930717958eac6f" +uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +version = "3.4.0" + +[[deps.AdvancedMH]] +deps = ["AbstractMCMC", "Distributions", "Random", "Requires"] +git-tree-sha1 = "d7a7dabeaef34e5106cdf6c2ac956e9e3f97f666" +uuid = "5b7e9947-ddc0-4b3f-9b55-0d8042f74170" +version = "0.6.8" + +[[deps.ArgCheck]] +git-tree-sha1 = "a3a402a35a2f7e0b87828ccabbd5ebfbebe356b4" +uuid = "dce04be8-c92d-5529-be00-80e4d2c0e197" +version = "2.3.0" + +[[deps.ArgTools]] +uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" +version = "1.1.1" + +[[deps.ArrayInterface]] +deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"] +git-tree-sha1 = "d84c956c4c0548b4caf0e4e96cf5b6494b5b1529" +uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" +version = "3.1.32" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" + +[[deps.AxisAlgorithms]] +deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"] +git-tree-sha1 = "66771c8d21c8ff5e3a93379480a2307ac36863f7" +uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950" +version = "1.0.1" + +[[deps.AxisArrays]] +deps = ["Dates", "IntervalSets", "IterTools", "RangeArrays"] +git-tree-sha1 = "1dd4d9f5beebac0c03446918741b1a03dc5e5788" +uuid = "39de3d68-74b9-583c-8d2d-e117c070f3a9" +version = "0.4.6" + +[[deps.BangBang]] +deps = ["Compat", "ConstructionBase", "Future", "InitialValues", "LinearAlgebra", "Requires", "Setfield", "Tables", "ZygoteRules"] +git-tree-sha1 = "7fe6d92c4f281cf4ca6f2fba0ce7b299742da7ca" +uuid = "198e06fe-97b7-11e9-32a5-e1d131e6ad66" +version = "0.3.37" + +[[deps.Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[deps.Baselet]] +git-tree-sha1 = "aebf55e6d7795e02ca500a689d326ac979aaf89e" +uuid = "9718e550-a3fa-408a-8086-8db961cd8217" +version = "0.1.1" + +[[deps.BenchmarkTools]] +deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"] +git-tree-sha1 = "d9a9701b899b30332bbcb3e1679c41cce81fb0e8" +uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +version = "1.3.2" + +[[deps.BitFlags]] +git-tree-sha1 = "43b1a4a8f797c1cddadf60499a8a077d4af2cd2d" +uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35" +version = "0.1.7" + +[[deps.Bzip2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2" +uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" +version = "1.0.8+0" + +[[deps.Cairo_jll]] +deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"] +git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2" +uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a" +version = "1.16.1+1" + +[[deps.CalibrateEmulateSample]] +deps = ["AbstractMCMC", "AdvancedMH", "Conda", "Distributions", "DocStringExtensions", "EnsembleKalmanProcesses", "GaussianProcesses", "LinearAlgebra", "MCMCChains", "Pkg", "Printf", "PyCall", "Random", "RandomFeatures", "ScikitLearn", "StableRNGs", "Statistics", "StatsBase"] +path = "../.." +uuid = "95e48a1f-0bec-4818-9538-3db4340308e3" +version = "0.2.0" + +[[deps.ChainRulesCore]] +deps = ["Compat", "LinearAlgebra", "SparseArrays"] +git-tree-sha1 = "e7ff6cadf743c098e08fca25c91103ee4303c9bb" +uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +version = "1.15.6" + +[[deps.ChangesOfVariables]] +deps = ["ChainRulesCore", "LinearAlgebra", "Test"] +git-tree-sha1 = "38f7a08f19d8810338d4f5085211c7dfa5d5bdd8" +uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" +version = "0.1.4" + +[[deps.CodecBzip2]] +deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"] +git-tree-sha1 = "2e62a725210ce3c3c2e1a3080190e7ca491f18d7" +uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd" +version = "0.7.2" + +[[deps.CodecZlib]] +deps = ["TranscodingStreams", "Zlib_jll"] +git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da" +uuid = "944b1d66-785c-5afd-91f1-9de20f533193" +version = "0.7.0" + +[[deps.ColorSchemes]] +deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Random", "SnoopPrecompile"] +git-tree-sha1 = "aa3edc8f8dea6cbfa176ee12f7c2fc82f0608ed3" +uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4" +version = "3.20.0" + +[[deps.ColorTypes]] +deps = ["FixedPointNumbers", "Random"] +git-tree-sha1 = "eb7f0f8307f71fac7c606984ea5fb2817275d6e4" +uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" +version = "0.11.4" + +[[deps.ColorVectorSpace]] +deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"] +git-tree-sha1 = "600cc5508d66b78aae350f7accdb58763ac18589" +uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4" +version = "0.9.10" + +[[deps.Colors]] +deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] +git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a" +uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" +version = "0.12.10" + +[[deps.CommonSubexpressions]] +deps = ["MacroTools", "Test"] +git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" +uuid = "bbf7d656-a473-5ed7-a52c-81e309532950" +version = "0.3.0" + +[[deps.Compat]] +deps = ["Dates", "LinearAlgebra", "UUIDs"] +git-tree-sha1 = "00a2cccc7f098ff3b66806862d275ca3db9e6e5a" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "4.5.0" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "0.5.2+0" + +[[deps.CompositionsBase]] +git-tree-sha1 = "455419f7e328a1a2493cabc6428d79e951349769" +uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b" +version = "0.1.1" + +[[deps.Conda]] +deps = ["Downloads", "JSON", "VersionParsing"] +git-tree-sha1 = "6e47d11ea2776bc5627421d59cdcc1296c058071" +uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d" +version = "1.7.0" + +[[deps.ConsoleProgressMonitor]] +deps = ["Logging", "ProgressMeter"] +git-tree-sha1 = "3ab7b2136722890b9af903859afcf457fa3059e8" +uuid = "88cd18e8-d9cc-4ea6-8889-5259c0d15c8b" +version = "0.1.2" + +[[deps.ConstructionBase]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "fb21ddd70a051d882a1686a5a550990bbe371a95" +uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" +version = "1.4.1" + +[[deps.Contour]] +git-tree-sha1 = "d05d9e7b7aedff4e5b51a029dced05cfb6125781" +uuid = "d38c429a-6771-53c6-b99e-75d170b6e991" +version = "0.6.2" + +[[deps.Convex]] +deps = ["AbstractTrees", "BenchmarkTools", "LDLFactorizations", "LinearAlgebra", "MathOptInterface", "OrderedCollections", "SparseArrays", "Test"] +git-tree-sha1 = "c90364e06afb0da76e72728fc758b82857ed14e2" +uuid = "f65535da-76fb-5f13-bab9-19810c17039a" +version = "0.15.2" + +[[deps.Crayons]] +git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15" +uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" +version = "4.1.1" + +[[deps.DataAPI]] +git-tree-sha1 = "e8119c1a33d267e16108be441a287a6981ba1630" +uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" +version = "1.14.0" + +[[deps.DataFrames]] +deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Random", "Reexport", "SnoopPrecompile", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"] +git-tree-sha1 = "d4f69885afa5e6149d0cab3818491565cf41446d" +uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +version = "1.4.4" + +[[deps.DataStructures]] +deps = ["Compat", "InteractiveUtils", "OrderedCollections"] +git-tree-sha1 = "d1fff3a548102f48987a52a2e0d114fa97d730f0" +uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +version = "0.18.13" + +[[deps.DataValueInterfaces]] +git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" +uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464" +version = "1.0.0" + +[[deps.Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[deps.DefineSingletons]] +git-tree-sha1 = "0fba8b706d0178b4dc7fd44a96a92382c9065c2c" +uuid = "244e2a9f-e319-4986-a169-4d1fe445cd52" +version = "0.1.2" + +[[deps.DelimitedFiles]] +deps = ["Mmap"] +uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" + +[[deps.DensityInterface]] +deps = ["InverseFunctions", "Test"] +git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b" +uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d" +version = "0.4.0" + +[[deps.DiffResults]] +deps = ["StaticArraysCore"] +git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621" +uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" +version = "1.1.0" + +[[deps.DiffRules]] +deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"] +git-tree-sha1 = "c5b6685d53f933c11404a3ae9822afe30d522494" +uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" +version = "1.12.2" + +[[deps.Distances]] +deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"] +git-tree-sha1 = "3258d0659f812acde79e8a74b11f17ac06d0ca04" +uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" +version = "0.10.7" + +[[deps.Distributed]] +deps = ["Random", "Serialization", "Sockets"] +uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" + +[[deps.Distributions]] +deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"] +git-tree-sha1 = "a7756d098cbabec6b3ac44f369f74915e8cfd70a" +uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" +version = "0.25.79" + +[[deps.DocStringExtensions]] +deps = ["LibGit2"] +git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" +uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +version = "0.9.3" + +[[deps.Downloads]] +deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] +uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +version = "1.6.0" + +[[deps.ElasticArrays]] +deps = ["Adapt"] +git-tree-sha1 = "e1c40d78de68e9a2be565f0202693a158ec9ad85" +uuid = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4" +version = "1.2.11" + +[[deps.ElasticPDMats]] +deps = ["LinearAlgebra", "MacroTools", "PDMats"] +git-tree-sha1 = "5157c93fe9431a041e4cd84265dfce3d53a52323" +uuid = "2904ab23-551e-5aed-883f-487f97af5226" +version = "0.2.2" + +[[deps.EnsembleKalmanProcesses]] +deps = ["Convex", "Distributions", "DocStringExtensions", "LinearAlgebra", "MathOptInterface", "Optim", "QuadGK", "Random", "SCS", "SparseArrays", "Statistics", "StatsBase", "TOML"] +git-tree-sha1 = "b2945201dc88e865f8766bd09344eeeb9358d7de" +uuid = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d" +version = "0.13.1" + +[[deps.Expat_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "bad72f730e9e91c08d9427d5e8db95478a3c323d" +uuid = "2e619515-83b5-522b-bb60-26c02a35a201" +version = "2.4.8+0" + +[[deps.FFMPEG]] +deps = ["FFMPEG_jll"] +git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8" +uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" +version = "0.4.1" + +[[deps.FFMPEG_jll]] +deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Pkg", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"] +git-tree-sha1 = "74faea50c1d007c85837327f6775bea60b5492dd" +uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5" +version = "4.4.2+2" + +[[deps.FFTW]] +deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"] +git-tree-sha1 = "90630efff0894f8142308e334473eba54c433549" +uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +version = "1.5.0" + +[[deps.FFTW_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "c6033cc3892d0ef5bb9cd29b7f2f0331ea5184ea" +uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a" +version = "3.3.10+0" + +[[deps.FastGaussQuadrature]] +deps = ["LinearAlgebra", "SpecialFunctions", "StaticArrays"] +git-tree-sha1 = "58d83dd5a78a36205bdfddb82b1bb67682e64487" +uuid = "442a2c76-b920-505d-bb47-c5924d526838" +version = "0.4.9" + +[[deps.FileIO]] +deps = ["Pkg", "Requires", "UUIDs"] +git-tree-sha1 = "7be5f99f7d15578798f338f5433b6c432ea8037b" +uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +version = "1.16.0" + +[[deps.FileWatching]] +uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" + +[[deps.FillArrays]] +deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"] +git-tree-sha1 = "9a0472ec2f5409db243160a8b030f94c380167a3" +uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" +version = "0.13.6" + +[[deps.FiniteDiff]] +deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"] +git-tree-sha1 = "ec299fdc8f49ae450807b0cb1d161c6b76fd2b60" +uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" +version = "2.10.1" + +[[deps.FixedPointNumbers]] +deps = ["Statistics"] +git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc" +uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" +version = "0.8.4" + +[[deps.Fontconfig_jll]] +deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"] +git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03" +uuid = "a3f928ae-7b40-5064-980b-68af3947d34b" +version = "2.13.93+0" + +[[deps.Formatting]] +deps = ["Printf"] +git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8" +uuid = "59287772-0a20-5a39-b81b-1366585eb4c0" +version = "0.4.2" + +[[deps.ForwardDiff]] +deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"] +git-tree-sha1 = "a69dd6db8a809f78846ff259298678f0d6212180" +uuid = "f6369f11-7733-5829-9624-2563aa707210" +version = "0.10.34" + +[[deps.FreeType2_jll]] +deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] +git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9" +uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7" +version = "2.10.4+0" + +[[deps.FriBidi_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91" +uuid = "559328eb-81f9-559d-9380-de523a88c83c" +version = "1.0.10+0" + +[[deps.Future]] +deps = ["Random"] +uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" + +[[deps.GLFW_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"] +git-tree-sha1 = "d972031d28c8c8d9d7b41a536ad7bb0c2579caca" +uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89" +version = "3.3.8+0" + +[[deps.GR]] +deps = ["Artifacts", "Base64", "DelimitedFiles", "Downloads", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Preferences", "Printf", "Random", "Serialization", "Sockets", "TOML", "Tar", "Test", "UUIDs", "p7zip_jll"] +git-tree-sha1 = "387d2b8b3ca57b791633f0993b31d8cb43ea3292" +uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" +version = "0.71.3" + +[[deps.GR_jll]] +deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"] +git-tree-sha1 = "5982b5e20f97bff955e9a2343a14da96a746cd8c" +uuid = "d2c73de3-f751-5644-a686-071e5b155ba9" +version = "0.71.3+0" + +[[deps.GaussianProcesses]] +deps = ["Distances", "Distributions", "ElasticArrays", "ElasticPDMats", "FastGaussQuadrature", "ForwardDiff", "LinearAlgebra", "Optim", "PDMats", "Printf", "ProgressMeter", "Random", "RecipesBase", "ScikitLearnBase", "SpecialFunctions", "StaticArrays", "Statistics", "StatsFuns"] +git-tree-sha1 = "31749ff6868caf6dd50902eec652a724071dbed3" +uuid = "891a1506-143c-57d2-908e-e1f8e92e6de9" +version = "0.12.5" + +[[deps.Gettext_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"] +git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046" +uuid = "78b55507-aeef-58d4-861c-77aaff3498b1" +version = "0.21.0+0" + +[[deps.Glib_jll]] +deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Pkg", "Zlib_jll"] +git-tree-sha1 = "d3b3624125c1474292d0d8ed0f65554ac37ddb23" +uuid = "7746bdde-850d-59dc-9ae8-88ece973131d" +version = "2.74.0+2" + +[[deps.Graphite2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011" +uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472" +version = "1.3.14+0" + +[[deps.Grisu]] +git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2" +uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe" +version = "1.0.2" + +[[deps.HTTP]] +deps = ["Base64", "CodecZlib", "Dates", "IniFile", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"] +git-tree-sha1 = "fd9861adba6b9ae4b42582032d0936d456c8602d" +uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" +version = "1.6.3" + +[[deps.HarfBuzz_jll]] +deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"] +git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3" +uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566" +version = "2.8.1+1" + +[[deps.IfElse]] +git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1" +uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" +version = "0.1.1" + +[[deps.IniFile]] +git-tree-sha1 = "f550e6e32074c939295eb5ea6de31849ac2c9625" +uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f" +version = "0.5.1" + +[[deps.InitialValues]] +git-tree-sha1 = "4da0f88e9a39111c2fa3add390ab15f3a44f3ca3" +uuid = "22cec73e-a1b8-11e9-2c92-598750a2cf9c" +version = "0.3.1" + +[[deps.IntelOpenMP_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "d979e54b71da82f3a65b62553da4fc3d18c9004c" +uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" +version = "2018.0.3+2" + +[[deps.InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[deps.Interpolations]] +deps = ["Adapt", "AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"] +git-tree-sha1 = "721ec2cf720536ad005cb38f50dbba7b02419a15" +uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +version = "0.14.7" + +[[deps.IntervalSets]] +deps = ["Dates", "Random", "Statistics"] +git-tree-sha1 = "16c0cc91853084cb5f58a78bd209513900206ce6" +uuid = "8197267c-284f-5f27-9208-e0e47529a953" +version = "0.7.4" + +[[deps.InverseFunctions]] +deps = ["Test"] +git-tree-sha1 = "49510dfcb407e572524ba94aeae2fced1f3feb0f" +uuid = "3587e190-3f89-42d0-90ee-14403ec27112" +version = "0.1.8" + +[[deps.InvertedIndices]] +git-tree-sha1 = "82aec7a3dd64f4d9584659dc0b62ef7db2ef3e19" +uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f" +version = "1.2.0" + +[[deps.IrrationalConstants]] +git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151" +uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" +version = "0.1.1" + +[[deps.IterTools]] +git-tree-sha1 = "fa6287a4469f5e048d763df38279ee729fbd44e5" +uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e" +version = "1.4.0" + +[[deps.IteratorInterfaceExtensions]] +git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856" +uuid = "82899510-4779-5014-852e-03e436cf321d" +version = "1.0.0" + +[[deps.JLD2]] +deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "Pkg", "Printf", "Reexport", "TranscodingStreams", "UUIDs"] +git-tree-sha1 = "ec8a9c9f0ecb1c687e34c1fda2699de4d054672a" +uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +version = "0.4.29" + +[[deps.JLFzf]] +deps = ["Pipe", "REPL", "Random", "fzf_jll"] +git-tree-sha1 = "f377670cda23b6b7c1c0b3893e37451c5c1a2185" +uuid = "1019f520-868f-41f5-a6de-eb00f4b6a39c" +version = "0.1.5" + +[[deps.JLLWrappers]] +deps = ["Preferences"] +git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.4.1" + +[[deps.JSON]] +deps = ["Dates", "Mmap", "Parsers", "Unicode"] +git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e" +uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +version = "0.21.3" + +[[deps.JpegTurbo_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "b53380851c6e6664204efb2e62cd24fa5c47e4ba" +uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8" +version = "2.1.2+0" + +[[deps.KernelDensity]] +deps = ["Distributions", "DocStringExtensions", "FFTW", "Interpolations", "StatsBase"] +git-tree-sha1 = "9816b296736292a80b9a3200eb7fbb57aaa3917a" +uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b" +version = "0.6.5" + +[[deps.LAME_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c" +uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d" +version = "3.100.1+0" + +[[deps.LDLFactorizations]] +deps = ["AMD", "LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "743544bcdba7b4ad744bfd5d062c977a9df553a7" +uuid = "40e66cde-538c-5869-a4ad-c39174c6795b" +version = "0.9.0" + +[[deps.LERC_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "bf36f528eec6634efc60d7ec062008f171071434" +uuid = "88015f11-f218-50d7-93a8-a6af411a945d" +version = "3.0.0+1" + +[[deps.LZO_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6" +uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac" +version = "2.10.1+0" + +[[deps.LaTeXStrings]] +git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996" +uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +version = "1.3.0" + +[[deps.Latexify]] +deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "OrderedCollections", "Printf", "Requires"] +git-tree-sha1 = "2422f47b34d4b127720a18f86fa7b1aa2e141f29" +uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" +version = "0.15.18" + +[[deps.LazyArtifacts]] +deps = ["Artifacts", "Pkg"] +uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" + +[[deps.LeftChildRightSiblingTrees]] +deps = ["AbstractTrees"] +git-tree-sha1 = "fb6803dafae4a5d62ea5cab204b1e657d9737e7f" +uuid = "1d6d02ad-be62-4b6b-8a6d-2f90e265016e" +version = "0.2.0" + +[[deps.LibCURL]] +deps = ["LibCURL_jll", "MozillaCACerts_jll"] +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" +version = "0.6.3" + +[[deps.LibCURL_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] +uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" +version = "7.84.0+0" + +[[deps.LibGit2]] +deps = ["Base64", "NetworkOptions", "Printf", "SHA"] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" + +[[deps.LibSSH2_jll]] +deps = ["Artifacts", "Libdl", "MbedTLS_jll"] +uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" +version = "1.10.2+0" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[deps.Libffi_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290" +uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490" +version = "3.2.2+1" + +[[deps.Libgcrypt_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"] +git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae" +uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4" +version = "1.8.7+0" + +[[deps.Libglvnd_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"] +git-tree-sha1 = "6f73d1dd803986947b2c750138528a999a6c7733" +uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29" +version = "1.6.0+0" + +[[deps.Libgpg_error_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9" +uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8" +version = "1.42.0+0" + +[[deps.Libiconv_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "c7cb1f5d892775ba13767a87c7ada0b980ea0a71" +uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" +version = "1.16.1+2" + +[[deps.Libmount_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73" +uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9" +version = "2.35.0+0" + +[[deps.Libtiff_jll]] +deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"] +git-tree-sha1 = "3eb79b0ca5764d4799c06699573fd8f533259713" +uuid = "89763e89-9b03-5906-acba-b20f662cd828" +version = "4.4.0+0" + +[[deps.Libuuid_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066" +uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700" +version = "2.36.0+0" + +[[deps.LineSearches]] +deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"] +git-tree-sha1 = "7bbea35cec17305fc70a0e5b4641477dc0789d9d" +uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" +version = "7.2.0" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[deps.LogDensityProblems]] +deps = ["ArgCheck", "DocStringExtensions", "Random", "UnPack"] +git-tree-sha1 = "05fdf369bd52030212a7c730645478dc5bcd67b6" +uuid = "6fdf6af0-433a-55f7-b3ed-c6c6e0b8df7c" +version = "2.1.0" + +[[deps.LogExpFunctions]] +deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"] +git-tree-sha1 = "946607f84feb96220f480e0422d3484c49c00239" +uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +version = "0.3.19" + +[[deps.Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[deps.LoggingExtras]] +deps = ["Dates", "Logging"] +git-tree-sha1 = "5d4d2d9904227b8bd66386c1138cf4d5ffa826bf" +uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36" +version = "0.4.9" + +[[deps.MCMCChains]] +deps = ["AbstractMCMC", "AxisArrays", "Dates", "Distributions", "Formatting", "IteratorInterfaceExtensions", "KernelDensity", "LinearAlgebra", "MCMCDiagnosticTools", "MLJModelInterface", "NaturalSort", "OrderedCollections", "PrettyTables", "Random", "RecipesBase", "Serialization", "Statistics", "StatsBase", "StatsFuns", "TableTraits", "Tables"] +git-tree-sha1 = "8577f598f0991e7bdf9472a8a59938b1e253b493" +uuid = "c7f686f2-ff18-58e9-bc7b-31028e88f75d" +version = "5.6.1" + +[[deps.MCMCDiagnosticTools]] +deps = ["AbstractFFTs", "DataAPI", "DataStructures", "Distributions", "LinearAlgebra", "MLJModelInterface", "Random", "SpecialFunctions", "Statistics", "StatsBase", "Tables"] +git-tree-sha1 = "d1737c39191aa26f42a64e320de313f1d1fd74b1" +uuid = "be115224-59cd-429b-ad48-344e309966f0" +version = "0.2.1" + +[[deps.MKL_jll]] +deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"] +git-tree-sha1 = "2ce8695e1e699b68702c03402672a69f54b8aca9" +uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" +version = "2022.2.0+0" + +[[deps.MLJModelInterface]] +deps = ["Random", "ScientificTypesBase", "StatisticalTraits"] +git-tree-sha1 = "c8b7e632d6754a5e36c0d94a4b466a5ba3a30128" +uuid = "e80e1ace-859a-464e-9ed9-23947d8ae3ea" +version = "1.8.0" + +[[deps.MacroTools]] +deps = ["Markdown", "Random"] +git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2" +uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +version = "0.5.10" + +[[deps.Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[deps.MathOptInterface]] +deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "Printf", "SparseArrays", "SpecialFunctions", "Test", "Unicode"] +git-tree-sha1 = "f21d9bb78cea58d9b535b9d9f2a1e601112ef0ed" +uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +version = "1.11.2" + +[[deps.MbedTLS]] +deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "Random", "Sockets"] +git-tree-sha1 = "03a9b9718f5682ecb107ac9f7308991db4ce395b" +uuid = "739be429-bea8-5141-9913-cc70e7f3736d" +version = "1.1.7" + +[[deps.MbedTLS_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" +version = "2.28.0+0" + +[[deps.Measures]] +git-tree-sha1 = "c13304c81eec1ed3af7fc20e75fb6b26092a1102" +uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e" +version = "0.3.2" + +[[deps.MicroCollections]] +deps = ["BangBang", "InitialValues", "Setfield"] +git-tree-sha1 = "4d5917a26ca33c66c8e5ca3247bd163624d35493" +uuid = "128add7d-3638-4c79-886c-908ea0c25c34" +version = "0.1.3" + +[[deps.Missings]] +deps = ["DataAPI"] +git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272" +uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" +version = "1.1.0" + +[[deps.Mmap]] +uuid = "a63ad114-7e13-5084-954f-fe012c677804" + +[[deps.MozillaCACerts_jll]] +uuid = "14a3606d-f60d-562e-9121-12d972cd8159" +version = "2022.2.1" + +[[deps.MutableArithmetics]] +deps = ["LinearAlgebra", "SparseArrays", "Test"] +git-tree-sha1 = "aa532179d4a643d4bd9f328589ca01fa20a0d197" +uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +version = "1.1.0" + +[[deps.NLSolversBase]] +deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"] +git-tree-sha1 = "a0b464d183da839699f4c79e7606d9d186ec172c" +uuid = "d41bc354-129a-5804-8e4c-c37616107c6c" +version = "7.8.3" + +[[deps.NaNMath]] +deps = ["OpenLibm_jll"] +git-tree-sha1 = "a7c3d1da1189a1c2fe843a3bfa04d18d20eb3211" +uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +version = "1.0.1" + +[[deps.NaturalSort]] +git-tree-sha1 = "eda490d06b9f7c00752ee81cfa451efe55521e21" +uuid = "c020b1a1-e9b0-503a-9c33-f039bfc54a85" +version = "1.0.0" + +[[deps.NetworkOptions]] +uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" +version = "1.2.0" + +[[deps.OffsetArrays]] +deps = ["Adapt"] +git-tree-sha1 = "f71d8950b724e9ff6110fc948dff5a329f901d64" +uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +version = "1.12.8" + +[[deps.Ogg_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "887579a3eb005446d514ab7aeac5d1d027658b8f" +uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051" +version = "1.3.5+1" + +[[deps.OpenBLAS32_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "9c6c2ed4b7acd2137b878eb96c68e63b76199d0f" +uuid = "656ef2d0-ae68-5445-9ca0-591084a874a2" +version = "0.3.17+0" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.20+0" + +[[deps.OpenLibm_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "05823500-19ac-5b8b-9628-191a04bc5112" +version = "0.8.1+0" + +[[deps.OpenSSL]] +deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "OpenSSL_jll", "Sockets"] +git-tree-sha1 = "df6830e37943c7aaa10023471ca47fb3065cc3c4" +uuid = "4d8831e6-92b7-49fb-bdf8-b643e874388c" +version = "1.3.2" + +[[deps.OpenSSL_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "f6e9dba33f9f2c44e08a020b0caf6903be540004" +uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" +version = "1.1.19+0" + +[[deps.OpenSpecFun_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" +uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" +version = "0.5.5+0" + +[[deps.Optim]] +deps = ["Compat", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"] +git-tree-sha1 = "1903afc76b7d01719d9c30d3c7d501b61db96721" +uuid = "429524aa-4258-5aef-a3af-852621145aeb" +version = "1.7.4" + +[[deps.Opus_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720" +uuid = "91d4177d-7536-5919-b921-800302f37372" +version = "1.3.2+0" + +[[deps.OrderedCollections]] +git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c" +uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +version = "1.4.1" + +[[deps.PCRE2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15" +version = "10.40.0+0" + +[[deps.PDMats]] +deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse", "Test"] +git-tree-sha1 = "95a4038d1011dfdbde7cecd2ad0ac411e53ab1bc" +uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" +version = "0.10.1" + +[[deps.Parameters]] +deps = ["OrderedCollections", "UnPack"] +git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe" +uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a" +version = "0.12.3" + +[[deps.Parsers]] +deps = ["Dates", "SnoopPrecompile"] +git-tree-sha1 = "6466e524967496866901a78fca3f2e9ea445a559" +uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" +version = "2.5.2" + +[[deps.Pipe]] +git-tree-sha1 = "6842804e7867b115ca9de748a0cf6b364523c16d" +uuid = "b98c9c47-44ae-5843-9183-064241ee97a0" +version = "1.3.0" + +[[deps.Pixman_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29" +uuid = "30392449-352a-5448-841d-b1acce4e97dc" +version = "0.40.1+0" + +[[deps.Pkg]] +deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] +uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +version = "1.8.0" + +[[deps.PlotThemes]] +deps = ["PlotUtils", "Statistics"] +git-tree-sha1 = "1f03a2d339f42dca4a4da149c7e15e9b896ad899" +uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a" +version = "3.1.0" + +[[deps.PlotUtils]] +deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "SnoopPrecompile", "Statistics"] +git-tree-sha1 = "5b7690dd212e026bbab1860016a6601cb077ab66" +uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043" +version = "1.3.2" + +[[deps.Plots]] +deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "JLFzf", "JSON", "LaTeXStrings", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "Preferences", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "RelocatableFolders", "Requires", "Scratch", "Showoff", "SnoopPrecompile", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "Unzip"] +git-tree-sha1 = "02ecc6a3427e7edfff1cebcf66c1f93dd77760ec" +uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +version = "1.38.1" + +[[deps.PooledArrays]] +deps = ["DataAPI", "Future"] +git-tree-sha1 = "a6062fe4063cdafe78f4a0a81cfffb89721b30e7" +uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720" +version = "1.4.2" + +[[deps.PositiveFactorizations]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "17275485f373e6673f7e7f97051f703ed5b15b20" +uuid = "85a6dd25-e78a-55b7-8502-1745935b8125" +version = "0.2.4" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.3.0" + +[[deps.PrettyTables]] +deps = ["Crayons", "Formatting", "LaTeXStrings", "Markdown", "Reexport", "StringManipulation", "Tables"] +git-tree-sha1 = "96f6db03ab535bdb901300f88335257b0018689d" +uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" +version = "2.2.2" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[deps.Profile]] +deps = ["Printf"] +uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79" + +[[deps.ProgressLogging]] +deps = ["Logging", "SHA", "UUIDs"] +git-tree-sha1 = "80d919dee55b9c50e8d9e2da5eeafff3fe58b539" +uuid = "33c8b6b6-d38a-422a-b730-caa89a2f386c" +version = "0.1.4" + +[[deps.ProgressMeter]] +deps = ["Distributed", "Printf"] +git-tree-sha1 = "d7a7aef8f8f2d537104f170139553b14dfe39fe9" +uuid = "92933f4c-e287-5a05-a399-4b506db050ca" +version = "1.7.2" + +[[deps.PyCall]] +deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"] +git-tree-sha1 = "b32c4b415f41f10c671cba02ae3275027dea8892" +uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" +version = "1.95.0" + +[[deps.PyPlot]] +deps = ["Colors", "LaTeXStrings", "PyCall", "Sockets", "Test", "VersionParsing"] +git-tree-sha1 = "f9d953684d4d21e947cb6d642db18853d43cb027" +uuid = "d330b81b-6aea-500a-939a-2ce795aea3ee" +version = "2.11.0" + +[[deps.Qt5Base_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"] +git-tree-sha1 = "0c03844e2231e12fda4d0086fd7cbe4098ee8dc5" +uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1" +version = "5.15.3+2" + +[[deps.QuadGK]] +deps = ["DataStructures", "LinearAlgebra"] +git-tree-sha1 = "97aa253e65b784fd13e83774cadc95b38011d734" +uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +version = "2.6.0" + +[[deps.REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + +[[deps.Random]] +deps = ["SHA", "Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[deps.RandomFeatures]] +deps = ["Distributions", "DocStringExtensions", "EnsembleKalmanProcesses", "LinearAlgebra", "Random", "SpecialFunctions", "Statistics", "StatsBase", "Tullio"] +git-tree-sha1 = "73119392cba800e1f433dd056fcb2d948932c7ab" +uuid = "36c3bae2-c0c3-419d-b3b4-eebadd35c5e5" +version = "0.2.1" + +[[deps.RangeArrays]] +git-tree-sha1 = "b9039e93773ddcfc828f12aadf7115b4b4d225f5" +uuid = "b3c3ace0-ae52-54e7-9d0b-2c1406fd6b9d" +version = "0.3.2" + +[[deps.Ratios]] +deps = ["Requires"] +git-tree-sha1 = "dc84268fe0e3335a62e315a3a7cf2afa7178a734" +uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439" +version = "0.4.3" + +[[deps.RecipesBase]] +deps = ["SnoopPrecompile"] +git-tree-sha1 = "261dddd3b862bd2c940cf6ca4d1c8fe593e457c8" +uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" +version = "1.3.3" + +[[deps.RecipesPipeline]] +deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase", "SnoopPrecompile"] +git-tree-sha1 = "e974477be88cb5e3040009f3767611bc6357846f" +uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c" +version = "0.6.11" + +[[deps.Reexport]] +git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "1.2.2" + +[[deps.RelocatableFolders]] +deps = ["SHA", "Scratch"] +git-tree-sha1 = "90bc7a7c96410424509e4263e277e43250c05691" +uuid = "05181044-ff0b-4ac5-8273-598c1e38db00" +version = "1.0.0" + +[[deps.Requires]] +deps = ["UUIDs"] +git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" +uuid = "ae029012-a4dd-5104-9daa-d747884805df" +version = "1.3.0" + +[[deps.Rmath]] +deps = ["Random", "Rmath_jll"] +git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f" +uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa" +version = "0.7.0" + +[[deps.Rmath_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7" +uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f" +version = "0.3.0+0" + +[[deps.SCS]] +deps = ["MathOptInterface", "Requires", "SCS_GPU_jll", "SCS_MKL_jll", "SCS_jll", "SparseArrays"] +git-tree-sha1 = "2e3ca40559ecaed6ffe9410b06aabcc1e087215d" +uuid = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" +version = "1.1.3" + +[[deps.SCS_GPU_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "OpenBLAS32_jll", "Pkg"] +git-tree-sha1 = "2b3799ff650d0530a19c2a3bd4b158a4f3e4581a" +uuid = "af6e375f-46ec-5fa0-b791-491b0dfa44a4" +version = "3.2.1+0" + +[[deps.SCS_MKL_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "MKL_jll", "Pkg"] +git-tree-sha1 = "a4923177e60fdb7f802e1a42a73d0af400eea163" +uuid = "3f2553a9-4106-52be-b7dd-865123654657" +version = "3.2.2+0" + +[[deps.SCS_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "OpenBLAS32_jll", "Pkg"] +git-tree-sha1 = "5544538910047c7522908cf87bb0c884a7afff92" +uuid = "f4f2fc5b-1d94-523c-97ea-2ab488bedf4b" +version = "3.2.1+0" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.ScientificTypesBase]] +git-tree-sha1 = "a8e18eb383b5ecf1b5e6fc237eb39255044fd92b" +uuid = "30f210dd-8aff-4c5f-94ba-8e64358c1161" +version = "3.0.0" + +[[deps.ScikitLearn]] +deps = ["Compat", "Conda", "DataFrames", "Distributed", "IterTools", "LinearAlgebra", "MacroTools", "Parameters", "Printf", "PyCall", "Random", "ScikitLearnBase", "SparseArrays", "StatsBase", "VersionParsing"] +git-tree-sha1 = "32c302f5e6bd7a0c6a44668bfe5133c2ff076191" +uuid = "3646fa90-6ef7-5e7e-9f22-8aca16db6324" +version = "0.6.6" + +[[deps.ScikitLearnBase]] +deps = ["LinearAlgebra", "Random", "Statistics"] +git-tree-sha1 = "7877e55c1523a4b336b433da39c8e8c08d2f221f" +uuid = "6e75b9c4-186b-50bd-896f-2d2496a4843e" +version = "0.5.0" + +[[deps.Scratch]] +deps = ["Dates"] +git-tree-sha1 = "f94f779c94e58bf9ea243e77a37e16d9de9126bd" +uuid = "6c6a2e73-6563-6170-7368-637461726353" +version = "1.1.1" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[deps.Setfield]] +deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"] +git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac" +uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46" +version = "1.1.1" + +[[deps.SharedArrays]] +deps = ["Distributed", "Mmap", "Random", "Serialization"] +uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" + +[[deps.Showoff]] +deps = ["Dates", "Grisu"] +git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de" +uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f" +version = "1.0.3" + +[[deps.SimpleBufferStream]] +git-tree-sha1 = "874e8867b33a00e784c8a7e4b60afe9e037b74e1" +uuid = "777ac1f9-54b0-4bf8-805c-2214025038e7" +version = "1.1.0" + +[[deps.SnoopPrecompile]] +git-tree-sha1 = "f604441450a3c0569830946e5b33b78c928e1a85" +uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c" +version = "1.0.1" + +[[deps.Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" + +[[deps.SortingAlgorithms]] +deps = ["DataStructures"] +git-tree-sha1 = "a4ada03f999bd01b3a25dcaa30b2d929fe537e00" +uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" +version = "1.1.0" + +[[deps.SparseArrays]] +deps = ["LinearAlgebra", "Random"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[deps.SpecialFunctions]] +deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] +git-tree-sha1 = "d75bda01f8c31ebb72df80a46c88b25d1c79c56d" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "2.1.7" + +[[deps.SplittablesBase]] +deps = ["Setfield", "Test"] +git-tree-sha1 = "e08a62abc517eb79667d0a29dc08a3b589516bb5" +uuid = "171d559e-b47b-412a-8079-5efa626c420e" +version = "0.1.15" + +[[deps.StableRNGs]] +deps = ["Random", "Test"] +git-tree-sha1 = "3be7d49667040add7ee151fefaf1f8c04c8c8276" +uuid = "860ef19b-820b-49d6-a774-d7a799459cd3" +version = "1.0.0" + +[[deps.Static]] +deps = ["IfElse"] +git-tree-sha1 = "a8f30abc7c64a39d389680b74e749cf33f872a70" +uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" +version = "0.3.3" + +[[deps.StaticArrays]] +deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"] +git-tree-sha1 = "6954a456979f23d05085727adb17c4551c19ecd1" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "1.5.12" + +[[deps.StaticArraysCore]] +git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a" +uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +version = "1.4.0" + +[[deps.StatisticalTraits]] +deps = ["ScientificTypesBase"] +git-tree-sha1 = "30b9236691858e13f167ce829490a68e1a597782" +uuid = "64bff920-2084-43da-a3e6-9bb72801c0c9" +version = "3.2.0" + +[[deps.Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[[deps.StatsAPI]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "f9af7f195fb13589dd2e2d57fdb401717d2eb1f6" +uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" +version = "1.5.0" + +[[deps.StatsBase]] +deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] +git-tree-sha1 = "d1bf48bfcc554a3761a133fe3a9bb01488e06916" +uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +version = "0.33.21" + +[[deps.StatsFuns]] +deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"] +git-tree-sha1 = "5950925ff997ed6fb3e985dcce8eb1ba42a0bbe7" +uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" +version = "0.9.18" + +[[deps.StringManipulation]] +git-tree-sha1 = "46da2434b41f41ac3594ee9816ce5541c6096123" +uuid = "892a3eda-7b42-436c-8928-eab12a02cf0e" +version = "0.3.0" + +[[deps.SuiteSparse]] +deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] +uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.0" + +[[deps.TableTraits]] +deps = ["IteratorInterfaceExtensions"] +git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39" +uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c" +version = "1.0.1" + +[[deps.Tables]] +deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"] +git-tree-sha1 = "c79322d36826aa2f4fd8ecfa96ddb47b174ac78d" +uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +version = "1.10.0" + +[[deps.Tar]] +deps = ["ArgTools", "SHA"] +uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" +version = "1.10.0" + +[[deps.TensorCore]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6" +uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50" +version = "0.1.1" + +[[deps.TerminalLoggers]] +deps = ["LeftChildRightSiblingTrees", "Logging", "Markdown", "Printf", "ProgressLogging", "UUIDs"] +git-tree-sha1 = "f53e34e784ae771eb9ccde4d72e578aa453d0554" +uuid = "5d786b92-1e48-4d6f-9151-6b4477ca9bed" +version = "0.1.6" + +[[deps.Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[deps.TranscodingStreams]] +deps = ["Random", "Test"] +git-tree-sha1 = "94f38103c984f89cf77c402f2a68dbd870f8165f" +uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" +version = "0.9.11" + +[[deps.Transducers]] +deps = ["Adapt", "ArgCheck", "BangBang", "Baselet", "CompositionsBase", "DefineSingletons", "Distributed", "InitialValues", "Logging", "Markdown", "MicroCollections", "Requires", "Setfield", "SplittablesBase", "Tables"] +git-tree-sha1 = "c42fa452a60f022e9e087823b47e5a5f8adc53d5" +uuid = "28d57a85-8fef-5791-bfe6-a80928e7c999" +version = "0.4.75" + +[[deps.Tullio]] +deps = ["ChainRulesCore", "DiffRules", "LinearAlgebra", "Requires"] +git-tree-sha1 = "7871a39eac745697ee512a87eeff06a048a7905b" +uuid = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" +version = "0.3.5" + +[[deps.URIs]] +git-tree-sha1 = "ac00576f90d8a259f2c9d823e91d1de3fd44d348" +uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4" +version = "1.4.1" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[[deps.UnPack]] +git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b" +uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" +version = "1.0.2" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" + +[[deps.UnicodeFun]] +deps = ["REPL"] +git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf" +uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1" +version = "0.4.1" + +[[deps.Unzip]] +git-tree-sha1 = "ca0969166a028236229f63514992fc073799bb78" +uuid = "41fe7b60-77ed-43a1-b4f0-825fd5a5650d" +version = "0.2.0" + +[[deps.VersionParsing]] +git-tree-sha1 = "58d6e80b4ee071f5efd07fda82cb9fbe17200868" +uuid = "81def892-9a0e-5fdd-b105-ffc91e053289" +version = "1.3.0" + +[[deps.Wayland_jll]] +deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"] +git-tree-sha1 = "ed8d92d9774b077c53e1da50fd81a36af3744c1c" +uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89" +version = "1.21.0+0" + +[[deps.Wayland_protocols_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "4528479aa01ee1b3b4cd0e6faef0e04cf16466da" +uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91" +version = "1.25.0+0" + +[[deps.WoodburyMatrices]] +deps = ["LinearAlgebra", "SparseArrays"] +git-tree-sha1 = "de67fa59e33ad156a590055375a30b23c40299d3" +uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6" +version = "0.5.5" + +[[deps.XML2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"] +git-tree-sha1 = "93c41695bc1c08c46c5899f4fe06d6ead504bb73" +uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" +version = "2.10.3+0" + +[[deps.XSLT_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"] +git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a" +uuid = "aed1982a-8fda-507f-9586-7b0439959a61" +version = "1.1.34+0" + +[[deps.Xorg_libX11_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"] +git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527" +uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc" +version = "1.6.9+4" + +[[deps.Xorg_libXau_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e" +uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec" +version = "1.0.9+4" + +[[deps.Xorg_libXcursor_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"] +git-tree-sha1 = "12e0eb3bc634fa2080c1c37fccf56f7c22989afd" +uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724" +version = "1.2.0+4" + +[[deps.Xorg_libXdmcp_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4" +uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05" +version = "1.1.3+4" + +[[deps.Xorg_libXext_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"] +git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3" +uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3" +version = "1.3.4+4" + +[[deps.Xorg_libXfixes_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"] +git-tree-sha1 = "0e0dc7431e7a0587559f9294aeec269471c991a4" +uuid = "d091e8ba-531a-589c-9de9-94069b037ed8" +version = "5.0.3+4" + +[[deps.Xorg_libXi_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"] +git-tree-sha1 = "89b52bc2160aadc84d707093930ef0bffa641246" +uuid = "a51aa0fd-4e3c-5386-b890-e753decda492" +version = "1.7.10+4" + +[[deps.Xorg_libXinerama_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll"] +git-tree-sha1 = "26be8b1c342929259317d8b9f7b53bf2bb73b123" +uuid = "d1454406-59df-5ea1-beac-c340f2130bc3" +version = "1.1.4+4" + +[[deps.Xorg_libXrandr_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll"] +git-tree-sha1 = "34cea83cb726fb58f325887bf0612c6b3fb17631" +uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484" +version = "1.5.2+4" + +[[deps.Xorg_libXrender_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"] +git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96" +uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa" +version = "0.9.10+4" + +[[deps.Xorg_libpthread_stubs_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb" +uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74" +version = "0.1.0+3" + +[[deps.Xorg_libxcb_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"] +git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6" +uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b" +version = "1.13.0+3" + +[[deps.Xorg_libxkbfile_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"] +git-tree-sha1 = "926af861744212db0eb001d9e40b5d16292080b2" +uuid = "cc61e674-0454-545c-8b26-ed2c68acab7a" +version = "1.1.0+4" + +[[deps.Xorg_xcb_util_image_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"] +git-tree-sha1 = "0fab0a40349ba1cba2c1da699243396ff8e94b97" +uuid = "12413925-8142-5f55-bb0e-6d7ca50bb09b" +version = "0.4.0+1" + +[[deps.Xorg_xcb_util_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll"] +git-tree-sha1 = "e7fd7b2881fa2eaa72717420894d3938177862d1" +uuid = "2def613f-5ad1-5310-b15b-b15d46f528f5" +version = "0.4.0+1" + +[[deps.Xorg_xcb_util_keysyms_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"] +git-tree-sha1 = "d1151e2c45a544f32441a567d1690e701ec89b00" +uuid = "975044d2-76e6-5fbe-bf08-97ce7c6574c7" +version = "0.4.0+1" + +[[deps.Xorg_xcb_util_renderutil_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"] +git-tree-sha1 = "dfd7a8f38d4613b6a575253b3174dd991ca6183e" +uuid = "0d47668e-0667-5a69-a72c-f761630bfb7e" +version = "0.3.9+1" + +[[deps.Xorg_xcb_util_wm_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"] +git-tree-sha1 = "e78d10aab01a4a154142c5006ed44fd9e8e31b67" +uuid = "c22f9ab0-d5fe-5066-847c-f4bb1cd4e361" +version = "0.4.1+1" + +[[deps.Xorg_xkbcomp_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxkbfile_jll"] +git-tree-sha1 = "4bcbf660f6c2e714f87e960a171b119d06ee163b" +uuid = "35661453-b289-5fab-8a00-3d9160c6a3a4" +version = "1.4.2+4" + +[[deps.Xorg_xkeyboard_config_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xkbcomp_jll"] +git-tree-sha1 = "5c8424f8a67c3f2209646d4425f3d415fee5931d" +uuid = "33bec58e-1273-512f-9401-5d533626f822" +version = "2.27.0+4" + +[[deps.Xorg_xtrans_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845" +uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10" +version = "1.4.0+3" + +[[deps.Zlib_jll]] +deps = ["Libdl"] +uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +version = "1.2.12+3" + +[[deps.Zstd_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "e45044cd873ded54b6a5bac0eb5c971392cf1927" +uuid = "3161d3a3-bdf6-5164-811a-617609db77b4" +version = "1.5.2+0" + +[[deps.ZygoteRules]] +deps = ["MacroTools"] +git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0" +uuid = "700de1a5-db45-46bc-99cf-38207098b444" +version = "0.2.2" + +[[deps.fzf_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "868e669ccb12ba16eaf50cb2957ee2ff61261c56" +uuid = "214eeab7-80f7-51ab-84ad-2988db7cef09" +version = "0.29.0+0" + +[[deps.libaom_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "3a2ea60308f0996d26f1e5354e10c24e9ef905d4" +uuid = "a4ae2306-e953-59d6-aa16-d00cac43593b" +version = "3.4.0+0" + +[[deps.libass_jll]] +deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] +git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47" +uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0" +version = "0.15.1+0" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl", "OpenBLAS_jll"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.1.1+0" + +[[deps.libfdk_aac_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55" +uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280" +version = "2.0.2+0" + +[[deps.libpng_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] +git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c" +uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f" +version = "1.6.38+0" + +[[deps.libvorbis_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"] +git-tree-sha1 = "b910cb81ef3fe6e78bf6acee440bda86fd6ae00c" +uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a" +version = "1.3.7+1" + +[[deps.nghttp2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" +version = "1.48.0+0" + +[[deps.p7zip_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" +version = "17.4.0+0" + +[[deps.x264_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2" +uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a" +version = "2021.5.5+0" + +[[deps.x265_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9" +uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76" +version = "3.5.0+0" + +[[deps.xkbcommon_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"] +git-tree-sha1 = "9ebfc140cc56e8c2156a15ceac2f0302e327ac0a" +uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd" +version = "1.4.1+0" diff --git a/examples/GCM/Project.toml b/examples/GCM/Project.toml new file mode 100644 index 000000000..7ae63f5a8 --- /dev/null +++ b/examples/GCM/Project.toml @@ -0,0 +1,17 @@ +[deps] +CalibrateEmulateSample = "95e48a1f-0bec-4818-9538-3db4340308e3" +Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +EnsembleKalmanProcesses = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[compat] +FiniteDiff = "~2.10" +julia = "~1.6" diff --git a/examples/GCM/data_from_eki_inflateyonly_100.jld2 b/examples/GCM/data_from_eki_inflateyonly_100.jld2 new file mode 100644 index 000000000..1eae5106d Binary files /dev/null and b/examples/GCM/data_from_eki_inflateyonly_100.jld2 differ diff --git a/examples/GCM/emulate_sample_script.jl b/examples/GCM/emulate_sample_script.jl new file mode 100644 index 000000000..b6cad6658 --- /dev/null +++ b/examples/GCM/emulate_sample_script.jl @@ -0,0 +1,218 @@ +# Script to run Emulation and Sampling on data from GCM + +# Import modules +using Distributions # probability distributions and associated functions +using LinearAlgebra +using Plots +using Random +using JLD2 +# CES +using CalibrateEmulateSample.Emulators +using CalibrateEmulateSample.MarkovChainMonteCarlo +using CalibrateEmulateSample.ParameterDistributions +using CalibrateEmulateSample.DataContainers + +@time begin + + rng_seed = 2413798 + Random.seed!(rng_seed) + + expname = "vrf-svd-nondiag-logdet-1-32" + # expname="vrf-nondiag-tik-33-64" + # expname="vrf-nondiag-tik-65-96" + #emulator_type ∈ ["GPR","ScalarRFR","VectorRFR-svd-diag","VectorRFR-svd-nondiag", "VectorRFR-nondiag"] + emulator_type = "VectorRFR-svd-nondiag" + + # Output figure save directory + example_directory = @__DIR__ + println(example_directory) + figure_save_directory = joinpath(example_directory, "output") + data_save_directory = joinpath(example_directory, "output") + if !isdir(figure_save_directory) + mkdir(figure_save_directory) + end + if !isdir(data_save_directory) + mkdir(data_save_directory) + end + + # Load data from file + datafile = "data_from_eki_inflateyonly_100.jld2" + inputs = load(datafile)["inputs"] #100 x 2 x 10 + outputs = load(datafile)["outputs"] #100 x 96 x 10 + truth = load(datafile)["truth"] # 96 + obs_noise_cov = load(datafile)["obs_noise_cov"] # 96 x 96 + + #take only first 400 points + iter_mask = 1:4 + data_mask = 1:32 + # data_mask= 33:64 + # data_mask= 65:96 + data_mask = 1:96 + + inputs = inputs[:, :, iter_mask] + outputs = outputs[:, data_mask, iter_mask] + obs_noise_cov = obs_noise_cov[data_mask, data_mask] + truth = truth[data_mask] + + # priorfile = "priors.jld2" + # prior = load(priorfile)["prior"] + + # derived quantities + N_ens, input_dim, N_iter = size(inputs) + output_dim = size(outputs, 2) + + stacked_inputs = reshape(permutedims(inputs, (1, 3, 2)), (N_ens * N_iter, input_dim)) + stacked_outputs = reshape(permutedims(outputs, (1, 3, 2)), (N_ens * N_iter, output_dim)) + input_output_pairs = PairedDataContainer(stacked_inputs, stacked_outputs, data_are_columns = false) #data are rows + normalized = true + + # setup random features + + + if emulator_type == "VectorRFR-svd-nondiag" || emulator_type == "VectorRFR-nondiag" + if emulator_type == "VectorRFR-svd-nondiag" + println("Running Vector RF model - using SVD and assuming non-diagonal variance ") + elseif emulator_type == "VectorRFR-nondiag" + println("Running Vector RF model - without SVD and assuming non-diagonal variance ") + end + + n_features = 40 * Int(floor(5 * sqrt(N_ens * N_iter))) + batch_sizes = Dict("train" => 500, "test" => 500, "feature" => 500) + println("build RF with $(N_ens*N_iter) training points and $(n_features) random features.") + eki_options_override = Dict("tikhonov" => 0) + + + mlt = VectorRandomFeatureInterface( + n_features, + input_dim, + output_dim, + batch_sizes = batch_sizes, + optimizer_options = eki_options_override, + ) + + elseif emulator_type == "VectorRFR-svd-diag" + + println("Running Vector RF model - using SVD and assuming diagonal variance") + n_features = 20 * Int(floor(5 * sqrt(N_ens * N_iter))) + batch_sizes = Dict("train" => 500, "test" => 500, "feature" => 500) + println("build RF with $(N_ens*N_iter) training points and $(n_features) random features.") + + mlt = VectorRandomFeatureInterface( + n_features, + input_dim, + output_dim, + batch_sizes = batch_sizes, + diagonalize_output = true, + ) + + elseif emulator_type == "ScalarRFR" + println("Running Scalar RF model") + n_features = 5 * Int(floor(5 * sqrt(N_ens * N_iter))) + batch_sizes = Dict("train" => 500, "test" => 500, "feature" => 500) + mlt = ScalarRandomFeatureInterface(n_features, input_dim, batch_sizes = batch_sizes) + + else + emulator_type == "GPR" + println("Running Gaussian Process model") + gppackage = SKLJL() + mlt = GaussianProcess(gppackage, noise_learn = false) + + end + + if emulator_type == "VectorRFR-nondiag" + emulator = Emulator( + mlt, + input_output_pairs; + obs_noise_cov = obs_noise_cov, + normalize_inputs = normalized, + decorrelate = false, + ) + else + emulator = Emulator(mlt, input_output_pairs; obs_noise_cov = obs_noise_cov, normalize_inputs = normalized) + + end + optimize_hyperparameters!(emulator) + # + # predict at some validation points + # + validate_id = ["phys", "mean", "rand"] + + for vid in validate_id + if vid == "phys" + new_input = [log(0.7 / 0.3) log(7200)]' # physical parameter value (at truth) + elseif vid == "mean" + new_input = [log(0.5 / (1 - 0.5)) log(43200)]' # mean-of-prior parameter value ("near-ish" truth) + elseif vid == "rand" + new_input = [log(0.31735951644387783 / (1 - 0.31735951644387783)) log(90632.50269636544)]' # random parameter value ("far" from truth + end + + pred_mean, pred_cov = predict(emulator, new_input, transform_to_real = true) + pred_sd = sqrt.([max(10 * eps(), pred_cov[1][i, i]) for i in 1:size(pred_cov[1], 1)]) + + + # NB pred_cov is a vector of matrices + tobj = load("truthobj_" * vid * "param.jld2")["truthobj"] + t_mean = tobj.mean[data_mask] + t_cov = tobj.cov[data_mask, data_mask] + + println("prediction error at truth for " * vid * " case:") + println(" mean: ", norm(t_mean - pred_mean)) + println(" cov: ", norm(t_cov - pred_cov[1])) + + save( + joinpath(data_save_directory, vid * "_" * expname * "_results.jld2"), + "pred_mean", + pred_mean, + "pred_cov", + pred_cov, + "pred_sd", + pred_sd, + ) + save(joinpath(data_save_directory, vid * "_" * expname * "_truth.jld2"), "true_mean", t_mean, "true_cov", t_cov) + end + + plot_input = [log(0.7 / 0.3) log(7200)]' # physical parameter value (at truth) + plot_mean, plot_cov = predict(emulator, plot_input, transform_to_real = true) + plot_sd = sqrt.([max(10 * eps(), plot_cov[1][i, i]) for i in 1:size(plot_cov[1], 1)]) + + + ids = [1:32, 33:64, 65:96] + plotnames = ["rh", "pr", "ext"] + + for (id, pn) in zip(ids, plotnames) + if data_mask == 1:96 + plt = plot( + collect(id), + plot_mean[id], + show = true, + ribbon = [2 * plot_sd[id]; 2 * plot_sd[id]], + linewidth = 5, + size = (600, 600), + label = "", + ) + figpath = joinpath(figure_save_directory, "predict_" * expname * "_" * pn * "_at_truth.png") + savefig(figpath) + println("plot saved at " * figpath) + else + if data_mask == id + plot_mask = 1:length(data_mask) + plt = plot( + collect(id), + plot_mean[plot_mask], + show = true, + ribbon = [2 * plot_sd[plot_mask]; 2 * plot_sd[plot_mask]], + linewidth = 5, + size = (600, 600), + label = "", + ) + figpath = joinpath(figure_save_directory, "predict_" * expname * "_" * pn * "_at_truth.png") + savefig(figpath) + println("plot saved at " * figpath) + end + end + + end + + + +end # for @time diff --git a/examples/GCM/priors.jld2 b/examples/GCM/priors.jld2 new file mode 100644 index 000000000..860634326 Binary files /dev/null and b/examples/GCM/priors.jld2 differ diff --git a/examples/GCM/sbatch_emulate_sample_jl b/examples/GCM/sbatch_emulate_sample_jl new file mode 100644 index 000000000..6d2bada2a --- /dev/null +++ b/examples/GCM/sbatch_emulate_sample_jl @@ -0,0 +1,36 @@ +#!/bin/bash +#Submit this script with: sbatch thefilename +#SBATCH --time=12:00:00 # walltime +#SBATCH --ntasks-per-node=1 # number of processor cores (i.e. tasks) +#SBATCH --cpus-per-task=28 +#SBATCH --mem-per-cpu=8000 +#SBATCH -J "emulate_sample" # job name +#SBATCH --output=output/out_err/slurm_%j.out +#SBATCH --error=output/out_err/slurm_%j.err + +#general +set -euo pipefail #kill job if anything fails\ +#set -x # + +#modules (not automatically loaded with session) +module load julia/1.8.5 + +#julia package management + +export JULIA_PROJECT=@. +#precompiling is now done manually before +#julia -e 'using Pkg; Pkg.instantiate(); Pkg.API.precompile()' + +#run code +start=$(date +%s) + + +julia --project -t 28 emulate_sample_script.jl + +end=$(date +%s) +runtime=$((end-start)) +echo "********************" +echo "run time: $runtime" +echo "********************" + + diff --git a/examples/GCM/truthobj_meanparam.jld2 b/examples/GCM/truthobj_meanparam.jld2 new file mode 100644 index 000000000..cef3b8f9c Binary files /dev/null and b/examples/GCM/truthobj_meanparam.jld2 differ diff --git a/examples/GCM/truthobj_physparam.jld2 b/examples/GCM/truthobj_physparam.jld2 new file mode 100644 index 000000000..b2ab8e962 Binary files /dev/null and b/examples/GCM/truthobj_physparam.jld2 differ diff --git a/examples/GCM/truthobj_randparam.jld2 b/examples/GCM/truthobj_randparam.jld2 new file mode 100644 index 000000000..248ce29f1 Binary files /dev/null and b/examples/GCM/truthobj_randparam.jld2 differ diff --git a/examples/Lorenz/calibrate.jl b/examples/Lorenz/calibrate.jl index 79038276a..5e8974df6 100644 --- a/examples/Lorenz/calibrate.jl +++ b/examples/Lorenz/calibrate.jl @@ -39,7 +39,7 @@ function main() dynamics = 2 # Transient is 2 # Statistics integration length # This has to be less than 360 and 360 must be divisible by Ts_days - Ts_days = 90.0 # Integration length in days + Ts_days = 30.0 # Integration length in days # Stats type, which statistics to construct from the L96 system # 4 is a linear fit over a batch of length Ts_days # 5 is the mean over a batch of length Ts_days @@ -190,7 +190,7 @@ function main() # EKP parameters N_ens = 20 # number of ensemble members - N_iter = 5 # number of EKI iterations + N_iter = 6 # number of EKI iterations # initial parameters: N_params x N_ens initial_params = EKP.construct_initial_ensemble(priors, N_ens; rng_seed = rng_seed) diff --git a/examples/Lorenz/emulate_sample.jl b/examples/Lorenz/emulate_sample.jl index 46e2286ba..f1f99569f 100644 --- a/examples/Lorenz/emulate_sample.jl +++ b/examples/Lorenz/emulate_sample.jl @@ -1,6 +1,5 @@ # Import modules include(joinpath(@__DIR__, "..", "ci", "linkfig.jl")) -include(joinpath(@__DIR__, "GModel.jl")) # Contains Lorenz 96 source code # Import modules using Distributions # probability distributions and associated functions @@ -35,180 +34,212 @@ end function main() - ##### - # Should be loaded: - τc = 5.0 - F_true = 8.0 # Mean F - A_true = 2.5 # Transient F amplitude - ω_true = 2.0 * π / (360.0 / τc) # Frequency of the transient F - #### - - exp_name = "Lorenz_histogram_F$(F_true)_A$(A_true)-w$(ω_true)" - rng_seed = 44009 - Random.seed!(rng_seed) - # loading relevant data - homedir = pwd() - println(homedir) - figure_save_directory = joinpath(homedir, "output/") - data_save_directory = joinpath(homedir, "output/") - data_save_file = joinpath(data_save_directory, "calibrate_results.jld2") - - if !isfile(data_save_file) - throw( - ErrorException( - "data file $data_save_file not found. \n First run: \n > julia --project calibrate.jl \n and store results $data_save_file", - ), - ) - end - ekiobj = load(data_save_file)["eki"] - priors = load(data_save_file)["priors"] - truth_sample = load(data_save_file)["truth_sample"] - truth_sample_mean = load(data_save_file)["truth_sample_mean"] - truth_params_constrained = load(data_save_file)["truth_input_constrained"] #true parameters in constrained space - truth_params = transform_constrained_to_unconstrained(priors, truth_params_constrained) - Γy = ekiobj.obs_noise_cov - ### - ### Emulate: Gaussian Process Regression - ### - - # Emulate-sample settings - standardize = true - retained_svd_frac = 0.95 - - gppackage = Emulators.GPJL() - pred_type = Emulators.YType() - gauss_proc = GaussianProcess( - gppackage; - kernel = nothing, # use default squared exponential kernel - prediction_type = pred_type, - noise_learn = false, - ) - - # Standardize the output data - # Use median over all data since all data are the same type - truth_sample_norm = vcat(truth_sample...) - norm_factor = get_standardizing_factors(truth_sample_norm) - println(size(norm_factor)) - #norm_factor = vcat(norm_factor...) - norm_factor = fill(norm_factor, size(truth_sample)) - println("Standardization factors") - println(norm_factor) - - # Get training points from the EKP iteration number in the second input term - N_iter = 5 - input_output_pairs = Utilities.get_training_points(ekiobj, N_iter) - # Save data - @save joinpath(data_save_directory, "input_output_pairs.jld2") input_output_pairs - - normalized = true - emulator = Emulator( - gauss_proc, - input_output_pairs; - obs_noise_cov = Γy, - normalize_inputs = normalized, - standardize_outputs = standardize, - standardize_outputs_factors = norm_factor, - retained_svd_frac = retained_svd_frac, - ) - optimize_hyperparameters!(emulator) - - # Check how well the Gaussian Process regression predicts on the - # true parameters - #if retained_svd_frac==1.0 - y_mean, y_var = Emulators.predict(emulator, reshape(truth_params, :, 1), transform_to_real = true) - - println("GP prediction on true parameters: ") - println(vec(y_mean)) - println(" GP variance") - println(diag(y_var[1], 0)) - println("true data: ") - println(truth_sample_mean) # same, regardless of norm_factor - println("GP MSE: ") - println(mean((truth_sample_mean - vec(y_mean)) .^ 2)) - - #end - ### - ### Sample: Markov Chain Monte Carlo - ### - # initial values - u0 = vec(mean(get_inputs(input_output_pairs), dims = 2)) - println("initial parameters: ", u0) - - # First let's run a short chain to determine a good step size - truth_sample = truth_sample - mcmc = MCMCWrapper(RWMHSampling(), truth_sample, priors, emulator; init_params = u0) - new_step = optimize_stepsize(mcmc; init_stepsize = 0.1, N = 2000, discard_initial = 0) - - # Now begin the actual MCMC - println("Begin MCMC - with step size ", new_step) - chain = MarkovChainMonteCarlo.sample(mcmc, 100_000; stepsize = new_step, discard_initial = 2_000) - posterior = MarkovChainMonteCarlo.get_posterior(mcmc, chain) - - post_mean = mean(posterior) - post_cov = cov(posterior) - println("post_mean") - println(post_mean) - println("post_cov") - println(post_cov) - println("D util") - println(det(inv(post_cov))) - println(" ") - - - # Plot the posteriors together with the priors and the true parameter values (in the constrained space) - n_params = length(truth_params) - save_directory = joinpath(figure_save_directory) - - posterior_distribution_unconstrained_dict = get_distribution(posterior) #as it is a Samples, the distribution are samples - posterior_samples_unconstrained_dict = get_distribution(posterior) #as it is a Samples, the distribution are samples - param_names = get_name(posterior) - posterior_samples_unconstrained_arr = vcat([posterior_samples_unconstrained_dict[n] for n in param_names]...) - posterior_samples = transform_unconstrained_to_constrained(priors, posterior_samples_unconstrained_arr) - - for idx in 1:n_params - if idx == 1 - param = "F" - xbounds = [F_true - 1.0, F_true + 1.0] - xs = collect(xbounds[1]:((xbounds[2] - xbounds[1]) / 100):xbounds[2]) - elseif idx == 2 - param = "A" - xbounds = [A_true - 1.0, A_true + 1.0] - xs = collect(xbounds[1]:((xbounds[2] - xbounds[1]) / 100):xbounds[2]) - elseif idx == 3 - param = "ω" - xs = collect((ω_true - 0.2):0.01:(ω_true + 0.2)) - xbounds = [xs[1], xs[end]] - else - throw("not implemented") + cases = [ + "GP", # diagonalize, train scalar GP, assume diag inputs + "RF-scalar-diagin", # diagonalize, train scalar RF, assume diag inputs (most comparable to GP) + "RF-scalar", # diagonalize, train scalar RF, don't asume diag inputs + "RF-vector-svd-diag", + "RF-vector-svd-nondiag", + "RF-vector-nosvd-diag", + "RF-vector-nosvd-nondiag", + ] + + mask = 4:7 + #### CHOOSE YOUR CASE: + for case in cases[mask] + #case = cases[7] + println("case: ", case) + + # overrides = Dict("n_iteration" => 10) # "defaults" + overrides = Dict("n_iteration" => 10, "inflation" => 1e-4) # "large ens" + + ### + ##### + # Should be loaded: + τc = 5.0 + F_true = 8.0 # Mean F + A_true = 2.5 # Transient F amplitude + ω_true = 2.0 * π / (360.0 / τc) # Frequency of the transient F + #### + + exp_name = "Lorenz_histogram_F$(F_true)_A$(A_true)-w$(ω_true)" + rng_seed = 44009 + Random.seed!(rng_seed) + # loading relevant data + homedir = pwd() + println(homedir) + figure_save_directory = joinpath(homedir, "output/") + data_save_directory = joinpath(homedir, "output/") + data_save_file = joinpath(data_save_directory, "calibrate_results.jld2") + + if !isfile(data_save_file) + throw( + ErrorException( + "data file $data_save_file not found. \n First run: \n > julia --project calibrate.jl \n and store results $data_save_file", + ), + ) end - label = "true " * param - - histogram(posterior_samples[idx, :], bins = 100, normed = true, fill = :slategray, lab = "posterior") - prior_plot = get_distribution(mcmc.prior) + ekiobj = load(data_save_file)["eki"] + priors = load(data_save_file)["priors"] + truth_sample = load(data_save_file)["truth_sample"] + truth_sample_mean = load(data_save_file)["truth_sample_mean"] + truth_params_constrained = load(data_save_file)["truth_input_constrained"] #true parameters in constrained space + truth_params = transform_constrained_to_unconstrained(priors, truth_params_constrained) + Γy = ekiobj.obs_noise_cov + + n_params = length(truth_params) # "input dim" + output_dim = size(Γy, 1) + ### + ### Emulate: Gaussian Process Regression + ### + + # Emulate-sample settings + # choice of machine-learning tool + if case == "GP" + gppackage = Emulators.GPJL() + pred_type = Emulators.YType() + mlt = GaussianProcess( + gppackage; + kernel = nothing, # use default squared exponential kernel + prediction_type = pred_type, + noise_learn = false, + ) + elseif case ∈ ["RF-scalar", "RF-scalar-diagin"] + n_features = 300 + batch_sizes = Dict("train" => 100, "test" => 100, "feature" => 100) + diagonalize_input = case == "RF-scalar-diagin" ? true : false + mlt = ScalarRandomFeatureInterface( + n_features, + n_params, + diagonalize_input = diagonalize_input, + batch_sizes = batch_sizes, + ) + elseif case ∈ ["RF-vector-svd-diag", "RF-vector-nosvd-diag", "RF-vector-svd-nondiag", "RF-vector-nosvd-nondiag"] + # do we want to assume that the outputs are decorrelated in the machine-learning problem? + diagonalize_output = case ∈ ["RF-vector-svd-diag", "RF-vector-nosvd-diag"] ? true : false + n_features = 300 + batch_sizes = Dict("train" => 200, "test" => 200, "feature" => 200) + mlt = VectorRandomFeatureInterface( + n_features, + n_params, + output_dim, + batch_sizes = batch_sizes, + diagonalize_output = diagonalize_output, # assume outputs are decorrelated + optimizer_options = overrides, + ) + end - # This requires StatsPlots - xs_rows = permutedims(repeat(xs, 1, size(posterior_samples, 1)), (2, 1)) # This hack will enable us to apply the transformation using the full prior - xs_rows_unconstrained = transform_constrained_to_unconstrained(priors, xs_rows) - plot!( - xs, - pdf.(prior_plot[param_names[idx]], xs_rows_unconstrained[idx, :]), - w = 2.6, - color = :blue, - lab = "prior", + # Standardize the output data + # Use median over all data since all data are the same type + truth_sample_norm = vcat(truth_sample...) + norm_factor = get_standardizing_factors(truth_sample_norm) + println(size(norm_factor)) + #norm_factor = vcat(norm_factor...) + norm_factor = fill(norm_factor, size(truth_sample)) + println("Standardization factors") + println(norm_factor) + + # Get training points from the EKP iteration number in the second input term + N_iter = 6 # note we have 6 iterations stored + input_output_pairs = Utilities.get_training_points(ekiobj, N_iter - 1) # 1:N-1 + input_output_pairs_test = Utilities.get_training_points(ekiobj, N_iter:N_iter) # just "next" iteration used for testing + # Save data + @save joinpath(data_save_directory, "input_output_pairs.jld2") input_output_pairs + + standardize = true + retained_svd_frac = 1.0 + normalized = true + # do we want to use SVD to decorrelate outputs + decorrelate = case ∈ ["RF-vector-nosvd-diag", "RF-vector-nosvd-nondiag"] ? false : true + + emulator = Emulator( + mlt, + input_output_pairs; + obs_noise_cov = Γy, + normalize_inputs = normalized, + standardize_outputs = standardize, + standardize_outputs_factors = norm_factor, + retained_svd_frac = retained_svd_frac, + decorrelate = decorrelate, ) - #plot!(xs, mcmc.prior[idx].dist, w=2.6, color=:blue, lab="prior") - plot!([truth_params[idx]], seriestype = "vline", w = 2.6, lab = label) - plot!(xlims = xbounds) + optimize_hyperparameters!(emulator) + + # Check how well the Gaussian Process regression predicts on the + # true parameters + #if retained_svd_frac==1.0 + y_mean, y_var = Emulators.predict(emulator, reshape(truth_params, :, 1), transform_to_real = true) + y_mean_test, y_var_test = + Emulators.predict(emulator, get_inputs(input_output_pairs_test), transform_to_real = true) + + println("ML prediction on true parameters: ") + println(vec(y_mean)) + println("true data: ") + println(truth_sample_mean) # same, regardless of norm_factor + println(" ML predicted variance") + println(diag(y_var[1], 0)) + println("ML MSE (truth): ") + println(mean((truth_sample_mean - vec(y_mean)) .^ 2)) + println("ML MSE (next ensemble): ") + println(mean((get_outputs(input_output_pairs_test) - y_mean_test) .^ 2)) + + #end + ### + ### Sample: Markov Chain Monte Carlo + ### + # initial values + u0 = vec(mean(get_inputs(input_output_pairs), dims = 2)) + println("initial parameters: ", u0) + + # First let's run a short chain to determine a good step size + truth_sample = truth_sample + mcmc = MCMCWrapper(RWMHSampling(), truth_sample, priors, emulator; init_params = u0) + new_step = optimize_stepsize(mcmc; init_stepsize = 0.1, N = 2000, discard_initial = 0) + + # Now begin the actual MCMC + println("Begin MCMC - with step size ", new_step) + chain = MarkovChainMonteCarlo.sample(mcmc, 100_000; stepsize = new_step, discard_initial = 2_000) + posterior = MarkovChainMonteCarlo.get_posterior(mcmc, chain) + + post_mean = mean(posterior) + post_cov = cov(posterior) + println("post_mean") + println(post_mean) + println("post_cov") + println(post_cov) + println("D util") + println(det(inv(post_cov))) + println(" ") + + constrained_truth_params = transform_unconstrained_to_constrained(posterior, truth_params) + param_names = get_name(posterior) + + posterior_samples = vcat([get_distribution(posterior)[name] for name in get_name(posterior)]...) #samples are columns + constrained_posterior_samples = + mapslices(x -> transform_unconstrained_to_constrained(posterior, x), posterior_samples, dims = 1) + + gr(dpi = 300, size = (300, 300)) + p = cornerplot(permutedims(constrained_posterior_samples, (2, 1)), label = param_names, compact = true) + plot!(p.subplots[1], [constrained_truth_params[1]], seriestype = "vline", w = 1.5, c = :steelblue, ls = :dash) # vline on top histogram + plot!(p.subplots[3], [constrained_truth_params[2]], seriestype = "hline", w = 1.5, c = :steelblue, ls = :dash) # hline on right histogram + plot!(p.subplots[2], [constrained_truth_params[1]], seriestype = "vline", w = 1.5, c = :steelblue, ls = :dash) # v & h line on scatter. + plot!(p.subplots[2], [constrained_truth_params[2]], seriestype = "hline", w = 1.5, c = :steelblue, ls = :dash) + figpath = joinpath(figure_save_directory, "posterior_2d-" * case * ".png") + savefig(figpath) - title!(param) - figpath = joinpath(figure_save_directory, "posterior_$(param)_" * exp_name * ".png") - savefig(figpath) - linkfig(figpath) + # Save data + save( + joinpath(data_save_directory, "posterior.jld2"), + "posterior", + posterior, + "input_output_pairs", + input_output_pairs, + "truth_params", + truth_params, + ) end - end main() diff --git a/src/CalibrateEmulateSample.jl b/src/CalibrateEmulateSample.jl index 9f227d0ee..f5956b175 100644 --- a/src/CalibrateEmulateSample.jl +++ b/src/CalibrateEmulateSample.jl @@ -14,13 +14,14 @@ import EnsembleKalmanProcesses: EnsembleKalmanProcesses, ParameterDistributions, export EnsembleKalmanProcesses, ParameterDistributions, Observations, DataContainers -# No internal deps, heavy external deps -#include("GaussianProcessEmulator.jl") -include("Emulator.jl") # Internal deps, light external deps include("Utilities.jl") +# No internal deps, heavy external deps +#include("GaussianProcessEmulator.jl") +include("Emulator.jl") + # Internal deps, light external deps include("MarkovChainMonteCarlo.jl") diff --git a/src/Emulator.jl b/src/Emulator.jl index e4fa76aa8..e6b085ff4 100644 --- a/src/Emulator.jl +++ b/src/Emulator.jl @@ -7,6 +7,7 @@ using Statistics using Distributions using LinearAlgebra using DocStringExtensions +using Random export Emulator @@ -24,7 +25,7 @@ abstract type MachineLearningTool end # include the different <: ML models include("GaussianProcess.jl") #for GaussianProcess -# include("RandomFeature.jl") +include("RandomFeature.jl") # include("NeuralNetwork.jl") # etc. @@ -75,6 +76,8 @@ struct Emulator{FT <: AbstractFloat} retained_svd_frac::FT end +get_machine_learning_tool(emulator::Emulator) = emulator.machine_learning_tool + # Constructor for the Emulator Object function Emulator( machine_learning_tool::MachineLearningTool, @@ -83,6 +86,7 @@ function Emulator( normalize_inputs::Bool = true, standardize_outputs::Bool = false, standardize_outputs_factors::Union{AbstractVector{FT}, Nothing} = nothing, + decorrelate::Bool = true, retained_svd_frac::FT = 1.0, ) where {FT <: AbstractFloat} @@ -118,24 +122,35 @@ function Emulator( training_outputs = get_outputs(input_output_pairs) end - # [3.] Decorrelating the outputs [always performed] - - #Transform data if obs_noise_cov available - # (if obs_noise_cov==nothing, transformed_data is equal to data) - decorrelated_training_outputs, decomposition = - svd_transform(training_outputs, obs_noise_cov, retained_svd_frac = retained_svd_frac) + # [3.] Decorrelating the outputs, not performed for vector RF + if decorrelate + #Transform data if obs_noise_cov available + # (if obs_noise_cov==nothing, transformed_data is equal to data) + decorrelated_training_outputs, decomposition = + svd_transform(training_outputs, obs_noise_cov, retained_svd_frac = retained_svd_frac) + + # write new pairs structure + if retained_svd_frac < 1.0 + #note this changes the dimension of the outputs + training_pairs = PairedDataContainer(training_inputs, decorrelated_training_outputs) + input_dim, output_dim = size(training_pairs, 1) + else + training_pairs = PairedDataContainer(training_inputs, decorrelated_training_outputs) + end - # write new pairs structure - if retained_svd_frac < 1.0 - #note this changes the dimension of the outputs - training_pairs = PairedDataContainer(training_inputs, decorrelated_training_outputs) - input_dim, output_dim = size(training_pairs, 1) + # [4.] build an emulator + build_models!(machine_learning_tool, training_pairs) else - training_pairs = PairedDataContainer(training_inputs, decorrelated_training_outputs) + if decorrelate || !isa(machine_learning_tool, VectorRandomFeatureInterface) + throw(ArgumentError("$machine_learning_tool is incompatible with option Emulator(...,decorrelate = false)")) + end + decomposition = nothing + training_pairs = PairedDataContainer(training_inputs, training_outputs) + + # [4.] build an emulator - providing the noise covariance as a Tikhonov regularizer + build_models!(machine_learning_tool, training_pairs, regularization_matrix = obs_noise_cov) end - # [4.] build an emulator - build_models!(machine_learning_tool, training_pairs) return Emulator{FT}( machine_learning_tool, @@ -170,6 +185,7 @@ function predict( emulator::Emulator{FT}, new_inputs::AbstractMatrix{FT}; transform_to_real = false, + vector_rf_unstandardize = true, ) where {FT <: AbstractFloat} # Check if the size of new_inputs is consistent with the GP model's input # dimension. @@ -177,7 +193,7 @@ function predict( N_samples = size(new_inputs, 2) size(new_inputs, 1) == input_dim || - throw(ArgumentError("GP object and input observations do not have consistent dimensions")) + throw(ArgumentError("Emulator object and input observations do not have consistent dimensions")) # [1.] normalize normalized_new_inputs = normalize(emulator, new_inputs) @@ -186,28 +202,70 @@ function predict( ds_outputs, ds_output_var = predict(emulator.machine_learning_tool, normalized_new_inputs) # [3.] transform back to real coordinates or remain in decorrelated coordinates - if transform_to_real && emulator.decomposition === nothing - throw(ArgumentError("""Need SVD decomposition to transform back to original space, - but GaussianProcess.decomposition == nothing. - Try setting transform_to_real=false""")) - elseif transform_to_real && emulator.decomposition !== nothing - #transform back to real coords - cov becomes dense - s_outputs, s_output_cov = svd_reverse_transform_mean_cov(ds_outputs, ds_output_var, emulator.decomposition) - if output_dim == 1 - s_output_cov = [s_output_cov[i][1] for i in 1:N_samples] + if !isa(get_machine_learning_tool(emulator), VectorRandomFeatureInterface) + if transform_to_real && emulator.decomposition === nothing + throw(ArgumentError("""Need SVD decomposition to transform back to original space, + but GaussianProcess.decomposition == nothing. + Try setting transform_to_real=false""")) + elseif transform_to_real && emulator.decomposition !== nothing + #transform back to real coords - cov becomes dense + s_outputs, s_output_cov = svd_reverse_transform_mean_cov(ds_outputs, ds_output_var, emulator.decomposition) + if output_dim == 1 + s_output_cov = [s_output_cov[i][1] for i in 1:N_samples] + end + # [4.] unstandardize + return reverse_standardize(emulator, s_outputs, s_output_cov) + else + # remain in decorrelated, standardized coordinates (cov remains diagonal) + # Convert to vector of matrices to match the format + # when transform_to_real=true + ds_output_diagvar = vec([Diagonal(ds_output_var[:, j]) for j in 1:N_samples]) + if output_dim == 1 + ds_output_diagvar = [ds_output_diagvar[i][1] for i in 1:N_samples] + end + return ds_outputs, ds_output_diagvar end - # [4.] unstandardize - return reverse_standardize(emulator, s_outputs, s_output_cov) else - # remain in decorrelated, standardized coordinates (cov remains diagonal) - # Convert to vector of matrices to match the format - # when transform_to_real=true - ds_output_diagvar = vec([Diagonal(ds_output_var[:, j]) for j in 1:N_samples]) - if output_dim == 1 - ds_output_diagvar = [ds_output_diagvar[i][1] for i in 1:N_samples] + # create a vector of covariance matrices + ds_output_covvec = vec([ds_output_var[:, :, j] for j in 1:N_samples]) + #=if output_dim == 1 + ds_output_covvec = [reshape(ds_output_covvec[i][1, 1],1,1) for i in 1:N_samples] + end + =# + + if emulator.decomposition === nothing # without applying SVD + if vector_rf_unstandardize + # [4.] unstandardize + return reverse_standardize(emulator, ds_outputs, ds_output_covvec) + else + return ds_outputs, ds_output_covvec + end + + else #if we applied SVD + if transform_to_real # ...and want to return coordinates + s_outputs, s_output_cov = + svd_reverse_transform_mean_cov(ds_outputs, ds_output_covvec, emulator.decomposition) + if output_dim == 1 + s_output_cov = [s_output_cov[i][1] for i in 1:N_samples] + end + # [4.] unstandardize + return reverse_standardize(emulator, s_outputs, s_output_cov) + else #... and want to stay in decorrelated standardized coordinates + diag_out_flag = get_diagonalize_output(get_machine_learning_tool(emulator)) + if diag_out_flag + ds_output_diagvar = vec([Diagonal(ds_output_covvec[j]) for j in 1:N_samples]) # extracts diagonal + if output_dim == 1 + ds_output_diagvar = [ds_output_covvec[i][1, 1] for i in 1:N_samples] # extracts value + end + return ds_outputs, ds_output_diagvar + else + return ds_outputs, ds_output_covvec + end + end end - return ds_outputs, ds_output_diagvar + end + end # Normalization, Standardization, and Decorrelation @@ -375,6 +433,7 @@ $(DocStringExtensions.TYPEDSIGNATURES) Transform the mean and covariance back to the original (correlated) coordinate system - `μ` - predicted mean; size *output\\_dim* × *N\\_predicted\\_points*. - `σ2` - predicted variance; size *output\\_dim* × *N\\_predicted\\_points*. + - predicted covariance; size *N\\_predicted\\_points* array of size *output\\_dim* × *output\\_dim*. - `decomposition` - SVD decomposition of *obs\\_noise\\_cov*. Returns the transformed mean (size *output\\_dim* × *N\\_predicted\\_points*) and variance. @@ -386,12 +445,20 @@ each element is a matrix of size *output\\_dim* × *output\\_dim*. """ function svd_reverse_transform_mean_cov( μ::AbstractMatrix{FT}, - σ2::AbstractMatrix{FT}, + σ2::AbstractVector,# vector of output_dim x output_dim matrices decomposition::SVD, ) where {FT <: AbstractFloat} - @assert size(μ) == size(σ2) - output_dim, N_predicted_points = size(σ2) + N_predicted_points = length(σ2) + if !(eltype(σ2) <: AbstractMatrix) + throw( + ArgumentError( + "input σ2 must be a vector of eltype <: AbstractMatrix, instead eltype(σ2) = ", + string(eltype(σ2)), + ), + ) + end + # We created meanvGP = D_inv * Vt * mean_v so meanv = V * D * meanvGP sqrt_singular_values = Diagonal(sqrt.(decomposition.S)) transformed_μ = decomposition.V * sqrt_singular_values * μ @@ -400,13 +467,27 @@ function svd_reverse_transform_mean_cov( # Back transformation for j in 1:N_predicted_points - σ2_j = decomposition.V * sqrt_singular_values * Diagonal(σ2[:, j]) * sqrt_singular_values * decomposition.Vt + σ2_j = decomposition.V * sqrt_singular_values * σ2[j] * sqrt_singular_values * decomposition.Vt transformed_σ2[j] = Symmetric(σ2_j) end return transformed_μ, transformed_σ2 + +end + +function svd_reverse_transform_mean_cov( + μ::AbstractMatrix{FT}, + σ2::AbstractMatrix{FT}, + decomposition::SVD, +) where {FT <: AbstractFloat} + @assert size(μ) == size(σ2) + output_dim, N_predicted_points = size(σ2) + σ2_as_cov = vec([Diagonal(σ2[:, j]) for j in 1:N_predicted_points]) + + return svd_reverse_transform_mean_cov(μ, σ2_as_cov, decomposition) end + end diff --git a/src/MarkovChainMonteCarlo.jl b/src/MarkovChainMonteCarlo.jl index 1abb1134b..2e376b32e 100644 --- a/src/MarkovChainMonteCarlo.jl +++ b/src/MarkovChainMonteCarlo.jl @@ -165,7 +165,10 @@ function EmulatorPosteriorModel( # Vector of N_samples covariance matrices. For MH, N_samples is always 1, so we # have to reshape()/re-cast input/output; simpler to do here than add a # predict() method. - g, g_cov = Emulators.predict(em, reshape(θ, :, 1), transform_to_real = false) + g, g_cov = + Emulators.predict(em, reshape(θ, :, 1), transform_to_real = false, vector_rf_unstandardize = false) + #TODO vector_rf will always unstandardize, but other methods will not, so we require this additional flag. + if isa(g_cov[1], Real) return logpdf(MvNormal(obs_sample, g_cov[1] * I), vec(g)) + get_logpdf(prior, θ) else diff --git a/src/RandomFeature.jl b/src/RandomFeature.jl new file mode 100644 index 000000000..48bf696eb --- /dev/null +++ b/src/RandomFeature.jl @@ -0,0 +1,283 @@ + +using RandomFeatures +const RF = RandomFeatures +using EnsembleKalmanProcesses +const EKP = EnsembleKalmanProcesses +using ..ParameterDistributions +using ..Utilities +using StableRNGs +using Random + +abstract type RandomFeatureInterface <: MachineLearningTool end + + +#common functions +function flat_to_chol(x::AbstractVector) + choldim = Int(floor(sqrt(2 * length(x)))) + cholmat = zeros(choldim, choldim) + for i in 1:choldim + for j in 1:i + cholmat[i, j] = x[sum(0:(i - 1)) + j] + end + end + return cholmat +end + +function hyperparameters_from_flat( + x::AbstractVector, + input_dim::Int, + output_dim::Int, + diagonalize_input::Bool, + diagonalize_output::Bool, +) + # if dim = 1 then parameters are a 1x1 matrix + # if we diagonalize then parameters are diagonal entries + constant + # if we don't diagonalize then parameters are a cholesky product + constant on diagonal. + + #input space + if input_dim == 1 + in_max = 1 + U = x[in_max] * ones(1, 1) + elseif diagonalize_input + in_max = input_dim + 2 + U = convert(Matrix, x[in_max - 1] * (Diagonal(x[1:(in_max - 2)]) + x[in_max] * I)) + elseif !diagonalize_input + in_max = Int(input_dim * (input_dim + 1) / 2) + 2 + cholU = flat_to_chol(x[1:(in_max - 2)]) + U = x[in_max - 1] * (cholU * permutedims(cholU, (2, 1)) + x[in_max] * I) + end + + # output_space + out_min = in_max + 1 + if output_dim == 1 + V = x[end] * ones(1, 1) + elseif diagonalize_output + V = convert(Matrix, x[end - 1] * (Diagonal(x[out_min:(end - 2)]) + x[end] * I)) + elseif !diagonalize_output + cholV = flat_to_chol(x[out_min:(end - 2)]) + V = x[end - 1] * (cholV * permutedims(cholV, (2, 1)) + x[end] * I) + end + + return U, V + +end + +function hyperparameters_from_flat(x::AbstractVector, input_dim::Int, diagonalize_input::Bool) + output_dim = 1 + diagonalize_output = false + U, V = hyperparameters_from_flat(x, input_dim, output_dim, diagonalize_input, diagonalize_output) + # Note that the scalar setting has a slight inconsistency from the 1-D output vector case + # We correct here by rescaling U using the single hyperparameter in V + return V[1, 1] * U +end + +function calculate_n_hyperparameters(input_dim::Int, output_dim::Int, diagonalize_input::Bool, diagonalize_output::Bool) + + n_hp = 0 + # inputs: diagonal or cholesky + n_hp += diagonalize_input ? input_dim : Int(input_dim * (input_dim + 1) / 2) + n_hp += diagonalize_output ? output_dim : Int(output_dim * (output_dim + 1) / 2) + + # add two more for constant diagonal in each case. + n_hp += (input_dim > 1) ? 2 : 0 + n_hp += (output_dim > 1) ? 2 : 0 + return n_hp + +end + +function calculate_n_hyperparameters(input_dim::Int, diagonalize_input::Bool) + #scalar case consistent with vector case + output_dim = 1 + diagonalize_output = false + return calculate_n_hyperparameters(input_dim, output_dim, diagonalize_input, diagonalize_output) +end + + + +function calculate_mean_cov_and_coeffs( + rfi::RFI, + rng::AbstractRNG, + l::Union{Real, AbstractVecOrMat}, + regularization::Union{Matrix, UniformScaling, Diagonal}, + n_features::Int, + n_train::Int, + n_test::Int, + batch_sizes::Union{Dict, Nothing}, + io_pairs::PairedDataContainer; +) where {RFI <: RandomFeatureInterface} + + # split data into train/test + itrain = get_inputs(io_pairs)[:, 1:n_train] + otrain = get_outputs(io_pairs)[:, 1:n_train] + io_train_cost = PairedDataContainer(itrain, otrain) + itest = get_inputs(io_pairs)[:, (n_train + 1):end] + otest = get_outputs(io_pairs)[:, (n_train + 1):end] + input_dim = size(itrain, 1) + output_dim = size(otrain, 1) + + # build and fit the RF + rfm = RFM_from_hyperparameters(rfi, rng, l, regularization, n_features, batch_sizes, input_dim, output_dim) + decomp_type = "svd" + if decomp_type == "chol" + fitted_features = RF.Methods.fit(rfm, io_train_cost, decomposition_type = "cholesky") + else + fitted_features = RF.Methods.fit(rfm, io_train_cost, decomposition_type = "svd") + end + test_batch_size = RF.Methods.get_batch_size(rfm, "test") + batch_inputs = RF.Utilities.batch_generator(itest, test_batch_size, dims = 2) # input_dim x batch_size + + #we want to calc 1/var(y-mean)^2 + lambda/m * coeffs^2 in the end + pred_mean, pred_cov = RF.Methods.predict(rfm, fitted_features, DataContainer(itest)) + # sizes (output_dim x n_test), (output_dim x output_dim x n_test) + scaled_coeffs = sqrt(1 / n_features) * RF.Methods.get_coeffs(fitted_features) + + #we add the additional complexity term + if decomp_type == "chol" + 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 + complexity = sum(log, svd_singval) # note this is log(abs(det)) + end + return pred_mean, pred_cov, scaled_coeffs, complexity + +end + +function estimate_mean_and_coeffnorm_covariance( + rfi::RFI, + rng::AbstractRNG, + l::Union{Real, AbstractVecOrMat}, + regularization::Union{Matrix, UniformScaling, Diagonal}, + n_features::Int, + n_train::Int, + n_test::Int, + batch_sizes::Union{Dict, Nothing}, + io_pairs::PairedDataContainer, + n_samples::Int; + repeats::Int = 1, +) where {RFI <: RandomFeatureInterface} + output_dim = size(get_outputs(io_pairs), 1) + + means = zeros(n_test, n_samples, output_dim) + mean_of_covs = zeros(output_dim, output_dim, n_test) + complexity = zeros(1, n_samples) + coeffl2norm = zeros(1, n_samples) + println("estimate cov with " * string(n_samples * repeats) * " iterations...") + + for i in ProgressBar(1:n_samples) + for j in 1:repeats + + m, v, c, cplxty = calculate_mean_cov_and_coeffs( + rfi, + rng, + l, + regularization, + n_features, + n_train, + n_test, + batch_sizes, + io_pairs, + ) + # m output_dim x n_test + # v output_dim x output_dim x n_test + # c n_features + + # update vbles needed for cov + means[:, i, :] .+= m' ./ repeats + coeffl2norm[1, i] += sqrt(sum(abs2, c)) / repeats + complexity[1, i] += cplxty / repeats + + # update vbles needed for mean + mean_of_covs[:, :, :] .+= v ./ (repeats * n_samples) + + end + end + mean_of_covs = permutedims(mean_of_covs, (3, 1, 2)) + + approx_σ2 = zeros(n_test * output_dim, n_test * output_dim) + blockmeans = zeros(n_test * output_dim, n_samples) + for i in 1:n_test + id = ((i - 1) * output_dim + 1):(i * output_dim) + 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 + #symmeterize + approx_σ2 = 0.5 * (approx_σ2 + permutedims(approx_σ2, (2, 1))) + if !isposdef(approx_σ2) + println("approx_σ2 not posdef") + correction = abs(minimum(eigvals(approx_σ2))) + 1e8 * eps() + approx_σ2 += correction * I # to give a little help with pos def + end + + Γ = cov(vcat(blockmeans, coeffl2norm, complexity), dims = 2) + return Γ, approx_σ2 + +end + +function calculate_ensemble_mean_and_coeffnorm( + rfi::RFI, + rng::AbstractRNG, + lvecormat::AbstractVecOrMat, + regularization::Union{Matrix, UniformScaling, Diagonal}, + n_features::Int, + n_train::Int, + n_test::Int, + batch_sizes::Union{Dict, Nothing}, + io_pairs::PairedDataContainer; + repeats::Int = 1, +) where {RFI <: RandomFeatureInterface} + if isa(lvecormat, AbstractVector) + lmat = reshape(lvecormat, 1, :) + else + lmat = lvecormat + end + N_ens = size(lmat, 2) + output_dim = size(get_outputs(io_pairs), 1) + + means = zeros(n_test, N_ens, output_dim) + mean_of_covs = zeros(output_dim, output_dim, n_test) + complexity = zeros(1, N_ens) + coeffl2norm = zeros(1, N_ens) + println("calculating " * string(N_ens * repeats) * " ensemble members...") + + for i in ProgressBar(1:N_ens) + for j in collect(1:repeats) + l = lmat[:, i] + m, v, c, cplxty = calculate_mean_cov_and_coeffs( + rfi, + rng, + l, + regularization, + n_features, + n_train, + n_test, + batch_sizes, + io_pairs, + ) + # m output_dim x n_test + # v output_dim x output_dim x n_test + # c n_features + means[:, i, :] += m' ./ repeats + mean_of_covs[:, :, :] += v ./ (repeats * N_ens) + coeffl2norm[1, i] += sqrt(sum(c .^ 2)) / repeats + complexity[1, i] += cplxty / repeats + end + end + mean_of_covs = permutedims(mean_of_covs, (3, 1, 2)) + blockcovmat = zeros(n_test * output_dim, n_test * output_dim) + blockmeans = zeros(n_test * output_dim, N_ens) + for i in 1:n_test + id = ((i - 1) * output_dim + 1):(i * output_dim) + blockcovmat[id, id] = mean_of_covs[i, :, :] + blockmeans[id, :] = permutedims(means[i, :, :], (2, 1)) + end + #symmeterize + blockcovmat = 0.5 * (blockcovmat + permutedims(blockcovmat, (2, 1))) + + return vcat(blockmeans, coeffl2norm, complexity), blockcovmat + +end + + +include("ScalarRandomFeature.jl") +include("VectorRandomFeature.jl") diff --git a/src/ScalarRandomFeature.jl b/src/ScalarRandomFeature.jl new file mode 100644 index 000000000..15028570c --- /dev/null +++ b/src/ScalarRandomFeature.jl @@ -0,0 +1,263 @@ +export ScalarRandomFeatureInterface +export get_rfms, + get_fitted_features, + get_batch_sizes, + get_n_features, + get_input_dim, + get_diagonalize_input, + get_rng, + get_optimizer_options + + + +struct ScalarRandomFeatureInterface <: RandomFeatureInterface + "vector of `RandomFeatureMethod`s, contains the feature structure, batch-sizes and regularization" + rfms::Vector{RF.Methods.RandomFeatureMethod} + "vector of `Fit`s, containing the matrix decomposition and coefficients of RF when fitted to data" + fitted_features::Vector{RF.Methods.Fit} + "batch sizes" + batch_sizes::Union{Dict, Nothing} + "n_features" + n_features::Union{Int, Nothing} + "input dimension" + input_dim::Int + "choice of random number generator" + rng::AbstractRNG + "Assume inputs are decorrelated for ML" + diagonalize_input::Bool + "dictionary of options for hyperparameter optimizer" + optimizer_options::Dict +end + +get_rfms(srfi::ScalarRandomFeatureInterface) = srfi.rfms +get_fitted_features(srfi::ScalarRandomFeatureInterface) = srfi.fitted_features +get_batch_sizes(srfi::ScalarRandomFeatureInterface) = srfi.batch_sizes +get_n_features(srfi::ScalarRandomFeatureInterface) = srfi.n_features +get_input_dim(srfi::ScalarRandomFeatureInterface) = srfi.input_dim +get_rng(srfi::ScalarRandomFeatureInterface) = srfi.rng +get_diagonalize_input(srfi::ScalarRandomFeatureInterface) = srfi.diagonalize_input +get_optimizer_options(srfi::ScalarRandomFeatureInterface) = srfi.optimizer_options + +function ScalarRandomFeatureInterface( + n_features::Int, + input_dim::Int; + diagonalize_input::Bool = false, + batch_sizes::Union{Dict, Nothing} = nothing, + rng::AbstractRNG = Random.GLOBAL_RNG, + optimizer_options::Union{Dict, Nothing} = nothing, +) + # Initialize vector for GP models + rfms = Vector{RF.Methods.RandomFeatureMethod}(undef, 0) + fitted_features = Vector{RF.Methods.Fit}(undef, 0) + + n_hp = calculate_n_hyperparameters(input_dim, diagonalize_input) + prior = constrained_gaussian("MvNormal_parameters", 1.0, 1.0, 0.0, Inf, repeats = n_hp) + + # default optimizer settings + optimizer_opts = Dict( + "prior" => prior, #the hyperparameter_prior + "n_ensemble" => ndims(prior) + 1, #number of ensemble + "n_iteration" => 5, # number of eki iterations + "inflation" => 1e-4, # additive inflation ∈ [0,1] with 0 being no inflation + ) + + if !isnothing(optimizer_options) + for key in keys(optimizer_options) + optimizer_opts[key] = optimizer_options[key] + end + end + + return ScalarRandomFeatureInterface( + rfms, + fitted_features, + batch_sizes, + n_features, + input_dim, + rng, + diagonalize_input, + optimizer_opts, + ) +end + +function RFM_from_hyperparameters( + srfi::ScalarRandomFeatureInterface, + rng::AbstractRNG, + l::Union{Real, AbstractVecOrMat}, + regularization::Union{Matrix, UniformScaling, Diagonal}, # just a 1x1 matrix though + n_features::Int, + batch_sizes::Union{Dict, Nothing}, + input_dim::Int; +) + diagonalize_input = get_diagonalize_input(srfi) + + M = zeros(input_dim) #scalar output + U = hyperparameters_from_flat(l, input_dim, diagonalize_input) + if !isposdef(U) + println("U not posdef - adding to diagonal") + correction = abs(minimum(eigvals(U))) + 1e8 * eps() + U += correction * I + println("added $correction to regularize") + end + + dist = MvNormal(M, U) + pd = ParameterDistribution( + Dict( + "distribution" => Parameterized(dist), + "constraint" => repeat([no_constraint()], input_dim), + "name" => "xi", + ), + ) + feature_sampler = RF.Samplers.FeatureSampler(pd, rng = rng) + # Learn hyperparameters for different feature types + + sf = RF.Features.ScalarFourierFeature(n_features, feature_sampler) + if isnothing(batch_sizes) + return RF.Methods.RandomFeatureMethod(sf, regularization = regularization) + else + return RF.Methods.RandomFeatureMethod(sf, regularization = regularization, batch_sizes = batch_sizes) + end +end + +#removes vector-only input arguments +RFM_from_hyperparameters( + srfi::ScalarRandomFeatureInterface, + rng::AbstractRNG, + l::Union{Real, AbstractVecOrMat}, + regularization::Union{Matrix, UniformScaling}, # just a 1x1 matrix though + n_features::Int, + batch_sizes::Union{Dict, Nothing}, + input_dim::Int, + output_dim::Int; +) = RFM_from_hyperparameters(srfi, rng, l, regularization, n_features, batch_sizes, input_dim) + + +function build_models!( + srfi::ScalarRandomFeatureInterface, + input_output_pairs::PairedDataContainer{FT}, +) where {FT <: AbstractFloat} + # get inputs and outputs + input_values = get_inputs(input_output_pairs) + output_values = get_outputs(input_output_pairs) + n_rfms, n_data = size(output_values) + noise_sd = 1.0 + + input_dim = size(input_values, 1) + + rfms = get_rfms(srfi) + fitted_features = get_fitted_features(srfi) + n_features = get_n_features(srfi) + batch_sizes = get_batch_sizes(srfi) + rng = get_rng(srfi) + optimizer_options = get_optimizer_options(srfi) + + #regularization = I = 1.0 in scalar case + regularization = I + # Optimize features with EKP for each output dim + # [1.] Split data into test/train 80/20 + n_train = Int(floor(0.8 * n_data)) + n_test = n_data - n_train + n_features_opt = min(n_features, Int(floor(2 * n_test))) #we take a specified number of features for optimization. + idx_shuffle = randperm(rng, n_data) + + for i in 1:n_rfms + + + io_pairs_opt = PairedDataContainer( + input_values[:, idx_shuffle], + reshape(output_values[i, idx_shuffle], 1, size(output_values, 2)), + ) + + # [2.] Estimate covariance at mean value + prior = optimizer_options["prior"] + μ_hp = transform_unconstrained_to_constrained(prior, mean(prior)) + cov_samples = 2 * Int(floor(n_test)) + internal_Γ, approx_σ2 = estimate_mean_and_coeffnorm_covariance( + srfi, + rng, + μ_hp, + regularization, + n_features_opt, + n_train, + n_test, + batch_sizes, + io_pairs_opt, + cov_samples, + ) + + Γ = internal_Γ + Γ[1:n_test, 1:n_test] += approx_σ2 + Γ[(n_test + 1):end, (n_test + 1):end] += I + + # [3.] set up EKP optimization + n_ensemble = optimizer_options["n_ensemble"] + n_iteration = optimizer_options["n_iteration"] + + initial_params = construct_initial_ensemble(rng, prior, n_ensemble) + min_complexity = log(det(regularization)) + data = vcat(get_outputs(io_pairs_opt)[(n_train + 1):end], 0.0, min_complexity) + ekiobj = EKP.EnsembleKalmanProcess(initial_params, data, Γ, Inversion(), rng = rng) + err = zeros(n_iteration) + + # [4.] optimize with EKP + for i in 1:n_iteration + + #get parameters: + lvec = transform_unconstrained_to_constrained(prior, get_u_final(ekiobj)) + g_ens, _ = calculate_ensemble_mean_and_coeffnorm( + srfi, + rng, + lvec, + regularization, + n_features_opt, + n_train, + n_test, + batch_sizes, + io_pairs_opt, + ) + inflation = optimizer_options["inflation"] + if inflation > 0 + EKP.update_ensemble!(ekiobj, g_ens, additive_inflation = true, use_prior_cov = true, s = inflation) # small regularizing inflation + else + EKP.update_ensemble!(ekiobj, g_ens) # small regularizing inflation + end + err[i] = get_error(ekiobj)[end] #mean((params_true - mean(params_i,dims=2)).^2) + + end + + # [5.] extract optimal hyperparameters + hp_optimal = get_ϕ_mean_final(prior, ekiobj)[:] + println("RFM $i, learnt hyperparameters: $(hp_optimal)") + + io_pairs_i = PairedDataContainer(input_values, reshape(output_values[i, :], 1, size(output_values, 2))) + # Now, fit new RF model with the optimized hyperparameters + rfm_i = RFM_from_hyperparameters(srfi, rng, hp_optimal, regularization, n_features, batch_sizes, input_dim) + fitted_features_i = RF.Methods.fit(rfm_i, io_pairs_i, decomposition_type = "svd") #fit features + + push!(rfms, rfm_i) + push!(fitted_features, fitted_features_i) + end + +end + + +function optimize_hyperparameters!(srfi::ScalarRandomFeatureInterface, args...; kwargs...) + println("Random Features already trained. continuing...") +end + + +function predict(srfi::ScalarRandomFeatureInterface, new_inputs::AbstractMatrix{FT}) where {FT <: AbstractFloat} + M = length(get_rfms(srfi)) + N_samples = size(new_inputs, 2) + # Predicts columns of inputs: input_dim × N_samples + μ = zeros(M, N_samples) + σ2 = zeros(M, N_samples) + for i in 1:M + μ[i, :], σ2[i, :] = + RF.Methods.predict(get_rfms(srfi)[i], get_fitted_features(srfi)[i], DataContainer(new_inputs)) + end + + # add the noise contribution from the regularization + σ2[:, :] = σ2[:, :] .+ 1.0 + + return μ, σ2 +end diff --git a/src/Utilities.jl b/src/Utilities.jl index eb354842e..e9ccd36d9 100644 --- a/src/Utilities.jl +++ b/src/Utilities.jl @@ -14,7 +14,7 @@ export get_training_points export get_obs_sample export orig2zscore export zscore2orig - +export posdef_correct """ $(DocStringExtensions.TYPEDSIGNATURES) @@ -22,13 +22,20 @@ Extract the training points needed to train the Gaussian process regression. - `ekp` - EnsembleKalmanProcess holding the parameters and the data that were produced during the Ensemble Kalman (EK) process. -- `N_train_iterations` - Number of EK layers/iterations to train on. +- `train_iterations` - Number (or indices) EK layers/iterations to train on. """ -function get_training_points(ekp::EnsembleKalmanProcess{FT, IT, P}, N_train_iterations::IT) where {FT, IT, P} - - # Note u[end] does not have an equivalent g - iter_range = (get_N_iterations(ekp) - N_train_iterations + 1):get_N_iterations(ekp) +function get_training_points( + ekp::EnsembleKalmanProcess{FT, IT, P}, + train_iterations::Union{IT, AbstractVector{IT}}, +) where {FT, IT, P} + + if !isa(train_iterations, AbstractVector) + # Note u[end] does not have an equivalent g + iter_range = (get_N_iterations(ekp) - train_iterations + 1):get_N_iterations(ekp) + else + iter_range = train_iterations + end u_tp = [] g_tp = [] @@ -114,4 +121,18 @@ function zscore2orig(Z::AbstractMatrix{FT}, mean::AbstractVector{FT}, std::Abstr end +""" + posdef_correct(mat::AbstractMatrix; tol::Real=1e8*eps()) + +Makes square matrix `mat` positive definite, by symmetrizing and bounding the minimum eigenvalue below by `tol` +""" +function posdef_correct(mat::AbstractMatrix; tol::Real = 1e8 * eps()) + out = 0.5 * (mat + permutedims(mat, (2, 1))) #symmetrize + out += (abs(minimum(eigvals(out))) + tol) * I #add to diag + return out +end + + + + end # module diff --git a/src/VectorRandomFeature.jl b/src/VectorRandomFeature.jl new file mode 100644 index 000000000..7ade5c301 --- /dev/null +++ b/src/VectorRandomFeature.jl @@ -0,0 +1,382 @@ +using ProgressBars + +export VectorRandomFeatureInterface +export get_rfms, + get_fitted_features, + get_batch_sizes, + get_n_features, + get_input_dim, + get_output_dim, + get_rng, + get_diagonalize_input, + get_diagonalize_output, + get_optimizer_options + +struct VectorRandomFeatureInterface <: RandomFeatureInterface + "A vector of `RandomFeatureMethod`s, contains the feature structure, batch-sizes and regularization" + rfms::Vector{RF.Methods.RandomFeatureMethod} + "vector of `Fit`s, containing the matrix decomposition and coefficients of RF when fitted to data" + fitted_features::Vector{RF.Methods.Fit} + "batch sizes" + batch_sizes::Union{Dict, Nothing} + "number of features" + n_features::Union{Int, Nothing} + "input dimension" + input_dim::Int + "output_dimension" + output_dim::Int + "rng" + rng::AbstractRNG + "regularization" + regularization::Vector{Union{Matrix, UniformScaling, Diagonal}} + "Assume inputs are decorrelated for ML" + diagonalize_input::Bool + "Assume outputs are decorrelated for ML. Note: emulator(..., decorrelate=true) by default" + diagonalize_output::Bool + "dictionary of options for hyperparameter optimizer" + optimizer_options::Dict +end + +get_rfms(vrfi::VectorRandomFeatureInterface) = vrfi.rfms +get_fitted_features(vrfi::VectorRandomFeatureInterface) = vrfi.fitted_features +get_batch_sizes(vrfi::VectorRandomFeatureInterface) = vrfi.batch_sizes +get_n_features(vrfi::VectorRandomFeatureInterface) = vrfi.n_features +get_input_dim(vrfi::VectorRandomFeatureInterface) = vrfi.input_dim +get_output_dim(vrfi::VectorRandomFeatureInterface) = vrfi.output_dim +get_rng(vrfi::VectorRandomFeatureInterface) = vrfi.rng +get_regularization(vrfi::VectorRandomFeatureInterface) = vrfi.regularization +get_diagonalize_input(vrfi::VectorRandomFeatureInterface) = vrfi.diagonalize_input +get_diagonalize_output(vrfi::VectorRandomFeatureInterface) = vrfi.diagonalize_output +get_optimizer_options(vrfi::VectorRandomFeatureInterface) = vrfi.optimizer_options + + + +function VectorRandomFeatureInterface( + n_features::Int, + input_dim::Int, + output_dim::Int; + diagonalize_input::Bool = false, + diagonalize_output::Bool = false, + batch_sizes::Union{Dict, Nothing} = nothing, + rng::AbstractRNG = Random.GLOBAL_RNG, + optimizer_options::Union{Dict, Nothing} = nothing, +) + + # Initialize vector for RF model + rfms = Vector{RF.Methods.RandomFeatureMethod}(undef, 0) + fitted_features = Vector{RF.Methods.Fit}(undef, 0) + regularization = Vector{Union{Matrix, UniformScaling, Nothing}}(undef, 0) + + #Optimization Defaults + n_hp = calculate_n_hyperparameters(input_dim, output_dim, diagonalize_input, diagonalize_output) + prior = constrained_gaussian("MatrixNormal_parameters", 1.0, 1.0, 0.0, Inf, repeats = n_hp) + + optimizer_opts = Dict( + "prior" => prior, #the hyperparameter_prior + "n_ensemble" => ndims(prior) + 1, #number of ensemble + "n_iteration" => 5, # number of eki iterations + "tikhonov" => 0, # tikhonov regularization parameter if >0 + "inflation" => 1e-4, # additive inflation ∈ [0,1] with 0 being no inflation + ) + + if !isnothing(optimizer_options) + for key in keys(optimizer_options) + optimizer_opts[key] = optimizer_options[key] + end + end + + println("hyperparameter optimization with EKI configured with:", optimizer_opts) + + return VectorRandomFeatureInterface( + rfms, + fitted_features, + batch_sizes, + n_features, + input_dim, + output_dim, + rng, + regularization, + diagonalize_input, + diagonalize_output, + optimizer_opts, + ) +end + +function RFM_from_hyperparameters( + vrfi::VectorRandomFeatureInterface, + rng::AbstractRNG, + l::Union{Real, AbstractVecOrMat}, + regularization::Union{Matrix, UniformScaling, Diagonal}, + n_features::Int, + batch_sizes::Union{Dict, Nothing}, + input_dim::Int, + output_dim::Int; +) + diagonalize_input = get_diagonalize_input(vrfi) + diagonalize_output = get_diagonalize_output(vrfi) + + M = zeros(input_dim, output_dim) # n x p mean + U, V = hyperparameters_from_flat(l, input_dim, output_dim, diagonalize_input, diagonalize_output) + + if !isposdef(U) + println("U not posdef - correcting") + U = posdef_correct(U) + end + if !isposdef(V) + println("V not posdef - correcting") + V = posdef_correct(V) + end + + dist = MatrixNormal(M, U, V) + pd = ParameterDistribution( + Dict( + "distribution" => Parameterized(dist), + "constraint" => repeat([no_constraint()], input_dim * output_dim), + "name" => "xi", + ), + ) + feature_sampler = RF.Samplers.FeatureSampler(pd, output_dim, rng = copy(rng)) + vff = RF.Features.VectorFourierFeature(n_features, output_dim, feature_sampler) + + if isnothing(batch_sizes) + return RF.Methods.RandomFeatureMethod(vff, regularization = regularization) + else + return RF.Methods.RandomFeatureMethod(vff, regularization = regularization, batch_sizes = batch_sizes) + end +end + + + + +function build_models!( + vrfi::VectorRandomFeatureInterface, + input_output_pairs::PairedDataContainer{FT}; + regularization_matrix::Union{Matrix, Nothing} = nothing, +) where {FT <: AbstractFloat} + + # get inputs and outputs + input_values = get_inputs(input_output_pairs) + output_values = get_outputs(input_output_pairs) + n_rfms, n_data = size(output_values) + + input_dim = size(input_values, 1) + output_dim = size(output_values, 1) + # there are less hyperparameters when we diagonalize + diagonalize_input = get_diagonalize_input(vrfi) + diagonalize_output = get_diagonalize_output(vrfi) + n_hp = calculate_n_hyperparameters(input_dim, output_dim, diagonalize_input, diagonalize_output) + + rfms = get_rfms(vrfi) + fitted_features = get_fitted_features(vrfi) + n_features = get_n_features(vrfi) + batch_sizes = get_batch_sizes(vrfi) + optimizer_options = get_optimizer_options(vrfi) + prior = optimizer_options["prior"] + rng = get_rng(vrfi) + + if ndims(prior) > n_hp + # comes from having a truncated output_dimension + # TODO not really a truncation here, resetting to default + @info "truncating hyperparameter space to default in first $n_hp dimensions, due to truncation of output space" + prior = constrained_gaussian("cholesky_factors", 1.0, 1.0, 0.0, Inf, repeats = n_hp) + + end + + # Optimize feature cholesky factors with EKP + # [1.] Split data into test/train (e.g. 80/20) + n_train = Int(floor(0.8 * n_data)) # 20% split + n_test = n_data - n_train + if diagonalize_output + n_features_opt = min(n_features, Int(floor(4 * n_train))) #we take a specified number of features for optimization. MAGIC NUMBER + else + # note the n_features_opt should NOT exceed output_dim * n_train or the regularization is replaced with a diagonal approx. + n_features_opt = max(min(n_features, Int(floor(4 * sqrt(n_train * output_dim)))), Int(floor(1.9 * n_train)))#we take a specified number of features for optimization. MAGIC NUMBER + end + println( + "hyperparameter learning using: $n_train training points, $n_test validation points and $n_features_opt features", + ) + + # regularization_matrix = nothing when we use scaled SVD to decorrelate the space, + # in this setting, noise = I + if regularization_matrix === nothing + regularization = I + else + reg_mat = regularization_matrix + if !isposdef(regularization_matrix) + reg_mat = posdef_correct(regularization_matrix) + println("RF regularization matrix is not positive definite, correcting") + + else + # think of the regularization_matrix as the observational noise covariance, or a related quantity + regularization = exp((1 / output_dim) * sum(log.(eigvals(reg_mat)))) * I #i.e. det(M)^{1/output_dim} I + + #regularization = reg_mat #using the full p.d. tikhonov exp. EXPENSIVE, and challenge get complexity terms + end + end + + + idx_shuffle = randperm(rng, n_data) + + io_pairs_opt = PairedDataContainer( + input_values[:, idx_shuffle], + reshape(output_values[:, idx_shuffle], :, size(output_values, 2)), + ) + + # [2.] Estimate covariance at mean value + μ_hp = transform_unconstrained_to_constrained(prior, mean(prior)) + if diagonalize_output + n_sample_diag = 20 # diagonal case, less samples needed, MAGIC NUMBER + n_samples = max(2, n_sample_diag) + else + n_samples = (n_test * output_dim + 2) + end + internal_Γ, approx_σ2 = estimate_mean_and_coeffnorm_covariance( + vrfi, + rng, + μ_hp, # take mean values + regularization, + n_features_opt, + n_train, + n_test, + batch_sizes, + io_pairs_opt, + n_samples, + ) + + tikhonov_opt_val = optimizer_options["tikhonov"] + if tikhonov_opt_val == 0 + # Build the covariance + Γ = internal_Γ + Γ[1:(n_test * output_dim), 1:(n_test * output_dim)] += approx_σ2 + Γ[(n_test * output_dim + 1):end, (n_test * output_dim + 1):end] += I + + #in diag case we have data logdet = λ^m, in non diag case we have logdet(Λ^) to match the different reg matrices. + min_complexity = + isa(regularization, UniformScaling) ? n_features_opt * log(regularization.λ) : + n_features_opt / output_dim * 2 * sum(log.(diag(cholesky(regularization).L))) + + data = vcat(reshape(get_outputs(io_pairs_opt)[:, (n_train + 1):end], :, 1), 0.0, min_complexity) #flatten data + + elseif tikhonov_opt_val > 0 + # augment the state to add tikhonov + 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)] += approx_σ2 + Γ[(n_test * output_dim + 1):outsize, (n_test * output_dim + 1):outsize] += I + + Γ[(outsize + 1):end, (outsize + 1):end] = tikhonov_opt_val .* cov(prior) + + #TODO the min complexity here is not the correct object in the non-diagonal case + min_complexity = + isa(regularization, UniformScaling) ? n_features_opt * log(regularization.λ) : + n_features_opt / output_dim * 2 * sum(log.(diag(cholesky(regularization).L))) + + data = vcat( + reshape(get_outputs(io_pairs_opt)[:, (n_train + 1):end], :, 1), + 0.0, + min_complexity, + zeros(size(Γ, 1) - outsize, 1), + ) #flatten data with additional zeros + else + throw( + ArgumentError( + "Tikhonov parameter must be non-negative, instead received tikhonov_opt_val=$tikhonov_opt_val", + ), + ) + end + + # [3.] set up EKP optimization + n_ensemble = optimizer_options["n_ensemble"] # minimal ensemble size n_hp, + n_iteration = optimizer_options["n_iteration"] + + if !isposdef(Γ) + Γ = posdef_correct(Γ) + end + + initial_params = construct_initial_ensemble(rng, prior, n_ensemble) + + ekiobj = EKP.EnsembleKalmanProcess(initial_params, data[:], Γ, Inversion(), rng = rng) + err = zeros(n_iteration) + + # [4.] optimize with EKP + for i in 1:n_iteration + println("iteration $i") + #get parameters: + lvec = get_ϕ_final(prior, ekiobj) + g_ens, _ = calculate_ensemble_mean_and_coeffnorm( + vrfi, + rng, + lvec, + regularization, + n_features_opt, + n_train, + n_test, + batch_sizes, + io_pairs_opt, + ) + if tikhonov_opt_val > 0 + # augment with the computational parameters (u not ϕ) + uvecormat = get_u_final(ekiobj) + if isa(uvecormat, AbstractVector) + umat = reshape(uvecormat, 1, :) + else + umat = uvecormat + end + + g_ens = vcat(g_ens, umat) + end + inflation = optimizer_options["inflation"] + if inflation > 0 + EKP.update_ensemble!(ekiobj, g_ens, additive_inflation = true, use_prior_cov = true, s = inflation) # small regularizing inflation + else + EKP.update_ensemble!(ekiobj, g_ens) # small regularizing inflation + end + err[i] = get_error(ekiobj)[end] #mean((params_true - mean(params_i,dims=2)).^2) + + end + + # [5.] extract optimal hyperparameters + hp_optimal = get_ϕ_mean_final(prior, ekiobj)[:] + println("Vector RFM, learnt hyperparameters: $(hp_optimal)") + # Now, fit new RF model with the optimized hyperparameters + + rfm_tmp = + RFM_from_hyperparameters(vrfi, rng, hp_optimal, regularization, n_features, batch_sizes, input_dim, output_dim) + fitted_features_tmp = fit(rfm_tmp, input_output_pairs, decomposition_type = "svd") + + push!(rfms, rfm_tmp) + push!(fitted_features, fitted_features_tmp) + push!(get_regularization(vrfi), regularization) + +end + + +function optimize_hyperparameters!(vrfi::VectorRandomFeatureInterface, args...; kwargs...) + println("Random Features already trained. continuing...") +end + + +function predict(vrfi::VectorRandomFeatureInterface, new_inputs::AbstractMatrix{FT}) where {FT <: AbstractFloat} + input_dim = get_input_dim(vrfi) + output_dim = get_output_dim(vrfi) + rfm = get_rfms(vrfi)[1] + ff = get_fitted_features(vrfi)[1] + N_samples = size(new_inputs, 2) + # Predicts columns of inputs: input_dim × N_samples + μ, σ2 = RF.Methods.predict(rfm, ff, DataContainer(new_inputs)) + # sizes (output_dim x n_test), (output_dim x output_dim x n_test) + + # add the noise contribution from the regularization + # note this is because we are predicting the data here, not the latent function. + lambda = get_regularization(vrfi)[1] + for i in 1:N_samples + σ2[:, :, i] = 0.5 * (σ2[:, :, i] + permutedims(σ2[:, :, i], (2, 1))) + lambda + + if !isposdef(σ2[:, :, i]) + σ2[:, :, i] = posdef_correct(σ2[:, :, i]) + end + end + + return μ, σ2 +end diff --git a/test/GaussianProcess/runtests.jl b/test/GaussianProcess/runtests.jl index 53624cc56..5935f834a 100644 --- a/test/GaussianProcess/runtests.jl +++ b/test/GaussianProcess/runtests.jl @@ -214,7 +214,7 @@ using CalibrateEmulateSample.DataContainers @test μ4_noise_from_Σ[:, 3] ≈ [0.0, 0.0] atol = tol_mu @test μ4_noise_from_Σ[:, 4] ≈ [0.0, -2.0] atol = tol_mu - # check match between the means and variances (should be similar at least + # check match between the variances (should be similar at least) @test all(isapprox.(σ4²_noise_from_Σ, σ4²_noise_learnt, rtol = 2 * tol_mu)) diff --git a/test/RandomFeature/runtests.jl b/test/RandomFeature/runtests.jl new file mode 100644 index 000000000..b750d23d9 --- /dev/null +++ b/test/RandomFeature/runtests.jl @@ -0,0 +1,126 @@ +using Test +using Random +using RandomFeatures + +using CalibrateEmulateSample.Emulators +using CalibrateEmulateSample.DataContainers +using CalibrateEmulateSample.ParameterDistributions +using RandomFeatures + +seed = 10101010 + +@testset "RandomFeatures" begin + + @testset "ScalarRandomFeatureInterface" begin + + rng = Random.MersenneTwister(seed) + + input_dim = 2 + n_features = 200 + batch_sizes = Dict("train" => 100, "test" => 100, "feature" => 100) + diagonalize_input = false + #build interface + n_input_hp = Int(input_dim * (input_dim + 1) / 2) + 2 # see calculate_n_hyperparameters for details + n_output_hp = 1 + prior = constrained_gaussian("MvNormal_parameters", 1.0, 1.0, 0.0, Inf, repeats = n_input_hp + n_output_hp) + optimizer_options = + Dict("prior" => prior, "n_ensemble" => ndims(prior) + 1, "n_iteration" => 5, "inflation" => 1e-4) + + srfi = ScalarRandomFeatureInterface( + n_features, + input_dim, + batch_sizes = batch_sizes, + rng = rng, + diagonalize_input = diagonalize_input, + optimizer_options = optimizer_options, + ) + + @test isa(get_rfms(srfi), Vector{RandomFeatures.Methods.RandomFeatureMethod}) + @test isa(get_fitted_features(srfi), Vector{RandomFeatures.Methods.Fit}) + @test get_batch_sizes(srfi) == batch_sizes + @test get_n_features(srfi) == n_features + @test get_input_dim(srfi) == input_dim + @test get_rng(srfi) == rng + @test get_diagonalize_input(srfi) == diagonalize_input + @test get_optimizer_options(srfi) == optimizer_options + + # check defaults + srfi2 = ScalarRandomFeatureInterface(n_features, input_dim) + @test get_batch_sizes(srfi2) === nothing + @test get_rng(srfi2) == Random.GLOBAL_RNG + @test get_diagonalize_input(srfi2) == false + @test get_optimizer_options(srfi2) == optimizer_options # we just set the defaults above + + end + + @testset "VectorRandomFeatureInterface" begin + + rng = Random.MersenneTwister(seed) + + input_dim = 2 + output_dim = 3 + n_features = 200 + batch_sizes = Dict("train" => 100, "test" => 100, "feature" => 100) + diagonalize_input = false + diagonalize_output = false + + n_input_hp = Int(input_dim * (input_dim + 1) / 2) + 2 + n_output_hp = Int(output_dim * (output_dim + 1) / 2) + 2 + prior = constrained_gaussian("MatrixNormal_parameters", 1.0, 1.0, 0.0, Inf, repeats = n_input_hp + n_output_hp) + + optimizer_options = Dict( + "prior" => prior, + "n_ensemble" => ndims(prior) + 1, + "n_iteration" => 5, + "tikhonov" => 0, + "inflation" => 1e-4, + ) + + #build interfaces + vrfi = VectorRandomFeatureInterface( + n_features, + input_dim, + output_dim, + rng = rng, + batch_sizes = batch_sizes, + diagonalize_input = diagonalize_input, + diagonalize_output = diagonalize_output, + optimizer_options = optimizer_options, + ) + + @test isa(get_rfms(vrfi), Vector{RandomFeatures.Methods.RandomFeatureMethod}) + @test isa(get_fitted_features(vrfi), Vector{RandomFeatures.Methods.Fit}) + @test get_batch_sizes(vrfi) == batch_sizes + @test get_n_features(vrfi) == n_features + @test get_input_dim(vrfi) == input_dim + @test get_output_dim(vrfi) == output_dim + @test get_rng(vrfi) == rng + @test get_diagonalize_input(vrfi) == diagonalize_input + @test get_diagonalize_output(vrfi) == diagonalize_output + @test get_optimizer_options(vrfi) == optimizer_options + + + #check defaults + vrfi2 = VectorRandomFeatureInterface(n_features, input_dim, output_dim) + + @test get_batch_sizes(vrfi2) === nothing + @test get_rng(vrfi2) == Random.GLOBAL_RNG + @test get_diagonalize_input(vrfi2) == false + @test get_diagonalize_output(vrfi2) == false + @test get_optimizer_options(vrfi2) == optimizer_options # we just set the defaults above + + end + + + #= + # Training data + output_dim = 3 + n = 20 # number of training points + x = reshape(2.0 * π * rand(n), input_dim, n) # predictors/features: 2 x n + y = reshape(sin.(x) + 0.05 * randn(rng, n)', output_dim, n) # predictands/targets: 3 x n + iopairs = PairedDataContainer(x, y, data_are_columns = true) + =# + + + +end diff --git a/test/Utilities/runtests.jl b/test/Utilities/runtests.jl index bec9283ef..54e1ad7c9 100644 --- a/test/Utilities/runtests.jl +++ b/test/Utilities/runtests.jl @@ -60,4 +60,17 @@ using CalibrateEmulateSample.DataContainers @test get_inputs(training_points) ≈ initial_ensemble @test get_outputs(training_points) ≈ g_ens + #positive definiteness + mat = reshape(collect(-3:(-3 + 10^2 - 1)), 10, 10) + tol = 1e12 * eps() + @test !isposdef(mat) + + pdmat = posdef_correct(mat) + @test isposdef(pdmat) + @test minimum(eigvals(pdmat)) >= (1 - 1e-4) * 1e8 * eps() #eigvals is approximate so need a bit of give here + + pdmat2 = posdef_correct(mat, tol = tol) + @test isposdef(pdmat2) + @test minimum(eigvals(pdmat2)) >= (1 - 1e-4) * tol + end diff --git a/test/runtests.jl b/test/runtests.jl index c48b2aa9d..eeda04c7c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -24,7 +24,7 @@ end end end - for submodule in ["Emulator", "GaussianProcess", "MarkovChainMonteCarlo", "Utilities"] + for submodule in ["Emulator", "GaussianProcess", "RandomFeature", "MarkovChainMonteCarlo", "Utilities"] if all_tests || has_submodule(submodule) || "CalibrateEmulateSample" in ARGS include_test(submodule) end