Bugfix: cases with different time intervals

src/CalibrateCMP/Diagnostics.jl

@@ -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
@@ -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
@@ -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
@@ -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

src/CalibrateCMP/HelperFuncs.jl

@@ -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}
@@ -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)

src/CalibrateCMP/KiDUtils.jl

@@ -434,10 +434,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]

src/CalibrateCMP/ReferenceModels.jl

@@ -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

src/CalibrateCMP/ReferenceStats.jl

@@ -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(
@@ -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,
@@ -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
@@ -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)
@@ -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
@@ -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
         ]
@@ -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)
@@ -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
 
@@ -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
 
@@ -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
@@ -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
@@ -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]

test/unit_tests/calibration_pipeline/test_diagnostics.jl

@@ -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

test/unit_tests/calibration_pipeline/test_helper_funcs.jl

@@ -45,7 +45,7 @@ end
 
     #test
     @test length(nums) == 3
-    @test typeof(nums) == @NamedTuple{n_cases::Int, n_heights::Vector{Int}, n_times::Int}
+    @test typeof(nums) == @NamedTuple{n_cases::Int, n_heights::Vector{Int}, n_times::Vector{Int}}
     @test length(nums.n_heights) == length(config["observations"]["data_names"])
 end
 
@@ -57,11 +57,11 @@ end
         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))
 
     #action
     single_case_vec = KCP.get_case_i_vec(vec, 1, n_single_case)
-    fields = KCP.get_single_case_fields(single_case_vec, n_z, n_t)
+    fields = KCP.get_single_case_fields(single_case_vec, n_z, n_t[1])
 
     #test
     length(single_case_vec) == n_single_case