lorenz 2d statsplot

CliMA · Feb 25, 2023 · f533dda · f533dda
1 parent 5af6607
commit f533dda
Show file tree

Hide file tree

Showing 3 changed files with 43 additions and 150 deletions.
diff --git a/examples/Lorenz/Lorenz_example.jl b/examples/Lorenz/Lorenz_example.jl
@@ -287,7 +287,7 @@ end
 
 # Emulate-sample settings
 standardize = true
-retained_svd_frac = 0.95
+retained_svd_frac = 1.0
 
 gppackage = Emulators.GPJL()
 pred_type = Emulators.YType()
@@ -378,51 +378,17 @@ end
 n_params = length(true_values)
 y_folder = "qoi_" * string(lorenz_settings.stats_type) * "_"
 save_directory = figure_save_directory * y_folder
-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")
-    end
-
-    if log_normal == true
-        xs = log.(xs)
-        xbounds = log.(xbounds)
-    end
-
-    label = "true " * param
-    posterior_distributions = get_distribution(posterior)
-
-    histogram(
-        posterior_distributions[param_names[idx]][1, :],
-        bins = 100,
-        normed = true,
-        fill = :slategray,
-        lab = "posterior",
-    )
-    prior_plot = get_distribution(mcmc.prior)
-    # This requires StatsPlots
-    plot!(xs, prior_plot[param_names[idx]], w = 2.6, color = :blue, lab = "prior")
-    #plot!(xs, mcmc.prior[idx].dist, w=2.6, color=:blue, lab="prior")
-    plot!([true_values[idx]], seriestype = "vline", w = 2.6, lab = label)
-    plot!(xlims = xbounds)
-
-    title!(param)
-
-    figpath = joinpath(figure_save_directory, "posterior_$(param)_T_$(T)_w_$(ω_true).png")
-    savefig(figpath)
-    linkfig(figpath)
-end
 
+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,[true_values[1]], seriestype = "vline", w = 2.6, lab = param_names[1])
+#plot!(p,[true_values[2]], seriestype = "hline", w = 2.6, lab = param_names[2])
+figpath = joinpath(figure_save_directory, "posterior_samples_2d_gp.png")
+savefig(figpath)
+
 # Save data
 @save joinpath(data_save_directory, "input_output_pairs.jld2") input_output_pairs
diff --git a/examples/Lorenz/Lorenz_example_RF.jl b/examples/Lorenz/Lorenz_example_RF.jl
@@ -286,19 +286,12 @@ end
 
 # Emulate-sample settings
 standardize = true
-retained_svd_frac = 0.95
+retained_svd_frac = 1.0
 
 # setup random features
 n_features = 300
 batch_sizes = Dict("train" => 100, "test" => 100, "feature" => 100)
-
-# hyperparameter prior
-μ_l = 5.0
-σ_l = 10.0
-
-prior_lengthscale = constrained_gaussian("lengthscale", μ_l, σ_l, 0.0, Inf, repeats = size(params_true, 1))
-srfi = ScalarRandomFeatureInterface(n_features, prior_lengthscale, batch_sizes = batch_sizes)
-
+srfi = ScalarRandomFeatureInterface(n_features, n_param, batch_sizes = batch_sizes)
 # Standardize the output data
 # Use median over all data since all data are the same type
 yt_norm = vcat(yt...)
@@ -379,51 +372,18 @@ end
 n_params = length(true_values)
 y_folder = "qoi_" * string(lorenz_settings.stats_type) * "_"
 save_directory = figure_save_directory * y_folder
-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")
-    end
 
-    if log_normal == true
-        xs = log.(xs)
-        xbounds = log.(xbounds)
-    end
+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)
 
-    label = "true " * param
-    posterior_distributions = get_distribution(posterior)
+gr(dpi = 300, size = (300,300))
+p = cornerplot(permutedims(constrained_posterior_samples, (2, 1)), label = param_names, compact = true)
+#plot!(p,[true_values[1]], seriestype = "vline", w = 2.6, lab = param_names[1])
+#plot!(p,[true_values[2]], seriestype = "hline", w = 2.6, lab = param_names[2])
+figpath = joinpath(figure_save_directory, "posterior_samples_2d_scalarRF.png")
+savefig(figpath)
 
