Add a minimal calibration example

This commit adds a simple calibration example, with detailed
explanations. We calibrate two parameters, g1 and g0, on
canopy transpiration data target, using EnsembleKalmanProcesses.jl.
We use a single site model, at the ozark site, and calibrate on
simulation output, with the goal of demonstrating the ability to
recover the original parameters and demonstrate a simple calibration
framework.

experiments/calibration/README.md

@@ -0,0 +1,9 @@
+# ClimaLand calibration
+
+This folder contains scripts to perform calibration on ClimaLand simulations.
+
+## Documentation
+
+To understand the dependencies of our calibration workflow, read the documentation of
+[ClimaCalibrate.jl](https://github.com/CliMA/ClimaCalibrate.jl) and [EnsembleKalmanProcesses.jl](https://github.com/CliMA/EnsembleKalmanProcesses.jl)
+

experiments/calibration/minimal_working_example/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+ClimaLand = "08f4d4ce-cf43-44bb-ad95-9d2d5f413532"
+ClimaLandSimulations = "348a0bd3-1299-4261-8002-d2cd97df6055"
+EnsembleKalmanProcesses = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d"
+GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"

experiments/calibration/minimal_working_example/minimal_working_example.jl

@@ -0,0 +1,96 @@
+# This script describes the simplest possible calibration of ClimaLand,
+# with step by step explanations.
+# If one line or comment confuses you, or you believe it could be simplified,
+# reach out to Alexis Renchon.
+
+# We calibrate a single site simulation, using EnsembleKalmanProcesses.jl.
+# 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 transpiration output given the parameters we want to calibrate.
+# 2. the "truth" target data, to calibrate on.
+# 3. the prior distribution of these parameters.
+
+# Let's go!
+
+# Packages
+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 transpiration output given the parameters to calibrate
+function ClimaLand_transpiration(params) # params is a 2 element Array
+    g1, g0 = params
+
+    sv = CLS.run_fluxnet(
+                         "US-MOz";
+                         params = CLS.ozark_default_params(; conductance = CLS.conductance_ozark(; g1 = g1, g0 = g0)),
+                        )[1];
+
+    inputs = CLS.make_inputs_df("US-MOz")[1];
+
+    simulation_output = CLS.make_output_df("US-MOz", sv, inputs);
+
+    transpiration = simulation_output.transpiration .* 1e9 # 1e9 because very small otherwise
+
+    return transpiration
+end
+
+# "truth" target data to calibrate on
+transpiration_target = ClimaLand_transpiration([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? Is this to be of same order of magnitude?
+prior_g0 = EKP.constrained_gaussian("g0", 0.00015, 0.0001, 0, Inf)
+prior = EKP.combine_distributions([prior_g1, prior_g0])
+
+# generate the initial ensemble and set up the ensemble Kalman inversion
+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, 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 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
+
+# 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])
+trange = 1:1:length(transpiration_target)
+
+l1 = lines!(ax, trange, ClimaLand_transpiration(theta_true), color = :black)
+[lines!(
+    ax,
+    trange,
+    ClimaLand_transpiration(EKP.get_ϕ(prior, ensemble_kalman_process, 1)[:, i]),
+    color = :red,
+   ) for i in 1:N_ensemble]
+[lines!(
+    ax,
+    trange,
+    ClimaLand_transpiration(final_ensemble[:, i]),
+    color = :blue,
+   ) for i in 1:N_ensemble]
+