diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp
index feda6c71f01e72ee8872f4b17253e62f7c7b87ee..0ec601138edd0f9ecc51cf959ea38b4369bf843a 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 9ff21ff208203a64bbcc43b615bd4e7ce0ac778a..61538595544ebb510d0f40a55e7f32f0c6c45e89 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 f392f9ad5885619219288e8fddfec4dda9754dbf..bc5bc45eae1edc00f6107d5cec6ceaac99ee8b61 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(