diff --git a/.gitignore b/.gitignore index 68692582..c5958958 100644 --- a/.gitignore +++ b/.gitignore @@ -39,4 +39,4 @@ Results/ *.code-workspace -data +data \ No newline at end of file diff --git a/Manifest.toml b/Manifest.toml index 9c3baa6d..ec7cb980 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -340,9 +340,9 @@ version = "0.21.3" [[JSON3]] deps = ["Dates", "Mmap", "Parsers", "StructTypes", "UUIDs"] -git-tree-sha1 = "7d58534ffb62cd947950b3aa9b993e63307a6125" +git-tree-sha1 = "175b6ff26cd0fa01dd60021ce76bbdefdf91e4a0" uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" -version = "1.9.2" +version = "1.9.3" [[KernelAbstractions]] deps = ["Adapt", "Cassette", "InteractiveUtils", "MacroTools", "SpecialFunctions", "StaticArrays", "UUIDs"] @@ -549,9 +549,9 @@ version = "0.11.6" [[Parsers]] deps = ["Dates"] -git-tree-sha1 = "13468f237353112a01b2d6b32f3d0f80219944aa" +git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.2.2" +version = "2.2.3" [[PencilArrays]] deps = ["ArrayInterface", "JSON3", "LinearAlgebra", "MPI", "OffsetArrays", "Random", "Reexport", "Requires", "StaticArrays", "StaticPermutations", "Strided", "TimerOutputs", "VersionParsing"] @@ -784,10 +784,10 @@ uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c" version = "1.0.1" [[Tables]] -deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"] -git-tree-sha1 = "bb1064c9a84c52e277f1096cf41434b675cd368b" +deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"] +git-tree-sha1 = "5ce79ce186cc678bbb5c5681ca3379d1ddae11a1" uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" -version = "1.6.1" +version = "1.7.0" [[Tar]] deps = ["ArgTools", "SHA"] diff --git a/examples/calibrate_to_LESbrary/calibrate_CATKE_to_LESbrary.jl b/examples/calibrate_to_LESbrary/calibrate_CATKE_to_LESbrary.jl deleted file mode 100644 index 51a263af..00000000 --- a/examples/calibrate_to_LESbrary/calibrate_CATKE_to_LESbrary.jl +++ /dev/null @@ -1,82 +0,0 @@ -pushfirst!(LOAD_PATH, joinpath(@__DIR__, "../..")) - -using Oceananigans -using Plots, LinearAlgebra, Distributions, JLD2 -using Oceananigans.Units -using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity -using OceanTurbulenceParameterEstimation - -include("utils/lesbrary_paths.jl") -include("utils/one_dimensional_ensemble_model.jl") -include("utils/parameters.jl") -include("utils/visualize_profile_predictions.jl") - -# two_day_suite_dir = "/Users/gregorywagner/Projects/OceanTurbulenceParameterEstimation/data/2DaySuite" -# four_day_suite_dir = "/Users/gregorywagner/Projects/OceanTurbulenceParameterEstimation/data/4DaySuite" -# six_day_suite_dir = "/Users/gregorywagner/Projects/OceanTurbulenceParameterEstimation/data/6DaySuite" - -# two_day_suite = TwoDaySuite(two_day_suite_dir) -# four_day_suite = FourDaySuite(four_day_suite_dir) -# six_day_suite = SixDaySuite(six_day_suite_dir) - -# calibration = InverseProblem(two_day_suite, parameters; relative_weights = relative_weight_options["all_but_e"], -# architecture = GPU(), ensemble_size = 10, Δt = 30.0) - -# validation = InverseProblem(four_day_suite, calibration; Nz = 32); - -##### -##### Set up ensemble model -##### - -## NEED TO IMPLEMENT COARSE-GRAINING -## -## - -lesbrary_directory = "/Users/adelinehillier/Desktop/dev/" - -observations = TwoDaySuite(lesbrary_directory; first_iteration = 13, last_iteration = nothing, normalize = ZScore, Nz = 128) - -parameter_set = CATKEParametersRiDependent -closure = closure_with_parameter_set(CATKEVerticalDiffusivity(Float64;), parameter_set) - -ensemble_model = OneDimensionalEnsembleModel(observations; - architecture = CPU(), - ensemble_size = 20, - closure = closure -) - -ensemble_simulation = Simulation(ensemble_model; Δt = 10seconds, stop_time = 2days) - -pop!(ensemble_simulation.diagnostics, :nan_checker) - -##### -##### Build free parameters -##### - -free_parameter_names = keys(parameter_set.defaults) -free_parameter_means = collect(values(parameter_set.defaults)) -priors = NamedTuple(pname => ConstrainedNormal(0.0, 1.0, bounds(pname) .* 0.5...) for pname in free_parameter_names) - -free_parameters = FreeParameters(priors) - -##### -##### Build the Inverse Problem -##### - -calibration = InverseProblem(observations, ensemble_simulation, free_parameters); - -# ##### -# ##### Calibrate -# ##### - -iterations = 5 -# eki = EnsembleKalmanInversion(calibration; noise_covariance = 1e-2) -# params = iterate!(eki; iterations = iterations) - -directory = -# visualize!(calibration, params; -# field_names = [:u, :v, :b, :e], -# directory = @__DIR__, -# filename = "perfect_model_visual_calibrated.png" -# ) -# @show params diff --git a/examples/calibrate_to_LESbrary/calibrate_catke_to_perfect_model.jl b/examples/calibrate_to_LESbrary/calibrate_catke_to_perfect_model.jl deleted file mode 100644 index b582a656..00000000 --- a/examples/calibrate_to_LESbrary/calibrate_catke_to_perfect_model.jl +++ /dev/null @@ -1,172 +0,0 @@ -pushfirst!(LOAD_PATH, joinpath(@__DIR__, "../..")) - -using OceanTurbulenceParameterEstimation, LinearAlgebra, CairoMakie -using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity -using OceanTurbulenceParameterEstimation.EnsembleKalmanInversions: NaNResampler, FullEnsembleDistribution - -include("utils/lesbrary_paths.jl") -include("utils/parameters.jl") -include("utils/visualize_profile_predictions.jl") - -examples_path = joinpath(pathof(OceanTurbulenceParameterEstimation), "../../examples") -include(joinpath(examples_path, "intro_to_inverse_problems.jl")) - -# Pick a parameter set defined in `./parameters.jl` -parameter_set = CATKEParametersRiDependent -closure = closure_with_parameters(CATKEVerticalDiffusivity(Float64;), parameter_set.settings) - -# Pick the secret `true_parameters` -true_parameters = (Cᵟu = 0.5, CᴷRiʷ = 1.0, Cᵂu★ = 2.0, CᵂwΔ = 1.0, Cᴷeʳ = 5.0, Cᵟc = 0.5, Cᴰ = 2.0, Cᴷc⁻ = 0.5, Cᴷe⁻ = 0.2, Cᴷcʳ = 3.0, Cᴸᵇ = 1.0, CᴷRiᶜ = 1.0, Cᴷuʳ = 4.0, Cᴷu⁻ = 0.8, Cᵟe = 0.5) -true_closure = closure_with_parameters(closure, true_parameters) - -# Generate and load synthetic observations -kwargs = (tracers = (:b, :e), Δt = 10.0, stop_time = 4hours) - -# Note: if an output file of the same name already exist, `generate_synthetic_observations` will return the existing path and skip re-generating the data. -observ_path1 = generate_synthetic_observations("perfect_model_observation1"; Qᵘ = 2e-4, Qᵇ = 4e-8, f₀ = 1e-4, closure = true_closure, kwargs...) -observ_path2 = generate_synthetic_observations("perfect_model_observation2"; Qᵘ = -1e-4, Qᵇ = 0, f₀ = 0, closure = true_closure, kwargs...) -observation1 = SyntheticObservations(observ_path1, field_names = (:b, :e, :u, :v), normalize = ZScore) -observation2 = SyntheticObservations(observ_path2, field_names = (:b, :e, :u, :v), normalize = ZScore) -observations = [observation1, observation2] - -# Set up ensemble model -ensemble_simulation, closure = build_ensemble_simulation(observations; Nensemble = 50) - -# Build free parameters -build_prior(name) = ConstrainedNormal(0.0, 1.0, bounds(name) .* 0.5...) -free_parameters = FreeParameters(named_tuple_map(names(parameter_set), build_prior)) - -# Pack everything into Inverse Problem `calibration` -calibration = InverseProblem(observations, ensemble_simulation, free_parameters, output_map = ConcatenatedOutputMap()); - -# Make sure the forward map evaluated on the true parameters matches the observation map -x = forward_map(calibration, true_parameters)[:, 1:1]; -y = observation_map(calibration); -@show x == y - -directory = "calibrate_catke_to_perfect_model/" -isdir(directory) || mkpath(directory) - -visualize!(calibration, true_parameters; - field_names = (:u, :v, :b, :e), - directory = directory, - filename = "perfect_model_visual_true_params.png" -) - -# Calibrate - -noise_covariance = Matrix(Diagonal(observation_map_variance_across_time(calibration)[1, :, 1] .+ 1e-5)) -resampler = NaNResampler(; abort_fraction=0.5, distribution=FullEnsembleDistribution()) - -eki = EnsembleKalmanInversion(calibration; noise_covariance = noise_covariance, - resampler = resampler) -params = iterate!(eki; iterations = 5) - -### -### Summary Plots -### - -using CairoMakie -using LinearAlgebra - -visualize!(calibration, params; - field_names = [:u, :v, :b, :e], - directory = directory, - filename = "perfect_model_visual_calibrated.png" -) -@show params - -# Last, we visualize the outputs of EKI calibration. - -### Parameter convergence plot - -# Vector of NamedTuples, ensemble mean at each iteration -ensemble_means = getproperty.(eki.iteration_summaries, :ensemble_mean) - -# N_param x N_iter matrix, ensemble covariance at each iteration -θθ_std_arr = sqrt.(hcat(diag.(getproperty.(eki.iteration_summaries, :ensemble_cov))...)) - -N_param, N_iter = size(θθ_std_arr) -iter_range = 0:(N_iter-1) -pnames = calibration.free_parameters.names - -n_cols = 3 -n_rows = Int(ceil(N_param / n_cols)) -ax_coords = [(i, j) for i = 1:n_rows, j = 1:n_cols] - -f = Figure(resolution = (500n_cols, 200n_rows)) -for (i, pname) in enumerate(pnames) - coords = ax_coords[i] - ax = Axis(f[coords...], - xlabel = "Iteration", - xticks = iter_range, - ylabel = string(pname)) - - ax.ylabelsize = 20 - - mean_values = [getproperty.(ensemble_means, pname)...] - lines!(ax, iter_range, mean_values) - band!(ax, iter_range, mean_values .+ θθ_std_arr[i, :], mean_values .- θθ_std_arr[i, :]) - hlines!(ax, [true_parameters[pname]], color = :red) -end -save(joinpath(directory, "perfect_model_calibration_parameter_convergence.png"), f); - -### Pairwise ensemble plots - -N_param, N_iter = size(θθ_std_arr) -for pname1 in pnames, pname2 in pnames - if pname1 != pname2 - - f = Figure() - axtop = Axis(f[1, 1]) - axmain = Axis(f[2, 1], - xlabel = string(pname1), - ylabel = string(pname2) - ) - axright = Axis(f[2, 2]) - scatters = [] - for iteration in [0, 1, 2, N_iter - 1] - ensemble = eki.iteration_summaries[iteration].parameters - ensemble = [[particle[pname1], particle[pname2]] for particle in ensemble] - ensemble = transpose(hcat(ensemble...)) # N_ensemble x 2 - push!(scatters, scatter!(axmain, ensemble)) - density!(axtop, ensemble[:, 1]) - density!(axright, ensemble[:, 2], direction = :y) - end - vlines!(axmain, [true_parameters[pname1]], color = :red) - vlines!(axtop, [true_parameters[pname1]], color = :red) - hlines!(axmain, [true_parameters[pname2]], color = :red, alpha = 0.6) - hlines!(axright, [true_parameters[pname2]], color = :red, alpha = 0.6) - colsize!(f.layout, 1, Fixed(300)) - colsize!(f.layout, 2, Fixed(200)) - rowsize!(f.layout, 1, Fixed(200)) - rowsize!(f.layout, 2, Fixed(300)) - Legend(f[1, 2], scatters, - ["Initial ensemble", "Iteration 1", "Iteration 2", "Iteration $N_iter"], - position = :lb) - hidedecorations!(axtop, grid = false) - hidedecorations!(axright, grid = false) - xlims!(axright, 0, 10) - ylims!(axtop, 0, 10) - save(joinpath(directory, "eki_$(pname1)_$(pname2).png"), f) - end -end - -# Compare EKI result to true values - -weight_distances = [norm(collect(ensemble_means[iter]) .- collect(true_parameters)) for iter = 0:N_iter-1] -output_distances = [mapslices(norm, (forward_map(calibration, [ensemble_means...])[:, 1:N_iter] .- y), dims = 1)...] - -f = Figure() -lines(f[1, 1], iter_range, weight_distances, color = :red, linewidth = 2, - axis = (title = "Parameter distance", - xlabel = "Iteration", - ylabel = "|θ̅ₙ - θ⋆|", - yscale = log10)) -lines(f[1, 2], iter_range, output_distances, color = :blue, linewidth = 2, - axis = (title = "Output distance", - xlabel = "Iteration", - ylabel = "|G(θ̅ₙ) - y|", - yscale = log10)) - -save(joinpath(directory, "error_convergence_summary.png"), f); diff --git a/examples/calibrate_to_LESbrary/calibrate_convadj_to_lesbrary.jl b/examples/calibrate_to_LESbrary/calibrate_convadj_to_lesbrary.jl deleted file mode 100644 index 6b284063..00000000 --- a/examples/calibrate_to_LESbrary/calibrate_convadj_to_lesbrary.jl +++ /dev/null @@ -1,246 +0,0 @@ -# Calibrate convective adjustment closure parameters to LESbrary 2-day "free_convection" simulation - -using OceanTurbulenceParameterEstimation, LinearAlgebra, CairoMakie -using Oceananigans.Units - -include("utils/lesbrary_paths.jl") -include("utils/one_dimensional_ensemble_model.jl") - -# Build an observation from "free convection" LESbrary simulation - -LESbrary_directory = "/Users/adelinehillier/Desktop/dev/2DaySuite/" - -suite = OrderedDict("2d_free_convection" => ( - filename = joinpath(LESbrary_directory, "free_convection/instantaneous_statistics.jld2"), - fields = (:b,))) - -observations = SyntheticObservationsBatch(suite; first_iteration = 13, last_iteration = nothing, normalize = ZScore, Nz = 32) - -closure = ConvectiveAdjustmentVerticalDiffusivity(; - convective_κz = 1.0, - background_κz = 1e-4 -) - -# Build an ensemble simulation based on observation - -ensemble_size = 30 -ensemble_model = OneDimensionalEnsembleModel(observations; - architecture = CPU(), - ensemble_size = ensemble_size, - closure = closure -) - -ensemble_simulation = Simulation(ensemble_model; Δt = 10seconds, stop_time = 4days) - -# Specify priors - -priors = ( - convective_κz = ConstrainedNormal(0.0, 1.0, 0.1, 1.0), - background_κz = ConstrainedNormal(0.0, 1.0, 0.0, 10e-4) -) - -free_parameters = FreeParameters(priors) - -# Specify an output map that tracks 3 uniformly spaced time steps, ignoring the initial condition -track_times = floor.(range(1, stop = length(observations[1].times), length = 3)) -popfirst!(track_times) -output_map = ConcatenatedOutputMap(track_times) - -# Build `InverseProblem` -calibration = InverseProblem(observations, ensemble_simulation, free_parameters; output_map = output_map) - -# Ensemble Kalman Inversion - -eki = EnsembleKalmanInversion(calibration; noise_covariance = 0.01) - -iterations = 10 -iterate!(eki; iterations = iterations) - -# Visualize the outputs of EKI calibration. Plots will be stored in `directory`. - -directory = "calibrate_convadj_to_lesbrary/$(iterations)_iters_$(ensemble_size)_particles/" -isdir(directory) || mkpath(directory) - -### Parameter convergence plot - -# Vector of NamedTuples, ensemble mean at each iteration -ensemble_means = getproperty.(eki.iteration_summaries, :ensemble_mean) - -# N_param x N_iter matrix, ensemble covariance at each iteration -θθ_std_arr = sqrt.(hcat(diag.(getproperty.(eki.iteration_summaries, :ensemble_cov))...)) - -N_param, N_iter = size(θθ_std_arr) -iter_range = 0:(N_iter-1) -pnames = calibration.free_parameters.names - -n_cols = 3 -n_rows = Int(ceil(N_param / n_cols)) -ax_coords = [(i, j) for i = 1:n_rows, j = 1:n_cols] - -f = Figure(resolution = (500n_cols, 200n_rows)) -for (i, pname) in enumerate(pnames) - coords = ax_coords[i] - ax = Axis(f[coords...], - xlabel = "Iteration", - xticks = iter_range, - ylabel = string(pname)) - - ax.ylabelsize = 20 - - mean_values = [getproperty.(ensemble_means, pname)...] - lines!(ax, iter_range, mean_values) - band!(ax, iter_range, mean_values .+ θθ_std_arr[i, :], mean_values .- θθ_std_arr[i, :]) -end - -save(joinpath(directory, "conv_adj_to_LESbrary_parameter_convergence.pdf"), f); - -### Pairwise ensemble plots - -N_param, N_iter = size(θθ_std_arr) -for pname1 in pnames, pname2 in pnames - if pname1 != pname2 - - f = Figure() - axtop = Axis(f[1, 1]) - axmain = Axis(f[2, 1], - xlabel = string(pname1), - ylabel = string(pname2) - ) - axright = Axis(f[2, 2]) - scatters = [] - for iteration in [0, 1, 2, N_iter - 1] - ensemble = eki.iteration_summaries[iteration].parameters - ensemble = [[particle[pname1], particle[pname2]] for particle in ensemble] - ensemble = transpose(hcat(ensemble...)) # N_ensemble x 2 - push!(scatters, scatter!(axmain, ensemble)) - density!(axtop, ensemble[:, 1]) - density!(axright, ensemble[:, 2], direction = :y) - end - colsize!(f.layout, 1, Fixed(300)) - colsize!(f.layout, 2, Fixed(200)) - rowsize!(f.layout, 1, Fixed(200)) - rowsize!(f.layout, 2, Fixed(300)) - Legend(f[1, 2], scatters, - ["Initial ensemble", "Iteration 1", "Iteration 2", "Iteration $N_iter"], - position = :lb) - hidedecorations!(axtop, grid = false) - hidedecorations!(axright, grid = false) - xlims!(axright, 0, 10) - ylims!(axtop, 0, 10) - save(joinpath(directory, "conv_adj_to_LESbrary_eki_$(pname1)_$(pname2).pdf"), f) - end -end - -# Compare EKI result to true values - -y = observation_map(calibration) -output_distances = [mapslices(norm, (forward_map(calibration, [ensemble_means...])[:, 1:N_iter] .- y), dims = 1)...] - -f = Figure() -lines(f[1, 1], iter_range, output_distances, color = :blue, linewidth = 2, - axis = (title = "Output distance", - xlabel = "Iteration", - ylabel = "|G(θ̅ₙ) - y|", - yscale = log10)) -save(joinpath(directory, "conv_adj_to_LESbrary_error_convergence_summary.pdf"), f); - -include("./utils/visualize_profile_predictions.jl") -visualize!(calibration, ensemble_means[end]; - field_names = (:b,), - directory = directory, - filename = "realizations.pdf" -) - -θglobalmin = NamedTuple((:convective_κz => 0.275, :background_κz => 0.000275)) -visualize!(calibration, θglobalmin; - field_names = (:b,), - directory = directory, - filename = "realizations_θglobalmin.pdf" -) - -## Visualize loss landscape - -name = "Loss landscape" - -pvalues = Dict( - :convective_κz => collect(0.075:0.025:1.025), - :background_κz => collect(0e-4:0.25e-4:10e-4), -) - -ni = length(pvalues[:convective_κz]) -nj = length(pvalues[:background_κz]) - -params = hcat([[pvalues[:convective_κz][i], pvalues[:background_κz][j]] for i = 1:ni, j = 1:nj]...) -xc = params[1, :] -yc = params[2, :] - -# build an `InverseProblem` that can accommodate `ni*nj` ensemble members -ensemble_model = OneDimensionalEnsembleModel(observations; - architecture = CPU(), - ensemble_size = ni * nj, - closure = closure) -ensemble_simulation = Simulation(ensemble_model; Δt = 10seconds, stop_time = 6days) -calibration = InverseProblem(observations, ensemble_simulation, free_parameters) - -y = observation_map(calibration) - -using FileIO -# a = forward_map(calibration, params) .- y -# save("calibrate_convadj_to_lesbrary/loss_landscape.jld2", "a", a) - -a = load("calibrate_convadj_to_lesbrary/loss_landscape.jld2")["a"] -zc = [mapslices(norm, a, dims = 1)...] - -# 2D contour plot with EKI particles superimposed -begin - f = Figure() - ax1 = Axis(f[1, 1], - title = "EKI Particle Traversal Over Loss Landscape", - xlabel = "convective_κz", - ylabel = "background_κz") - - co = CairoMakie.contourf!(ax1, xc, yc, zc, levels = 50, colormap = :default) - - cvt(iter) = hcat(collect.(eki.iteration_summaries[iter].parameters)...) - diffc = cvt(2) .- cvt(1) - diff_mag = mapslices(norm, diffc, dims = 1) - # diffc ./= 2 - us = diffc[1, :] - vs = diffc[2, :] - xs = cvt(1)[1, :] - ys = cvt(1)[2, :] - - arrows!(xs, ys, us, vs, arrowsize = 10, lengthscale = 0.3, - arrowcolor = :yellow, linecolor = :yellow) - - am = argmin(zc) - minimizing_params = [xc[am] yc[am]] - - scatters = [scatter!(ax1, minimizing_params, marker = :x, markersize = 30)] - for (i, iteration) in enumerate([1, 2, iterations]) - ensemble = eki.iteration_summaries[iteration].parameters - ensemble = [[particle[:convective_κz], particle[:background_κz]] for particle in ensemble] - ensemble = transpose(hcat(ensemble...)) # N_ensemble x 2 - push!(scatters, scatter!(ax1, ensemble)) - end - Legend(f[1, 2], scatters, - ["Global minimum", "Initial ensemble", "Iteration 1", "Iteration $(iterations)"], - position = :lb) - - save(joinpath(directory, "loss_contour.pdf"), f) -end - -# 3D loss landscape -begin - f = Figure() - ax1 = Axis3(f[1, 1], - title = "Loss Landscape", - xlabel = "convective_κz", - ylabel = "background_κz", - zlabel = "MSE loss" - ) - - CairoMakie.surface!(ax1, xc, yc, zc, colorscheme = :thermal) - - save(joinpath(directory, "loss_landscape.png"), f) -end \ No newline at end of file diff --git a/examples/calibrate_to_LESbrary/future_calibrate_catke_to_lesbrary.jl b/examples/calibrate_to_LESbrary/future_calibrate_catke_to_lesbrary.jl deleted file mode 100644 index 11e8681d..00000000 --- a/examples/calibrate_to_LESbrary/future_calibrate_catke_to_lesbrary.jl +++ /dev/null @@ -1,69 +0,0 @@ -# In this example, we use EKI to tune the closure parameters of a HydrostaticFreeSurfaceModel -# with a CATKEBasedVerticalDiffusivity closure in order to align the predictions of the model -# to those of a high-resolution LES data generated in LESbrary.jl. Here `predictions` refers to the -# 1-D profiles of temperature, velocity, and turbulent kinetic energy horizontally averaged over a -# 3-D physical domain. - -pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..")) -pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", "projects", "OceanBoundaryLayerParameterizations", "src")) - -image_dir = joinpath(@__DIR__, "quick_calibrate") -mkpath(image_dir) - -using Oceananigans -using OceanTurbulenceParameterEstimation -using OceanBoundaryLayerParameterizations - -two_day_suite_dir = "/Users/gregorywagner/Projects/OceanTurbulenceParameterEstimation/data/2DaySuite" -four_day_suite_dir = "/Users/gregorywagner/Projects/OceanTurbulenceParameterEstimation/data/4DaySuite" -six_day_suite_dir = "/Users/gregorywagner/Projects/OceanTurbulenceParameterEstimation/data/6DaySuite" - -two_day_suite = TwoDaySuite(two_day_suite_dir) -four_day_suite = FourDaySuite(four_day_suite_dir) -six_day_suite = SixDaySuite(six_day_suite_dir) - -calibration = InverseProblem(two_day_suite, parameters; relative_weights = relative_weight_options["all_but_e"], - architecture = GPU(), ensemble_size = 10, Δt = 30.0) - -validation = InverseProblem(four_day_suite, calibration; Nz = 64); - -# Loss on default parameters -l0 = calibration() - -# Example parameters -θ = calibration.default_parameters - -# Loss on parameters θ. -# θ can be -# 1. a vector -# 2. a FreeParameters object -# 3. a vector of parameter vectors (one for each ensemble member) -# 4. or a vector of FreeParameter objects (one for each ensemble member) -# If (1) or (2), the ensemble members are redundant and the loss is computed for just the one parameter set. -lθ = calibration(θ) - -# Output files/figures -directory = joinpath(pwd(), "quick_calibrate") - -# Run the model forward and store the solution -output = model_time_series(calibration, θ) - -# Run the model forward with parameters θ and visualize the solution compared to the truth -visualize!(output; filename = "visualize_default_parameters.png") - -#= - -# Use EKI to calibrate the model parameters -calibration_algorithm = EKI(noise_level = 10^(-2.0), - N_iter = 15, - stds_within_bounds = 0.6, - informed_priors = false) - -best_parameters = calibrate(calibration; algorithm = calibration_algorithm) - -loss = calibration(best_parameters) - -=# - -# Runs `visualize!` and records a summary of the calibration results in a `result.txt` file. -visualize_and_save!(calibration, validation, θ, directory) diff --git a/examples/calibrate_to_LESbrary/utils/animate_profile_predictions.jl b/examples/calibrate_to_LESbrary/utils/animate_profile_predictions.jl deleted file mode 100644 index dcf7730d..00000000 --- a/examples/calibrate_to_LESbrary/utils/animate_profile_predictions.jl +++ /dev/null @@ -1,71 +0,0 @@ - -colors = Dict("u" => :blue, "v" => :green, "b" => :red, :e => :yellow) -x_lims = Dict("u" => (-0.3,0.4), "v" => (-0.3,0.1), "b" => (19.6,20.0), "e" => (-0.5,4)) - -function animate_LESbrary_suite(ce, directory; parameters=ce.default_parameters, targets=ce.validation.loss.loss.batch[1].loss.targets) - - truth = Dict() - model = Dict() - for myloss in ce.validation.loss.loss.batch - truth[myloss.data.name] = myloss.data - model[myloss.data.name] = model_time_series(parameters, myloss.model, myloss.data) - end - - function plot_(name, file, t) - start=13 - truth_field = getproperty(truth[file], Symbol(name)) - catke_field = getproperty(model[file], Symbol(name)) - Nz = truth_field[1].grid.Nz - z = znodes(truth_field[1])[1:Nz-1] - Nz_catke = catke_field[1].grid.Nz - z_catke = znodes(catke_field[1])[1:Nz_catke-1] - p = Plots.plot(legend=false, plot_titlefontsize=20, xlims=x_lims[name]) - for t = t - truth_profile = truth_field[t+13].data[1:Nz-1] - catke_profile = catke_field[t].data[1:Nz_catke-1] - plot!(truth_profile, z, color=colors[name], linewidth=10, la=0.3) - if !any(i -> isnan(i), catke_profile) - plot!(catke_profile, z_catke, color=colors[name], linewidth=3, linestyle=:solid) - # scatter!(catke_profile, z_catke, markersize=3, color=colors[name], linewidth=10) - end - if name=="u" - truth_profile = getproperty(truth[file], :v)[t+13].data[1:Nz-1] - catke_profile = getproperty(model[file], :v)[t].data[1:Nz_catke-1] - plot!(truth_profile, z, color=colors["v"], linewidth=10, la=0.3) - if !any(i -> isnan(i), catke_profile) - plot!(catke_profile, z_catke, color=colors["v"], linewidth=3, linestyle=:solid) - # scatter!(catke_profile, z_catke, markersize=3, color=colors["v"], linewidth=10) - end - end - end - p - end - - function stacked_(file, t) - u = plot_("u", file, t) - t = plot_("b", file, t) - layout=@layout[a; b] - p = Plots.plot(t, u, layout=layout) - plot!(tickfontsize=20, ylims=(-256,0), ticks=false) - plot!(widen=true, grid=false, framestyle=:none) - return p - end - - function LESbrary_suite_snapshot(t) - a = stacked_("free_convection", t) - b = stacked_("strong_wind", t) - c = stacked_("strong_wind_no_rotation", t) - d = stacked_("weak_wind_strong_cooling", t) - e = stacked_("strong_wind_weak_cooling", t) - layout = @layout [a b c d e] - p = Plots.plot(a, b, c, d, e, layout=layout, framestyle=:none) - plot!(bottom_margin=0*Plots.mm, size=(1800, 800)) - return p - end - - anim = @animate for t=targets - p = LESbrary_suite_snapshot(t) - end - - Plots.gif(anim, directory, fps=400) -end diff --git a/examples/calibrate_to_LESbrary/utils/legacy_data_time_serieses.jl b/examples/calibrate_to_LESbrary/utils/legacy_data_time_serieses.jl deleted file mode 100644 index b1954b0c..00000000 --- a/examples/calibrate_to_LESbrary/utils/legacy_data_time_serieses.jl +++ /dev/null @@ -1,51 +0,0 @@ -using JLD2 -using Oceananigans.Fields: location, XFaceField, YFaceField, ZFaceField, CenterField -using Oceananigans.Grids: Face, Center - -location_guide = Dict(:u => (Face, Center, Center), - :v => (Center, Face, Center), - :w => (Center, Center, Face)) - -# If the location isn't correct in `observations`, `interior` won't be computed correctly -# in InverseProblems.transpose_model_output -function infer_location(field_name) - if field_name in keys(location_guide) - return location_guide[field_name] - else - return (Center, Center, Center) - end -end - -function legacy_data_field_time_serieses(path, field_names, times) - - # Build a grid, assuming it's a 1D RectilinearGrid - file = jldopen(path) - - old_Nz = file["grid/Nz"] - Hz = file["grid/Hz"] - Lz = file["grid/Lz"] - - αg = file["buoyancy/gravitational_acceleration"] * file["buoyancy/equation_of_state/α"] - Qᵇ = file["parameters/boundary_condition_θ_top"] * αg - Qᵘ = file["parameters/boundary_condition_u_top"] - u_bottom = file["parameters/boundary_condition_u_bottom"] - b_bottom = file["parameters/boundary_condition_θ_bottom"] * αg - - close(file) - - b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ), bottom = GradientBoundaryCondition(b_bottom)) - u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ), bottom = GradientBoundaryCondition(u_bottom)) - - grid = RectilinearGrid(size = old_Nz, halo = Hz, z = (-Lz, 0), topology = (Flat, Flat, Bounded)) - - boundary_conditions = (; u = u_bcs, b = b_bcs) - - field_time_serieses = NamedTuple(name => FieldTimeSeries(path, string(name); - times = times, - location = infer_location(name), - grid = grid, - boundary_conditions = boundary_conditions) - for name in field_names) - - return field_time_serieses -end \ No newline at end of file diff --git a/examples/calibrate_to_LESbrary/utils/lesbrary_paths.jl b/examples/calibrate_to_LESbrary/utils/lesbrary_paths.jl deleted file mode 100644 index 7e2cd326..00000000 --- a/examples/calibrate_to_LESbrary/utils/lesbrary_paths.jl +++ /dev/null @@ -1,134 +0,0 @@ -using OrderedCollections -using Oceananigans.OutputReaders: FieldTimeSeries - -##### -##### The LESbrary (so to speak) -##### - -include("./legacy_data_time_serieses.jl") - -function get_times(path) - file = jldopen(path) - file_iterations = parse.(Int, keys(file["timeseries/t"])) - file_times = [file["timeseries/t/$i"] for i in file_iterations] - return file_times -end - -# https://engaging-web.mit.edu/~alir/lesbrary/6DaySuite/ -function SyntheticObservationsBatch(suite; first_iteration = 1, stride = 1, last_iteration = nothing, normalize = ZScore, Nz) - - observations = Vector{SyntheticObservations}() - - for case in values(suite) - path = case.filename - times = get_times(path) - last_iteration = isnothing(last_iteration) ? length(times) : last_iteration - times = times[first_iteration:stride:last_iteration] - field_names = case.fields - - field_time_serieses = legacy_data_field_time_serieses(path, field_names, times) - - observation = SyntheticObservations(path; field_names, normalize, times, field_time_serieses, grid_size = (1, 1, Nz)) - push!(observations, observation) - end - - return observations -end - -# https://engaging-web.mit.edu/~alir/lesbrary/2DaySuite/ -function TwoDaySuite(directory; first_iteration = 13, stride = 1, last_iteration = nothing, normalize = ZScore, Nz = 128) - directory = joinpath(directory, "2DaySuite") - suite = OrderedDict( - "2d_free_convection" => ( - filename = joinpath(directory, "free_convection/instantaneous_statistics.jld2"), - fields = (:b, :e)), - "2d_strong_wind" => ( - filename = joinpath(directory, "strong_wind/instantaneous_statistics.jld2"), - fields = (:b, :u, :v, :e)), - "2d_strong_wind_no_rotation" => ( - filename = joinpath(directory, "strong_wind_no_rotation/instantaneous_statistics.jld2"), - fields = (:b, :u, :e)), - "2d_strong_wind_weak_cooling" => ( - filename = joinpath(directory, "strong_wind_weak_cooling/instantaneous_statistics.jld2"), - fields = (:b, :u, :v, :e)), - "2d_weak_wind_strong_cooling" => ( - filename = joinpath(directory, "weak_wind_strong_cooling/instantaneous_statistics.jld2"), - fields = (:b, :u, :v, :e)), - ) - - return SyntheticObservationsBatch(suite; first_iteration, stride, last_iteration, normalize, Nz) -end - -# https://engaging-web.mit.edu/~alir/lesbrary/4DaySuite/ -function FourDaySuite(directory; first_iteration = 13, stride = 1, last_iteration = nothing, normalize = ZScore, Nz = 128) - directory = joinpath(directory, "4DaySuite") - suite = OrderedDict( - "4d_free_convection" => ( - filename = joinpath(directory, "free_convection/instantaneous_statistics.jld2"), - fields = (:b, :e)), - "4d_strong_wind" => ( - filename = joinpath(directory, "strong_wind/instantaneous_statistics.jld2"), - fields = (:b, :u, :v, :e)), - "4d_strong_wind_no_rotation" => ( - filename = joinpath(directory, "strong_wind_no_rotation/instantaneous_statistics.jld2"), - fields = (:b, :u, :e)), - "4d_strong_wind_weak_cooling" => ( - filename = joinpath(directory, "strong_wind_weak_cooling/instantaneous_statistics.jld2"), - fields = (:b, :u, :v, :e)), - "4d_weak_wind_strong_cooling" => ( - filename = joinpath(directory, "weak_wind_strong_cooling/instantaneous_statistics.jld2"), - fields = (:b, :u, :v, :e)), - ) - - return SyntheticObservationsBatch(suite; first_iteration, stride, last_iteration, normalize, Nz) -end - -function SixDaySuite(directory; first_iteration = 13, stride = 1, last_iteration = nothing, normalize = ZScore, Nz = 128) - directory = joinpath(directory, "6DaySuite") - suite = OrderedDict( - "6d_free_convection" => ( - filename = joinpath(directory, "free_convection/instantaneous_statistics.jld2"), - fields = (:b, :e)), - "6d_strong_wind" => ( - filename = joinpath(directory, "strong_wind/instantaneous_statistics.jld2"), - fields = (:b, :u, :v, :e)), - "6d_strong_wind_no_rotation" => ( - filename = joinpath(directory, "strong_wind_no_rotation/instantaneous_statistics.jld2"), - fields = (:b, :u, :e)), - "6d_strong_wind_weak_cooling" => ( - filename = joinpath(directory, "strong_wind_weak_cooling/instantaneous_statistics.jld2"), - fields = (:b, :u, :v, :e)), - "6d_weak_wind_strong_cooling" => ( - filename = joinpath(directory, "weak_wind_strong_cooling/instantaneous_statistics.jld2"), - fields = (:b, :u, :v, :e)), - ) - - return SyntheticObservationsBatch(suite; first_iteration, stride, last_iteration, normalize, Nz) -end - -function GeneralStrat(directory) - directory = joinpath(directory, "6DaySuite") - suite = OrderedDict( - "general_strat_4" => ( - filename = joinpath(directory, "general_strat_4/instantaneous_statistics.jld2"), - fields = (:b,), - first = 37, # cut out the first 6 hours - last = 288), - "general_strat_8" => ( - filename = joinpath(directory, "general_strat_8/instantaneous_statistics.jld2"), - fields = (:b,), - first = 13, - last = 648), - "general_strat_16" => ( - filename = joinpath(directory, "general_strat_16/instantaneous_statistics.jld2"), - fields = (:b,), - first = 13, - last = nothing), - "general_strat_32" => ( - filename = joinpath(directory, "general_strat_32/instantaneous_statistics.jld2"), - fields = (:b,), - first = 13, - last = nothing), - ) -end - diff --git a/examples/calibrate_to_LESbrary/utils/one_dimensional_ensemble_model.jl b/examples/calibrate_to_LESbrary/utils/one_dimensional_ensemble_model.jl deleted file mode 100644 index 9f1c03a9..00000000 --- a/examples/calibrate_to_LESbrary/utils/one_dimensional_ensemble_model.jl +++ /dev/null @@ -1,79 +0,0 @@ -using Oceananigans -using JLD2 -using Oceananigans.Architectures: arch_array -using Oceananigans: AbstractModel -using Oceananigans.Models.HydrostaticFreeSurfaceModels: HydrostaticFreeSurfaceModel -using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize - -function get_parameter(observation, parameter_path) - file = jldopen(observation.path) - parameter = file[parameter_path] - close(file) - return parameter -end - -vectorize(observation) = [observation] -vectorize(observations::Vector) = observations - -""" - OneDimensionalEnsembleModel(observations::SyntheticObservationsBatch; architecture = CPU(), ensemble_size = 50, kwargs...) - -Build an Oceananigans `HydrostaticFreeSurfaceModel` with many independent -columns. The model grid is given by the data in `observations`, and the -dynamics in each column is attached to its own `CATKEVerticalDiffusivity` closure stored in -an `(Nx, Ny)` Matrix of closures. `ensemble_size = 50` is the -desired count of ensemble members for calibration with Ensemble Kalman Inversion (EKI), and the -remaining keyword arguments `kwargs` define the default closure across all columns. - -In the "many columns" configuration, we run the model on a 3D grid with `(Flat, Flat, Bounded)` boundary -conditions so that many independent columns can be evolved at once with much of the computational overhead -split among the columns. The `Nx` rows of vertical columns are each reserved for an "ensemble" member -whose attached parameter value (updated at each iteration of EKI) sets the diffusivity closure -used to predict the model solution for the `Ny` physical scenarios described by the simulation-specific -`SyntheticObservations` objects in `observations`. -""" -function OneDimensionalEnsembleModel(observations; - architecture = CPU(), - ensemble_size = 50, - closure = CATKEVerticalDiffusivity(Float64; warning = false) -) - observations = vectorize(observations) - - Lz = get_parameter(observations[1], "grid/Lz") - g = get_parameter(observations[1], "buoyancy/gravitational_acceleration") - α = get_parameter(observations[1], "buoyancy/equation_of_state/α") - - Nx = ensemble_size - Ny = length(observations) - Nz = observations[1].grid.Nz - column_ensemble_size = ColumnEnsembleSize(Nz = Nz, ensemble = (Nx, Ny), Hz = 1) - - ensemble_grid = RectilinearGrid(size = column_ensemble_size, z = (-Lz, 0), topology = (Flat, Flat, Bounded)) - closure_ensemble = arch_array(architecture, [closure for i = 1:Nx, j = 1:Ny]) - - get_Qᵇ(observation) = get_parameter(observation, "parameters/boundary_condition_θ_top") * α * g - get_Qᵘ(observation) = get_parameter(observation, "parameters/boundary_condition_u_top") - get_f₀(observation) = get_parameter(observation, "coriolis/f") - get_dudz_bottom(observation) = get_parameter(observation, "parameters/boundary_condition_u_bottom") - get_dbdz_bottom(observation) = get_parameter(observation, "parameters/boundary_condition_θ_bottom") * α * g - - to_arch_array(f) = arch_array(architecture, [f(observations[j]) for i = 1:Nx, j = 1:Ny]) - - Qᵇ_ensemble = to_arch_array(get_Qᵇ) - Qᵘ_ensemble = to_arch_array(get_Qᵘ) - ensemble_u_bottom = to_arch_array(get_dudz_bottom) - ensemble_b_bottom = to_arch_array(get_dbdz_bottom) - - ensemble_b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ_ensemble), bottom = GradientBoundaryCondition(ensemble_b_bottom)) - ensemble_u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ_ensemble), bottom = GradientBoundaryCondition(ensemble_u_bottom)) - coriolis_ensemble = arch_array(architecture, [FPlane(f = get_f₀(observations[j])) for i = 1:Nx, j = 1:Ny]) - - ensemble_model = HydrostaticFreeSurfaceModel(grid = ensemble_grid, - tracers = (:b, :e), - buoyancy = BuoyancyTracer(), - boundary_conditions = (; u = ensemble_u_bcs, b = ensemble_b_bcs), - coriolis = coriolis_ensemble, - closure = closure_ensemble) - - return ensemble_model -end \ No newline at end of file diff --git a/examples/calibrate_to_LESbrary/utils/parameters.jl b/examples/calibrate_to_LESbrary/utils/parameters.jl deleted file mode 100644 index f61df463..00000000 --- a/examples/calibrate_to_LESbrary/utils/parameters.jl +++ /dev/null @@ -1,74 +0,0 @@ -using LaTeXStrings - -import OceanTurbulenceParameterEstimation.TurbulenceClosureParameters: closure_with_parameters - -parameter_guide = Dict(:Cᴰ => (name = "Dissipation parameter (TKE equation)", latex = L"C^D", - default = 2.9079, bounds = (0.0, 10.0)), :Cᴸᵇ => (name = "Mixing length parameter", latex = L"C^{\ell}_b", - default = 1.1612, bounds = (0.0, 10.0)), :Cᵟu => (name = "Ratio of mixing length to grid spacing", latex = L"C^{\delta}u", - default = 0.5, bounds = (0.0, 3.0)), :Cᵟc => (name = "Ratio of mixing length to grid spacing", latex = L"C^{\delta}c", - default = 0.5, bounds = (0.0, 3.0)), :Cᵟe => (name = "Ratio of mixing length to grid spacing", latex = L"C^{\delta}e", - default = 0.5, bounds = (0.0, 3.0)), :Cᵂu★ => (name = "Mixing length parameter", latex = L"C^W_{u\star}", - default = 3.6188, bounds = (0.0, 20.0)), :CᵂwΔ => (name = "Mixing length parameter", latex = L"C^Ww\Delta", - default = 1.3052, bounds = (0.0, 5.0)), :CᴷRiʷ => (name = "Stability function parameter", latex = L"C^KRi^w", - default = 0.7213, bounds = (0.0, 5.0)), :CᴷRiᶜ => (name = "Stability function parameter", latex = L"C^KRi^c", - default = 0.7588, bounds = (-1.5, 4.0)), :Cᴷu⁻ => (name = "Velocity diffusivity LB", latex = L"C^Ku^-", - default = 0.1513, bounds = (0.0, 2.0)), :Cᴷuʳ => (name = "Velocity diffusivity (UB-LB)/LB", latex = L"C^Ku^r", - default = 3.8493, bounds = (0.0, 50.0)), :Cᴷc⁻ => (name = "Velocity diffusivity LB", latex = L"C^Kc^-", - default = 0.3977, bounds = (0.0, 5.0)), :Cᴷcʳ => (name = "Velocity diffusivity (UB-LB)/LB", latex = L"C^Kc^r", - default = 3.4601, bounds = (0.0, 50.0)), :Cᴷe⁻ => (name = "Velocity diffusivity LB", latex = L"C^Ke^-", - default = 0.1334, bounds = (0.0, 3.0)), :Cᴷeʳ => (name = "Velocity diffusivity (UB-LB)/LB", latex = L"C^Ke^r", - default = 8.1806, bounds = (0.0, 50.0)), :Cᴬu => (name = "Convective adjustment velocity parameter", latex = L"C^A_U", - default = 0.0057, bounds = (0.0, 0.2)), :Cᴬc => (name = "Convective adjustment tracer parameter", latex = L"C^A_C", - default = 0.6706, bounds = (0.0, 2.0)), :Cᴬe => (name = "Convective adjustment TKE parameter", latex = L"C^A_E", - default = 0.2717, bounds = (0.0, 2.0)), -) - -bounds(name) = parameter_guide[name].bounds -default(name) = parameter_guide[name].default - -function named_tuple_map(names, f) - names = Tuple(names) - return NamedTuple{names}(f.(names)) -end - -""" - ParameterSet() - -Parameter set containing the names `names` of parameters, and a -NamedTuple `settings` mapping names of "background" parameters -to their fixed values to be maintained throughout the calibration. -""" -struct ParameterSet{N,S} - names :: N - settings :: S -end - -""" - ParameterSet(names::Set; nullify = Set()) - -Construct a `ParameterSet` containing all of the information necessary -to build a closure with the specified default parameters and settings, -given a Set `names` of the parameter names to be tuned, and a Set `nullify` -of parameters to be set to zero. -""" -function ParameterSet(names::Set; nullify = Set()) - zero_set = named_tuple_map(nullify, name -> 0.0) - bkgd_set = named_tuple_map(keys(parameter_guide), name -> default(name)) - settings = merge(bkgd_set, zero_set) # order matters: `zero_set` overrides `bkgd_set` - return ParameterSet(Tuple(names), settings) -end - -names(ps::ParameterSet) = ps.names - -### -### Define some convenient parameter sets based on the present CATKE formulation -### - -required_params = Set([:Cᵟu, :Cᵟc, :Cᵟe, :Cᴰ, :Cᴸᵇ, :Cᵂu★, :CᵂwΔ, :Cᴷu⁻, :Cᴷc⁻, :Cᴷe⁻]) -ri_depen_params = Set([:CᴷRiʷ, :CᴷRiᶜ, :Cᴷuʳ, :Cᴷcʳ, :Cᴷeʳ]) -conv_adj_params = Set([:Cᴬu, :Cᴬc, :Cᴬe]) - -CATKEParametersRiDependent = ParameterSet(union(required_params, ri_depen_params); nullify = conv_adj_params) -CATKEParametersRiIndependent = ParameterSet(union(required_params); nullify = union(conv_adj_params, ri_depen_params)) -CATKEParametersRiDependentConvectiveAdjustment = ParameterSet(union(required_params, conv_adj_params, ri_depen_params)) -CATKEParametersRiIndependentConvectiveAdjustment = ParameterSet(union(required_params, conv_adj_params); nullify = ri_depen_params) diff --git a/examples/calibrate_to_LESbrary/utils/visualize_profile_predictions.jl b/examples/calibrate_to_LESbrary/utils/visualize_profile_predictions.jl deleted file mode 100644 index 8d9aca15..00000000 --- a/examples/calibrate_to_LESbrary/utils/visualize_profile_predictions.jl +++ /dev/null @@ -1,165 +0,0 @@ -using OceanTurbulenceParameterEstimation.InverseProblems: vectorize, forward_run!, transpose_model_output -using CairoMakie - -include("visualize_profile_predictions_utils.jl") - -""" - visualize!(ip::InverseProblem, parameters; - field_names = [:u, :v, :b, :e], - directory = pwd(), - filename = "realizations.png" - ) - - For visualizing 1-dimensional time series predictions. -""" -function visualize!(ip::InverseProblem, parameters; - field_names = [:u, :v, :b, :e], - directory = pwd(), - filename = "realizations.png" - ) - - isdir(directory) || makedir(directory) - - model = ip.simulation.model - - n_fields = length(field_names) - - observations = vectorize(ip.observations) - - forward_run!(ip, parameters) - - # Vector of SyntheticObservations objects, one for each observation - predictions = transpose_model_output(ip.time_series_collector, observations) - - fig = Figure(resolution = (200*(length(field_names)+1), 200*(length(ip.observations)+1)), font = "CMU Serif") - colors = [:black, :red, :blue] - - function empty_plot!(fig_position) - ax = fig_position = Axis(fig_position) - hidedecorations!(ax) - hidespines!(ax, :t, :b, :l, :r) - end - - calibration.observations[1].field_time_serieses.b[1].grid - - for (oi, observation) in enumerate(observations) - - i = oi + 1 - prediction = predictions[oi] - - targets = observation.times - snapshots = round.(Int, range(1, length(targets), length=3)) - - Qᵇ = model.tracers.b.boundary_conditions.top.condition[1,oi] - Qᵘ = model.velocities.u.boundary_conditions.top.condition[1,oi] - fv = model.coriolis[1,oi].f - - empty_plot!(fig[i,1]) - text!(fig[i,1], "Qᵇ = $(tostring(Qᵇ)) m⁻¹s⁻³\nQᵘ = $(tostring(Qᵘ)) m⁻¹s⁻²\nf = $(tostring(fv)) s⁻¹", - position = (0, 0), - align = (:center, :center), - textsize = 15, - justification = :left) - - for (j, field_name) in enumerate(field_names) - - middle = j > 1 && j < n_fields - remove_spines = j == 1 ? (:t, :r) : j == n_fields ? (:t, :l) : (:t, :l, :r) - axis_position = j == n_fields ? (ylabelposition=:right, yaxisposition=:right) : NamedTuple() - - j += 1 # reserve the first column for row labels - - info = field_guide[field_name] - - grid = observation.grid - - z = field_name ∈ [:u, :v] ? grid.zᵃᵃᶠ[1:grid.Nz] : grid.zᵃᵃᶜ[1:grid.Nz] - - to_plot = field_name ∈ keys(observation.field_time_serieses) - - if to_plot - - # field data for each time step - obsn = getproperty(observation.field_time_serieses, field_name) - pred = getproperty(prediction.field_time_serieses, field_name) - - ax = Axis(fig[i,j]; xlabelpadding=0, xtickalign=1, ytickalign=1, - merge(axis_position, info.axis_args)...) - - hidespines!(ax, remove_spines...) - - middle && hideydecorations!(ax, grid=false) - - lins = [] - for (color_index, t) in enumerate(snapshots) - l = lines!([obsn[1,1,1:grid.Nz,t] .* info.scaling ...], z; color = (colors[color_index], 0.4)) - push!(lins, l) - l = lines!([pred[1,1,1:grid.Nz,t] .* info.scaling ...], z; color = (colors[color_index], 1.0), linestyle = :dash) - push!(lins, l) - end - - times = @. round((targets[snapshots] - targets[snapshots[1]]) / 86400, sigdigits=2) - - legendlabel(time) = ["LES, t = $time days", "Model, t = $time days"] - Legend(fig[1,2:3], lins, vcat([legendlabel(time) for time in times]...), nbanks=2) - lins = [] - else - - empty_plot!(fig[i,j]) - end - end - end - - save(joinpath(directory, filename), fig, px_per_unit = 2.0) - return nothing -end - -# visualize!(ip::InverseProblem, parameters; kwargs...) = visualize!(model_time_series(ip, parameters); kwargs...) - -# function visualize_and_save!(calibration, validation, parameters, directory; fields=[:u, :v, :b, :e]) -# isdir(directory) || makedir(directory) - -# path = joinpath(directory, "results.txt") -# o = open_output_file(path) -# write(o, "Training relative weights: $(calibration.relative_weights) \n") -# write(o, "Validation relative weights: $(validation.relative_weights) \n") -# write(o, "Training default parameters: $(validation.default_parameters) \n") -# write(o, "Validation default parameters: $(validation.default_parameters) \n") - -# write(o, "------------ \n \n") -# default_parameters = calibration.default_parameters -# train_loss_default = calibration(default_parameters) -# valid_loss_default = validation(default_parameters) -# write(o, "Default parameters: $(default_parameters) \nLoss on training: $(train_loss_default) \nLoss on validation: $(valid_loss_default) \n------------ \n \n") - -# train_loss = calibration(parameters) -# valid_loss = validation(parameters) -# write(o, "Parameters: $(parameters) \nLoss on training: $(train_loss) \nLoss on validation: $(valid_loss) \n------------ \n \n") - -# write(o, "Training loss reduction: $(train_loss/train_loss_default) \n") -# write(o, "Validation loss reduction: $(valid_loss/valid_loss_default) \n") -# close(o) - -# parameters = calibration.parameters.ParametersToOptimize(parameters) - -# for inverse_problem in [calibration, validation] - -# all_data = inverse_problem.observations -# simulation = inverse_problem.simulation -# set!(simulation.model, parameters) - -# for data_length in Set(length.(getproperty.(all_data, :t))) - -# observations = [d for d in all_data if length(d.t) == data_length] -# days = observations[1].t[end]/86400 - -# new_ip = InverseProblem() - -# visualize!(simulation, observations, parameters; -# fields = fields, -# filename = joinpath(directory, "$(days)_day_simulations.png")) -# end -# end - -# return nothing -# end diff --git a/examples/calibrate_to_LESbrary/utils/visualize_profile_predictions_utils.jl b/examples/calibrate_to_LESbrary/utils/visualize_profile_predictions_utils.jl deleted file mode 100644 index 5ad3b1e1..00000000 --- a/examples/calibrate_to_LESbrary/utils/visualize_profile_predictions_utils.jl +++ /dev/null @@ -1,47 +0,0 @@ -function open_output_file(directory) - isdir(directory) || mkpath(directory) - file = joinpath(directory, "output.txt") - touch(file) - o = open(file, "w") - return o -end - -function writeout(o, name, loss, params) - param_vect = [params...] - loss_value = loss(params) - write(o, "----------- \n") - write(o, "$(name) \n") - write(o, "Parameters: $(param_vect) \n") - write(o, "Loss: $(loss_value) \n") - saveplot(params, name, loss) -end - -field_guide = Dict( - :u => ( - axis_args = (ylabel="z (m)", xlabel="U velocity (dm/s)"), - scaling = 1e1, - ), - - :v => ( - axis_args = (xlabel="V velocity (dm/s)",), - scaling = 1e1, - ), - - :b => ( - axis_args = (xlabel="Buoyancy (mN/kg)",), - scaling = 1e3, - ), - - :e => ( - axis_args = (ylabel="z (m)", xlabel="TKE (cm²/s²)"), - scaling = 1e4, - ) -) - -function tostring(num) - num == 0 && return "0" - om = Int(floor(log10(abs(num)))) - num /= 10.0^om - num = num%1 ≈ 0 ? Int(num) : round(num; digits=2) - return "$(num)e$om" -end diff --git a/projects/OceanBoundaryLayerParameterizations/Manifest.toml b/projects/OceanBoundaryLayerParameterizations/Manifest.toml deleted file mode 100644 index f45eecff..00000000 --- a/projects/OceanBoundaryLayerParameterizations/Manifest.toml +++ /dev/null @@ -1,2 +0,0 @@ -# This file is machine-generated - editing it directly is not advised - diff --git a/projects/OceanBoundaryLayerParameterizations/Project.toml b/projects/OceanBoundaryLayerParameterizations/Project.toml deleted file mode 100644 index b49b2907..00000000 --- a/projects/OceanBoundaryLayerParameterizations/Project.toml +++ /dev/null @@ -1,4 +0,0 @@ -name = "OceanBoundaryLayerParameterizations" -uuid = "f2a495e6-5045-4354-8fbd-e5ecaaea5da5" -authors = ["Adeline Hillier "] -version = "0.1.0" diff --git a/projects/OceanBoundaryLayerParameterizations/src/EKI_hyperparameter_search.jl b/projects/OceanBoundaryLayerParameterizations/src/EKI_hyperparameter_search.jl deleted file mode 100644 index 578a25ae..00000000 --- a/projects/OceanBoundaryLayerParameterizations/src/EKI_hyperparameter_search.jl +++ /dev/null @@ -1,112 +0,0 @@ - -# -# Hyperparameter search -# - -function loss_reduction(calibration, validation, initial_parameters, kwargs) - params, losses, mean_vars = eki(calibration, initial_parameters; kwargs...) - valid_loss_start = validation(initial_parameters) - valid_loss_final = validation(params) - - @info "parameters: $params" - @info "final loss / initial loss: $(losses[end] / losses[1])" - - train_loss_reduction = losses[end] / losses[1] - valid_loss_reduction = valid_loss_final / valid_loss_start - return train_loss_reduction, valid_loss_reduction -end - -## Stds within bounds - -function plot_stds_within_bounds(calibration, validation, initial_parameters, directory; xrange=-3:0.25:5) - loss_reductions = Dict() - val_loss_reductions = Dict() - for stds_within_bounds in xrange - train_loss_reduction, val_loss_reduction = loss_reduction(calibration, validation, initial_parameters, (stds_within_bounds = stds_within_bounds,)) - loss_reductions[stds_within_bounds] = train_loss_reduction - val_loss_reductions[stds_within_bounds] = val_loss_reduction - end - p = Plots.plot(title="Loss reduction versus prior std's within bounds (n_σ)", ylabel="Loss reduction (Final / Initial)", xlabel="prior std's spanned by bound width (n_σ)", legend=false, lw=3) - plot!(loss_reductions, label="training", color=:purple, lw=4) - plot!(val_loss_reductions, label="validation", color=:blue, lw=4) - Plots.savefig(p, joinpath(directory, "stds_within_bounds.pdf")) - println("loss-minimizing stds within bounds: $(argmin(val_loss_reductions))") -end - -## Prior Variance -function plot_prior_variance(calibration, validation, initial_parameters, directory; xrange=0.1:0.1:1.0) - loss_reductions = Dict() - val_loss_reductions = Dict() - for variance = xrange - train_loss_reduction, val_loss_reduction = loss_reduction(calibration, validation, initial_parameters, (stds_within_bounds = variance,)) - loss_reductions[variance] = train_loss_reduction - val_loss_reductions[variance] = val_loss_reduction - end - p = Plots.plot(title="Loss reduction versus prior variance", ylabel="Loss reduction (Final / Initial)", xlabel="Prior variance") - plot!(loss_reductions, label="training", lw=4, color=:purple) - plot!(val_loss_reductions, label="validation", lw=4, color=:blue) - Plots.savefig(p, joinpath(directory, "variance.pdf")) - v = argmin(var_val_loss_reductions) - println("loss-minimizing variance: $(v)") - return v -end - -## Number of ensemble members - -function plot_num_ensemble_members(calibration, validation, initial_parameters, directory; xrange=1:5:30) - loss_reductions = Dict() - val_loss_reductions = Dict() - for N_ens = xrange - train_loss_reduction, val_loss_reduction = loss_reduction(calibration, validation, initial_parameters, (N_ens = N_ens,)) - loss_reductions[N_ens] = train_loss_reduction - val_loss_reductions[N_ens] = val_loss_reduction - end - p = Plots.plot(title="Loss reduction versus N_ens", xlabel="N_ens", ylabel="Loss reduction (Final / Initial)", legend=false, lw=3) - plot!(loss_reductions, label="training") - plot!(val_loss_reductions, label="validation") - Plots.savefig(p, joinpath(directory, "N_ens.pdf")) -end - -## Observation noise level - -function plot_observation_noise_level(calibration, validation, initial_parameters, directory; xrange=-2.0:0.2:3.0) - loss_reductions = Dict() - val_loss_reductions = Dict() - for log_noise_level = xrange - train_loss_reduction, val_loss_reduction = loss_reduction(calibration, validation, initial_parameters, (noise_level = 10.0^log_noise_level,)) - loss_reductions[log_noise_level] = train_loss_reduction - val_loss_reductions[log_noise_level] = val_loss_reduction - end - p = Plots.plot(title="Loss Reduction versus Observation Noise Level", xlabel="log₁₀(Observation noise level)", ylabel="Loss reduction (Final / Initial)", legend=:topleft) - plot!(loss_reductions, label="training", lw=4, color=:purple) - plot!(val_loss_reductions, label="validation", lw=4, color=:blue) - Plots.savefig(p, joinpath(directory, "obs_noise_level.pdf")) - nl = argmin(val_loss_reductions) - println("loss-minimizing obs noise level: $(nl)") - return 10^nl -end - -function plot_prior_variance_and_obs_noise_level(calibration, validation, initial_parameters, directory; vrange=0.40:0.025:0.90, nlrange=-2.5:0.1:0.5) - Γθs = collect(vrange) - Γys = 10 .^ collect(nlrange) - losses = zeros((length(Γθs), length(Γys))) - counter = 1 - countermax = length(Γθs)*length(Γys) - for i in 1:length(Γθs) - for j in 1:length(Γys) - println("progress $(counter)/$(countermax)") - Γθ = Γθs[i] - Γy = Γys[j] - train_loss_reduction, val_loss_reduction = loss_reduction(calibration, validation, initial_parameters, (stds_within_bounds=Γθ, noise_level=Γy, N_iter=10, N_ens=10, informed_priors=false)) - losses[i, j] = val_loss_reduction - counter += 1 - end - end - p = Plots.heatmap(Γys, Γθs, losses, xlabel=L"\Gamma_y", ylabel=L"\Gamma_\theta", size=(250,250), yscale=:log10) - Plots.savefig(p, joinpath(directory, "GammaHeatmap.pdf")) - v = Γθs[argmin(losses)[1]] - nl = Γys[argmin(losses)[2]] - println("loss-minimizing Γθ: $(v)") - println("loss-minimizing log10(Γy): $(log10(nl))") - return v, nl -end diff --git a/projects/OceanBoundaryLayerParameterizations/src/OceanBoundaryLayerParameterizations.jl b/projects/OceanBoundaryLayerParameterizations/src/OceanBoundaryLayerParameterizations.jl deleted file mode 100644 index b0cd6557..00000000 --- a/projects/OceanBoundaryLayerParameterizations/src/OceanBoundaryLayerParameterizations.jl +++ /dev/null @@ -1,32 +0,0 @@ -module OceanBoundaryLayerParameterizations - -pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", "..")) - -using OceanTurbulenceParameterEstimation -using OceanTurbulenceParameterEstimation.Models: set! -using CairoMakie, LaTeXStrings, OrderedCollections -using Oceananigans, FileIO - -export - # lesbrary_paths.jl - TwoDaySuite, FourDaySuite, SixDaySuite, GeneralStrat, - - # one_dimensional_ensemble_model.jl - OneDimensionalEnsembleModel, - - # EKI_hyperparameter_search.jl - plot_stds_within_bounds, plot_prior_variance, plot_num_ensemble_members, - plot_observation_noise_level, plot_prior_variance_and_obs_noise_level, - - # visualize_profile_predictions.jl - visualize!, - visualize_and_save! - -include("utils/lesbrary_paths.jl") -include("utils/one_dimensional_ensemble_model.jl") -include("EKI_hyperparameter_search.jl") -include("visualize_profile_predictions_utils.jl") -include("utils/visualize_profile_predictions.jl") - - -end # module diff --git a/projects/OceanBoundaryLayerParameterizations/src/visualize_profile_predictions_utils.jl b/projects/OceanBoundaryLayerParameterizations/src/visualize_profile_predictions_utils.jl deleted file mode 100644 index 31fd6dbf..00000000 --- a/projects/OceanBoundaryLayerParameterizations/src/visualize_profile_predictions_utils.jl +++ /dev/null @@ -1,47 +0,0 @@ -function open_output_file(directory) - isdir(directory) || mkpath(directory) - file = joinpath(directory, "output.txt") - touch(file) - o = open(file, "w") - return o -end - -function writeout(o, name, loss, params) - param_vect = [params...] - loss_value = loss(params) - write(o, "----------- \n") - write(o, "$(name) \n") - write(o, "Parameters: $(param_vect) \n") - write(o, "Loss: $(loss_value) \n") - saveplot(params, name, loss) -end - -field_guide = Dict( - :u => ( - axis_args = (ylabel="z (m)", xlabel="U velocity (dm/s)"), - scaling = 1e1, - ), - - :v => ( - axis_args = (xlabel="V velocity (dm/s)",), - scaling = 1e1, - ), - - :b => ( - axis_args = (xlabel="Buoyancy (cN/kg)",), - scaling = 1e2, - ), - - :e => ( - axis_args = (ylabel="z (m)", xlabel="TKE (cm²/s²)"), - scaling = 1e4, - ) -) - -function tostring(num) - num == 0 && return "0" - om = Int(floor(log10(abs(num)))) - num /= 10.0^om - num = num%1 ≈ 0 ? Int(num) : round(num; digits=2) - return "$(num)e$om" -end diff --git a/src/EnsembleKalmanInversions.jl b/src/EnsembleKalmanInversions.jl index c3796c20..6b43cd24 100644 --- a/src/EnsembleKalmanInversions.jl +++ b/src/EnsembleKalmanInversions.jl @@ -205,12 +205,13 @@ Return function iterate!(eki::EnsembleKalmanInversion; iterations = 1, convergence_rate = eki.convergence_rate, - show_progress = true) + show_progress = true, + adaptive_step_parameters = adaptive_step_parameters) iterator = show_progress ? ProgressBar(1:iterations) : 1:iterations for _ in iterator - eki.unconstrained_parameters = step_parameters(eki, convergence_rate) + eki.unconstrained_parameters = step_parameters(eki, convergence_rate, adaptive_step_parameters) eki.iteration += 1 # Forward map @@ -238,7 +239,7 @@ end # it's not adaptive adaptive_step_parameters(::Nothing, Xⁿ, Gⁿ, y, Γy, process) = step_parameters(X, G, y, Γy, process) -function step_parameters(eki::EnsembleKalmanInversion, convergence_rate) +function step_parameters(eki::EnsembleKalmanInversion, convergence_rate, adaptive_step_parameters) process = eki.ensemble_kalman_process y = eki.mapped_observations Γy = eki.noise_covariance diff --git a/src/EnsembleSimulations.jl b/src/EnsembleSimulations.jl index 01d9432e..3e453605 100644 --- a/src/EnsembleSimulations.jl +++ b/src/EnsembleSimulations.jl @@ -64,4 +64,3 @@ function ensemble_column_model_simulation(observations; end end # module - diff --git a/src/InverseProblems.jl b/src/InverseProblems.jl index 61b204bd..0051ad7b 100644 --- a/src/InverseProblems.jl +++ b/src/InverseProblems.jl @@ -34,6 +34,35 @@ using Oceananigans.Grids: Flat, Bounded, using Oceananigans.Models.HydrostaticFreeSurfaceModels: SingleColumnGrid, YZSliceGrid, ColumnEnsembleSize +import ..Transformations: normalize! + +##### +##### Output maps (maps from simulation output to observation space) +##### + +output_map_type(fp) = output_map_str(fp) + +""" + struct ConcatenatedOutputMap end + +Forward map transformation of simulation output to the concatenated +vectors of the simulation output. +""" +struct ConcatenatedOutputMap end + +output_map_str(::ConcatenatedOutputMap) = "ConcatenatedOutputMap" + +""" + struct ConcatenatedVectorNormMap() + +Forward map transformation of simulation output to a scalar by +taking a naive `norm` of the difference between concatenated vectors of the +observations and simulation output. +""" +struct ConcatenatedVectorNormMap end + +output_map_str(::ConcatenatedVectorNormMap) = "ConcatenatedVectorNormMap" + ##### ##### InverseProblems ##### @@ -144,9 +173,12 @@ Nensemble(ip::InverseProblem) = Nensemble(ip.simulation.model.grid) """ observation_map(ip::InverseProblem) -Transform and return `ip.observations` appropriately for `ip.output_map`. +Transform and return `ip.observations` appropriate for `ip.output_map`. """ -observation_map(ip::InverseProblem) = transform_time_series(ip.output_map, ip.observations) +observation_map(ip::InverseProblem) = observation_map(ip.output_map, ip.observations) + +observation_map(map::ConcatenatedOutputMap, observations) = transform_time_series(map, observations) +observation_map(map::ConcatenatedVectorNormMap, observations) = hcat(0.0) """ forward_run!(ip, parameters) @@ -218,11 +250,6 @@ end ##### ConcatenatedOutputMap ##### -# Need docstrings -struct ConcatenatedOutputMap end - -observation_map(map::ConcatenatedOutputMap, observations) = transform_time_series(map, observations) - """ transform_time_series(::ConcatenatedOutputMap, observation::SyntheticObservations) @@ -269,6 +296,18 @@ function transform_forward_map_output(map::ConcatenatedOutputMap, return transform_time_series(map, transposed_forward_map_output) end +function transform_output(output_map::ConcatenatedVectorNormMap, + observations::Union{SyntheticObservations, Vector{<:SyntheticObservations}}, + time_series_collector) + + concat_map = ConcatenatedOutputMap() + fwd_map = transform_output(concat_map, observations, time_series_collector) + obs_map = transform_time_series(concat_map, observations) + + diffn = fwd_map .- obs_map + return sqrt.(mapslices(norm, diffn; dims = 1)) +end + vectorize(observation) = [observation] vectorize(observations::Vector) = observations diff --git a/src/ParameterEstimocean.jl b/src/ParameterEstimocean.jl index b8b1d2b0..d668b1dd 100644 --- a/src/ParameterEstimocean.jl +++ b/src/ParameterEstimocean.jl @@ -16,6 +16,7 @@ export observation_map_variance_across_time, ensemble_column_model_simulation, ConcatenatedOutputMap, + ConcatenatedVectorNormMap, eki, lognormal, ScaledLogitNormal, @@ -23,8 +24,7 @@ export EnsembleKalmanInversion, Resampler, FullEnsembleDistribution, - SuccessfulEnsembleDistribution, - ConstrainedNormal + SuccessfulEnsembleDistribution include("Utils.jl") include("Transformations.jl") diff --git a/src/iteration_summary.jl b/src/iteration_summary.jl index 026b524e..9df46191 100644 --- a/src/iteration_summary.jl +++ b/src/iteration_summary.jl @@ -1,12 +1,52 @@ -struct IterationSummary{P, M, C, V, E} +using ..Parameters: transform_to_unconstrained + +struct IterationSummary{P, M, C, V, E, O} parameters :: P # constrained ensemble_mean :: M # constrained ensemble_cov :: C # constrained ensemble_var :: V mean_square_errors :: E + objective_values :: O iteration :: Int end +""" + eki_objective(eki, θ, G) + +Given forward map `G` and parameters `θ`, return a tuple `(Φ₁, Φ₂)` +of terms in the EKI regularized objective function, where + + Φ = (1/2)*(Φ₁ + Φ₂) + +Φ₁ measures output misfit `|| Γy^(-¹/₂) * (y .- G(θ)) ||²` and +Φ₂ measures prior misfit `|| Γθ^(-¹/₂) * (θ .- μθ) ||²`, where `y` is the observation +map, `G(θ)` is the forward map, `Γy` is the observation noise covariance, `Γθ` is +the prior covariance, and `μθ` represents the prior means. Note that `Γ^(-1/2) = +inv(sqrt(Γ))`. The keyword argument `constrained` is `true` if the input `θ` +represents constrained parameters. +""" +function eki_objective(eki, θ::AbstractVector, G::AbstractVector; constrained = false) + y = eki.mapped_observations + Γy = eki.noise_covariance + + fp = eki.inverse_problem.free_parameters + priors = fp.priors + unconstrained_priors = [unconstrained_prior(priors[name]) for name in fp.names] + μθ = getproperty.(unconstrained_priors, :μ) + Γθ = diagm( getproperty.(unconstrained_priors, :σ).^2 ) + + if constrained + θ = [transform_to_unconstrained(priors[name], θ[i]) + for (i, name) in enumerate(keys(priors))] + end + + # Φ₁ = || Γy^(-½) * (y - G) ||² + Φ₁ = norm(inv(sqrt(Γy)) * (y .- G))^2 + # Φ₂ = || Γθ^(-½) * (θ - μθ) ||² + Φ₂ = norm(inv(sqrt(Γθ)) * (θ .- μθ))^2 + return (Φ₁, Φ₂) +end + """ IterationSummary(eki, X, forward_map_output=nothing) @@ -25,20 +65,24 @@ function IterationSummary(eki, X, forward_map_output=nothing) constrained_parameters = transform_to_constrained(priors, X) + G = forward_map_output if !isnothing(forward_map_output) Nobs, Nens= size(forward_map_output) y = eki.mapped_observations - G = forward_map_output mean_square_errors = [mapreduce((x, y) -> (x - y)^2, +, y, view(G, :, k)) / Nobs for k = 1:Nens] else mean_square_errors = nothing end + # Vector of (Φ₁, Φ₂) pairs, one for each ensemble member at the current iteration + objective_values = [eki_objective(eki, X[:, j], G[:, j]) for j in 1:size(G, 2)] + return IterationSummary(constrained_parameters, constrained_ensemble_mean, constrained_ensemble_covariance, constrained_ensemble_variance, mean_square_errors, + objective_values, eki.iteration) end @@ -89,4 +133,3 @@ particle_str(particle, error, parameters) = @sprintf("% 11s particle: ", particle) * string(param_str.(values(parameters))...) * @sprintf("error = %.6e", error) -