From 9fc694bd558237d37d9c877331e2939e2c61ac03 Mon Sep 17 00:00:00 2001 From: Thibault Hallouin <thibault.hallouin@inrae.fr> Date: Tue, 27 Dec 2022 12:13:10 +0100 Subject: [PATCH] add missing opening/closing curly brackets in for/if-else structures --- .../evalhyd/detail/determinist/evaluator.hpp | 69 +++++++---- include/evalhyd/detail/masks.hpp | 108 ++++++++++++++++-- .../evalhyd/detail/probabilist/evaluator.hpp | 42 ++++--- include/evalhyd/detail/uncertainty.hpp | 32 ++++-- include/evalhyd/detail/utils.hpp | 20 ++++ include/evalhyd/evald.hpp | 54 ++++++++- include/evalhyd/evalp.hpp | 74 +++++++++++- tests/test_determinist.cpp | 6 +- tests/test_probabilist.cpp | 3 +- 9 files changed, 341 insertions(+), 67 deletions(-) diff --git a/include/evalhyd/detail/determinist/evaluator.hpp b/include/evalhyd/detail/determinist/evaluator.hpp index bd2b745..71a93c1 100644 --- a/include/evalhyd/detail/determinist/evaluator.hpp +++ b/include/evalhyd/detail/determinist/evaluator.hpp @@ -74,7 +74,8 @@ namespace evalhyd t_msk = xt::ones<bool>({n_msk, n_srs, n_tim}); xt::view(t_msk, xt::all()) = xt::view(msk, xt::all(), xt::newaxis(), xt::all()); - for (std::size_t m = 0; m < n_msk; m++) { + for (std::size_t m = 0; m < n_msk; m++) + { xt::view(t_msk, m) = xt::where(obs_nan | prd_nan, false, xt::view(t_msk, m)); @@ -116,13 +117,15 @@ namespace evalhyd void Evaluator<D2, B2>::calc_mean_obs() { mean_obs = xt::zeros<double>(mean_dims); - for (std::size_t m = 0; m < n_msk; m++) { + for (std::size_t m = 0; m < n_msk; m++) + { // apply the mask // (using NaN workaround until reducers work on masked_view) auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN); // compute variable one sample at a time - for (std::size_t e = 0; e < n_exp; e++) { + for (std::size_t e = 0; e < n_exp; e++) + { // apply the bootstrap sampling auto obs = xt::view(obs_masked, xt::all(), b_exp[e]); xt::view(mean_obs, m, e) = @@ -143,13 +146,15 @@ namespace evalhyd void Evaluator<D2, B2>::calc_mean_prd() { mean_prd = xt::zeros<double>(mean_dims); - for (std::size_t m = 0; m < n_msk; m++) { + for (std::size_t m = 0; m < n_msk; m++) + { // apply the mask // (using NaN workaround until reducers work on masked_view) auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN); // compute variable one sample at a time - for (std::size_t e = 0; e < n_exp; e++) { + for (std::size_t e = 0; e < n_exp; e++) + { // apply the bootstrap sampling auto prd = xt::view(prd_masked, xt::all(), b_exp[e]); xt::view(mean_prd, m, e) = @@ -190,12 +195,14 @@ namespace evalhyd void Evaluator<D2, B2>::calc_quad_obs() { quad_obs = xt::zeros<double>(inter_dims); - for (std::size_t m = 0; m < n_msk; m++) { + for (std::size_t m = 0; m < n_msk; m++) + { // apply the mask // (using NaN workaround until reducers work on masked_view) auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN); - for (std::size_t e = 0; e < n_exp; e++) { + for (std::size_t e = 0; e < n_exp; e++) + { xt::view(quad_obs, m, e) = xt::square( obs_masked - xt::view(mean_obs, m, e) ); @@ -218,12 +225,14 @@ namespace evalhyd void Evaluator<D2, B2>::calc_quad_prd() { quad_prd = xt::zeros<double>(inter_dims); - for (std::size_t m = 0; m < n_msk; m++) { + for (std::size_t m = 0; m < n_msk; m++) + { // apply the mask // (using NaN workaround until reducers work on masked_view) auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN); - for (std::size_t e = 0; e < n_exp; e++) { + for (std::size_t e = 0; e < n_exp; e++) + { xt::view(quad_prd, m, e) = xt::square( prd_masked - xt::view(mean_prd, m, e) ); @@ -260,14 +269,16 @@ namespace evalhyd // calculate error in timing and dynamics $r_{pearson}$ // (Pearson's correlation coefficient) r_pearson = xt::zeros<double>(inner_dims); - for (std::size_t m = 0; m < n_msk; m++) { + for (std::size_t m = 0; m < n_msk; m++) + { // apply the mask // (using NaN workaround until reducers work on masked_view) auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN); auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN); // compute variable one sample at a time - for (std::size_t e = 0; e < n_exp; e++) { + for (std::size_t e = 0; e < n_exp; e++) + { auto prd = xt::view(prd_masked, xt::all(), b_exp[e]); auto obs = xt::view(obs_masked, xt::all(), b_exp[e]); auto r_num = xt::nansum( @@ -309,14 +320,16 @@ namespace evalhyd { // calculate error in spread of flow $alpha$ alpha = xt::zeros<double>(inner_dims); - for (std::size_t m = 0; m < n_msk; m++) { + for (std::size_t m = 0; m < n_msk; m++) + { // apply the mask // (using NaN workaround until reducers work on masked_view) auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN); auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN); // compute variable one sample at a time - for (std::size_t e = 0; e < n_exp; e++) { + for (std::size_t e = 0; e < n_exp; e++) + { auto prd = xt::view(prd_masked, xt::all(), b_exp[e]); auto obs = xt::view(obs_masked, xt::all(), b_exp[e]); xt::view(alpha, m, e) = @@ -342,14 +355,16 @@ namespace evalhyd { // calculate $bias$ bias = xt::zeros<double>(inner_dims); - for (std::size_t m = 0; m < n_msk; m++) { + for (std::size_t m = 0; m < n_msk; m++) + { // apply the mask // (using NaN workaround until reducers work on masked_view) auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN); auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN); // compute variable one sample at a time - for (std::size_t e = 0; e < n_exp; e++) { + for (std::size_t e = 0; e < n_exp; e++) + { auto prd = xt::view(prd_masked, xt::all(), b_exp[e]); auto obs = xt::view(obs_masked, xt::all(), b_exp[e]); xt::view(bias, m, e) = @@ -371,13 +386,15 @@ namespace evalhyd { // compute RMSE RMSE = xt::zeros<double>(final_dims); - for (std::size_t m = 0; m < n_msk; m++) { + for (std::size_t m = 0; m < n_msk; m++) + { // apply the mask // (using NaN workaround until reducers work on masked_view) auto quad_err_masked = xt::where(xt::view(t_msk, m), quad_err, NAN); // compute variable one sample at a time - for (std::size_t e = 0; e < n_exp; e++) { + for (std::size_t e = 0; e < n_exp; e++) + { auto err2 = xt::view(quad_err_masked, xt::all(), b_exp[e]); xt::view(RMSE, xt::all(), m, e) = xt::sqrt(xt::nanmean(err2, -1)); } @@ -399,13 +416,15 @@ namespace evalhyd void Evaluator<D2, B2>::calc_NSE() { NSE = xt::zeros<double>(final_dims); - for (std::size_t m = 0; m < n_msk; m++) { + for (std::size_t m = 0; m < n_msk; m++) + { // apply the mask // (using NaN workaround until reducers work on masked_view) auto quad_err_masked = xt::where(xt::view(t_msk, m), quad_err, NAN); // compute variable one sample at a time - for (std::size_t e = 0; e < n_exp; e++) { + for (std::size_t e = 0; e < n_exp; e++) + { // compute squared errors operands auto err2 = xt::view(quad_err_masked, xt::all(), b_exp[e]); xt::xtensor<double, 1> f_num = @@ -438,8 +457,10 @@ namespace evalhyd void Evaluator<D2, B2>::calc_KGE() { KGE = xt::zeros<double>(final_dims); - for (std::size_t m = 0; m < n_msk; m++) { - for (std::size_t e = 0; e < n_exp; e++) { + for (std::size_t m = 0; m < n_msk; m++) + { + for (std::size_t e = 0; e < n_exp; e++) + { // compute KGE xt::view(KGE, xt::all(), m, e) = 1 - xt::sqrt( xt::square(xt::view(r_pearson, m, e) - 1) @@ -474,8 +495,10 @@ namespace evalhyd void Evaluator<D2, B2>::calc_KGEPRIME() { KGEPRIME = xt::zeros<double>(final_dims); - for (std::size_t m = 0; m < n_msk; m++) { - for (std::size_t e = 0; e < n_exp; e++) { + for (std::size_t m = 0; m < n_msk; m++) + { + for (std::size_t e = 0; e < n_exp; e++) + { // calculate error in spread of flow $gamma$ auto gamma = xt::view(alpha, m, e) * (xt::view(mean_obs, m, e, xt::all(), 0) diff --git a/include/evalhyd/detail/masks.hpp b/include/evalhyd/detail/masks.hpp index d4a1c6b..46313e5 100644 --- a/include/evalhyd/detail/masks.hpp +++ b/include/evalhyd/detail/masks.hpp @@ -62,32 +62,42 @@ namespace evalhyd std::set<std::string> supported_op = {"<", ">", "<=", ">=", "!=", "=="}; if (mt[1].str().empty()) + { throw std::runtime_error( "missing operator for streamflow masking condition" ); + } else if (supported_op.find(mt[1]) != supported_op.end()) { if ((mt[2].str() == "median") || (mt[2].str() == "mean") || (mt[2].str() == "quantile")) + { conditions.push_back({mt[1].str(), mt[2].str(), mt[3].str()}); + } else + { // it is a simple numerical value, swap last two conditions.push_back({mt[1].str(), mt[3].str(), mt[2].str()}); + } } else + { throw std::runtime_error( "invalid operator for streamflow masking " "condition: " + mt[1].str() ); + } } // check that a maximum of two conditions were provided if (conditions.size() > 2) + { throw std::runtime_error( "no more than two streamflow masking conditions " "can be provided" ); + } subset[var] = conditions; } @@ -119,11 +129,15 @@ namespace evalhyd // check whether it is all indices, a range of indices, or an index if (m[1] == ":") + { // it is all indices (i.e. t{:}) so keep everything condition.emplace_back(); + } else if (m[2].str().empty()) + { // it is an index (i.e. t{#}) condition.push_back({m[1].str()}); + } else { // it is a range of indices (i.e. t{#:#}) @@ -160,9 +174,11 @@ namespace evalhyd // check if conditions were found in parsing if (subset.empty()) + { throw std::runtime_error( "no valid condition found to generate mask(s)" ); + } // initialise a boolean expression for the masks xt::xtensor<bool, 1> t_msk = xt::zeros<bool>(q_obs.shape()); @@ -179,25 +195,33 @@ namespace evalhyd { // preprocess streamflow depending on kind auto get_q = [&]() { - if (var == "q_obs") { + if (var == "q_obs") + { return q_obs; } - else if (var == "q_prd_median") { + else if (var == "q_prd_median") + { if (q_prd.size() < 1) + { throw std::runtime_error( "condition on streamflow predictions " "not allowed for generating masks" ); + } xt::xtensor<double, 1> q_prd_median = xt::median(q_prd, 0); return q_prd_median; } - else { // i.e. (var == "q_prd_mean") + else + { + // i.e. (var == "q_prd_mean") if (q_prd.size() < 1) + { throw std::runtime_error( "condition on streamflow predictions " "not allowed for generating masks" ); + } xt::xtensor<double, 1> q_prd_mean = xt::mean(q_prd, 0); return q_prd_mean; @@ -206,16 +230,26 @@ namespace evalhyd auto q = get_q(); // define lambda function to precompute mean/median/quantile - auto get_val = [&](const std::string& str, const std::string& num) { - if (str.empty()) - // it is a simple numerical value + auto get_val = + [&](const std::string& str, const std::string& num) + { + if (str.empty()) // it is a simple numerical value + { return std::stod(num); + } else if (str == "median") + { return xt::median(q); + } else if (str == "mean") + { return xt::mean(q)(); + } else // (str == "quantile") + { return eh::maths::quantile(q, std::stod(num)); + } + }; // preprocess conditions to identify special cases @@ -238,8 +272,13 @@ namespace evalhyd if ((opr2 == ">") || (opr2 == ">=")) { if (val2 > val1) + { without = true; - else { within = true; } + } + else + { + within = true; + } } } else if ((opr1 == ">") || (opr1 == ">=")) @@ -247,8 +286,13 @@ namespace evalhyd if ((opr2 == "<") || (opr2 == "<=")) { if (val2 > val1) + { within = true; - else { without = true; } + } + else + { + without = true; + } } } } @@ -258,58 +302,90 @@ namespace evalhyd if (within) { if ((opr1 == "<") && (opr2 == ">")) + { t_msk = xt::where((q < val1) & (q > val2), 1, t_msk); + } else if ((opr1 == "<=") && (opr2 == ">")) + { t_msk = xt::where((q <= val1) & (q > val2), 1, t_msk); + } else if ((opr1 == "<") && (opr2 == ">=")) + { t_msk = xt::where((q < val1) & (q >= val2), 1, t_msk); + } else if ((opr1 == "<=") && (opr2 == ">=")) + { t_msk = xt::where((q <= val1) & (q >= val2), 1, t_msk); + } if ((opr2 == "<") && (opr1 == ">")) + { t_msk = xt::where((q < val2) & (q > val1), 1, t_msk); + } else if ((opr2 == "<=") && (opr1 == ">")) + { t_msk = xt::where((q <= val2) & (q > val1), 1, t_msk); + } else if ((opr2 == "<") && (opr1 == ">=")) + { t_msk = xt::where((q < val2) & (q >= val1), 1, t_msk); + } else if ((opr2 == "<=") && (opr1 == ">=")) + { t_msk = xt::where((q <= val2) & (q >= val1), 1, t_msk); + } } else if (without) { if ((opr1 == "<") && (opr2 == ">")) + { t_msk = xt::where((q < val1) | (q > val2), 1, t_msk); + } else if ((opr1 == "<=") && (opr2 == ">")) + { t_msk = xt::where((q <= val1) | (q > val2), 1, t_msk); + } else if ((opr1 == "<") && (opr2 == ">=")) + { t_msk = xt::where((q < val1) | (q >= val2), 1, t_msk); + } else if ((opr1 == "<=") && (opr2 == ">=")) + { t_msk = xt::where((q <= val1) & (q >= val2), 1, t_msk); + } if ((opr2 == "<") && (opr1 == ">")) + { t_msk = xt::where((q < val2) | (q > val1), 1, t_msk); + } else if ((opr2 == "<=") && (opr1 == ">")) + { t_msk = xt::where((q <= val2) | (q > val1), 1, t_msk); + } else if ((opr2 == "<") && (opr1 == ">=")) + { t_msk = xt::where((q < val2) | (q >= val1), 1, t_msk); + } else if ((opr2 == "<=") && (opr1 == ">=")) + { t_msk = xt::where((q <= val2) | (q >= val1), 1, t_msk); + } } else { @@ -320,40 +396,54 @@ namespace evalhyd // apply masking condition to given subset if (opr == "<") + { t_msk = xt::where( q < val, 1, t_msk ); + } else if (opr == ">") + { t_msk = xt::where( q > val, 1, t_msk ); + } else if (opr == "<=") + { t_msk = xt::where( q <= val, 1, t_msk ); + } else if (opr == ">=") + { t_msk = xt::where( q >= val, 1, t_msk ); + } else if (opr == "==") + { t_msk = xt::where( xt::equal(q, val), 1, t_msk ); + } else if (opr == "!=") + { t_msk = xt::where( xt::not_equal(q, val), 1, t_msk ); + } } } } - // condition on time index + // condition on time index else if (var == "t") { for (const auto & sequence : cond) { if (sequence.empty()) + { // i.e. t{:} xt::view(t_msk, xt::all()) = 1; + } else { // convert string indices to integer indices diff --git a/include/evalhyd/detail/probabilist/evaluator.hpp b/include/evalhyd/detail/probabilist/evaluator.hpp index 7574b19..ca209d0 100644 --- a/include/evalhyd/detail/probabilist/evaluator.hpp +++ b/include/evalhyd/detail/probabilist/evaluator.hpp @@ -22,31 +22,31 @@ namespace evalhyd { private: using view1d_xtensor2d_double_type = decltype( - xt::view( - std::declval<const D2&>(), - std::declval<std::size_t>(), - xt::all() - ) + xt::view( + std::declval<const D2&>(), + std::declval<std::size_t>(), + xt::all() + ) ); using view2d_xtensor4d_double_type = decltype( - xt::view( - std::declval<const D4&>(), - std::declval<std::size_t>(), - std::declval<std::size_t>(), - xt::all(), - xt::all() - ) + xt::view( + std::declval<const D4&>(), + std::declval<std::size_t>(), + std::declval<std::size_t>(), + xt::all(), + xt::all() + ) ); using view2d_xtensor4d_bool_type = decltype( - xt::view( - std::declval<const B4&>(), - std::declval<std::size_t>(), - std::declval<std::size_t>(), - xt::all(), - xt::all() - ) + xt::view( + std::declval<const B4&>(), + std::declval<std::size_t>(), + std::declval<std::size_t>(), + xt::all(), + xt::all() + ) ); // members for input data @@ -81,7 +81,9 @@ namespace evalhyd // initialise a mask if none provided // (corresponding to no temporal subset) if (msk.size() < 1) + { t_msk = xt::ones<bool>({std::size_t {1}, q_obs.size()}); + } // determine size for recurring dimensions n = q_obs.size(); @@ -96,9 +98,11 @@ namespace evalhyd auto sum_nan = xt::eval(xt::sum(prd_nan, -1)); if (xt::amin(sum_nan) != xt::amax(sum_nan)) + { throw std::runtime_error( "predictions across members feature non-equal lengths" ); + } auto msk_nan = xt::where(obs_nan | xt::row(prd_nan, 0))[0]; diff --git a/include/evalhyd/detail/uncertainty.hpp b/include/evalhyd/detail/uncertainty.hpp index e55d046..9599b44 100644 --- a/include/evalhyd/detail/uncertainty.hpp +++ b/include/evalhyd/detail/uncertainty.hpp @@ -39,7 +39,8 @@ namespace evalhyd std::tm tm = {}; std::istringstream ss(str); ss >> std::get_time(&tm, "%Y-%m-%d %H:%M:%S"); - if (ss.fail()) { + if (ss.fail()) + { throw std::runtime_error("datetime string parsing failed"); } tm.tm_year += 400; // add 400y to avoid dates prior 1970 @@ -59,11 +60,14 @@ namespace evalhyd // check constant time interval auto ti = x_timepoints[1] - x_timepoints[0]; for (std::size_t t = 1; t < x_timepoints.size() - 1; t++) - if (x_timepoints[t + 1] - x_timepoints[t] != ti) { + { + if (x_timepoints[t + 1] - x_timepoints[t] != ti) + { throw std::runtime_error( "time interval not constant across datetimes" ); } + } // identify start and end years for period int start_yr = v_tm.front().tm_year + 1900; @@ -73,7 +77,8 @@ namespace evalhyd std::tm start_hy = v_tm.front(); xt::xtensor<int, 1> year_blocks = xt::zeros<int>({v_tm.size()}); - for (int y = start_yr; y < end_yr; y++) { + for (int y = start_yr; y < end_yr; y++) + { // define window for year blocks start_hy.tm_year = y - 1900; auto start = std::chrono::system_clock::from_time_t( @@ -89,7 +94,8 @@ namespace evalhyd // check that year is complete (without a rigorous leap year check) int n_days = xt::sum(wdw)(); - if ((n_days != 365) && (n_days != 366)) { + if ((n_days != 365) && (n_days != 366)) + { throw std::runtime_error( "year starting in " + std::to_string(y) + " is incomplete" @@ -101,7 +107,8 @@ namespace evalhyd } // check that time series ends on the last day of a year block - if (year_blocks(year_blocks.size() - 1) == 0) { + if (year_blocks(year_blocks.size() - 1) == 0) + { throw std::runtime_error( "final day of final year not equal to first day of " "first year minus one time step" @@ -116,7 +123,8 @@ namespace evalhyd std::vector<xt::xkeep_slice<int>> samples; // compute metrics for each sample - for (int s = 0; s < n_samples; s++) { + for (int s = 0; s < n_samples; s++) + { // select bootstrapped years auto exp = xt::view(experiment, s); @@ -128,7 +136,8 @@ namespace evalhyd ); xt::xtensor<int, 1> idx = xt::concatenate(xt::xtuple(i0, i1), 0); - for (std::size_t p = 2; p < exp.size(); p++) { + for (std::size_t p = 2; p < exp.size(); p++) + { auto i = xt::flatten_indices( xt::argwhere(xt::equal(year_blocks, exp(p))) ); @@ -146,7 +155,8 @@ namespace evalhyd // TODO: wait for xt::quantile to be available // or implement it internally for n-dim expressions // summary 2: series of quantiles across samples - if (summary == 2) { + if (summary == 2) + { // // adjust last axis size (from samples to number of statistics) // auto s = values.shape(); // s.pop_back(); @@ -161,7 +171,8 @@ namespace evalhyd // TODO: wait for xt::stddev to be fixed for rtensor // or implement it internally for n-dim expressions // summary 1: mean and standard deviation across samples - else if (summary == 1) { + else if (summary == 1) + { // // adjust last axis size (from samples to number of statistics) // auto s = values.shape(); // s.pop_back(); @@ -175,7 +186,8 @@ namespace evalhyd return values; } // summary 0: raw (keep all samples) - else { + else + { return values; } } diff --git a/include/evalhyd/detail/utils.hpp b/include/evalhyd/detail/utils.hpp index 369fd4f..1d11b4b 100644 --- a/include/evalhyd/detail/utils.hpp +++ b/include/evalhyd/detail/utils.hpp @@ -41,11 +41,13 @@ namespace evalhyd { // add elements to pre-computation set for (const auto& element : elements[metric]) + { if (found_elements.find(element) == found_elements.end()) { found_elements.insert(element); required_elements.push_back(element); } + } // add metric dependencies to pre-computation set if (dependencies.find(metric) != dependencies.end()) @@ -59,11 +61,13 @@ namespace evalhyd } // add dependency elements to pre-computation set for (const auto& element : elements[dependency]) + { if (found_elements.find(element) == found_elements.end()) { found_elements.insert(element); required_elements.push_back(element); } + } } } } @@ -104,30 +108,42 @@ namespace evalhyd { // check n_samples if (bootstrap.find("n_samples") == bootstrap.end()) + { throw std::runtime_error( "number of samples missing for bootstrap" ); + } // check len_sample if (bootstrap.find("len_sample") == bootstrap.end()) + { throw std::runtime_error( "length of sample missing for bootstrap" ); + } // check summary if (bootstrap.find("summary") == bootstrap.end()) + { throw std::runtime_error( "summary missing for bootstrap" ); + } auto s = bootstrap.find("summary")->second; // TODO: change upper bound when mean+stddev and quantiles implemented if ((s < 0) || (s > 0)) + { throw std::runtime_error( "invalid value for bootstrap summary" ); + } // set seed if (bootstrap.find("seed") == bootstrap.end()) + { xt::random::seed(time(nullptr)); + } else + { xt::random::seed(bootstrap.find("seed")->second); + } } namespace evalp @@ -151,11 +167,15 @@ namespace evalhyd { if (std::find(threshold_metrics.begin(), threshold_metrics.end(), metric) != threshold_metrics.end()) + { if (thresholds.size() < 1) + { throw std::runtime_error( "missing thresholds *q_thr* required to " "compute " + metric ); + } + } } } } diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp index 625ae23..4cd63a2 100644 --- a/include/evalhyd/evald.hpp +++ b/include/evalhyd/evald.hpp @@ -180,27 +180,40 @@ namespace evalhyd // check that data dimensions are compatible // > time if (q_obs_.shape(1) != q_prd_.shape(1)) + { throw std::runtime_error( "observations and predictions feature different " "temporal lengths" ); + } if (t_msk_.size() > 0) + { if (q_obs_.shape(1) != t_msk_.shape(1)) + { throw std::runtime_error( "observations and masks feature different " "temporal lengths" ); + } + } if (!dts.empty()) + { if (q_obs_.shape(1) != dts.size()) + { throw std::runtime_error( "observations and datetimes feature different " "temporal lengths" ); + } + } + // > series if (q_obs_.shape(0) != 1) + { throw std::runtime_error( "observations contain more than one time series" ); + } // retrieve dimensions std::size_t n_tim = q_obs_.shape(1); @@ -209,26 +222,33 @@ namespace evalhyd // initialise a mask if none provided // (corresponding to no temporal subset) - auto gen_msk = [&]() { + auto gen_msk = [&]() + { // if t_msk provided, it takes priority if (t_msk_.size() > 0) + { return t_msk_; + } // else if m_cdt provided, use them to generate t_msk else if (m_cdt.size() > 0) { B2 c_msk = xt::zeros<bool>({n_msk, n_tim}); for (std::size_t m = 0; m < n_msk; m++) + { xt::view(c_msk, m) = masks::generate_mask_from_conditions( m_cdt[0], xt::view(q_obs_, 0), q_prd_ ); + } return c_msk; } // if neither t_msk nor m_cdt provided, generate dummy mask else + { return B2({xt::ones<bool>({std::size_t{1}, n_tim})}); + } }; auto msk = gen_msk(); @@ -247,16 +267,20 @@ namespace evalhyd else if ( transform == "inv" ) { if ( epsilon == -9 ) + { // determine an epsilon value to avoid zero divide epsilon = xt::mean(q_obs_)() * 0.01; + } return D2(1. / (q + epsilon)); } else if ( transform == "log" ) { if ( epsilon == -9 ) + { // determine an epsilon value to avoid log zero epsilon = xt::mean(q_obs_)() * 0.01; + } return D2(xt::log(q + epsilon)); } @@ -265,8 +289,10 @@ namespace evalhyd if ( exponent < 0 ) { if ( epsilon == -9 ) + { // determine an epsilon value to avoid zero divide epsilon = xt::mean(q_obs_)() * 0.01; + } return D2(xt::pow(q + epsilon, exponent)); } @@ -293,9 +319,11 @@ namespace evalhyd if ((n_samples != -9) && (len_sample != -9)) { if (dts.empty()) + { throw std::runtime_error( "bootstrap requested but datetimes not provided" ); + } exp = uncertainty::bootstrap( dts, n_samples, len_sample @@ -337,21 +365,37 @@ namespace evalhyd for ( const auto& element : req_elt ) { if ( element == "mean_obs" ) + { evaluator.calc_mean_obs(); + } else if ( element == "mean_prd" ) + { evaluator.calc_mean_prd(); + } else if ( element == "quad_err" ) + { evaluator.calc_quad_err(); + } else if ( element == "quad_obs" ) + { evaluator.calc_quad_obs(); + } else if ( element == "quad_prd" ) + { evaluator.calc_quad_prd(); + } else if ( element == "r_pearson" ) + { evaluator.calc_r_pearson(); + } else if ( element == "alpha" ) + { evaluator.calc_alpha(); + } else if ( element == "bias" ) + { evaluator.calc_bias(); + } } // pre-compute required dep @@ -371,28 +415,36 @@ namespace evalhyd { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) + { evaluator.calc_RMSE(); + } r.emplace_back(uncertainty::summarise(evaluator.RMSE, summary)); } else if ( metric == "NSE" ) { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) + { evaluator.calc_NSE(); + } r.emplace_back(uncertainty::summarise(evaluator.NSE, summary)); } else if ( metric == "KGE" ) { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) + { evaluator.calc_KGE(); + } r.emplace_back(uncertainty::summarise(evaluator.KGE, summary)); } else if ( metric == "KGEPRIME" ) { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) + { evaluator.calc_KGEPRIME(); + } r.emplace_back(uncertainty::summarise(evaluator.KGEPRIME, summary)); } } diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp index 7a044c3..c41277a 100644 --- a/include/evalhyd/evalp.hpp +++ b/include/evalhyd/evalp.hpp @@ -159,47 +159,75 @@ namespace evalhyd // check that data dimensions are compatible // > time if (q_obs_.shape(1) != q_prd_.shape(3)) + { throw std::runtime_error( "observations and predictions feature different " "temporal lengths" ); + } if (t_msk_.size() > 0) + { if (q_obs_.shape(1) != t_msk_.shape(3)) + { throw std::runtime_error( "observations and masks feature different " "temporal lengths" ); + } + } if (!dts.empty()) + { if (q_obs_.shape(1) != dts.size()) + { throw std::runtime_error( "observations and datetimes feature different " "temporal lengths" ); + } + } + // > leadtimes if (t_msk_.size() > 0) + { if (q_prd_.shape(1) != t_msk_.shape(1)) + { throw std::runtime_error( "predictions and temporal masks feature different " "numbers of lead times" ); + } + } + // > sites if (q_obs_.shape(0) != q_prd_.shape(0)) + { throw std::runtime_error( "observations and predictions feature different " "numbers of sites" ); + } if (t_msk_.size() > 0) + { if (q_obs_.shape(0) != t_msk_.shape(0)) + { throw std::runtime_error( "observations and temporal masks feature different " "numbers of sites" ); + } + } + if (m_cdt.size() > 0) + { if (q_obs_.shape(0) != m_cdt.shape(0)) + { throw std::runtime_error( "observations and masking conditions feature different " "numbers of sites" ); + } + } + // retrieve dimensions std::size_t n_sit = q_prd_.shape(0); @@ -246,18 +274,28 @@ namespace evalhyd utils::find_requirements(metrics, elt, dep, req_elt, req_dep); // generate masks from conditions if provided - auto gen_msk = [&]() { + auto gen_msk = [&]() + { B4 c_msk = xt::zeros<bool>({n_sit, n_ltm, n_msk, n_tim}); if (m_cdt.size() > 0) + { for (std::size_t s = 0; s < n_sit; s++) + { for (std::size_t l = 0; l < n_ltm; l++) + { for (std::size_t m = 0; m < n_msk; m++) + { xt::view(c_msk, s, l, m) = masks::generate_mask_from_conditions( xt::view(m_cdt, s, m), xt::view(q_obs_, s), xt::view(q_prd_, s, l) ); + } + } + } + } + return c_msk; }; const B4 c_msk = gen_msk(); @@ -269,9 +307,11 @@ namespace evalhyd if ((n_samples != -9) && (len_sample != -9)) { if (dts.empty()) + { throw std::runtime_error( "bootstrap requested but datetimes not provided" ); + } b_exp = uncertainty::bootstrap(dts, n_samples, len_sample); } @@ -286,12 +326,15 @@ namespace evalhyd // initialise data structure for outputs std::vector<xt::xarray<double>> r; for (const auto& metric : metrics) + { r.emplace_back(xt::zeros<double>(dim[metric])); + } auto summary = bootstrap.find("summary")->second; // compute variables one site at a time to minimise memory imprint for (std::size_t s = 0; s < n_sit; s++) + { // compute variables one lead time at a time to minimise memory imprint for (std::size_t l = 0; l < n_ltm; l++) { @@ -314,24 +357,38 @@ namespace evalhyd for (const auto& element : req_elt) { if ( element == "o_k" ) + { evaluator.calc_o_k(); + } else if ( element == "bar_o" ) + { evaluator.calc_bar_o(); + } else if ( element == "y_k" ) + { evaluator.calc_y_k(); + } else if ( element == "q_qnt" ) + { evaluator.calc_q_qnt(); + } } // pre-compute required dep for (const auto& dependency : req_dep) { if ( dependency == "bs" ) + { evaluator.calc_bs(); + } else if ( dependency == "qs" ) + { evaluator.calc_qs(); + } else if ( dependency == "crps" ) + { evaluator.calc_crps(); + } } // retrieve or compute requested metrics @@ -343,7 +400,9 @@ namespace evalhyd { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) + { evaluator.calc_BS(); + } // (sites, lead times, subsets, samples, thresholds) xt::view(r[m], s, l, xt::all(), xt::all(), xt::all()) = uncertainty::summarise(evaluator.BS, summary); @@ -352,7 +411,9 @@ namespace evalhyd { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) + { evaluator.calc_BSS(); + } // (sites, lead times, subsets, samples, thresholds) xt::view(r[m], s, l, xt::all(), xt::all(), xt::all()) = uncertainty::summarise(evaluator.BSS, summary); @@ -361,7 +422,9 @@ namespace evalhyd { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) + { evaluator.calc_BS_CRD(); + } // (sites, lead times, subsets, samples, thresholds, components) xt::view(r[m], s, l, xt::all(), xt::all(), xt::all(), xt::all()) = uncertainty::summarise(evaluator.BS_CRD, summary); @@ -370,7 +433,10 @@ namespace evalhyd { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) + { evaluator.calc_BS_LBD(); + } + // (sites, lead times, subsets, samples, thresholds, components) xt::view(r[m], s, l, xt::all(), xt::all(), xt::all(), xt::all()) = uncertainty::summarise(evaluator.BS_LBD, summary); @@ -379,7 +445,9 @@ namespace evalhyd { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) + { evaluator.calc_QS(); + } // (sites, lead times, subsets, samples, quantiles) xt::view(r[m], s, l, xt::all(), xt::all(), xt::all()) = uncertainty::summarise(evaluator.QS, summary); @@ -388,14 +456,16 @@ namespace evalhyd { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) + { evaluator.calc_CRPS(); + } // (sites, lead times, subsets, samples) xt::view(r[m], s, l, xt::all(), xt::all()) = uncertainty::summarise(evaluator.CRPS, summary); } } } - + } return r; }; } diff --git a/tests/test_determinist.cpp b/tests/test_determinist.cpp index 5239076..d5f725e 100644 --- a/tests/test_determinist.cpp +++ b/tests/test_determinist.cpp @@ -301,8 +301,10 @@ TEST(DeterministTests, TestMissingData) std::vector<xt::xarray<double>> metrics_nan = evalhyd::evald<xt::xtensor<double, 2>, xt::xtensor<bool, 2>>(observed, predicted, metrics); - for (std::size_t m = 0; m < metrics.size(); m++) { - for (std::size_t p = 0; p < predicted.shape(0); p++) { + for (std::size_t m = 0; m < metrics.size(); m++) + { + for (std::size_t p = 0; p < predicted.shape(0); p++) + { // compute metrics on subset of observations and predictions (i.e. // eliminating region with NaN in observations or predictions manually) xt::xtensor<double, 1> obs = diff --git a/tests/test_probabilist.cpp b/tests/test_probabilist.cpp index c258b0c..342cbf0 100644 --- a/tests/test_probabilist.cpp +++ b/tests/test_probabilist.cpp @@ -377,7 +377,8 @@ TEST(ProbabilistTests, TestMissingData) ); // check that numerical results are identical - for (std::size_t m = 0; m < metrics.size(); m++) { + for (std::size_t m = 0; m < metrics.size(); m++) + { // for leadtime 1 EXPECT_TRUE( xt::allclose( -- GitLab