Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Bugfix: cases with different time intervals #190

Merged
merged 4 commits into from
Aug 22, 2024
Merged

Bugfix: cases with different time intervals #190

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
3940cc7
time_bug
crocicc Aug 22, 2024
9ad6f15
changing n_times from int to a vector of ints; this is necessary for …
sajjadazimi Aug 22, 2024
03cc520
minor
sajjadazimi Aug 22, 2024
177d0d8
minor
sajjadazimi Aug 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 12 additions & 11 deletions src/CalibrateCMP/Diagnostics.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,18 +14,18 @@ function lwp(vec::Vector{FT}, config) where {FT <: Real}

(; n_cases, n_heights, n_times) = get_numbers_from_config(config)

lwp = zeros(n_times, n_cases)
ind_ρ = findall(x -> x == "rho", config["observations"]["data_names"])[1]
ind_ql = findall(x -> x == "ql", config["observations"]["data_names"])[1]
dz = (config["model"]["z_max"] - config["model"]["z_min"]) / n_heights[ind_ρ]
n_single_case = sum(n_heights) * n_times

lwp = []
for i in 1:n_cases
_v = get_case_i_vec(vec, i, n_single_case)
_fields = get_single_case_fields(_v, n_heights, n_times)
_fields = get_single_case_fields(_v, n_heights, n_times[i])
ρ = _fields[ind_ρ]
ql = _fields[ind_ql]
lwp[:, i] = sum(ql .* ρ, dims = 1) .* dz
push!(lwp, sum(ql .* ρ, dims = 1) .* dz)
end
return lwp
end
Expand All @@ -47,18 +47,18 @@ function rwp(vec::Vector{FT}, config) where {FT <: Real}

(; n_cases, n_heights, n_times) = get_numbers_from_config(config)

rwp = zeros(n_times, n_cases)
ind_ρ = findall(x -> x == "rho", config["observations"]["data_names"])[1]
ind_qr = findall(x -> x == "qr", config["observations"]["data_names"])[1]
dz = (config["model"]["z_max"] - config["model"]["z_min"]) / n_heights[ind_ρ]
n_single_case = sum(n_heights) * n_times

