diff --git a/experiments/calibration/minimal_working_example/minimal_working_example.jl b/experiments/calibration/minimal_working_example/minimal_working_example.jl index 31f248c6c8..88033b9859 100644 --- a/experiments/calibration/minimal_working_example/minimal_working_example.jl +++ b/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)