perfect model ice nuc calibs + tests updated

papers/ice_nucleation_2024/calibration.jl

@@ -11,15 +11,17 @@ import CloudMicrophysics.Parameters as CMP
 import Thermodynamics as TD
 import StatsBase as SB
 
-# format: off
+#! format: off
 # definition of the ODE problem for parcel model
 include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
 include(joinpath(pkgdir(CM), "papers", "ice_nucleation_2024", "calibration_setup.jl"))
 
 # Define model which wraps around parcel and overwrites calibrated parameters
 function run_model(p, coefficients, FT, IC, end_sim; calibration = false)
     # grabbing parameters
-    ABDINM_m_calibrated, ABDINM_c_calibrated, ABIFM_m_calibrated, ABIFM_c_calibrated, ABHOM_m_calibrated, ABHOM_c_calibrated = coefficients
+    ABDINM_m_calibrated, ABDINM_c_calibrated,
+    ABIFM_m_calibrated, ABIFM_c_calibrated,
+    ABHOM_m_calibrated, ABHOM_c_calibrated = coefficients
     (; prescribed_thermodynamics, t_profile, T_profile, P_profile) = p
     (; const_dt, w, t_max, ips) = p
     (; aerosol_act, aerosol, r_nuc) = p
@@ -103,7 +105,7 @@ function run_model(p, coefficients, FT, IC, end_sim; calibration = false)
     # solve ODE
     local sol = run_parcel(IC, FT(0), t_max, params)
     if calibration == true
-        return sol[9, end - end_sim:end] ./ (IC[7] + IC[8] + IC[9])
+        return sol[9, (end - end_sim):end] ./ (IC[7] + IC[8] + IC[9])
     elseif calibration == false
         return sol
     end
@@ -197,7 +199,7 @@ function calibrate_J_parameters_EKI(FT, IN_mode, params, IC, y_truth, end_sim, 
     initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ensemble)
     EKI_obj = EKP.EnsembleKalmanProcess(
         initial_ensemble,
-        y_truth[end - end_sim:end],
+        y_truth[(end - end_sim):end],
         Γ,
         EKP.Inversion();
         rng = rng,
@@ -323,8 +325,11 @@ function ensemble_means(ϕ_n_values, N_iterations, N_ensemble)
             Distributions.mean(ϕ_n_values[iter][5, i] for i in 1:N_ensemble)
         ABHOM_c_mean[iter] =
             Distributions.mean(ϕ_n_values[iter][6, i] for i in 1:N_ensemble)
-
     end
 
-    return [ABDINM_m_mean, ABDINM_c_mean, ABIFM_m_mean, ABIFM_c_mean, ABHOM_m_mean, ABHOM_c_mean]
+    return [
+        ABDINM_m_mean, ABDINM_c_mean,
+        ABIFM_m_mean, ABIFM_c_mean,
+        ABHOM_m_mean, ABHOM_c_mean
+    ]
 end

papers/ice_nucleation_2024/calibration_setup.jl

@@ -135,17 +135,15 @@ function perf_model_IC(FT, IN_mode)
     return [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)]
 end
 
-function perf_model_pseudo_data(FT, IN_mode, params, IC)
+function perf_model_pseudo_data(FT, IN_mode, params, IC, end_sim)
     n_samples = 10
 
-    if IN_mode == "ABDINM"
-        coeff_true = [FT(27.551), FT(-2.2209)]
-    elseif IN_mode == "ABIFM"
-        coeff_true = [FT(54.58834), FT(-10.54758)]
-    elseif IN_mode == "ABHOM"
-        coeff_true = [FT(255.927125), FT(-68.553283)]
-    end
-    sol_ICNC = run_calibrated_model(FT, IN_mode, coeff_true, params, IC)[9, :]
+    coeff_true = [
+        FT(27.551), FT(-2.2209),        # ABDINM
+        FT(54.58834), FT(-10.54758),    # ABIFM
+        FT(255.927125), FT(-68.553283), # ABHOM
+    ]
+    sol_ICNC = run_model(params, coeff_true, FT, IC, end_sim)[9, :]
     G_truth = sol_ICNC ./ (IC[7] + IC[8] + IC[9])
     dim_output = length(G_truth)
 

