WIP

experiments/calibration/minimal_working_example/minimal_working_example.jl

@@ -4,12 +4,12 @@
 # reach out to Alexis Renchon.
 
 # We calibrate a single site simulation, using EnsembleKalmanProcesses.jl.
-# We calibrate the parameter(s) g1, to better match latent heat flux (lhf) target data.
+# We calibrate the parameters g1 and g0, to better match latent heat flux (lhf) target data.
 
 # To perform this calibration we will need:
-# 1. a function returning our model lhf output given the parameter(s) we want to calibrate.
+# 1. a function returning our model lhf output given the parameters we want to calibrate.
 # 2. the "truth" target data, to calibrate on.
-# 3. the prior distribution of these parameter(s).
+# 3. the prior distribution of these parameters.
 
 # Let's go!
 
@@ -19,10 +19,12 @@ import EnsembleKalmanProcesses as EKP # to perform the calibration
 import Random # we need it for ...
 
 # a function returning our model lhf output given the params to calibrate
-function ClimaLand_lhf(g1) # Should this return a vector or a ::OS?
+function ClimaLand_lhf(params) # Should this return a vector or a ::OS?
+    g1, g0 = params
+
     sv = CLS.run_fluxnet(
     "US-MOz";
-    params = CLS.ozark_default_params(; conductance = CLS.conductance_ozark(; g1 = g1)),
+    params = CLS.ozark_default_params(; conductance = CLS.conductance_ozark(; g1 = g1, g0 = g0)),
     )[1];
 
     inputs = CLS.make_inputs_df("US-MOz")[1];
@@ -31,38 +33,41 @@ function ClimaLand_lhf(g1) # Should this return a vector or a ::OS?
 
     lhf = simulation_output.lhf
 
-    lhf_obs = EKP.Observation(lhf, EKP.I(length(lhf)), "lhf") # target data needs to be of type ::OS. Not sure about I!!
-    return lhf_obs
+    # lhf_obs = EKP.Observation(lhf, EKP.I(length(lhf)), "lhf") # target data needs to be of type ::OS. Not sure about I!!
+    return lhf
 end
 
 # "truth" target data to calibrate on
-lhf_target = ClimaLand_lhf(141.0) # our default for Ozark, in sqrt(Pa). This is equal to 4.46 sqrt(kPa).
+lhf_target = ClimaLand_lhf([141.0, 0.0001]) # our default for Ozark, g1 in sqrt(Pa). This is equal to 4.46 sqrt(kPa).
 
 # parameters prior
 prior_g1 = EKP.constrained_gaussian("g1", 221, 100, 0, Inf) # mean of 221 sqrt(Pa) = 7 sqrt(kPa), std of 100 (3 kPa)
-# Ollie: why does it return μ=5.3, σ=0.4? Are the values transformed?
+# Ollie: why does it return μ=5.3, σ=0.4? Are the values transformed? Is this to be of same order of magnitude?
+prior_g0 = EKP.constrained_gaussian("g0", 0.001, 0.01, 0, Inf)
+prior = EKP.combine_distributions([prior_g1, prior_g0])
 
 # generate the initial ensemble and set up the ensemble Kalman inversion
-N_ensemble = 5 # Ollie: is this parameter set at each iteration?
-N_iterations = 5 # Ollie: is this always equal to N_ensemble?
+N_ensemble = 5
+N_iterations = 5
 rng_seed = 41 # Ollie: why 41?
 rng = Random.MersenneTwister(rng_seed) # Ollie: pseudo-random number?
+Γ = 0.1 * EKP.I # What is this? Do I need to add it to target data?
 
-initial_ensemble = EKP.construct_initial_ensemble(rng, prior_g1, N_ensemble) # Ollie: what does this do? 5 random values of g1?
+initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ensemble) # Ollie: what does this do? 5 random values of g1?
 
-ensemble_kalman_process = EKP.EnsembleKalmanProcess(initial_ensemble, lhf_target, EKP.Inversion(); rng = rng) # Ollie: observation_series (lhf_target) needs to be ::OS? how?, do we need Γ? Do we need rng?
+ensemble_kalman_process = EKP.EnsembleKalmanProcess(initial_ensemble, lhf_target, Γ, EKP.Inversion(); rng = rng) # Ollie: do we need Γ? Do we need rng?
 
 # We are now ready to carry out the inversion. At each iteration, we get the ensemble from the last iteration, apply
-# ClimaLand_lhf(g1) to each ensemble member, and apply the Kalman update to the ensemble.
+# ClimaLand_lhf(params) to each ensemble member, and apply the Kalman update to the ensemble.
 
-for i in 1:N_iterations
-    params_i = EKP.get_ϕ_final(prior_g1, ensemble_kalman_process) # 5 values of g1
+for i in 1:N_iterations # This will run the model 5*5 = 25 times ( ? )
+    params_i = EKP.get_ϕ_final(prior, ensemble_kalman_process) # 5 values of params
 
-    ClimaLand_ens = hcat([ClimaLand_lhf(params_i[i]) for i in 1:N_ensemble]...) # adjusted for 1 param
+    ClimaLand_ens = hcat([ClimaLand_lhf(params_i[:, i]) for i in 1:N_ensemble]...)
 
     EKP.update_ensemble!(ensemble_kalman_process, ClimaLand_ens)
-end # This runs the simulation a bunch of times, but eventually crashed...
+end # EKP.update_ensemble!() crashes... ERROR: LinearAlgebra.SingularException(3)
 
 # Done! Here are the parameters:
 
-final_ensemble = EKP.get_ϕ_final(prior_g1, ensemble_kalman_process)
+final_ensemble = EKP.get_ϕ_final(prior, ensemble_kalman_process)