diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp
index a3af81b522bd7fd5247b88febbc67dfac56613e9..248071dcdffbc5f58044813a1ef54848a57be116 100644
--- a/include/evalhyd/evald.hpp
+++ b/include/evalhyd/evald.hpp
@@ -35,14 +35,17 @@ namespace evalhyd
/// :Parameters:
///
/// q_obs: ``xt::xexpression<A>``
- /// Streamflow observations. Time steps with missing observation
+ /// Streamflow observations. Time steps with missing observations
/// must be assigned `NAN` values. Those time steps will be ignored
/// both in the observations and the predictions before the
/// *metrics* are computed.
/// shape: ({1+, ...}, time)
///
/// q_prd: ``xt::xexpression<A>``
- /// Streamflow predictions.
+ /// Streamflow predictions. Time steps with missing predictions
+ /// must be assigned `NAN` values. Those time steps will be ignored
+ /// both in the observations and the predictions before the
+ /// *metrics* are computed.
/// shape: ({1+, ...}, time)
///
/// metrics: ``std::vector<std::string>``
diff --git a/src/determinist/evaluator.hpp b/src/determinist/evaluator.hpp
index 95f034643a401d6e9adf9d239fe2e8fea613320c..124b9f06d3bc566a8bf7b761277b2426d0d77e9e 100644
--- a/src/determinist/evaluator.hpp
+++ b/src/determinist/evaluator.hpp
@@ -26,10 +26,20 @@ namespace evalhyd
template <class A>
class Evaluator
{
+ // define type for function to be used to "mask" (i.e. yield NaNs) the
+ // time steps for which the observations/predictions are missing data
+ typedef xt::xfunction<xt::detail::conditional_ternary,
+ xt::xfunction<xt::math::isnan_fun, const A&>,
+ xt::xscalar<float>,
+ const A&> nan_function;
+
private:
// members for input data
- const A& q_obs;
- const A& q_prd;
+ const A& _q_obs;
+ const A& _q_prd;
+ nan_function q_obs;
+ nan_function q_prd;
+
// members for computational elements
A mean_obs;
A mean_prd;
@@ -43,17 +53,33 @@ namespace evalhyd
// constructor method
Evaluator(const xt::xexpression<A>& obs,
const xt::xexpression<A>& prd) :
- q_obs{obs.derived_cast()},
- q_prd{prd.derived_cast()}
+ _q_obs{obs.derived_cast()},
+ _q_prd{prd.derived_cast()},
+ // "mask" (i.e. use NaN) observations where predictions are
+ // missing (i.e. NaN)
+ q_obs{xt::where(xt::isnan(prd.derived_cast()),
+ NAN, obs.derived_cast())},
+ // "mask" (i.e. use NaN) predictions where observations are
+ // missing (i.e. NaN)
+ q_prd{xt::where(xt::isnan(obs.derived_cast()),
+ NAN, prd.derived_cast())}
{
// check that data dimensions are compatible
- if (q_prd.dimension() != q_obs.dimension())
+ if (_q_prd.dimension() != _q_obs.dimension())
{
throw std::runtime_error(
"number of dimensions of 'sim' and 'obs' "
"must be identical"
);
}
+ else if (_q_prd.shape(_q_prd.dimension() - 1)
+ != _q_obs.shape( _q_obs.dimension() - 1))
+ {
+ throw std::runtime_error(
+ "length of time series of 'sim' and 'obs' "
+ "must be identical"
+ );
+ }
};
// members for evaluation metrics
@@ -103,10 +129,7 @@ namespace evalhyd
template <class A>
void Evaluator<A>::calc_mean_prd()
{
- mean_prd = xt::nanmean(
- xt::where(xt::isnan(q_obs), NAN, q_prd),
- -1, xt::keep_dims
- );
+ mean_prd = xt::nanmean(q_prd, -1, xt::keep_dims);
}
// Compute the quadratic error between observations and predictions.
@@ -157,9 +180,7 @@ namespace evalhyd
template <class A>
void Evaluator<A>::calc_quad_prd()
{
- quad_prd = xt::square(
- xt::where(xt::isnan(q_obs), NAN, q_prd) - mean_prd
- );
+ quad_prd = xt::square(q_prd - mean_prd);
}
// Compute the Pearson correlation coefficient.
@@ -191,8 +212,7 @@ namespace evalhyd
// calculate error in timing and dynamics $r_{pearson}$
// (Pearson's correlation coefficient)
auto r_num = xt::nansum(
- (xt::where(xt::isnan(q_obs), NAN, q_prd) - mean_prd)
- * (q_obs - mean_obs),
+ (q_prd - mean_prd) * (q_obs - mean_obs),
-1, xt::keep_dims
);
auto r_den = xt::sqrt(
@@ -218,8 +238,7 @@ namespace evalhyd
void Evaluator<A>::calc_bias()
{
// calculate $bias$
- bias = xt::nansum(xt::where(xt::isnan(q_obs), NAN, q_prd),
- -1, xt::keep_dims)
+ bias = xt::nansum(q_prd, -1, xt::keep_dims)
/ xt::nansum(q_obs, -1, xt::keep_dims);
}
@@ -287,8 +306,7 @@ namespace evalhyd
void Evaluator<A>::calc_KGE()
{
// calculate error in spread of flow $alpha$
- auto q_prd_ = xt::where(xt::isnan(q_obs), NAN, q_prd);
- auto alpha = detail::std_dev(q_prd_, mean_prd)
+ auto alpha = detail::std_dev(q_prd, mean_prd)
/ detail::std_dev(q_obs, mean_obs);
// compute KGE
@@ -326,8 +344,7 @@ namespace evalhyd
void Evaluator<A>::calc_KGEPRIME()
{
// calculate error in spread of flow $gamma$
- auto q_prd_ = xt::where(xt::isnan(q_obs), NAN, q_prd);
- auto gamma = (detail::std_dev(q_prd_, mean_prd) / mean_prd)
+ auto gamma = (detail::std_dev(q_prd, mean_prd) / mean_prd)
/ (detail::std_dev(q_obs, mean_obs) / mean_obs);
// compute KGE
diff --git a/tests/test_determinist.cpp b/tests/test_determinist.cpp
index 8ef2fec6e58d2d7ff1d40d82da2cfbb688966265..89c213e56172fa0483de3a2af82a20efbbaa7405 100644
--- a/tests/test_determinist.cpp
+++ b/tests/test_determinist.cpp
@@ -189,6 +189,12 @@ TEST(DeterministTests, TestMissingData)
// add some missing observations artificially by assigning NaN values
xt::view(observed, xt::all(), xt::range(0, 20)) = NAN;
+ // add some missing predictions artificially by assigning NaN values
+ xt::view(observed, 0, xt::range(20, 23)) = NAN;
+ xt::view(observed, 1, xt::range(20, 26)) = NAN;
+ xt::view(observed, 2, xt::range(20, 29)) = NAN;
+ xt::view(observed, 3, xt::range(20, 32)) = NAN;
+ xt::view(observed, 4, xt::range(20, 35)) = NAN;
// compute metrics with observations containing NaN values
std::vector<xt::xarray<double>> metrics_nan =
@@ -196,19 +202,27 @@ TEST(DeterministTests, TestMissingData)
observed, predicted, metrics
);
- // compute metrics on subset of observations and predictions
- // (i.e. eliminating region with NaN in observations manually)
- 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_sbs =
- evalhyd::evald<xt::xtensor<double, 2>>(
- obs, prd, metrics
- );
-
for (int m = 0; m < metrics.size(); m++) {
- EXPECT_TRUE(
- xt::allclose(metrics_nan[m], metrics_sbs[m])
- ) << "Failure for (" << metrics[m] << ")";
+ for (int p = 0; p < predicted.shape(0); p++) {
+ // compute metrics on subset of observations and predictions
+ // (i.e. eliminating region with NaN in observations manually)
+ xt::xtensor<double, 1> obs =
+ xt::view(observed, 0, xt::range(3*(p+1), _));
+ xt::xtensor<double, 1> prd =
+ xt::view(predicted, p, xt::range(3*(p+1), _));
+
+ std::vector<xt::xarray<double>> metrics_sbs =
+ evalhyd::evald<xt::xtensor<double, 1>>(
+ obs, prd, {metrics[m]}
+ );
+
+ // compare to check results are the same
+ EXPECT_TRUE(
+ xt::allclose(
+ xt::view(metrics_nan[m], p),
+ metrics_sbs[0]
+ )
+ ) << "Failure for (" << metrics[m] << ")";
+ }
}
}
\ No newline at end of file