diff --git a/experiments/calibration/minimal_working_example/minimal_working_example.jl b/experiments/calibration/minimal_working_example/minimal_working_example.jl index 6fc9508196..5a61a22060 100644 --- a/experiments/calibration/minimal_working_example/minimal_working_example.jl +++ b/experiments/calibration/minimal_working_example/minimal_working_example.jl @@ -35,6 +35,7 @@ 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_Obs = EKP.Observation(lhf_target, EKP.I(length(lhf_target)), "lhf") # target data needs to be of type ::OS # 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) @@ -48,19 +49,19 @@ rng = Random.MersenneTwister(rng_seed) # Ollie: pseudo-random number? initial_ensemble = EKP.construct_initial_ensemble(rng, prior_g1, 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_Obs, EKP.Inversion(); rng = rng) # Ollie: observation_series (lhf_target) needs to be ::OS? how?, 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. for i in 1:N_iterations - params_i = get_ϕ_final(prior_g1, ensemble_kalman_process) + params_i = EKP.get_ϕ_final(prior_g1, ensemble_kalman_process) # 5 values of g1 - ClimaLand_ens = hcat([ClimaLand_lhf(params_i[:, i]) for i in 1:N_ensemble]...) + ClimaLand_ens = hcat([ClimaLand_lhf(params_i[i]) for i in 1:N_ensemble]...) # adjusted for 1 param EKP.update_ensemble!(ensemble_kalman_process, ClimaLand_ens) -end +end # This runs the simulation a bunch of times, but eventually crashed... ERROR: LinearAlgebra.SingularException(2) # Done! Here are the parameters: -final_ensemble = get_ϕ_final(prior_g1, ensemble_kalman_process) +final_ensemble = EKP.get_ϕ_final(prior_g1, ensemble_kalman_process)