Skip to content

Commit

Permalink
almost...
Browse files Browse the repository at this point in the history
  • Loading branch information
AlexisRenchon committed Sep 11, 2024
1 parent 250bce1 commit 5bd6ce6
Showing 1 changed file with 6 additions and 5 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)

0 comments on commit 5bd6ce6

Please sign in to comment.