wip

experiments/calibration/minimal_working_example/minimal_working_example.jl

@@ -4,10 +4,10 @@
 # reach out to Alexis Renchon.
 
 # We calibrate a single site simulation, using EnsembleKalmanProcesses.jl.
-# We calibrate the parameters g1 and g0, to better match latent heat flux (lhf) target data.
+# We calibrate the parameters g1 and g0, to better match canopy transpiration target data.
 
 # To perform this calibration we will need:
-# 1. a function returning our model lhf output given the parameters we want to calibrate.
+# 1. a function returning our model transpiration output given the parameters we want to calibrate.
 # 2. the "truth" target data, to calibrate on.
 # 3. the prior distribution of these parameters.
 
@@ -18,8 +18,8 @@ import ClimaLandSimulations.Fluxnet as CLS # to run the model
 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_transpiration(params) # Should this return a vector or a ::OS?
+# a function returning our model transpiration output given the parameters to calibrate
+function ClimaLand_transpiration(params) # params is a 2 element Array
     g1, g0 = params
 
     sv = CLS.run_fluxnet(
@@ -31,9 +31,8 @@ function ClimaLand_transpiration(params) # Should this return a vector or a ::OS
 
     simulation_output = CLS.make_output_df("US-MOz", sv, inputs);
 
-    transpiration = simulation_output.transpiration .* 1e9
+    transpiration = simulation_output.transpiration .* 1e9 # 1e9 because very small otherwise
 
-    # lhf_obs = EKP.Observation(lhf, EKP.I(length(lhf)), "lhf") # target data needs to be of type ::OS. Not sure about I!!
     return transpiration
 end
 
@@ -53,50 +52,45 @@ 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, 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 params, normmalized?
 
 ensemble_kalman_process = EKP.EnsembleKalmanProcess(initial_ensemble, transpiration_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_transpiration(params) to each ensemble member, and apply the Kalman update to the ensemble.
 
 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
+    params_i = EKP.get_ϕ_final(prior, ensemble_kalman_process) # 5 random values of params from distribution?
 
     ClimaLand_ens = hcat([ClimaLand_transpiration(params_i[:, i]) for i in 1:N_ensemble]...)
 
     EKP.update_ensemble!(ensemble_kalman_process, ClimaLand_ens)
-end # EKP.update_ensemble!() crashes... ERROR: LinearAlgebra.SingularException(3)
+end
 
 # Done! Here are the parameters:
 
 final_ensemble = EKP.get_ϕ_final(prior, ensemble_kalman_process)
 
 
-
-
-
-
 # Plot
 using GLMakie
 GLMakie.activate!()
 theta_true = [141.0, 0.0001]
 fig = Figure()
-ax = Axis(fig[1,1])
+ax = Axis(fig[1, 1])
 trange = 1:1:length(transpiration_target)
 
 l1 = lines!(ax, trange, ClimaLand_transpiration(theta_true), color = :black)
-l2 = lines!(
+[lines!(
     ax,
     trange,
-    ClimaLand_transpiration(EKP.get_ϕ(prior, ensemble_kalman_process, 1)[:, 1]),
+    ClimaLand_transpiration(EKP.get_ϕ(prior, ensemble_kalman_process, 1)[:, i]),
     color = :red,
-)
-l3 = lines!(
+   ) for i in 1:N_ensemble]
+[lines!(
     ax,
     trange,
-    ClimaLand_transpiration(final_ensemble[:, 1]),
+    ClimaLand_transpiration(final_ensemble[:, i]),
     color = :blue,
-)
-
+   ) for i in 1:N_ensemble]