-    histogram(
-        posterior_distributions[param_names[idx]][1, :],
-        bins = 100,
-        normed = true,
-        fill = :slategray,
-        lab = "posterior",
-    )
-    prior_plot = get_distribution(mcmc.prior)
-    # This requires StatsPlots
-    plot!(xs, prior_plot[param_names[idx]], w = 2.6, color = :blue, lab = "prior")
-    #plot!(xs, mcmc.prior[idx].dist, w=2.6, color=:blue, lab="prior")
-    plot!([true_values[idx]], seriestype = "vline", w = 2.6, lab = label)
-    plot!(xlims = xbounds)
-
-    title!(param)
-
-    figpath = joinpath(figure_save_directory, "posterior_$(param)_T_$(T)_w_$(ω_true).png")
-    savefig(figpath)
-    linkfig(figpath)
-end
 
 # Save data
 @save joinpath(data_save_directory, "input_output_pairs.jld2") input_output_pairs
diff --git a/examples/Lorenz/Lorenz_example_vector_RF.jl b/examples/Lorenz/Lorenz_example_vector_RF.jl
@@ -287,19 +287,18 @@ end
 
 # Emulate-sample settings
 standardize = true
-retained_svd_frac = 0.95
+retained_svd_frac = 1.0
 
 # setup random features
 n_features = 200
 batch_sizes = Dict("train" => 200, "test" => 200, "feature" => 200)
-
-# hyperparameter prior
-μ_l = 5.0
-σ_l = 10.0
-
-#prior_lengthscale = constrained_gaussian("lengthscale", μ_l, σ_l, 0.0, Inf, repeats = size(params_true, 1))
-vrfi = VectorRandomFeatureInterface(n_features, input_dim, output_dim, batch_sizes = batch_sizes)
-
+vrfi = VectorRandomFeatureInterface(
+    n_features,
+    input_dim,
+    output_dim,
+    batch_sizes = batch_sizes,
+    diagonalize_output=false, # assume outputs are decorrelated
+)
 
 # Standardize the output data
 # Use median over all data since all data are the same type
@@ -323,6 +322,7 @@ emulator = Emulator(
     standardize_outputs = standardize,
     standardize_outputs_factors = norm_factor,
     retained_svd_frac = retained_svd_frac,
+    decorrelate=true, # apply SVD to data? 
 )
 optimize_hyperparameters!(emulator)
 
@@ -381,51 +381,18 @@ end
 n_params = length(true_values)
 y_folder = "qoi_" * string(lorenz_settings.stats_type) * "_"
 save_directory = figure_save_directory * y_folder
-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")
-    end
 
-    if log_normal == true
-        xs = log.(xs)
-        xbounds = log.(xbounds)
-    end
-
-    label = "true " * param
-    posterior_distributions = get_distribution(posterior)
-
-    histogram(
-        posterior_distributions[param_names[idx]][1, :],
-        bins = 100,
-        normed = true,
-        fill = :slategray,
-        lab = "posterior",
-    )
-    prior_plot = get_distribution(mcmc.prior)
-    # This requires StatsPlots
-    plot!(xs, prior_plot[param_names[idx]], w = 2.6, color = :blue, lab = "prior")
-    #plot!(xs, mcmc.prior[idx].dist, w=2.6, color=:blue, lab="prior")
-    plot!([true_values[idx]], seriestype = "vline", w = 2.6, lab = label)
-    plot!(xlims = xbounds)
-
-    title!(param)
-
-    figpath = joinpath(figure_save_directory, "posterior_$(param)_T_$(T)_w_$(ω_true).png")
-    savefig(figpath)
-    linkfig(figpath)
-end
+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,[true_values[1]], seriestype = "vline", w = 2.6, lab = param_names[1])
+#plot!(p,[true_values[2]], seriestype = "hline", w = 2.6, lab = param_names[2])
+figpath = joinpath(figure_save_directory, "posterior_samples_2d.png")
+savefig(figpath)
+
 
 # Save data
 @save joinpath(data_save_directory, "input_output_pairs.jld2") input_output_pairs