diff --git a/CMakeLists.txt b/CMakeLists.txt index bff6e9365d385b61e9db6d2ba080fac0f5b8f0c3..4e0cfc77977f1211ee849b10e1359d84e393a5cf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,7 +7,7 @@ cmake_minimum_required(VERSION 3.15) project( EvalHyd LANGUAGES CXX - VERSION 0.1.0 + VERSION 0.1.1 DESCRIPTION "Utility to evaluate streamflow predictions" ) diff --git a/changelog.rst b/changelog.rst index bbc1cb24d2e30e758d3858d056455d4972dd8336..d47ce60b53d322c48ef57e5f4bb0af60165a7a96 100644 --- a/changelog.rst +++ b/changelog.rst @@ -6,6 +6,31 @@ Yet to be versioned and released. Only available from *dev* branch until then. + +v0.1.1 +------ + +Released on 2023-06-16. + +.. rubric:: Scope changes + +* remove `"WSS"` and `"AWI"` as probabilistic evaluation metric because of the + arbitrary nature of using the sample climatology as reference/benchmark + (`CPP#7 <https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-cpp/-/issues/7>`_) + +.. rubric:: Bug fixes + +* fix irrelevant rank check failure when passing arrays (and not tensors) + (`CPP#6 <https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-cpp/-/issues/6>`_, + `R#4 <https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-r/-/issues/4>`_) +* fix crashing conditional masking functionality at runtime with certain + compilers on Windows (in particular when building with Rtools for R>4.1.3) due + to the presence of square brackets in the regular expressions that do not seem + supported + (`R#6 <https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-r/-/issues/6>`_) +* fix bug when using calendar year as block for the bootstrapping functionality + (`CPP#8 <https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-cpp/-/issues/8>`_) + v0.1.0 ------ diff --git a/include/evalhyd/detail/determinist/efficiencies.hpp b/include/evalhyd/detail/determinist/efficiencies.hpp index df893670f3f70e4b900903e9417d368f5a3c0ce4..b5d230249aae761f65d8b102a877631a66cad75a 100644 --- a/include/evalhyd/detail/determinist/efficiencies.hpp +++ b/include/evalhyd/detail/determinist/efficiencies.hpp @@ -809,4 +809,4 @@ namespace evalhyd } } -#endif //EVALHYD_DETERMINIST_EFFICIENCIES_HPP \ No newline at end of file +#endif //EVALHYD_DETERMINIST_EFFICIENCIES_HPP diff --git a/include/evalhyd/detail/determinist/errors.hpp b/include/evalhyd/detail/determinist/errors.hpp index 646bbd4ad0776c7c97c4cb4f5a0bbbb7922125a3..4c778e20f362f2a4ee33cb7dbe629ff946c79032 100644 --- a/include/evalhyd/detail/determinist/errors.hpp +++ b/include/evalhyd/detail/determinist/errors.hpp @@ -478,4 +478,4 @@ namespace evalhyd } } -#endif //EVALHYD_DETERMINIST_ERRORS_HPP \ No newline at end of file +#endif //EVALHYD_DETERMINIST_ERRORS_HPP diff --git a/include/evalhyd/detail/masks.hpp b/include/evalhyd/detail/masks.hpp index d1d503e257a826f6448f199c42265daa874fcd65..fe40d078f6cd4836b63fd72499a96c56ffbb9382 100644 --- a/include/evalhyd/detail/masks.hpp +++ b/include/evalhyd/detail/masks.hpp @@ -35,7 +35,10 @@ namespace evalhyd // observed or predicted (median or mean for probabilist) streamflow // e.g. q{>9.} q{<9} q{>=99.0} q{<=99} q{>9,<99} q{==9} q{!=9} std::regex exp_q ( - R"((q_obs|q_prd_median|q_prd_mean)\{(((<|>|<=|>=|==|!=)(mean,?|median,?|qtl[0-1]\.[0-9]+,?|[0-9]+\.?[0-9]*,?))+)\})" + R"((q_obs|q_prd_median|q_prd_mean)\{(((<|>|<=|>=|==|!=)(mean,?|median,?|qtl(0|1)\.(0|1|2|3|4|5|6|7|8|9)+,?|(0|1|2|3|4|5|6|7|8|9)+\.?(0|1|2|3|4|5|6|7|8|9)*,?))+)\})" + // NOTE: this should be `R"((q_obs|q_prd_median|q_prd_mean)\{(((<|>|<=|>=|==|!=)(mean,?|median,?|qtl[0-1]\.[0-9]+,?|[0-9]+\.?[0-9]*,?))+)\})"` + // but there is a bug in the building chain for R packages + // https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-r/-/issues/6 ); for (std::sregex_iterator i = @@ -51,7 +54,12 @@ namespace evalhyd std::vector<std::vector<std::string>> conditions; // pattern supported to specify masking conditions based on streamflow - std::regex ex (R"((<|>|<=|>=|==|!=)(mean|median|qtl[0-1]\.[0-9]+|[0-9]+\.?[0-9]*))"); + std::regex ex ( + R"((<|>|<=|>=|==|!=)(mean|median|qtl(0|1)\.(0|1|2|3|4|5|6|7|8|9)+|(0|1|2|3|4|5|6|7|8|9)+\.?(0|1|2|3|4|5|6|7|8|9)*))" + // NOTE: this should be `R"((<|>|<=|>=|==|!=)(mean|median|qtl[0-1]\.[0-9]+|[0-9]+\.?[0-9]*))"` + // but there is a bug in the building chain for R packages + // https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-r/-/issues/6 + ); for (std::sregex_iterator j = std::sregex_iterator(str.begin(), str.end(), ex); @@ -92,7 +100,12 @@ namespace evalhyd // pattern supported to specify conditions to generate masks on time index // e.g. t{0:10} t{0:10,20:30} t{0,1,2,3} t{0:10,30,40,50} t{:} - std::regex exp_t (R"((t)\{(:|([0-9]+:[0-9]+,?|[0-9]+,?)+)\})"); + std::regex exp_t ( + R"((t)\{(:|((0|1|2|3|4|5|6|7|8|9)+:(0|1|2|3|4|5|6|7|8|9)+,?|(0|1|2|3|4|5|6|7|8|9)+,?)+)\})" + // NOTE: this should be `R"((t)\{(:|([0-9]+:[0-9]+,?|[0-9]+,?)+)\})"` + // but there is a bug in the building chain for R packages + // https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-r/-/issues/6 + ); for (std::sregex_iterator i = std::sregex_iterator(msk_str.begin(), msk_str.end(), exp_t); @@ -114,7 +127,12 @@ namespace evalhyd else { // pattern supported to specify masking conditions based on time index - std::regex e (R"([0-9]+:[0-9]+|[0-9]+)"); + std::regex e ( + R"((0|1|2|3|4|5|6|7|8|9)+:(0|1|2|3|4|5|6|7|8|9)+|(0|1|2|3|4|5|6|7|8|9)+)" + // NOTE: this should be `R"([0-9]+:[0-9]+|[0-9]+)"` + // but there is a bug in the building chain for R packages + // https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-r/-/issues/6 + ); for (std::sregex_iterator j = std::sregex_iterator(s.begin(), s.end(), e); diff --git a/include/evalhyd/detail/probabilist/cdf.hpp b/include/evalhyd/detail/probabilist/cdf.hpp index e0cf0b146c8a9d65554bc43e416075c1b4d83c1a..1b02c2cbbb993172f6d78666537aa036ae22111f 100644 --- a/include/evalhyd/detail/probabilist/cdf.hpp +++ b/include/evalhyd/detail/probabilist/cdf.hpp @@ -198,4 +198,4 @@ namespace evalhyd } } -#endif //EVALHYD_PROBABILIST_CDF_HPP \ No newline at end of file +#endif //EVALHYD_PROBABILIST_CDF_HPP diff --git a/include/evalhyd/detail/probabilist/contingency.hpp b/include/evalhyd/detail/probabilist/contingency.hpp index 0c3aa0ecd0fd669558eb7b69d0e602aca2ac6fad..a4685e4aec83615ebad1c612d3fbcdb72b42f97b 100644 --- a/include/evalhyd/detail/probabilist/contingency.hpp +++ b/include/evalhyd/detail/probabilist/contingency.hpp @@ -541,4 +541,4 @@ namespace evalhyd } } -#endif //EVALHYD_PROBABILIST_CONTINGENCY_HPP \ No newline at end of file +#endif //EVALHYD_PROBABILIST_CONTINGENCY_HPP diff --git a/include/evalhyd/detail/probabilist/evaluator.hpp b/include/evalhyd/detail/probabilist/evaluator.hpp index 6293ea29e546df8f9f3756608d213aa4e5773a9e..976145c58817f529c2786a246665c7eac55c5eb9 100644 --- a/include/evalhyd/detail/probabilist/evaluator.hpp +++ b/include/evalhyd/detail/probabilist/evaluator.hpp @@ -76,7 +76,6 @@ namespace evalhyd xtl::xoptional<xt::xtensor<double, 5>, bool> itv_bnds; xtl::xoptional<xt::xtensor<double, 4>, bool> obs_in_itv; xtl::xoptional<xt::xtensor<double, 4>, bool> itv_width; - xtl::xoptional<xt::xtensor<double, 6>, bool> clim_bnds; // members for intermediate evaluation metrics // (i.e. before the reduction along the temporal axis) @@ -124,9 +123,7 @@ namespace evalhyd xtl::xoptional<xt::xtensor<double, 5>, bool> CR; xtl::xoptional<xt::xtensor<double, 5>, bool> AW; xtl::xoptional<xt::xtensor<double, 5>, bool> AWN; - xtl::xoptional<xt::xtensor<double, 5>, bool> AWI; xtl::xoptional<xt::xtensor<double, 5>, bool> WS; - xtl::xoptional<xt::xtensor<double, 5>, bool> WSS; // > Multi-variate xtl::xoptional<xt::xtensor<double, 4>, bool> ES; @@ -367,19 +364,6 @@ namespace evalhyd return itv_width.value(); }; - - xt::xtensor<double, 6> get_clim_bnds() - { - if (!clim_bnds.has_value()) - { - clim_bnds = elements::calc_clim_bnds( - q_obs, get_c_lvl(), t_msk, b_exp, - n_sit, n_ldt, n_itv, n_msk, n_exp - ); - } - return clim_bnds.value(); - }; - // methods to compute intermediate metrics xt::xtensor<double, 4> get_bs() { @@ -798,17 +782,6 @@ namespace evalhyd return AWN.value(); }; - xt::xtensor<double, 5> get_AWI() - { - if (!AWI.has_value()) - { - AWI = metrics::calc_AWI( - get_AW(), get_clim_bnds() - ); - } - return AWI.value(); - }; - xt::xtensor<double, 5> get_WS() { if (!WS.has_value()) @@ -821,18 +794,6 @@ namespace evalhyd return WS.value(); }; - xt::xtensor<double, 5> get_WSS() - { - if (!WSS.has_value()) - { - WSS = metrics::calc_WSS( - q_obs, get_c_lvl(), get_clim_bnds(), get_WS(), - t_msk, b_exp, n_sit, n_ldt, n_itv, n_msk, n_exp - ); - } - return WSS.value(); - }; - xt::xtensor<double, 4> get_ES() { if (!ES.has_value()) diff --git a/include/evalhyd/detail/probabilist/intervals.hpp b/include/evalhyd/detail/probabilist/intervals.hpp index c36a235838ead80eec395b4ccd342684217b5afc..adbd487b592ab4903210072e9ae01836ba41e0ac 100644 --- a/include/evalhyd/detail/probabilist/intervals.hpp +++ b/include/evalhyd/detail/probabilist/intervals.hpp @@ -123,119 +123,6 @@ namespace evalhyd return (u - l); } - - /// Compute the bounds of the climatology corresponding to the - /// confidence levels. - /// - /// \param q_obs - /// Streamflow predictions. - /// shape: (sites, time) - /// \param c_lvl - /// Confidence levels for the predictive intervals. - /// shape: (intervals,) - /// \param t_msk - /// Temporal subsets of the whole record. - /// shape: (sites, lead times, subsets, time) - /// \param b_exp - /// Boostrap samples. - /// shape: (samples, time slice) - /// \param n_sit - /// Number of sites. - /// \param n_ldt - /// Number of lead times. - /// \param n_itv - /// Number of predictive intervals. - /// \param n_msk - /// Number of temporal subsets. - /// \param n_exp - /// Number of bootstrap samples. - /// \return - /// Climatology bounds. - /// shape: (sites, lead times, subsets, samples, intervals, bounds) - template <class XD2> - inline xt::xtensor<double, 6> calc_clim_bnds( - const XD2& q_obs, - const xt::xtensor<double, 1>& c_lvl, - const xt::xtensor<bool, 4>& t_msk, - const std::vector<xt::xkeep_slice<int>>& b_exp, - std::size_t n_sit, - std::size_t n_ldt, - std::size_t n_itv, - std::size_t n_msk, - std::size_t n_exp - ) - { - // compute "climatology" average width - xt::xtensor<double, 6> clim_bnds = - xt::zeros<double>({n_sit, n_ldt, n_msk, n_exp, - n_itv, std::size_t {2}}); - - // determine quantiles forming the predictive intervals - // from the confidence levels - xt::xtensor<double, 2> quantiles = - xt::zeros<double>({n_itv, std::size_t {2}}); - xt::col(quantiles, 0) = 0.5 - c_lvl / 200.; - xt::col(quantiles, 1) = 0.5 + c_lvl / 200.; - - // compute variable one mask at a time to minimise memory imprint - for (std::size_t m = 0; m < n_msk; m++) - { - // apply the mask - // (using NaN workaround until reducers work on masked_view) - auto q_obs_masked = xt::where( - xt::view(t_msk, xt::all(), xt::all(), m, xt::all()), - xt::view(q_obs, xt::all(), xt::newaxis(), xt::all()), - NAN - ); - - // compute variable one sample at a time - for (std::size_t e = 0; e < n_exp; e++) - { - // apply the bootstrap sampling - auto q_obs_masked_sampled = - xt::view(q_obs_masked, xt::all(), xt::all(), b_exp[e]); - - // compute "climatology" interval - for (std::size_t s = 0; s < n_sit; s++) - { - for (std::size_t l = 0; l < n_ldt; l++) - { - for (std::size_t i = 0; i < n_itv; i++) - { - auto obs = xt::view(q_obs_masked_sampled, - s, l, xt::all()); - auto obs_filtered = - xt::filter(obs, !xt::isnan(obs)); - - if (obs_filtered.size() > 0) - { - // lower bound - xt::view(clim_bnds, s, l, m, e, i, 0) = - xt::quantile( - obs_filtered, - {quantiles(i, 0)} - )(); - // upper bound - xt::view(clim_bnds, s, l, m, e, i, 1) = - xt::quantile( - obs_filtered, - {quantiles(i, 1)} - )(); - } - else - { - xt::view(clim_bnds, s, l, m, e, i, xt::all()) = - NAN; - } - - } - } - } - } - } - - return clim_bnds; - } } namespace intermediate @@ -487,34 +374,6 @@ namespace evalhyd - std::numeric_limits<double>::infinity()); } - /// Compute the Average Width Index (AWI). - /// - /// \param AW - /// Average widths. - /// shape: (sites, lead times, subsets, samples, intervals) - /// \param clim_bnds - /// Climatology bounds. - /// shape: (sites, lead times, subsets, samples, intervals, bounds) - /// \return - /// Average width indices. - /// shape: (sites, lead times, subsets, samples, intervals) - inline xt::xtensor<double, 5> calc_AWI( - const xt::xtensor<double, 5>& AW, - const xt::xtensor<double, 6>& clim_bnds - ) - { - // compute "climatology" average width - auto AW_clim = - xt::view(clim_bnds, xt::all(), xt::all(), xt::all(), - xt::all(), xt::all(), 1) - - xt::view(clim_bnds, xt::all(), xt::all(), xt::all(), - xt::all(), xt::all(), 0); - - return xt::where(AW_clim > 0, - 1 - (AW / AW_clim), - - std::numeric_limits<double>::infinity()); - } - /// Compute the Winkler scores (WS), also known as interval score. /// /// \param ws @@ -554,95 +413,6 @@ namespace evalhyd ws, t_msk, b_exp, n_sit, n_ldt, n_itv, n_msk, n_exp ); } - - /// Compute the Winkler skill scores (WSS). - /// - /// \param q_obs - /// Streamflow predictions. - /// shape: (sites, time) - /// \param c_lvl - /// Confidence levels for the predictive intervals. - /// shape: (intervals,) - /// \param clim_bnds - /// Climatology bounds. - /// shape: (sites, lead times, subsets, samples, intervals, bounds) - /// \param WS - /// Winkler scores. - /// shape: (sites, lead times, subsets, samples, intervals) - /// \param t_msk - /// Temporal subsets of the whole record. - /// shape: (sites, lead times, subsets, time) - /// \param b_exp - /// Boostrap samples. - /// shape: (samples, time slice) - /// \param n_sit - /// Number of sites. - /// \param n_ldt - /// Number of lead times. - /// \param n_itv - /// Number of predictive intervals. - /// \param n_msk - /// Number of temporal subsets. - /// \param n_exp - /// Number of bootstrap samples. - /// \return - /// Winkler skill scores. - /// shape: (sites, lead times, subsets, samples, intervals) - template <class XD2> - inline xt::xtensor<double, 5> calc_WSS( - const XD2& q_obs, - const xt::xtensor<double, 1>& c_lvl, - const xt::xtensor<double, 6>& clim_bnds, - const xt::xtensor<double, 5>& WS, - const xt::xtensor<bool, 4>& t_msk, - const std::vector<xt::xkeep_slice<int>>& b_exp, - std::size_t n_sit, - std::size_t n_ldt, - std::size_t n_itv, - std::size_t n_msk, - std::size_t n_exp - ) - { - // compute "climatology" Winkler score - xt::xtensor<double, 5> WS_clim = - xt::zeros<double>({n_sit, n_ldt, n_msk, n_exp, n_itv}); - - for (std::size_t l = 0; l < n_ldt; l++) - { - for (std::size_t m = 0; m < n_msk; m++) - { - for (std::size_t e = 0; e < n_exp; e++) - { - auto ws_clim = intermediate::calc_ws( - q_obs, c_lvl, - xt::view(clim_bnds, xt::all(), l, xt::newaxis(), - m, e, xt::all(), xt::all(), - xt::newaxis()) - ); - - // apply the mask - // (using NaN workaround until reducers work on masked_view) - auto ws_clim_masked = xt::where( - xt::view(t_msk, xt::all(), l, m, xt::newaxis(), xt::all()), - xt::view(ws_clim, xt::all(), l, xt::all(), xt::all()), - NAN - ); - - // apply the bootstrap sampling - auto ws_clim_masked_sampled = - xt::view(ws_clim_masked, xt::all(), xt::all(), b_exp[e]); - - xt::view(WS_clim, xt::all(), l, m, e, xt::all()) = - xt::nanmean(ws_clim_masked_sampled, -1); - } - } - } - - // compute the Winkler skill score - return xt::where(WS_clim > 0, - 1 - (WS / WS_clim), - - std::numeric_limits<double>::infinity()); - } } } } diff --git a/include/evalhyd/detail/probabilist/ranks.hpp b/include/evalhyd/detail/probabilist/ranks.hpp index ee3029006470993397eb619cfbda6aebeff8eb0f..ea68e6d0aa623e441ce1262b1e27c579fd244c19 100644 --- a/include/evalhyd/detail/probabilist/ranks.hpp +++ b/include/evalhyd/detail/probabilist/ranks.hpp @@ -449,4 +449,4 @@ namespace evalhyd } } -#endif //EVALHYD_PROBABILIST_RANKS_HPP \ No newline at end of file +#endif //EVALHYD_PROBABILIST_RANKS_HPP diff --git a/include/evalhyd/detail/uncertainty.hpp b/include/evalhyd/detail/uncertainty.hpp index 25b13af0957e7552f52fe648564f2abc67be0054..af55fc2425c19b68bdaba70488d9305b60478825 100644 --- a/include/evalhyd/detail/uncertainty.hpp +++ b/include/evalhyd/detail/uncertainty.hpp @@ -76,7 +76,16 @@ namespace evalhyd int start_yr = v_tm.front().tm_year + 1900; int end_yr = v_tm.back().tm_year + 1900; - // assume start of year block as start of time series + // deal with special case with a start on 1st of January + // (note: use string rather than *tm_yday* member of time_point + // because *tm_yday* is not set when using `std::get_time`) + if (datetimes[0].substr(5, 5) == "01-01") + { + // add one year to make sure last year is included in loop + end_yr += 1; + } + + // take start of year block as start of time series std::tm start_hy = v_tm.front(); xt::xtensor<int, 1> year_blocks = xt::zeros<int>({v_tm.size()}); @@ -100,7 +109,7 @@ namespace evalhyd if ((n_days != 365) && (n_days != 366)) { throw std::runtime_error( - "year starting in " + std::to_string(y) + "year starting in " + std::to_string(y - 400) + " is incomplete" ); } diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp index 19cbdbcecde37ea6f5e82ce593c273916d04bae4..9ddeeadf1ddb42ba5caed1bc54748ce24b351be2 100644 --- a/include/evalhyd/evald.hpp +++ b/include/evalhyd/evald.hpp @@ -210,19 +210,25 @@ namespace evalhyd xtl::missing<const std::vector<std::string>>() ) { - // check ranks of tensors - if (xt::get_rank<XD2>::value != 2) + // check ranks of expressions if they are tensors + if (xt::get_rank<XD2>::value != SIZE_MAX) { - throw std::runtime_error( - "observations and/or predictions and/or thresholds " - "are not two-dimensional" - ); + if (xt::get_rank<XD2>::value != 2) + { + throw std::runtime_error( + "observations and/or predictions and/or thresholds " + "are not two-dimensional" + ); + } } - if (xt::get_rank<XB3>::value != 3) + if (xt::get_rank<XB3>::value != SIZE_MAX) { - throw std::runtime_error( - "temporal masks are not three-dimensional" - ); + if (xt::get_rank<XB3>::value != 3) + { + throw std::runtime_error( + "temporal masks are not three-dimensional" + ); + } } // retrieve real types of the expressions @@ -366,7 +372,7 @@ namespace evalhyd } else { - return XB3({}); + return XB3(xt::zeros<bool>({0, 0, 0})); } }; const XB3 c_msk = gen_msk(); diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp index d7e1d4279985fec564b21f41b21043edf22768ad..c04412097f925f675cb5750037086c421dceda9e 100644 --- a/include/evalhyd/evalp.hpp +++ b/include/evalhyd/evalp.hpp @@ -166,7 +166,7 @@ namespace evalhyd /// /// .. code-block:: c++ /// - /// evalhyd::evalp(obs, prd, {"CRPS"}); + /// evalhyd::evalp(obs, prd, {"CRPS_FROM_QS"}); /// /// \endrst template <class XD2, class XD4, class XB4 = xt::xtensor<bool, 4>, @@ -190,24 +190,35 @@ namespace evalhyd xtl::missing<const std::vector<std::string>>() ) { - // check ranks of tensors - if (xt::get_rank<XD2>::value != 2) + // check ranks of expressions if they are tensors + if (xt::get_rank<XD2>::value != SIZE_MAX) { - throw std::runtime_error( - "observations and/or thresholds are not two-dimensional" - ); + if (xt::get_rank<XD2>::value != 2) + { + throw std::runtime_error( + "observations and/or thresholds are not two-dimensional" + ); + } } - if (xt::get_rank<XD4>::value != 4) + + if (xt::get_rank<XD4>::value != SIZE_MAX) { - throw std::runtime_error( - "predictions are not four-dimensional" - ); + if (xt::get_rank<XD4>::value != 4) + { + throw std::runtime_error( + "predictions are not four-dimensional" + ); + } } - if (xt::get_rank<XB4>::value != 4) + + if (xt::get_rank<XB4>::value != SIZE_MAX) { - throw std::runtime_error( - "temporal masks are not four-dimensional" - ); + if (xt::get_rank<XB4>::value != 4) + { + throw std::runtime_error( + "temporal masks are not four-dimensional" + ); + } } // retrieve real types of the expressions @@ -229,7 +240,7 @@ namespace evalhyd "QS", "CRPS_FROM_QS", "POD", "POFD", "FAR", "CSI", "ROCSS", "RANK_HIST", "DS", "AS", - "CR", "AW", "AWN", "AWI", "WS", "WSS", + "CR", "AW", "AWN", "WS", "ES"} ); @@ -368,7 +379,7 @@ namespace evalhyd } else { - return XB4({}); + return XB4(xt::zeros<bool>({0, 0, 0, 0})); } }; const XB4 c_msk = gen_msk(); @@ -536,24 +547,12 @@ namespace evalhyd uncertainty::summarise_p(evaluator.get_AWN(), summary) ); } - else if ( metric == "AWI" ) - { - r.emplace_back( - uncertainty::summarise_p(evaluator.get_AWI(), summary) - ); - } else if ( metric == "WS" ) { r.emplace_back( uncertainty::summarise_p(evaluator.get_WS(), summary) ); } - else if ( metric == "WSS" ) - { - r.emplace_back( - uncertainty::summarise_p(evaluator.get_WSS(), summary) - ); - } else if ( metric == "ES" ) { r.emplace_back( diff --git a/tests/expected/evalp/AWI.csv b/tests/expected/evalp/AWI.csv deleted file mode 100644 index 8b3b7e7f9598b23911dc06f440d71572b99b48e5..0000000000000000000000000000000000000000 --- a/tests/expected/evalp/AWI.csv +++ /dev/null @@ -1 +0,0 @@ -0.9821120161733,0.9880951944476 diff --git a/tests/expected/evalp/WSS.csv b/tests/expected/evalp/WSS.csv deleted file mode 100644 index fd929dabd8f58354be587455dcbd4402d6a9f83d..0000000000000000000000000000000000000000 --- a/tests/expected/evalp/WSS.csv +++ /dev/null @@ -1 +0,0 @@ -0.6621887740287,0.4360388849930 diff --git a/tests/test_probabilist.cpp b/tests/test_probabilist.cpp index 65b16cdde5c4e0b30280bed9a96a8cfea6f36ba9..f4e53f71cc694abffe77c5d5287af2fa6124c4c9 100644 --- a/tests/test_probabilist.cpp +++ b/tests/test_probabilist.cpp @@ -35,7 +35,7 @@ std::vector<std::string> all_metrics_p = { "QS", "CRPS_FROM_QS", "POD", "POFD", "FAR", "CSI", "ROCSS", "RANK_HIST", "DS", "AS", - "CR", "AW", "AWN", "AWI", "WS", "WSS", + "CR", "AW", "AWN", "WS", "ES" }; @@ -263,7 +263,7 @@ TEST(ProbabilistTests, TestIntervals) auto expected = load_expected_p(); // compute scores - std::vector<std::string> metrics = {"CR", "AW", "AWN", "AWI", "WS", "WSS"}; + std::vector<std::string> metrics = {"CR", "AW", "AWN", "WS"}; std::vector<xt::xarray<double>> results = evalhyd::evalp( @@ -271,7 +271,7 @@ TEST(ProbabilistTests, TestIntervals) xt::eval(xt::view(observed, xt::newaxis(), xt::all())), // shape: (sites [1], lead times [1], members [m], time [t]) xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())), - {"CR", "AW", "AWN", "AWI", "WS", "WSS"}, + {"CR", "AW", "AWN", "WS"}, xt::xtensor<double, 2>({}), "", // events {30., 80.} // c_lvl