diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp
index c3cc7697fd41e40ccb2a621ad9af23613f31d59a..44bd3fff6da8fc373db81f71f823ed141ff55e2f 100644
--- a/include/evalhyd/evald.hpp
+++ b/include/evalhyd/evald.hpp
@@ -89,6 +89,14 @@ namespace evalhyd
/// `Pushpalatha et al. (2012)
/// <https://doi.org/10.1016/j.jhydrol.2011.11.055>`_.
///
+ /// t_msk: ``xt::xexpression<A>``, optional
+ /// Mask(s) used to generate temporal subsets of the whole streamflow
+ /// time series (where True/False is used for the time steps to
+ /// include/discard in a given subset). Masks must feature the same
+ /// number of dimensions as observations and predictions, and it must
+ /// broadcastable with both of them.
+ /// shape: ({... ,} time)
+ ///
/// :Returns:
///
/// ``std::vector<xt::xarray<double>>``
@@ -127,29 +135,47 @@ namespace evalhyd
const xt::xexpression<A>& q_obs,
const xt::xexpression<A>& q_prd,
const std::vector<std::string>& metrics,
- const std::string& transform="none",
- const double exponent=1,
- double epsilon=-9
+ const std::string& transform = "none",
+ const double exponent = 1,
+ double epsilon = -9,
+ const xt::xexpression<A>& t_msk = A()
)
{
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)
+ A msk;
+ if (t_msk_.size() < 1)
+ msk = xt::ones<bool>(q_obs_.shape());
+ else
+ msk = std::move(t_msk.derived_cast());
- // check that observations and predictions dimensions are compatible
- if (q_prd_.dimension() != q_obs_.dimension())
+ // check that observations, predictions, and masks dimensions are compatible
+ 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())
+ {
+ throw std::runtime_error(
+ "observations and masks feature different numbers "
+ "of dimensions"
+ );
+ }
auto obs_shp = q_obs_.shape();
auto prd_shp = q_prd_.shape();
+ auto msk_shp = msk.shape();
if (obs_shp != prd_shp)
{
- // check observations and predictions are broadcastable
+ // check observations, predictions, and masks are broadcastable
for (int a = 0; a < obs_shp.size(); a++)
{
if (obs_shp[a] != prd_shp[a])
@@ -158,6 +184,16 @@ namespace evalhyd
"observations and predictions are not "
"broadcastable"
);
+ 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[a] != msk_shp[a])
+ if ((prd_shp[a] != 1) and (msk_shp[a] != 1))
+ throw std::runtime_error(
+ "predictions and masks are not broadcastable"
+ );
}
}
@@ -224,7 +260,7 @@ namespace evalhyd
);
// instantiate determinist evaluator
- eh::determinist::Evaluator<A> evaluator(obs, prd);
+ eh::determinist::Evaluator<A> 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 c83591d8e1c4cda9692a71bb5d8df7708a75298b..2baaa812eb11ae9f9b02c39fb248f55da8bdf062 100644
--- a/src/determinist/evaluator.hpp
+++ b/src/determinist/evaluator.hpp
@@ -28,10 +28,22 @@ namespace evalhyd
{
// 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;
+ typedef xt::xfunction<
+ xt::detail::conditional_ternary,
+ xt::xfunction<
+ xt::detail::bitwise_or,
+ xt::xfunction<
+ xt::math::isnan_fun,
+ const A&
+ >,
+ xt::xfunction<
+ xt::detail::logical_not,
+ const A&
+ >
+ >,
+ xt::xscalar<float>,
+ const A&
+ > nan_function;
private:
// members for input data
@@ -52,17 +64,26 @@ namespace evalhyd
public:
// constructor method
Evaluator(const xt::xexpression<A>& obs,
- const xt::xexpression<A>& prd) :
+ const xt::xexpression<A>& prd,
+ const xt::xexpression<A>& msk) :
_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())},
+ q_obs{
+ xt::where(
+ xt::isnan(prd.derived_cast()) | !msk.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())}
+ q_prd{
+ xt::where(
+ xt::isnan(obs.derived_cast()) | !msk.derived_cast(),
+ NAN, prd.derived_cast()
+ )
+ }
{};
// members for evaluation metrics
diff --git a/tests/test_determinist.cpp b/tests/test_determinist.cpp
index 1980c4d648e38ad23b9dae7dc9f63e32dcc13f5f..f392f9ad5885619219288e8fddfec4dda9754dbf 100644
--- a/tests/test_determinist.cpp
+++ b/tests/test_determinist.cpp
@@ -170,6 +170,45 @@ TEST(DeterministTests, TestTransform)
}
+TEST(DeterministTests, TestMasks)
+{
+ // read in data
+ std::ifstream ifs;
+ ifs.open("./data/q_obs.csv");
+ xt::xtensor<double, 2> observed = xt::load_csv<int>(ifs);
+ ifs.close();
+
+ ifs.open("./data/q_prd.csv");
+ xt::xtensor<double, 2> predicted = xt::load_csv<double>(ifs);
+ ifs.close();
+
+ // generate temporal subset by dropping 20 first time steps
+ xt::xtensor<double, 2> masks =
+ xt::ones<bool>({std::size_t {1}, std::size_t {observed.size()}});
+ xt::view(masks, 0, xt::range(0, 20)) = 0;
+
+ // compute scores using masks to subset whole record
+ std::vector<std::string> metrics =
+ {"RMSE", "NSE", "KGE", "KGEPRIME"};
+
+ std::vector<xt::xarray<double>> metrics_masked =
+ evalhyd::evald(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);
+
+ // check results are identical
+ for (int m = 0; m < metrics.size(); m++)
+ {
+ EXPECT_TRUE(xt::allclose(metrics_masked[m], metrics_subset[m]))
+ << "Failure for (" << metrics[m] << ")";
+ }
+}
+
TEST(DeterministTests, TestMissingData)
{
std::vector<std::string> metrics =