diff --git a/include/evalhyd/determinist.hpp b/include/evalhyd/determinist.hpp index 19dce0afb95ecb5a12fb2c403dab2f7e44b42bd9..06ad762ac4e786d0a27b36b42e170ca369ece6f2 100644 --- a/include/evalhyd/determinist.hpp +++ b/include/evalhyd/determinist.hpp @@ -4,6 +4,7 @@ #include <unordered_map> #include <vector> #include <xtensor/xexpression.hpp> +#include <xtensor/xarray.hpp> #include "utils.hpp" #include "determinist/evaluator.hpp" @@ -26,7 +27,7 @@ namespace evalhyd /// \return /// Vector of 1D array of metrics for each time series. template <class A> - std::vector<A> evald( + std::vector<xt::xarray<double>> evald( const xt::xexpression<A>& q_obs, const xt::xexpression<A>& q_prd, const std::vector<std::string>& metrics @@ -73,7 +74,7 @@ namespace evalhyd } // retrieve or compute requested metrics - std::vector<A> r; + std::vector<xt::xarray<double>> r; for ( const auto& metric : metrics ) { @@ -82,7 +83,7 @@ namespace evalhyd if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) evaluator.calc_NSE(); - r.emplace_back(evaluator.NSE); + r.emplace_back(xt::squeeze(evaluator.NSE)); } } diff --git a/include/evalhyd/probabilist.hpp b/include/evalhyd/probabilist.hpp index 9dbfe7d013268dd457fd96d82dbec43bbd0a7e8c..8c401695e383a54925857f89eb0b845bb0c9d6fd 100644 --- a/include/evalhyd/probabilist.hpp +++ b/include/evalhyd/probabilist.hpp @@ -5,6 +5,7 @@ #include <unordered_map> #include <vector> #include <xtensor/xtensor.hpp> +#include <xtensor/xarray.hpp> #include <xtensor/xview.hpp> #include "utils.hpp" @@ -39,7 +40,7 @@ namespace evalhyd /// shape: (1+, time) /// \return /// Vector of 2D array of metrics for each threshold. - std::vector<xt::xtensor<double, 3>> evalp( + std::vector<xt::xarray<double>> evalp( const xt::xtensor<double, 2>& q_obs, const xt::xtensor<double, 2>& q_prd, const std::vector<std::string>& metrics, @@ -107,7 +108,7 @@ namespace evalhyd } // retrieve or compute requested metrics - std::vector<xt::xtensor<double, 3>> r; + std::vector<xt::xarray<double>> r; for (const auto& metric : metrics) { @@ -116,42 +117,42 @@ namespace evalhyd if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) evaluator.calc_BS(); - r.emplace_back(evaluator.BS); + r.emplace_back(xt::squeeze(evaluator.BS)); } else if ( metric == "BSS" ) { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) evaluator.calc_BSS(); - r.emplace_back(evaluator.BSS); + r.emplace_back(xt::squeeze(evaluator.BSS)); } else if ( metric == "BS_CRD" ) { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) evaluator.calc_BS_CRD(); - r.emplace_back(evaluator.BS_CRD); + r.emplace_back(xt::squeeze(evaluator.BS_CRD)); } else if ( metric == "BS_LBD" ) { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) evaluator.calc_BS_LBD(); - r.emplace_back(evaluator.BS_LBD); + r.emplace_back(xt::squeeze(evaluator.BS_LBD)); } else if ( metric == "QS" ) { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) evaluator.calc_QS(); - r.emplace_back(evaluator.QS); + r.emplace_back(xt::squeeze(evaluator.QS)); } else if ( metric == "CRPS" ) { if (std::find(req_dep.begin(), req_dep.end(), metric) == req_dep.end()) evaluator.calc_CRPS(); - r.emplace_back(evaluator.CRPS); + r.emplace_back(xt::squeeze(evaluator.CRPS)); } } diff --git a/include/evalhyd/probabilist/evaluator.h b/include/evalhyd/probabilist/evaluator.h index be6908904adda60bf189084bcfd2d5184bea5ebe..26f0223430efc75759aa13eaa3461a4044e86fa8 100644 --- a/include/evalhyd/probabilist/evaluator.h +++ b/include/evalhyd/probabilist/evaluator.h @@ -55,12 +55,12 @@ namespace evalhyd xt::xtensor<double, 2> crps; // members for evaluation metrics - xt::xtensor<double, 3> BS; + xt::xtensor<double, 2> BS; xt::xtensor<double, 3> BS_CRD; xt::xtensor<double, 3> BS_LBD; - xt::xtensor<double, 3> BSS; - xt::xtensor<double, 3> QS; - xt::xtensor<double, 3> CRPS; + xt::xtensor<double, 2> BSS; + xt::xtensor<double, 2> QS; + xt::xtensor<double, 1> CRPS; // methods to compute derived data void calc_q_qnt(); diff --git a/include/evalhyd/probabilist/evaluator_brier.cpp b/include/evalhyd/probabilist/evaluator_brier.cpp index 10c2e338bab8c308790a6376b5a6c1ae80f6c98a..dcc8b7e17a6d63d20f9e4a7e8c47cde04b7db6cb 100644 --- a/include/evalhyd/probabilist/evaluator_brier.cpp +++ b/include/evalhyd/probabilist/evaluator_brier.cpp @@ -50,7 +50,7 @@ namespace evalhyd { // initialise output variable // shape: (subsets, thresholds, 1) - BS = xt::zeros<double>({n_msk, n_thr, std::size_t {1}}); + BS = xt::zeros<double>({n_msk, n_thr}); // compute variable one mask at a time to minimise memory imprint for (int m = 0; m < n_msk; m++) @@ -61,8 +61,8 @@ namespace evalhyd // calculate the mean over the time steps // $BS = \frac{1}{n} \sum_{k=1}^{n} (o_k - y_k)^2$ - xt::view(BS, m, xt::all(), xt::all()) = - xt::nanmean(bs_masked, -1, xt::keep_dims); + xt::view(BS, m, xt::all()) = + xt::nanmean(bs_masked, -1); } } @@ -306,7 +306,7 @@ namespace evalhyd { // declare and initialise output variable // shape: (subsets, thresholds, components) - BSS = xt::zeros<double>({n_msk, n_thr, std::size_t {1}}); + BSS = xt::zeros<double>({n_msk, n_thr}); // compute variable one mask at a time to minimise memory imprint for (int m = 0; m < n_msk; m++) @@ -332,8 +332,8 @@ namespace evalhyd // compute Brier skill score(s) // $BSS = \frac{1}{n} \sum_{k=1}^{n} 1 - \frac{bs}{bs_{ref}} - xt::view(BSS, m, xt::all(), xt::all()) = - xt::nanmean(1 - (bs_masked / bs_ref), -1, xt::keep_dims); + xt::view(BSS, m, xt::all()) = + xt::nanmean(1 - (bs_masked / bs_ref), -1); } } } diff --git a/include/evalhyd/probabilist/evaluator_quantiles.cpp b/include/evalhyd/probabilist/evaluator_quantiles.cpp index 42bb1b76fb1a533f49867b27bee7e41c3e784879..726d14d8f0a0a8ac0bcd2f36c485bb302e1146af 100644 --- a/include/evalhyd/probabilist/evaluator_quantiles.cpp +++ b/include/evalhyd/probabilist/evaluator_quantiles.cpp @@ -55,7 +55,7 @@ namespace evalhyd { // initialise output variable // shape: (subsets, quantiles, 1) - QS = xt::zeros<double>({n_msk, n_mbr, std::size_t {1}}); + QS = xt::zeros<double>({n_msk, n_mbr}); // compute variable one mask at a time to minimise memory imprint for (int m = 0; m < n_msk; m++) @@ -66,8 +66,8 @@ namespace evalhyd // calculate the mean over the time steps // $QS = \frac{1}{n} \sum_{k=1}^{n} qs$ - xt::view(QS, m, xt::all(), xt::all()) = - xt::nanmean(qs_masked, -1, xt::keep_dims); + xt::view(QS, m, xt::all()) = + xt::nanmean(qs_masked, -1); } } @@ -111,7 +111,7 @@ namespace evalhyd { // initialise output variable // shape: (subsets, thresholds, 1) - CRPS = xt::zeros<double>({n_msk, std::size_t {1}, std::size_t {1}}); + CRPS = xt::zeros<double>({n_msk}); // compute variable one mask at a time to minimise memory imprint for (int m = 0; m < n_msk; m++) @@ -122,8 +122,8 @@ namespace evalhyd // calculate the mean over the time steps // $CRPS = \frac{1}{n} \sum_{k=1}^{n} crps$ - xt::view(CRPS, m, xt::all(), xt::all()) = - xt::nanmean(crps_masked, -1, xt::keep_dims); + xt::view(CRPS, m) = + xt::nanmean(crps_masked, -1); } } } diff --git a/tests/test_determinist.cpp b/tests/test_determinist.cpp index c89553e57cff3a67463e1d684bc644562c4abf79..076d2072faa7f48f7056e779e4ed60d8209ab1ed 100644 --- a/tests/test_determinist.cpp +++ b/tests/test_determinist.cpp @@ -3,6 +3,7 @@ #include <vector> #include <gtest/gtest.h> #include <xtensor/xtensor.hpp> +#include <xtensor/xarray.hpp> #include <xtensor/xview.hpp> #include <xtensor/xmanipulation.hpp> #include <xtensor/xcsv.hpp> @@ -27,19 +28,19 @@ TEST(DeterministTests, TestNSE) { xt::xtensor<double, 1> predicted_1d = xt::row(predicted_2d, 0); // compute scores (both with 2D and 1D tensors) - std::vector<xt::xtensor<double, 2>> metrics_2d = + std::vector<xt::xarray<double>> metrics_2d = evalhyd::evald<xt::xtensor<double, 2>>( observed_2d, predicted_2d, {"NSE"} ); - std::vector<xt::xtensor<double, 1>> metrics_1d = + std::vector<xt::xarray<double>> metrics_1d = evalhyd::evald<xt::xtensor<double, 1>>( observed_1d, predicted_1d, {"NSE"} ); // check results (both with 2D and 1D tensors) - xt::xtensor<double, 2> nse_2d = - {{0.71891219}, {0.7190249} , {0.71835777}, {0.71810361}, {0.71776748}}; + xt::xtensor<double, 1> nse_2d = + {0.71891219, 0.7190249, 0.71835777, 0.71810361, 0.71776748}; EXPECT_TRUE(xt::allclose(metrics_2d[0], nse_2d)); xt::xtensor<double, 1> nse_1d = {0.71891219}; diff --git a/tests/test_probabilist.cpp b/tests/test_probabilist.cpp index 7bc1734953d9aba500179fdc6610e4e5229d4498..33a4009a767a149b6195995e32d3d08f3a1558bd 100644 --- a/tests/test_probabilist.cpp +++ b/tests/test_probabilist.cpp @@ -3,6 +3,7 @@ #include <vector> #include <gtest/gtest.h> #include <xtensor/xtensor.hpp> +#include <xtensor/xarray.hpp> #include <xtensor/xmanipulation.hpp> #include <xtensor/xcsv.hpp> @@ -22,7 +23,7 @@ TEST(ProbabilistTests, TestBrier) { // compute scores xt::xtensor<double, 1> thresholds = {690, 534, 445}; - std::vector<xt::xtensor<double, 3>> metrics = + std::vector<xt::xarray<double>> metrics = evalhyd::evalp( observed, predicted, {"BS", "BSS", "BS_CRD", "BS_LBD"}, thresholds @@ -30,30 +31,26 @@ TEST(ProbabilistTests, TestBrier) { // check results // Brier scores - xt::xtensor<double, 3> bs = - {{{0.10615136}, - {0.07395622}, - {0.08669186}}}; + xt::xtensor<double, 1> bs = + {0.10615136, 0.07395622, 0.08669186}; EXPECT_TRUE(xt::allclose(metrics[0], bs)); // Brier skill scores - xt::xtensor<double, 3> bss = - {{{0.5705594}, - {0.6661165}, - {0.5635126}}}; + xt::xtensor<double, 1> bss = + {0.5705594, 0.6661165, 0.5635126}; EXPECT_TRUE(xt::allclose(metrics[1], bss)); // Brier calibration-refinement decompositions - xt::xtensor<double, 3> bs_crd = - {{{0.011411758, 0.1524456, 0.2471852}, - {0.005532413, 0.1530793, 0.2215031}, - {0.010139431, 0.1220601, 0.1986125}}}; + xt::xtensor<double, 2> bs_crd = + {{0.011411758, 0.1524456, 0.2471852}, + {0.005532413, 0.1530793, 0.2215031}, + {0.010139431, 0.1220601, 0.1986125}}; EXPECT_TRUE(xt::allclose(metrics[2], bs_crd)); // Brier likelihood-base rate decompositions - xt::xtensor<double, 3> bs_lbd = - {{{0.012159881, 0.1506234, 0.2446149}, - {0.008031746, 0.1473869, 0.2133114}, - {0.017191279, 0.1048221, 0.1743227}}}; + xt::xtensor<double, 2> bs_lbd = + {{0.012159881, 0.1506234, 0.2446149}, + {0.008031746, 0.1473869, 0.2133114}, + {0.017191279, 0.1048221, 0.1743227}}; EXPECT_TRUE(xt::allclose(metrics[3], bs_lbd)); }