rwp = []
for i in 1:n_cases
_v = get_case_i_vec(vec, i, n_single_case)
_fields = get_single_case_fields(_v, n_heights, n_times)
_fields = get_single_case_fields(_v, n_heights, n_times[i])
ρ = _fields[ind_ρ]
qr = _fields[ind_qr]
rwp[:, i] = sum(qr .* ρ, dims = 1) .* dz
push!(rwp, sum(qr .* ρ, dims = 1) .* dz)
end
return rwp
end
Expand All @@ -81,7 +81,6 @@ function rainrate(vec::Vector{FT}, config; height::FT = FT(0)) where {FT <: Real

(; n_cases, n_heights, n_times) = get_numbers_from_config(config)

rainrate = zeros(n_times, n_cases)
ind_rr = findall(x -> x == "rainrate", config["observations"]["data_names"])[1]
dz = (config["model"]["z_max"] - config["model"]["z_min"]) / n_heights[ind_rr]
n_single_case = sum(n_heights) * n_times
Expand All @@ -91,12 +90,14 @@ function rainrate(vec::Vector{FT}, config; height::FT = FT(0)) where {FT <: Real
a1 = ((ind_z2 - 0.5) * dz - height) / dz
a2 = (height - (ind_z1 - 0.5) * dz) / dz

rainrate = []
for i in 1:n_cases
_v = get_case_i_vec(vec, i, n_single_case)
_fields = get_single_case_fields(_v, n_heights, n_times)
_rainrate = _fields[ind_rr]
rainrate[:, i] = a1 .* _rainrate[ind_z1, :] + a2 .* _rainrate[ind_z2, :]
_fields = get_single_case_fields(_v, n_heights, n_times[i])
_rainrate = a1 .* _fields[ind_rr][ind_z1, :] + a2 .* _fields[ind_rr][ind_z2, :]
_rainrate[findall(x -> x < 0, _rainrate)] .= FT(0)
push!(rainrate, _rainrate)
end
rainrate[findall(x -> x < 0, rainrate)] .= FT(0)

return rainrate
end
17 changes: 13 additions & 4 deletions src/CalibrateCMP/HelperFuncs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,15 +50,24 @@ function get_numbers_from_config(config::Dict)
n_heights =
config["model"]["filter"]["apply"] ? config["model"]["filter"]["nz_filtered"] :
config["model"]["filter"]["nz_unfiltered"]
n_times =
_n_times_default =
config["model"]["filter"]["apply"] ? length(config["model"]["t_calib"]) - 1 : length(config["model"]["t_calib"])
n_times = Int[]
for case in config["observations"]["cases"]
if :t_cal in collect(keys(case))
_n_times_case_i = config["model"]["filter"]["apply"] ? length(case.t_cal) - 1 : length(case.t_cal)
else
_n_times_case_i = _n_times_default
end
n_times = [n_times; _n_times_case_i]
end

@assert length(n_heights) == length(config["observations"]["data_names"])
return (; n_cases, n_heights, n_times)
end

function get_case_i_vec(vec::Vector{FT}, i::Int, n_single_case::Int) where {FT <: Real}
return vec[((i - 1) * n_single_case + 1):(i * n_single_case)]
function get_case_i_vec(vec::Vector{FT}, i::Int, n_single_case::Vector{Int}) where {FT <: Real}
return vec[(sum(n_single_case[1:(i - 1)]) + 1):sum(n_single_case[1:i])]
end

function get_single_case_fields(vec_single_case::Vector{FT}, n_heights::Vector{Int}, n_times::Int) where {FT <: Real}
Expand Down Expand Up @@ -147,7 +156,7 @@ function compute_error_metrics(
config::Dict,
obs::Matrix{FT},
) where {FT <: Real}
n_cases = n_cases = length(config["observations"]["cases"])
n_cases = length(config["observations"]["cases"])

ref_stats_list = make_ref_stats_list(obs, config["statistics"], get_numbers_from_config(config)...)
model_error = zeros(3)
Expand Down
10 changes: 5 additions & 5 deletions src/CalibrateCMP/IOUtils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -206,8 +206,8 @@ function compare_model_and_obs_contours(
for case_num in 1:n_cases
_v_model = get_case_i_vec(G, case_num, n_single_case)
_v_obs = get_case_i_vec(obs, case_num, n_single_case)
_fields_model = get_single_case_fields(_v_model, n_heights, n_times)
_fields_obs = get_single_case_fields(_v_obs, n_heights, n_times)
_fields_model = get_single_case_fields(_v_model, n_heights, n_times[case_num])
_fields_obs = get_single_case_fields(_v_obs, n_heights, n_times[case_num])

p = Array{Plots.Plot}(undef, n_variables * 2)

Expand Down Expand Up @@ -360,9 +360,9 @@ function plot_box_results(
_v_model = get_case_i_vec(G, case_num, n_single_case)
_v_obs_mean = get_case_i_vec(obs_mean, case_num, n_single_case)
_v_obs_std = get_case_i_vec(obs_std, case_num, n_single_case)
_fields_model = get_single_case_fields(_v_model, n_heights, n_times)
_fields_obs_mean = get_single_case_fields(_v_obs_mean, n_heights, n_times)
_fields_obs_std = get_single_case_fields(_v_obs_std, n_heights, n_times)
_fields_model = get_single_case_fields(_v_model, n_heights, n_times[case_num])
_fields_obs_mean = get_single_case_fields(_v_obs_mean, n_heights, n_times[case_num])
_fields_obs_std = get_single_case_fields(_v_obs_std, n_heights, n_times[case_num])

p = Array{Plots.Plot}(undef, n_variables)

Expand Down
15 changes: 11 additions & 4 deletions src/CalibrateCMP/KiDUtils.jl
View file
Open in desktop
sajjadazimi marked this conversation as resolved.
Show resolved Hide resolved
Original file line number Diff line number Diff line change
Expand Up @@ -43,14 +43,22 @@ function run_KiD_multiple_cases(u::Array{FT, 1}, u_names::Array{String, 1}, conf
config["model"]["w1"] = case.w1
config["model"]["p0"] = case.p0
config["model"]["Nd"] = case.Nd
if "t_cal" in collect(keys(case))
if :t_cal in collect(keys(case))
config["model"]["filter"] = make_filter_props(
config["model"]["filter"]["nz_unfiltered"],
case.t_cal;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here, I think we need to add

        else
            config["model"]["filter"] = make_filter_props(
                config["model"]["filter"]["nz_unfiltered"],
                config["model"]["t_calib"];
                apply = config["model"]["filter"]["apply"],
                nz_per_filtered_cell = config["model"]["filter"]["nz_per_filtered_cell"],
                nt_per_filtered_cell = config["model"]["filter"]["nt_per_filtered_cell"],
            )  
        end

If we check with an example in which the 1st case has the t_cal flag, while the second case doesn't have the t_cal flag, then for the case without the flag it keeps the last filter (the one created for the case with t_cal), not the original one.

Let me know if you agree and if it makes sense to you!

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doesn't it make a new filter for every case in this way (not considering the default one anymore)?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree and it makes total sense. Thank you!

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can check with different combinations of cases with and without t_cal, but apart from that I think it looks good!

apply = config["model"]["filter"]["apply"],
nz_per_filtered_cell = config["model"]["filter"]["nz_per_filtered_cell"],
nt_per_filtered_cell = config["model"]["filter"]["nt_per_filtered_cell"],
)
else
config["model"]["filter"] = make_filter_props(
config["model"]["filter"]["nz_unfiltered"],
config["model"]["t_calib"];
apply = config["model"]["filter"]["apply"],
nz_per_filtered_cell = config["model"]["filter"]["nz_per_filtered_cell"],
nt_per_filtered_cell = config["model"]["filter"]["nt_per_filtered_cell"],
)
end

ode_sol, aux, precip = run_KiD(u, u_names, config["model"])
Expand Down Expand Up @@ -434,10 +442,9 @@ function test_model(u::Array{FT, 1}, u_names::Array{String, 1}, config::Dict, ca
b = u[findall(name -> name == "b", u_names)]

n_single_case = sum(n_heights) * n_times
x = range(0.0, 1.0, n_single_case)

outputs = Float64[]
for case in config["observations"]["cases"][case_numbers]
for (i, case) in enumerate(config["observations"]["cases"][case_numbers])
x = range(0.0, 1.0, n_single_case[case_numbers[i]])
p = case.power
y = a .* x .^ p .+ b
outputs = [outputs; y]
Expand Down
9 changes: 5 additions & 4 deletions src/CalibrateCMP/ReferenceModels.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,11 +41,12 @@ function get_obs_matrix(
for (i, case) in enumerate(cases)
_dir::String = case.dir

if "t_cal" in collect(keys(case))
times = case.t_cal
if :t_cal in collect(keys(case))
_times = collect(case.t_cal)
else
_times = times
end
_data_matrix_single_case::Matrix{FT} =
get_obs_matrix(_dir, variables, heights, times; apply_filter = apply_filter)
_data_matrix_single_case = get_obs_matrix(_dir, variables, heights, _times; apply_filter = apply_filter)

if i == 1
_data_matrix = _data_matrix_single_case
Expand Down
51 changes: 28 additions & 23 deletions src/CalibrateCMP/ReferenceStats.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ Base.@kwdef struct ReferenceStatistics{FT <: Real}
centered::Bool
case_numbers::Vector{Int}
n_heights::Vector{Int}
n_times::Int
n_times::Vector{Int}


function ReferenceStatistics(
Expand All @@ -28,7 +28,7 @@ Base.@kwdef struct ReferenceStatistics{FT <: Real}
centered::Bool,
case_numbers::Vector{Int},
n_heights::Vector{Int},
n_times::Int,
n_times::Vector{Int},
)
return new{FT}(
y,
Expand Down Expand Up @@ -59,7 +59,7 @@ Base.@kwdef struct ReferenceStatistics{FT <: Real}
@assert sum(n_heights) * n_times == size(obs)[1]
case_numbers = [case_number]
_n_cases = 1
(var_mean_max, var_std_max) = find_mean_and_std_maximum(obs, _n_cases, n_heights, n_times)
(var_mean_max, var_std_max) = find_mean_and_std_maximum(obs, _n_cases, n_heights, [n_times])
if stats_config["normalization"] == "mean_normalized"
y_norm = var_mean_max
centered = false
Expand All @@ -73,7 +73,7 @@ Base.@kwdef struct ReferenceStatistics{FT <: Real}
weights =
"weights" in keys(stats_config) ? reshape(stats_config["weights"], 1, _n_variables) : ones(1, _n_variables)
y_norm = y_norm ./ weights
obs_normalized = normalize_obs(obs, y_norm, _n_cases, n_heights, n_times, centered = centered)
obs_normalized = normalize_obs(obs, y_norm, _n_cases, n_heights, [n_times], centered = centered)
obs_mean = mean(obs, dims = 2)[:]
y_full = mean(obs_normalized, dims = 2)[:]
Γ_full = cov(obs_normalized, dims = 2)
Expand Down Expand Up @@ -109,25 +109,30 @@ Base.@kwdef struct ReferenceStatistics{FT <: Real}
centered,
case_numbers,
n_heights,
n_times,
[n_times],
)
end

end

function find_mean_and_std_maximum(y::Matrix{FT}, n_cases::Int, n_heights::Vector{Int}, n_times::Int) where {FT <: Real}
function find_mean_and_std_maximum(
y::Matrix{FT},
n_cases::Int,
n_heights::Vector{Int},
n_times::Vector{Int},
) where {FT <: Real}

_n_variables = length(n_heights)
_n_single_case_sim = sum(n_heights) * n_times
@assert size(y)[1] == (n_cases * _n_single_case_sim)
@assert size(y)[1] == sum(_n_single_case_sim)

var_mean_max = ones(n_cases, _n_variables)
var_std_max = ones(n_cases, _n_variables)
for i in 1:n_cases
_y_single_case = y[((i - 1) * _n_single_case_sim + 1):(i * _n_single_case_sim), :]
_y_single_case = y[(sum(_n_single_case_sim[1:(i - 1)]) + 1):sum(_n_single_case_sim[1:i]), :]

_y_mean_vector_single_case = mean(_y_single_case, dims = 2)
_y_mean_matrix::Matrix{FT} = reshape(_y_mean_vector_single_case, sum(n_heights), n_times)
_y_mean_matrix::Matrix{FT} = reshape(_y_mean_vector_single_case, sum(n_heights), n_times[i])
_var_mean_max = [
maximum(abs.(_y_mean_matrix[(sum(n_heights[1:(j - 1)]) + 1):sum(n_heights[1:j]), :])) for
j in 1:_n_variables
Expand All @@ -136,7 +141,7 @@ function find_mean_and_std_maximum(y::Matrix{FT}, n_cases::Int, n_heights::Vecto
var_mean_max[i, :] = _var_mean_max

_y_std_vector_single_case = std(_y_single_case, dims = 2)
_y_std_matrix::Matrix{FT} = reshape(_y_std_vector_single_case, sum(n_heights), n_times)
_y_std_matrix::Matrix{FT} = reshape(_y_std_vector_single_case, sum(n_heights), n_times[i])
_var_std_max = [
maximum(abs.(_y_std_matrix[(sum(n_heights[1:(j - 1)]) + 1):sum(n_heights[1:j]), :])) for j in 1:_n_variables
]
Expand All @@ -152,13 +157,13 @@ function normalize_obs(
ynorm::Matrix{FT},
n_cases::Int,
n_heights::Vector{Int},
n_times::Int;
n_times::Vector{Int};
centered::Bool = true,
) where {FT <: Real}

_n_variables = length(n_heights)
_n_single_case_sim = sum(n_heights) * n_times
@assert size(y)[1] == (n_cases * _n_single_case_sim)
@assert size(y)[1] == sum(_n_single_case_sim)
@assert size(ynorm) == (n_cases, _n_variables)

_y_normalized = centered ? y .- mean(y, dims = 2) : copy(y)
Expand All @@ -169,8 +174,8 @@ function normalize_obs(
for j in 1:_n_variables
_norm_vec_extended_single_time[(sum(n_heights[1:(j - 1)]) + 1):sum(n_heights[1:j])] .= _norm_vec[j]
end
_norm_vec_extended = reshape(ones(sum(n_heights), n_times) .* _norm_vec_extended_single_time, :, 1)
_row_indeces = ((i - 1) * _n_single_case_sim + 1):(i * _n_single_case_sim)
_norm_vec_extended = reshape(ones(sum(n_heights), n_times[i]) .* _norm_vec_extended_single_time, :, 1)
_row_indeces = (sum(_n_single_case_sim[1:(i - 1)]) + 1):sum(_n_single_case_sim[1:i])
_y_normalized[_row_indeces, :] = _y_normalized[_row_indeces, :] ./ _norm_vec_extended
end

Expand Down Expand Up @@ -198,8 +203,8 @@ function normalize_sim(
for j in 1:_n_variables
_norm_vec_extended_single_time[(sum(_n_heights[1:(j - 1)]) + 1):sum(_n_heights[1:j])] .= _norm_vec[j]
end
_norm_vec_extended = reshape(ones(sum(_n_heights), _n_times) .* _norm_vec_extended_single_time, :, 1)
_indeces = ((i - 1) * _n_single_case_sim + 1):(i * _n_single_case_sim)
_norm_vec_extended = reshape(ones(sum(_n_heights), _n_times[i]) .* _norm_vec_extended_single_time, :, 1)
_indeces = (sum(_n_single_case_sim[1:(i - 1)]) + 1):sum(_n_single_case_sim[1:i])
_sim_vec_normalized[_indeces] = _sim_vec_normalized[_indeces] ./ _norm_vec_extended
end

Expand Down Expand Up @@ -227,8 +232,8 @@ function unnormalize_sim(
for j in 1:_n_variables
_norm_vec_extended_single_time[(sum(_n_heights[1:(j - 1)]) + 1):sum(_n_heights[1:j])] .= _norm_vec[j]
end
_norm_vec_extended = reshape(ones(sum(_n_heights), _n_times) .* _norm_vec_extended_single_time, :, 1)
_indeces = ((i - 1) * _n_single_case_sim + 1):(i * _n_single_case_sim)
_norm_vec_extended = reshape(ones(sum(_n_heights), _n_times[i]) .* _norm_vec_extended_single_time, :, 1)
_indeces = (sum(_n_single_case_sim[1:(i - 1)]) + 1):sum(_n_single_case_sim[1:i])
_sim_vec_unnormalized[_indeces] = _sim_vec_unnormalized[_indeces] .* _norm_vec_extended
end
if ref_stats.centered
Expand Down Expand Up @@ -270,19 +275,19 @@ function make_ref_stats_list(
stats_config::Dict,
n_cases::Int,
n_heights::Vector{Int},
n_times::Int,
n_times::Vector{Int},
) where {FT <: Real}

_n_single_case_sim = sum(n_heights) * n_times
@assert n_cases * _n_single_case_sim == size(obs)[1]
@assert sum(_n_single_case_sim) == size(obs)[1]
ref_stats_list = ReferenceStatistics[]
for i in 1:n_cases
_single_case_row_indices = ((i - 1) * _n_single_case_sim + 1):(i * _n_single_case_sim)
_single_case_row_indices = (sum(_n_single_case_sim[1:(i - 1)]) + 1):sum(_n_single_case_sim[1:i])
_obs_single_case = obs[_single_case_row_indices, :]

ref_stats_list = [
ref_stats_list
ReferenceStatistics(_obs_single_case, stats_config, i, n_heights, n_times)
ReferenceStatistics(_obs_single_case, stats_config, i, n_heights, n_times[i])
]
end
return ref_stats_list
Expand All @@ -308,7 +313,7 @@ function combine_ref_stats(ref_stats_list::Vector{ReferenceStatistics})
for ref_stats in ref_stats_list[2:end]
@assert centered == ref_stats.centered
@assert n_heights == ref_stats.n_heights
@assert n_times == ref_stats.n_times
n_times = [n_times; ref_stats.n_times]
case_numbers = [case_numbers; ref_stats.case_numbers]
y_norm = [y_norm; ref_stats.y_norm]
var_mean_max = [var_mean_max; ref_stats.var_mean_max]
Expand Down
19 changes: 9 additions & 10 deletions test/unit_tests/calibration_pipeline/test_diagnostics.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,24 +6,23 @@
KCP.make_filter_props(config["model"]["n_elem"] .* ones(Int, 4), config["model"]["t_calib"])
(n_c, n_z, n_t) = KCP.get_numbers_from_config(config)
n_single_case = sum(n_z) * n_t
vec = rand(n_c * n_single_case)
vec = rand(sum(n_single_case))
dz = (config["model"]["z_max"] - config["model"]["z_min"]) / n_z[1]
_lwp = zeros(n_t, n_c)
_rwp = zeros(n_t, n_c)
_rainrate = zeros(n_t, n_c)
_lwp = []
_rwp = []
_rainrate = []

for i in 1:n_c
v_ = KCP.get_case_i_vec(vec, i, n_single_case)
fields_ = KCP.get_single_case_fields(v_, n_z, n_t)
fields_ = KCP.get_single_case_fields(v_, n_z, n_t[i])
ρ = fields_[1]
ql = fields_[2]
qr = fields_[3]
_lwp[:, i] = sum(ql .* ρ, dims = 1) .* dz
_rwp[:, i] = sum(qr .* ρ, dims = 1) .* dz
_rr = fields_[4]
rainrate_0 = 1.5 * _rr[1, :] - 0.5 * _rr[2, :]
push!(_lwp, sum(ql .* ρ, dims = 1) .* dz)
push!(_rwp, sum(qr .* ρ, dims = 1) .* dz)
rainrate_0 = 1.5 * fields_[4][1, :] - 0.5 * fields_[4][2, :]
rainrate_0[findall(x -> x < 0, rainrate_0)] .= Float64(0)
_rainrate[:, i] = rainrate_0
push!(_rainrate, rainrate_0)
end

#action
Expand Down
Loading
Loading