From 7ded7f393a14dcf790591771b24c89486567ec11 Mon Sep 17 00:00:00 2001 From: Thibault Hallouin <thibault.hallouin@inrae.fr> Date: Tue, 30 Aug 2022 14:59:15 +0200 Subject: [PATCH] move away from fully generic API for evald The use of `xexpression` and template parameter was showing its limits as the data type could not be different between the different n-dim inputs (e.g. masks had to be as double rather than bool). Since there is no real need to be able to use `evald` on `xarray`, the genericity of its API is reduced so that it now expects inputs to be tensor only. The number of dimensions for the tensors becomes a template parameter. --- include/evalhyd/evald.hpp | 76 ++++++++++++------------ src/determinist/evaluator.hpp | 106 +++++++++++++++------------------- tests/test_determinist.cpp | 35 ++++------- 3 files changed, 95 insertions(+), 122 deletions(-) diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp index feda6c7..0ec6011 100644 --- a/include/evalhyd/evald.hpp +++ b/include/evalhyd/evald.hpp @@ -136,38 +136,34 @@ namespace evalhyd /// evalhyd::evald(obs, prd, {"NSE"}, "none", 1, -9, msk); /// /// \endrst - template <class C> + template <std::size_t N> std::vector<xt::xarray<double>> evald( - const xt::xexpression<C>& q_obs, - const xt::xexpression<C>& q_prd, + const xt::xtensor<double, N>& q_obs, + const xt::xtensor<double, N>& q_prd, const std::vector<std::string>& metrics, const std::string& transform = "none", const double exponent = 1, double epsilon = -9, - const xt::xexpression<C>& t_msk = C() + const xt::xtensor<bool, N>& t_msk = {} ) { - auto q_obs_ = q_obs.derived_cast(); - auto q_prd_ = q_prd.derived_cast(); - auto t_msk_ = t_msk.derived_cast(); - // initialise a mask if none provided // (corresponding to no temporal subset) - C msk; - if (t_msk_.size() < 1) - msk = xt::ones<bool>(q_obs_.shape()); + xt::xtensor<bool, N> msk; + if (t_msk.size() < 1) + msk = xt::ones<bool>(q_obs.shape()); else msk = std::move(t_msk.derived_cast()); // check that observations, predictions, and masks dimensions are compatible - if (q_obs_.dimension() != q_prd_.dimension()) + if (q_obs.dimension() != q_prd.dimension()) { throw std::runtime_error( "observations and predictions feature different numbers " "of dimensions" ); } - if (q_obs_.dimension() != msk.dimension()) + if (q_obs.dimension() != msk.dimension()) { throw std::runtime_error( "observations and masks feature different numbers " @@ -175,65 +171,67 @@ namespace evalhyd ); } - auto obs_shp = q_obs_.shape(); - auto prd_shp = q_prd_.shape(); + auto obs_shp = q_obs.shape(); + auto prd_shp = q_prd.shape(); auto msk_shp = msk.shape(); + // check observations, predictions, and masks are broadcastable if (obs_shp != prd_shp) - { - // check observations, predictions, and masks are broadcastable - for (int a = 0; a < obs_shp.size(); a++) - { + for (int a = 0; a < q_obs.dimension(); a++) if (obs_shp[a] != prd_shp[a]) if ((obs_shp[a] != 1) and (prd_shp[a] != 1)) throw std::runtime_error( "observations and predictions are not " "broadcastable" ); + + if (obs_shp != msk_shp) + for (int a = 0; a < q_obs.dimension(); a++) if (obs_shp[a] != msk_shp[a]) if ((obs_shp[a] != 1) and (msk_shp[a] != 1)) throw std::runtime_error( "observations and masks are not broadcastable" ); + + if (prd_shp != msk_shp) + for (int a = 0; a < q_prd.dimension(); a++) if (prd_shp[a] != msk_shp[a]) if ((prd_shp[a] != 1) and (msk_shp[a] != 1)) throw std::runtime_error( "predictions and masks are not broadcastable" ); - } - } // apply streamflow transformation if requested - C obs; - C prd; + xt::xtensor<double, N> obs; + xt::xtensor<double, N> prd; if ( transform == "none" || (transform == "pow" && exponent == 1)) { - obs = std::move(q_obs_); - prd = std::move(q_prd_); + obs = std::move(q_obs); + prd = std::move(q_prd); } else if ( transform == "sqrt" ) { - obs = xt::sqrt(q_obs_); - prd = xt::sqrt(q_prd_); + obs = xt::sqrt(q_obs); + prd = xt::sqrt(q_prd); } else if ( transform == "inv" ) { if ( epsilon == -9 ) // determine an epsilon value to avoid zero divide - epsilon = xt::mean(q_obs_)() * 0.01; + epsilon = xt::mean(q_obs)() * 0.01; - obs = 1. / (q_obs_ + epsilon); - prd = 1. / (q_prd_ + epsilon); + obs = 1. / (q_obs + epsilon); + prd = 1. / (q_prd + epsilon); } else if ( transform == "log" ) { if ( epsilon == -9 ) // determine an epsilon value to avoid log zero - epsilon = xt::mean(q_obs_)() * 0.01; + epsilon = xt::mean(q_obs)() * 0.01; - obs = xt::log(q_obs_ + epsilon); - prd = xt::log(q_prd_ + epsilon); + obs = xt::log(q_obs + epsilon); + prd = xt::log(q_prd + epsilon); } else if ( transform == "pow" ) { @@ -241,15 +239,15 @@ namespace evalhyd { if ( epsilon == -9 ) // determine an epsilon value to avoid zero divide - epsilon = xt::mean(q_obs_)() * 0.01; + epsilon = xt::mean(q_obs)() * 0.01; - obs = xt::pow(q_obs_ + epsilon, exponent); - prd = xt::pow(q_prd_ + epsilon, exponent); + obs = xt::pow(q_obs + epsilon, exponent); + prd = xt::pow(q_prd + epsilon, exponent); } else { - obs = xt::pow(q_obs_, exponent); - prd = xt::pow(q_prd_, exponent); + obs = xt::pow(q_obs, exponent); + prd = xt::pow(q_prd, exponent); } } else @@ -266,7 +264,7 @@ namespace evalhyd ); // instantiate determinist evaluator - eh::determinist::Evaluator<C> evaluator(obs, prd, msk); + eh::determinist::Evaluator<N> evaluator(obs, prd, msk); // declare maps for memoisation purposes std::unordered_map<std::string, std::vector<std::string>> elt; diff --git a/src/determinist/evaluator.hpp b/src/determinist/evaluator.hpp index 9ff21ff..6153859 100644 --- a/src/determinist/evaluator.hpp +++ b/src/determinist/evaluator.hpp @@ -23,7 +23,7 @@ namespace evalhyd namespace determinist { - template <class C> + template <std::size_t N> class Evaluator { // define type for function to be used to "mask" (i.e. yield NaNs) the @@ -34,63 +34,53 @@ namespace evalhyd xt::detail::bitwise_or, xt::xfunction< xt::math::isnan_fun, - const C& + const xt::xtensor<double, N>& >, xt::xfunction< xt::detail::logical_not, - const C& + const xt::xtensor<bool, N>& > >, xt::xscalar<float>, - const C& + const xt::xtensor<double, N>& > nan_function; private: // members for input data - const C& _q_obs; - const C& _q_prd; + const xt::xtensor<double, N>& _q_obs; + const xt::xtensor<double, N>& _q_prd; nan_function q_obs; nan_function q_prd; // members for computational elements - C mean_obs; - C mean_prd; - C quad_err; - C quad_obs; - C quad_prd; - C r_pearson; - C bias; + xt::xtensor<double, N> mean_obs; + xt::xtensor<double, N> mean_prd; + xt::xtensor<double, N> quad_err; + xt::xtensor<double, N> quad_obs; + xt::xtensor<double, N> quad_prd; + xt::xtensor<double, N> r_pearson; + xt::xtensor<double, N> bias; public: // constructor method - Evaluator(const xt::xexpression<C>& obs, - const xt::xexpression<C>& prd, - const xt::xexpression<C>& msk) : - _q_obs{obs.derived_cast()}, - _q_prd{prd.derived_cast()}, + Evaluator(const xt::xtensor<double, N>& obs, + const xt::xtensor<double, N>& prd, + const xt::xtensor<bool, N>& msk) : + _q_obs{obs}, + _q_prd{prd}, // "mask" (i.e. use NaN) observations where predictions are // missing (i.e. NaN) - q_obs{ - xt::where( - xt::isnan(prd.derived_cast()) | !msk.derived_cast(), - NAN, obs.derived_cast() - ) - }, + q_obs{xt::where(xt::isnan(prd) | !msk, NAN, obs)}, // "mask" (i.e. use NaN) predictions where observations are // missing (i.e. NaN) - q_prd{ - xt::where( - xt::isnan(obs.derived_cast()) | !msk.derived_cast(), - NAN, prd.derived_cast() - ) - } + q_prd{xt::where(xt::isnan(obs) | !msk, NAN, prd)} {}; // members for evaluation metrics - C RMSE; - C NSE; - C KGE; - C KGEPRIME; + xt::xtensor<double, N> RMSE; + xt::xtensor<double, N> NSE; + xt::xtensor<double, N> KGE; + xt::xtensor<double, N> KGEPRIME; // methods to compute elements void calc_mean_obs(); @@ -116,8 +106,8 @@ namespace evalhyd // \assign mean_obs: // Mean observed streamflow. // shape: ({... ,} 1) - template <class C> - void Evaluator<C>::calc_mean_obs() + template <std::size_t N> + void Evaluator<N>::calc_mean_obs() { mean_obs = xt::nanmean(q_obs, -1, xt::keep_dims); } @@ -130,8 +120,8 @@ namespace evalhyd // \assign mean_prd: // Mean predicted streamflow. // shape: ({... ,} 1) - template <class C> - void Evaluator<C>::calc_mean_prd() + template <std::size_t N> + void Evaluator<N>::calc_mean_prd() { mean_prd = xt::nanmean(q_prd, -1, xt::keep_dims); } @@ -147,8 +137,8 @@ namespace evalhyd // \assign quad_err: // Quadratic errors between observations and predictions. // shape: ({... ,} time) - template <class C> - void Evaluator<C>::calc_quad_err() + template <std::size_t N> + void Evaluator<N>::calc_quad_err() { quad_err = xt::square(q_obs - q_prd); } @@ -164,8 +154,8 @@ namespace evalhyd // \assign quad_obs: // Quadratic errors between observations and mean observation. // shape: ({... ,} time) - template <class C> - void Evaluator<C>::calc_quad_obs() + template <std::size_t N> + void Evaluator<N>::calc_quad_obs() { quad_obs = xt::square(q_obs - mean_obs); } @@ -181,8 +171,8 @@ namespace evalhyd // \assign quad_prd: // Quadratic errors between predictions and mean prediction. // shape: ({... ,} time) - template <class C> - void Evaluator<C>::calc_quad_prd() + template <std::size_t N> + void Evaluator<N>::calc_quad_prd() { quad_prd = xt::square(q_prd - mean_prd); } @@ -210,8 +200,8 @@ namespace evalhyd // \assign r_pearson: // Pearson correlation coefficients. // shape: ({... ,} time) - template <class C> - void Evaluator<C>::calc_r_pearson() + template <std::size_t N> + void Evaluator<N>::calc_r_pearson() { // calculate error in timing and dynamics $r_{pearson}$ // (Pearson's correlation coefficient) @@ -238,8 +228,8 @@ namespace evalhyd // \assign bias: // Biases. // shape: ({... ,} 1) - template <class C> - void Evaluator<C>::calc_bias() + template <std::size_t N> + void Evaluator<N>::calc_bias() { // calculate $bias$ bias = xt::nansum(q_prd, -1, xt::keep_dims) @@ -254,8 +244,8 @@ namespace evalhyd // \assign RMSE: // Root-mean-square errors. // shape: ({... ,} 1) - template <class C> - void Evaluator<C>::calc_RMSE() + template <std::size_t N> + void Evaluator<N>::calc_RMSE() { // compute RMSE RMSE = xt::sqrt(xt::nanmean(quad_err, -1, xt::keep_dims)); @@ -272,12 +262,12 @@ namespace evalhyd // \assign NSE: // Nash-Sutcliffe efficiencies. // shape: ({... ,} 1) - template <class C> - void Evaluator<C>::calc_NSE() + template <std::size_t N> + void Evaluator<N>::calc_NSE() { // compute squared errors operands - C f_num = xt::nansum(quad_err, -1, xt::keep_dims); - C f_den = xt::nansum(quad_obs, -1, xt::keep_dims); + xt::xtensor<double, N> f_num = xt::nansum(quad_err, -1, xt::keep_dims); + xt::xtensor<double, N> f_den = xt::nansum(quad_obs, -1, xt::keep_dims); // compute NSE NSE = 1 - (f_num / f_den); @@ -306,8 +296,8 @@ namespace evalhyd // \assign KGE: // Kling-Gupta efficiencies. // shape: ({... ,} 1) - template <class C> - void Evaluator<C>::calc_KGE() + template <std::size_t N> + void Evaluator<N>::calc_KGE() { // calculate error in spread of flow $alpha$ auto alpha = detail::std_dev(q_prd, mean_prd) @@ -344,8 +334,8 @@ namespace evalhyd // \assign KGEPRIME: // Modified Kling-Gupta efficiencies. // shape: ({... ,} 1) - template <class C> - void Evaluator<C>::calc_KGEPRIME() + template <std::size_t N> + void Evaluator<N>::calc_KGEPRIME() { // calculate error in spread of flow $gamma$ auto gamma = (detail::std_dev(q_prd, mean_prd) / mean_prd) diff --git a/tests/test_determinist.cpp b/tests/test_determinist.cpp index f392f9a..bc5bc45 100644 --- a/tests/test_determinist.cpp +++ b/tests/test_determinist.cpp @@ -26,7 +26,7 @@ TEST(DeterministTests, TestMetrics2D) // compute scores (with 2D tensors) std::vector<xt::xarray<double>> metrics = - evalhyd::evald<xt::xtensor<double, 2>>( + evalhyd::evald<2>( observed, predicted, {"RMSE", "NSE", "KGE", "KGEPRIME"} ); @@ -80,7 +80,7 @@ TEST(DeterministTests, TestMetrics1D) // compute scores (with 1D tensors) std::vector<xt::xarray<double>> metrics = - evalhyd::evald<xt::xtensor<double, 1>>( + evalhyd::evald<1>( observed, predicted, {"RMSE", "NSE", "KGE", "KGEPRIME"} ); @@ -114,9 +114,7 @@ TEST(DeterministTests, TestTransform) // compute and check results on square-rooted streamflow series std::vector<xt::xarray<double>> metrics = - evalhyd::evald<xt::xtensor<double, 2>>( - observed, predicted, {"NSE"}, "sqrt" - ); + evalhyd::evald<2>(observed, predicted, {"NSE"}, "sqrt"); xt::xtensor<double, 2> nse_sqrt = {{0.882817}, @@ -127,10 +125,7 @@ TEST(DeterministTests, TestTransform) EXPECT_TRUE(xt::allclose(metrics[0], nse_sqrt)); // compute and check results on inverted streamflow series - metrics = - evalhyd::evald<xt::xtensor<double, 2>>( - observed, predicted, {"NSE"}, "inv" - ); + metrics = evalhyd::evald<2>(observed, predicted, {"NSE"}, "inv"); xt::xtensor<double, 2> nse_inv = {{0.737323}, @@ -141,10 +136,7 @@ TEST(DeterministTests, TestTransform) EXPECT_TRUE(xt::allclose(metrics[0], nse_inv)); // compute and check results on square-rooted streamflow series - metrics = - evalhyd::evald<xt::xtensor<double, 2>>( - observed, predicted, {"NSE"}, "log" - ); + metrics = evalhyd::evald<2>(observed, predicted, {"NSE"}, "log"); xt::xtensor<double, 2> nse_log = {{0.893344}, @@ -155,10 +147,7 @@ TEST(DeterministTests, TestTransform) EXPECT_TRUE(xt::allclose(metrics[0], nse_log)); // compute and check results on power-transformed streamflow series - metrics = - evalhyd::evald<xt::xtensor<double, 2>>( - observed, predicted, {"NSE"}, "pow", 0.2 - ); + metrics = evalhyd::evald<2>(observed, predicted, {"NSE"}, "pow", 0.2); xt::xtensor<double, 2> nse_pow = {{0.899207}, @@ -192,14 +181,14 @@ TEST(DeterministTests, TestMasks) {"RMSE", "NSE", "KGE", "KGEPRIME"}; std::vector<xt::xarray<double>> metrics_masked = - evalhyd::evald(observed, predicted, metrics, "none", 1, -9, masks); + evalhyd::evald<2>(observed, predicted, metrics, "none", 1, -9, masks); // compute scores on pre-computed subset of whole record xt::xtensor<double, 2> obs = xt::view(observed, xt::all(), xt::range(20, _)); xt::xtensor<double, 2> prd = xt::view(predicted, xt::all(), xt::range(20, _)); std::vector<xt::xarray<double>> metrics_subset = - evalhyd::evald(obs, prd, metrics); + evalhyd::evald<2>(obs, prd, metrics); // check results are identical for (int m = 0; m < metrics.size(); m++) @@ -237,9 +226,7 @@ TEST(DeterministTests, TestMissingData) // compute metrics with observations containing NaN values std::vector<xt::xarray<double>> metrics_nan = - evalhyd::evald<xt::xtensor<double, 2>>( - observed, predicted, metrics - ); + evalhyd::evald<2>(observed, predicted, metrics); for (int m = 0; m < metrics.size(); m++) { for (int p = 0; p < predicted.shape(0); p++) { @@ -251,9 +238,7 @@ TEST(DeterministTests, TestMissingData) xt::view(predicted, p, xt::range(20+(3*(p+1)), _)); std::vector<xt::xarray<double>> metrics_sbs = - evalhyd::evald<xt::xtensor<double, 1>>( - obs, prd, {metrics[m]} - ); + evalhyd::evald<1>(obs, prd, {metrics[m]}); // compare to check results are the same EXPECT_TRUE( -- GitLab