papers/ice_nucleation_2024/perfect_calib.jl

@@ -11,12 +11,12 @@ IN_mode_list = ["ABDINM", "ABIFM", "ABHOM"]
 for IN_mode in IN_mode_list
     params = perf_model_params(FT, IN_mode)
     IC = perf_model_IC(FT, IN_mode)
+    end_sim = 25
 
-    pseudo_data = perf_model_pseudo_data(FT, IN_mode, params, IC)
+    pseudo_data = perf_model_pseudo_data(FT, IN_mode, params, IC, end_sim)
     Γ = pseudo_data[2]
     y_truth = pseudo_data[1]
     coeff_true = pseudo_data[3]
-    end_sim = 25
 
     output = calibrate_J_parameters_EKI(
         FT,
@@ -31,29 +31,68 @@ for IN_mode in IN_mode_list
 
     iterations = collect(1:size(output[2])[1])
     calibrated_parameters = output[1]
-    m_mean = []
-    m_mean = ensemble_means(
+    ϕ_n_values = output[2]
+    ABDINM_m_mean = []
+    ABDINM_c_mean = []
+    ABIFM_m_mean = []
+    ABIFM_c_mean = []
+    ABHOM_m_mean = []
+    ABHOM_c_mean = []
+
+    ABDINM_m_mean = ensemble_means(
         output[2],
-        size(output[2])[1],
-        size(output[2][1])[2],
+        size(ϕ_n_values)[1],
+        size(ϕ_n_values[1])[2],
     )[1]
-    c_mean = []
-    c_mean = ensemble_means(
+    ABDINM_c_mean = ensemble_means(
         output[2],
-        size(output[2])[1],
-        size(output[2][1])[2],
+        size(ϕ_n_values)[1],
+        size(ϕ_n_values[1])[2],
     )[2]
+    ABIFM_m_mean = ensemble_means(
+        output[2],
+        size(ϕ_n_values)[1],
+        size(ϕ_n_values[1])[2],
+    )[3]
+    ABIFM_c_mean = ensemble_means(
+        output[2],
+        size(ϕ_n_values)[1],
+        size(ϕ_n_values[1])[2],
+    )[4]
+    ABHOM_m_mean = ensemble_means(
+        output[2],
+        size(ϕ_n_values)[1],
+        size(ϕ_n_values[1])[2],
+    )[5]
+    ABHOM_c_mean = ensemble_means(
+        output[2],
+        size(ϕ_n_values)[1],
+        size(ϕ_n_values[1])[2],
+    )[6]
 
     # Plotting calibrated parameters per iteration
-    fig = MK.Figure(size = (800, 600))
-    ax1 = MK.Axis(fig[1, 1], ylabel = "m coefficient [-]", xlabel = "iteration number", title = IN_mode)
-    ax2 = MK.Axis(fig[2, 1], ylabel = "c coefficient [-]", xlabel = "iteration number")
-
-    MK.lines!(ax1, iterations, m_mean,                                 label = "ensemble mean", color = :orange)
-    MK.lines!(ax1, iterations, zeros(length(m_mean)) .+ coeff_true[1], label = "default value")
-    MK.lines!(ax2, iterations, c_mean,                                 label = "ensemble mean", color = :orange)
-    MK.lines!(ax2, iterations, zeros(length(c_mean)) .+ coeff_true[2], label = "default value")
+    fig = MK.Figure(size = (1100, 900))
+    ax1 = MK.Axis(fig[1, 1], ylabel = "ABDINM m coefficient [-]", title = IN_mode)
+    ax2 = MK.Axis(fig[1, 2], ylabel = "ABDINM c coefficient [-]", xlabel = "iteration #")
+    ax3 = MK.Axis(fig[2, 1], ylabel = "ABIFM m coefficient [-]", title = IN_mode)
+    ax4 = MK.Axis(fig[2, 2], ylabel = "ABIFM c coefficient [-]", xlabel = "iteration #")
+    ax5 = MK.Axis(fig[3, 1], ylabel = "ABHOM m coefficient [-]", title = IN_mode)
+    ax6 = MK.Axis(fig[3, 2], ylabel = "ABHOM c coefficient [-]", xlabel = "iteration #")
 
+
+    MK.lines!(ax1, iterations, ABDINM_m_mean,                                 label = "ensemble mean", color = :orange)
+    MK.lines!(ax1, iterations, zeros(length(ABDINM_m_mean)) .+ coeff_true[1], label = "default value")
+    MK.lines!(ax2, iterations, ABDINM_c_mean,                                 label = "ensemble mean", color = :orange)
+    MK.lines!(ax2, iterations, zeros(length(ABDINM_c_mean)) .+ coeff_true[2], label = "default value")
+    MK.lines!(ax3, iterations, ABIFM_m_mean,                                 label = "ensemble mean", color = :orange)
+    MK.lines!(ax3, iterations, zeros(length(ABIFM_m_mean)) .+ coeff_true[3], label = "default value")
+    MK.lines!(ax4, iterations, ABIFM_c_mean,                                 label = "ensemble mean", color = :orange)
+    MK.lines!(ax4, iterations, zeros(length(ABIFM_c_mean)) .+ coeff_true[4], label = "default value")
+    MK.lines!(ax5, iterations, ABHOM_m_mean,                                 label = "ensemble mean", color = :orange)
+    MK.lines!(ax5, iterations, zeros(length(ABHOM_m_mean)) .+ coeff_true[5], label = "default value")
+    MK.lines!(ax6, iterations, ABHOM_c_mean,                                 label = "ensemble mean", color = :orange)
+    MK.lines!(ax6, iterations, zeros(length(ABHOM_c_mean)) .+ coeff_true[6], label = "default value")
+
     MK.axislegend(ax1, framevisible = true, labelsize = 12, position = :rc)
     MK.axislegend(ax2, framevisible = true, labelsize = 12)
 
@@ -72,11 +111,11 @@ for IN_mode in IN_mode_list
     fig2 = MK.Figure(size = (800, 600))
     ax3 = MK.Axis(fig2[1, 1], ylabel = "Frozen Fraction [-]", xlabel = "time [s]", title = IN_mode)
 
-    soln = run_calibrated_model(FT, IN_mode, calibrated_parameters, params, IC)
-    soln_dflt = run_calibrated_model(FT, IN_mode, coeff_true, params, IC)
+    soln = run_model(params, calibrated_parameters, FT, IC, end_sim)
+    soln_dflt = run_model(params, coeff_true, FT, IC, end_sim)
 
     MK.lines!(ax3, soln.t, soln[9,:]./ (IC[7] + IC[8] + IC[9]), label = "calibrated")
-    MK.lines!(ax3, soln_dflt.t, soln_dflt[9,:]./ (IC[7] + IC[8] + IC[9]), label = "default", color = :orange)
+    MK.lines!(ax3, soln_dflt.t, soln_dflt[9,:]./ (IC[7] + IC[8] + IC[9]), label = "default", color = :orange, linestyle = :dot)
     plot_name = "perfect_calibration_ICNC_$mode_label.svg"
     MK.save(plot_name, fig2)
 end

test/ice_nucleation_calibration.jl

@@ -8,12 +8,12 @@ include(joinpath(pkgdir(CM), "papers", "ice_nucleation_2024", "calibration.jl"))
 function test_J_calibration(FT, IN_mode)
     params = perf_model_params(FT, IN_mode)
     IC = perf_model_IC(FT, IN_mode)
+    end_sim = 25
 
-    pseudo_data = perf_model_pseudo_data(FT, IN_mode, params, IC)
+    pseudo_data = perf_model_pseudo_data(FT, IN_mode, params, IC, end_sim)
     Γ = pseudo_data[2]
     y_truth = pseudo_data[1]
     coeff_true = pseudo_data[3]
-    end_sim = 25
 
     EKI_output = calibrate_J_parameters_EKI(
         FT,
@@ -39,10 +39,10 @@ function test_J_calibration(FT, IN_mode)
     EKI_calibrated_parameters = EKI_output[1]
     UKI_calibrated_parameters = UKI_output[1]
     EKI_calibrated_soln =
-        run_calibrated_model(FT, IN_mode, EKI_calibrated_parameters, params, IC)
+        run_model(params, EKI_calibrated_parameters, FT, IC, end_sim)
     UKI_calibrated_soln =
-        run_calibrated_model(FT, IN_mode, UKI_calibrated_parameters, params, IC)
-    true_soln = run_calibrated_model(FT, IN_mode, coeff_true, params, IC)
+        run_model(params, UKI_calibrated_parameters, FT, IC, end_sim)
+    true_soln = run_model(params, coeff_true, FT, IC, end_sim)
 
     TT.@testset "EKI Perfect Model Calibrations on AIDA" begin
         # test that coeffs are close to "true" values and that end ICNC are similar
@@ -52,13 +52,13 @@ function test_J_calibration(FT, IN_mode)
             TT.@test EKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol =
                 FT(0.3)
         elseif IN_mode == "ABIFM"
-            TT.@test EKI_calibrated_parameters[1] ≈ coeff_true[1] rtol = FT(0.3)
-            TT.@test EKI_calibrated_parameters[2] ≈ coeff_true[2] rtol = FT(0.3)
+            TT.@test EKI_calibrated_parameters[3] ≈ coeff_true[3] rtol = FT(0.3)
+            TT.@test EKI_calibrated_parameters[4] ≈ coeff_true[4] rtol = FT(0.3)
             TT.@test EKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol =
                 FT(0.3)
         elseif IN_mode == "ABHOM"
-            TT.@test EKI_calibrated_parameters[1] ≈ coeff_true[1] rtol = FT(0.3)
-            TT.@test EKI_calibrated_parameters[2] ≈ coeff_true[2] rtol = FT(0.3)
+            TT.@test EKI_calibrated_parameters[5] ≈ coeff_true[5] rtol = FT(0.3)
+            TT.@test EKI_calibrated_parameters[6] ≈ coeff_true[6] rtol = FT(0.3)
             TT.@test EKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol =
                 FT(0.3)
         end
@@ -72,13 +72,13 @@ function test_J_calibration(FT, IN_mode)
             TT.@test UKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol =
                 FT(0.3)
         elseif IN_mode == "ABIFM"
-            TT.@test UKI_calibrated_parameters[1] ≈ coeff_true[1] rtol = FT(0.3)
-            TT.@test UKI_calibrated_parameters[2] ≈ coeff_true[2] rtol = FT(0.3)
+            TT.@test UKI_calibrated_parameters[3] ≈ coeff_true[3] rtol = FT(0.3)
+            TT.@test UKI_calibrated_parameters[4] ≈ coeff_true[4] rtol = FT(0.3)
             TT.@test UKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol =
                 FT(0.3)
         elseif IN_mode == "ABHOM"
-            TT.@test UKI_calibrated_parameters[1] ≈ coeff_true[1] rtol = FT(0.3)
-            TT.@test UKI_calibrated_parameters[2] ≈ coeff_true[2] rtol = FT(0.3)
+            TT.@test UKI_calibrated_parameters[5] ≈ coeff_true[5] rtol = FT(0.3)
+            TT.@test UKI_calibrated_parameters[6] ≈ coeff_true[6] rtol = FT(0.3)
             TT.@test UKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol =
                 FT(0.3)
         end