diff --git a/include/evalhyd/detail/probabilist/evaluator.hpp b/include/evalhyd/detail/probabilist/evaluator.hpp index 0e2ba0561dddc8e7b42c153146ab7c820824e225..c0c960158ada6581fb1d3f6f769b31f2fb8ffe87 100644 --- a/include/evalhyd/detail/probabilist/evaluator.hpp +++ b/include/evalhyd/detail/probabilist/evaluator.hpp @@ -35,6 +35,9 @@ namespace evalhyd XB4 t_msk; const std::vector<xt::xkeep_slice<int>>& b_exp; + // member for "reproducible randomness" + const long int random_seed; + // members for dimensions std::size_t n_sit; std::size_t n_ldt; @@ -255,7 +258,7 @@ namespace evalhyd if (!r_k.has_value()) { r_k = elements::calc_r_k( - q_obs, get_q_qnt(), n_mbr + q_obs, get_q_qnt(), n_mbr, random_seed ); } return r_k.value(); @@ -358,9 +361,11 @@ namespace evalhyd const XD2& thr, xtl::xoptional<const std::string&, bool> events, const XB4& msk, - const std::vector<xt::xkeep_slice<int>>& exp) : + const std::vector<xt::xkeep_slice<int>>& exp, + const long int seed) : q_obs{obs}, q_prd{prd}, - _q_thr{thr}, _events{events}, t_msk(msk), b_exp(exp) + _q_thr{thr}, _events{events}, t_msk(msk), b_exp(exp), + random_seed{seed} { // initialise a mask if none provided // (corresponding to no temporal subset) diff --git a/include/evalhyd/detail/probabilist/ranks.hpp b/include/evalhyd/detail/probabilist/ranks.hpp index a03f4d00604ee9e8e939aacaa58ed81a7374f536..02082338ee6e983a492e97d293f520de227a0b8f 100644 --- a/include/evalhyd/detail/probabilist/ranks.hpp +++ b/include/evalhyd/detail/probabilist/ranks.hpp @@ -29,6 +29,8 @@ namespace evalhyd /// shape: (sites, lead times, quantiles, time) /// \param n_mbr /// Number of ensemble members. + /// \param seed + /// Seed to be used by random generator. /// \return /// Ranks of streamflow observations. /// shape: (sites, lead times, time) @@ -36,7 +38,8 @@ namespace evalhyd inline xt::xtensor<double, 3> calc_r_k( const XD2& q_obs, const XD4& q_qnt, - std::size_t n_mbr + std::size_t n_mbr, + long int seed ) { xt::xtensor<double, 3> ranks = xt::zeros<double>( @@ -96,6 +99,7 @@ namespace evalhyd ); // for ties, take random rank between min and max + xt::random::seed(seed); xt::view(ranks, xt::all()) = xt::where( !xt::isnan(min_ranks), min_ranks diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp index 4bd0f955d5afafa24742e81d9dbb85bc209705f8..59977c80086e021b1714fff620ad7e84928c56aa 100644 --- a/include/evalhyd/evalp.hpp +++ b/include/evalhyd/evalp.hpp @@ -373,7 +373,8 @@ namespace evalhyd probabilist::Evaluator<XD2, XD4, XB4> evaluator( q_obs_, q_prd_, q_thr_, events, t_msk_.size() > 0 ? t_msk_: (m_cdt.size() > 0 ? c_msk : t_msk_), - b_exp + b_exp, + random_seed ); // initialise data structure for outputs diff --git a/tests/test_probabilist.cpp b/tests/test_probabilist.cpp index 8a54c25261ffc8b119046fd1ef34b5feb3dc0beb..f8a1100c263dc7b7d817d023fc18e7b04e384691 100644 --- a/tests/test_probabilist.cpp +++ b/tests/test_probabilist.cpp @@ -6,7 +6,10 @@ #include <vector> #include <tuple> #include <array> + #include <gtest/gtest.h> + +#include <xtl/xoptional.hpp> #include <xtensor/xtensor.hpp> #include <xtensor/xview.hpp> #include <xtensor/xsort.hpp> @@ -25,7 +28,8 @@ using namespace xt::placeholders; // required for `_` to work std::vector<std::string> all_metrics_p = { "BS", "BSS", "BS_CRD", "BS_LBD", "QS", "CRPS", - "POD", "POFD", "FAR", "CSI", "ROCSS" + "POD", "POFD", "FAR", "CSI", "ROCSS", + "RANK_DIAG", "DS", "AS" }; std::tuple<xt::xtensor<double, 1>, xt::xtensor<double, 2>> load_data_p() @@ -411,6 +415,97 @@ TEST(ProbabilistTests, TestContingency) ); } + +TEST(ProbabilistTests, TestRanks) +{ + // read in data + xt::xtensor<double, 1> observed; + xt::xtensor<double, 2> predicted; + std::tie(observed, predicted) = load_data_p(); + + // compute scores + std::vector<xt::xarray<double>> metrics = + evalhyd::evalp( + // shape: (sites [1], time [t]) + 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())), + {"RANK_DIAG", "DS", "AS"}, + xt::xtensor<double, 2>({}), + "high", // events + xt::xtensor<bool, 4>({}), // t_msk + {}, // m_cdt + xtl::missing<const std::unordered_map<std::string, int>>(), // bootstrap + {}, // dts + 7 // seed + ); + + // check results + // Rank diagram + xt::xtensor<double, 5> rank_diag; +#if EVALHYD_TESTING_OS == WINDOWS + rank_diag = {{{{{ 0.607717, 0. , 0. , 0. , 0. , 0.003215, + 0. , 0. , 0. , 0. , 0. , 0. , + 0.003215, 0. , 0.003215, 0. , 0. , 0. , + 0. , 0. , 0. , 0. , 0. , 0. , + 0. , 0. , 0. , 0. , 0. , 0.003215, + 0. , 0. , 0. , 0. , 0. , 0.006431, + 0. , 0. , 0.003215, 0.006431, 0. , 0. , + 0. , 0.003215, 0. , 0. , 0.003215, 0.003215, + 0.003215, 0. , 0.006431, 0.344051}}}}}; +#elif EVALHYD_TESTING_OS == MACOS + rank_diag = {{{{{ 0.607717, 0. , 0. , 0. , 0. , 0.003215, + 0. , 0. , 0. , 0. , 0. , 0. , + 0.003215, 0. , 0.003215, 0. , 0. , 0. , + 0. , 0. , 0. , 0. , 0. , 0. , + 0. , 0. , 0. , 0. , 0. , 0.003215, + 0. , 0. , 0. , 0. , 0. , 0.006431, + 0. , 0. , 0.003215, 0.006431, 0. , 0. , + 0. , 0.003215, 0. , 0. , 0.003215, 0.003215, + 0.003215, 0. , 0.006431, 0.344051}}}}}; +#elif EVALHYD_TESTING_OS == LINUX + rank_diag = {{{{{ 0.607717, 0. , 0. , 0. , 0. , 0.003215, + 0. , 0. , 0. , 0. , 0. , 0. , + 0.003215, 0. , 0.003215, 0. , 0. , 0. , + 0. , 0. , 0. , 0. , 0. , 0. , + 0. , 0. , 0. , 0. , 0. , 0.003215, + 0. , 0. , 0. , 0. , 0. , 0.006431, + 0. , 0. , 0.003215, 0.006431, 0. , 0. , + 0. , 0.003215, 0. , 0. , 0.003215, 0.003215, + 0.003215, 0. , 0.006431, 0.344051}}}}}; +#endif + EXPECT_TRUE( + xt::all(xt::isclose(metrics[0], rank_diag, 1e-04, 1e-06, true)) + ); + + // Delta scores + xt::xtensor<double, 4> ds; +#if EVALHYD_TESTING_OS == WINDOWS + ds = {{{{ 148.790164}}}}; +#elif EVALHYD_TESTING_OS == MACOS + ds = {{{{ 148.790164}}}}; +#elif EVALHYD_TESTING_OS == LINUX + ds = {{{{ 148.790164}}}}; +#endif + EXPECT_TRUE( + xt::all(xt::isclose(metrics[1], ds, 1e-04, 1e-07, true)) + ); + + + // Alpha scores + xt::xtensor<double, 4> as; +#if EVALHYD_TESTING_OS == WINDOWS + as = {{{{ 0.491481}}}}; +#elif EVALHYD_TESTING_OS == MACOS + as = {{{{ 0.491481}}}}; +#elif EVALHYD_TESTING_OS == LINUX + as = {{{{ 0.491481}}}}; +#endif + EXPECT_TRUE( + xt::all(xt::isclose(metrics[2], as, 1e-04, 1e-07, true)) + ); +} + TEST(ProbabilistTests, TestMasks) { // read in data @@ -455,7 +550,21 @@ TEST(ProbabilistTests, TestMasks) // check results are identical for (std::size_t m = 0; m < all_metrics_p.size(); m++) { - EXPECT_TRUE(xt::allclose(metrics_masked[m], metrics_subset[m])) + // --------------------------------------------------------------------- + // /!\ skip ranks-based metrics because it contains a random process + // for which setting the seed will not work because the time series + // lengths are different between "masked" and "subset", which + // results in different tensor shapes, and hence in different + // random ranks for ties + if ((all_metrics_p[m] == "RANK_DIAG") + || (all_metrics_p[m] == "DS") + || (all_metrics_p[m] == "AS")) + { + continue; + } + // --------------------------------------------------------------------- + + EXPECT_TRUE(xt::all(xt::isclose(metrics_masked[m], metrics_subset[m], 1e-04, 1e-07, true))) << "Failure for (" << all_metrics_p[m] << ")"; } } @@ -746,6 +855,20 @@ TEST(ProbabilistTests, TestBootstrap) // check results are identical for (std::size_t m = 0; m < all_metrics_p.size(); m++) { + // --------------------------------------------------------------------- + // /!\ skip ranks-based metrics because it contains a random process + // for which setting the seed will not work because the time series + // lengths are different between "bts" and "rep", which + // results in different tensor shapes, and hence in different + // random ranks for ties + if ((all_metrics_p[m] == "RANK_DIAG") + || (all_metrics_p[m] == "DS") + || (all_metrics_p[m] == "AS")) + { + continue; + } + // --------------------------------------------------------------------- + EXPECT_TRUE( xt::allclose( metrics_bts[m], metrics_rep[m]