diff --git a/include/evalhyd/detail/determinist/elements.hpp b/include/evalhyd/detail/determinist/elements.hpp
index 74007923355a95f16279450a2fed03d9c6c68969..fb8dc43fa7ee5dcfde87093cad0c960e97c2c5d5 100644
--- a/include/evalhyd/detail/determinist/elements.hpp
+++ b/include/evalhyd/detail/determinist/elements.hpp
@@ -24,7 +24,7 @@ namespace evalhyd
/// shape: (1, time)
/// \param t_msk
/// Temporal subsets of the whole record.
- /// shape: (subsets, series, time)
+ /// shape: (series, subsets, time)
/// \param b_exp
/// Boostrap samples.
/// shape: (samples, time slice)
@@ -54,7 +54,8 @@ namespace evalhyd
{
// apply the mask
// (using NaN workaround until reducers work on masked_view)
- auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+ auto obs_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ q_obs, NAN);
// compute variable one sample at a time
for (std::size_t e = 0; e < n_exp; e++)
@@ -76,7 +77,7 @@ namespace evalhyd
/// shape: (series, time)
/// \param t_msk
/// Temporal subsets of the whole record.
- /// shape: (subsets, series, time)
+ /// shape: (series, subsets, time)
/// \param b_exp
/// Boostrap samples.
/// shape: (samples, time slice)
@@ -106,7 +107,8 @@ namespace evalhyd
{
// apply the mask
// (using NaN workaround until reducers work on masked_view)
- auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
+ auto prd_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ q_prd, NAN);
// compute variable one sample at a time
for (std::size_t e = 0; e < n_exp; e++)
@@ -151,7 +153,7 @@ namespace evalhyd
/// shape: (subsets, samples, series, 1)
/// \param t_msk
/// Temporal subsets of the whole record.
- /// shape: (subsets, series, time)
+ /// shape: (series, subsets, time)
/// \param n_srs
/// Number of prediction series.
/// \param n_tim
@@ -181,7 +183,8 @@ namespace evalhyd
{
// apply the mask
// (using NaN workaround until reducers work on masked_view)
- auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+ auto obs_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ q_obs, NAN);
for (std::size_t e = 0; e < n_exp; e++)
{
@@ -204,7 +207,7 @@ namespace evalhyd
/// shape: (subsets, samples, series, 1)
/// \param t_msk
/// Temporal subsets of the whole record.
- /// shape: (subsets, series, time)
+ /// shape: (series, subsets, time)
/// \param n_srs
/// Number of prediction series.
/// \param n_tim
@@ -234,7 +237,8 @@ namespace evalhyd
{
// apply the mask
// (using NaN workaround until reducers work on masked_view)
- auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
+ auto prd_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ q_prd, NAN);
for (std::size_t e = 0; e < n_exp; e++)
{
@@ -269,7 +273,7 @@ namespace evalhyd
/// shape: (subsets, samples, series, time)
/// \param t_msk
/// Temporal subsets of the whole record.
- /// shape: (subsets, series, time)
+ /// shape: (series, subsets, time)
/// \param b_exp
/// Boostrap samples.
/// shape: (samples, time slice)
@@ -306,8 +310,10 @@ namespace evalhyd
{
// apply the mask
// (using NaN workaround until reducers work on masked_view)
- auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
- auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+ auto prd_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ q_prd, NAN);
+ auto obs_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ q_obs, NAN);
// compute variable one sample at a time
for (std::size_t e = 0; e < n_exp; e++)
@@ -349,7 +355,7 @@ namespace evalhyd
/// shape: (subsets, samples, series, 1)
/// \param t_msk
/// Temporal subsets of the whole record.
- /// shape: (subsets, series, time)
+ /// shape: (series, subsets, time)
/// \param b_exp
/// Boostrap samples.
/// shape: (samples, time slice)
@@ -383,8 +389,10 @@ namespace evalhyd
{
// apply the mask
// (using NaN workaround until reducers work on masked_view)
- auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
- auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+ auto prd_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ q_prd, NAN);
+ auto obs_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ q_obs, NAN);
// compute variable one sample at a time
for (std::size_t e = 0; e < n_exp; e++)
@@ -410,7 +418,7 @@ namespace evalhyd
/// shape: (series, time)
/// \param t_msk
/// Temporal subsets of the whole record.
- /// shape: (subsets, series, time)
+ /// shape: (series, subsets, time)
/// \param b_exp
/// Boostrap samples.
/// shape: (samples, time slice)
@@ -442,8 +450,10 @@ namespace evalhyd
{
// apply the mask
// (using NaN workaround until reducers work on masked_view)
- auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
- auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+ auto prd_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ q_prd, NAN);
+ auto obs_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ q_obs, NAN);
// compute variable one sample at a time
for (std::size_t e = 0; e < n_exp; e++)
diff --git a/include/evalhyd/detail/determinist/evaluator.hpp b/include/evalhyd/detail/determinist/evaluator.hpp
index e17bfd54fbe3a5b40ae8155d5a21488acf4d04ea..6314795e7096cfc29a69c988e5fcfdd0a5f64605 100644
--- a/include/evalhyd/detail/determinist/evaluator.hpp
+++ b/include/evalhyd/detail/determinist/evaluator.hpp
@@ -19,14 +19,15 @@ namespace evalhyd
{
namespace determinist
{
- template <class XD2, class XB2>
+ template <class XD2, class XB3>
class Evaluator
{
private:
// members for input data
const XD2& q_obs;
const XD2& q_prd;
- xt::xtensor<bool, 3> t_msk;
+ // members for optional input data
+ XB3 t_msk;
const std::vector<xt::xkeep_slice<int>>& b_exp;
// members for dimensions
@@ -56,7 +57,7 @@ namespace evalhyd
{
if (!mean_obs.has_value())
{
- mean_obs = elements::calc_mean_obs<XD2>(
+ mean_obs = elements::calc_mean_obs(
q_obs, t_msk, b_exp, n_srs, n_msk, n_exp
);
}
@@ -149,28 +150,34 @@ namespace evalhyd
// constructor method
Evaluator(const XD2& obs,
const XD2& prd,
- const XB2& msk,
+ const XB3& msk,
const std::vector<xt::xkeep_slice<int>>& exp) :
- q_obs{obs}, q_prd{prd}, b_exp{exp}
+ q_obs{obs}, q_prd{prd}, t_msk{msk}, b_exp{exp}
{
+ // initialise a mask if none provided
+ // (corresponding to no temporal subset)
+ if (msk.size() < 1)
+ {
+ t_msk = xt::ones<bool>(
+ {q_prd.shape(0), std::size_t {1}, q_prd.shape(1)}
+ );
+ }
+
// determine size for recurring dimensions
- n_tim = q_prd.shape(1);
n_srs = q_prd.shape(0);
- n_msk = msk.shape(0);
+ n_tim = q_prd.shape(1);
+ n_msk = msk.shape(1);
n_exp = b_exp.size();
// drop time steps where observations or predictions are NaN
- auto obs_nan = xt::isnan(q_obs);
- auto prd_nan = xt::isnan(q_prd);
-
- t_msk = xt::ones<bool>({n_msk, n_srs, n_tim});
- xt::view(t_msk, xt::all()) =
- xt::view(msk, xt::all(), xt::newaxis(), xt::all());
- for (std::size_t m = 0; m < n_msk; m++)
+ for (std::size_t s = 0; s < n_srs; s++)
{
- xt::view(t_msk, m) =
- xt::where(obs_nan || prd_nan,
- false, xt::view(t_msk, m));
+ auto obs_nan = xt::isnan(xt::view(q_obs, 0));
+ auto prd_nan = xt::isnan(xt::view(q_prd, s));
+
+ auto msk_nan = xt::where(obs_nan || prd_nan)[0];
+
+ xt::view(t_msk, s, xt::all(), xt::keep(msk_nan)) = false;
}
};
diff --git a/include/evalhyd/detail/determinist/metrics.hpp b/include/evalhyd/detail/determinist/metrics.hpp
index 8ee318f0cfadff8f4c4933c79c85f1fbf3c5a6d1..550c69561efd023d56f7ad820242134d95341560 100644
--- a/include/evalhyd/detail/determinist/metrics.hpp
+++ b/include/evalhyd/detail/determinist/metrics.hpp
@@ -23,7 +23,7 @@ namespace evalhyd
/// shape: (series, time)
/// \param t_msk
/// Temporal subsets of the whole record.
- /// shape: (subsets, series, time)
+ /// shape: (series, subsets, time)
/// \param b_exp
/// Boostrap samples.
/// shape: (samples, time slice)
@@ -53,7 +53,8 @@ namespace evalhyd
{
// apply the mask
// (using NaN workaround until reducers work on masked_view)
- auto quad_err_masked = xt::where(xt::view(t_msk, m), quad_err, NAN);
+ auto quad_err_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ quad_err, NAN);
// compute variable one sample at a time
for (std::size_t e = 0; e < n_exp; e++)
@@ -76,7 +77,7 @@ namespace evalhyd
/// shape: (subsets, samples, series, time)
/// \param t_msk
/// Temporal subsets of the whole record.
- /// shape: (subsets, series, time)
+ /// shape: (series, subsets, time)
/// \param b_exp
/// Boostrap samples.
/// shape: (samples, time slice)
@@ -106,7 +107,8 @@ namespace evalhyd
{
// apply the mask
// (using NaN workaround until reducers work on masked_view)
- auto quad_err_masked = xt::where(xt::view(t_msk, m), quad_err, NAN);
+ auto quad_err_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ quad_err, NAN);
// compute variable one sample at a time
for (std::size_t e = 0; e < n_exp; e++)
diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp
index 9bdd3658a697b0d739e63fccc08b186dcd73ed37..31bc5440410e44783d7af37c80795a2549f14866 100644
--- a/include/evalhyd/evald.hpp
+++ b/include/evalhyd/evald.hpp
@@ -32,9 +32,9 @@ namespace evalhyd
/// (e.g. ``xt::xtensor<double, 2>``, ``xt::pytensor<double, 2>``,
/// ``xt::rtensor<double, 2>``, etc.).
///
- /// XB2: Any 2-dimensional container class storing boolean elements
- /// (e.g. ``xt::xtensor<bool, 4>``, ``xt::pytensor<bool, 4>``,
- /// ``xt::rtensor<bool, 4>``, etc.).
+ /// XB3: Any 3-dimensional container class storing boolean elements
+ /// (e.g. ``xt::xtensor<bool, 3>``, ``xt::pytensor<bool, 3>``,
+ /// ``xt::rtensor<bool, 3>``, etc.).
///
/// :Parameters:
///
@@ -76,13 +76,11 @@ namespace evalhyd
/// epsilon, as recommended by `Pushpalatha et al. (2012)
/// <https://doi.org/10.1016/j.jhydrol.2011.11.055>`_.
///
- /// t_msk: ``XB2``, optional
+ /// t_msk: ``XB3``, optional
/// Mask used to temporally subset of the whole streamflow time series
/// (where True/False is used for the time steps to include/discard in
- /// the subset). If provided, masks must feature the same number of
- /// dimensions as observations and predictions, and it must
- /// broadcastable with both of them.
- /// shape: (subsets, time)
+ /// the subset).
+ /// shape: (series, subsets, time)
///
/// .. seealso:: :doc:`../../functionalities/temporal-masking`
///
@@ -153,12 +151,12 @@ namespace evalhyd
///
/// .. code-block:: c++
///
- /// xt::xtensor<double, 2> msk = {{ 1, 1, 0, 1, 0 }};
+ /// xt::xtensor<double, 3> msk = {{{ 1, 1, 0, 1, 0 }}};
///
/// evalhyd::evald(obs, prd, {"NSE"}, "none", 1, -9, msk);
///
/// \endrst
- template <class XD2, class XB2 = xt::xtensor<bool, 2>>
+ template <class XD2, class XB3 = xt::xtensor<bool, 3>>
std::vector<xt::xarray<double>> evald(
const xt::xexpression<XD2>& q_obs,
const xt::xexpression<XD2>& q_prd,
@@ -169,7 +167,7 @@ namespace evalhyd
xtl::missing<double>(),
xtl::xoptional<double, bool> epsilon =
xtl::missing<double>(),
- const xt::xexpression<XB2>& t_msk = XB2({}),
+ const xt::xexpression<XB3>& t_msk = XB3({}),
const xt::xtensor<std::array<char, 32>, 1>& m_cdt = {},
xtl::xoptional<const std::unordered_map<std::string, int>, bool> bootstrap =
xtl::missing<const std::unordered_map<std::string, int>>(),
@@ -183,10 +181,10 @@ namespace evalhyd
"observations and/or predictions are not two-dimensional"
);
}
- if (xt::get_rank<XB2>::value != 2)
+ if (xt::get_rank<XB3>::value != 3)
{
throw std::runtime_error(
- "temporal masks are not two-dimensional"
+ "temporal masks are not three-dimensional"
);
}
@@ -194,7 +192,7 @@ namespace evalhyd
const XD2& q_obs_ = q_obs.derived_cast();
const XD2& q_prd_ = q_prd.derived_cast();
- const XB2& t_msk_ = t_msk.derived_cast();
+ const XB3& t_msk_ = t_msk.derived_cast();
// check that the metrics to be computed are valid
utils::check_metrics(
@@ -219,7 +217,7 @@ namespace evalhyd
}
if (t_msk_.size() > 0)
{
- if (q_obs_.shape(1) != t_msk_.shape(1))
+ if (q_obs_.shape(1) != t_msk_.shape(2))
{
throw std::runtime_error(
"observations and masks feature different "
@@ -245,44 +243,51 @@ namespace evalhyd
"observations contain more than one time series"
);
}
+ if (t_msk_.size() > 0)
+ {
+ if (q_prd_.shape(0) != t_msk_.shape(0))
+ {
+ throw std::runtime_error(
+ "predictions and masks feature different "
+ "number of series"
+ );
+ }
+ }
// retrieve dimensions
- std::size_t n_tim = q_obs_.shape(1);
- std::size_t n_msk = t_msk_.size() > 0 ? t_msk_.shape(0) :
- (m_cdt.size() > 0 ? m_cdt.shape(0) : 1);
+ std::size_t n_tim = q_prd_.shape(1);
- // initialise a mask if none provided
- // (corresponding to no temporal subset)
+ // generate masks from conditions if provided
auto gen_msk = [&]()
{
- // if t_msk provided, it takes priority
- if (t_msk_.size() > 0)
+ if ((t_msk_.size() < 1) && (m_cdt.size() > 0))
{
- return t_msk_;
- }
- // else if m_cdt provided, use them to generate t_msk
- else if (m_cdt.size() > 0)
- {
- XB2 c_msk = xt::zeros<bool>({n_msk, n_tim});
+ std::size_t n_srs = q_prd_.shape(0);
+ std::size_t n_msk = m_cdt.shape(0);
- for (std::size_t m = 0; m < n_msk; m++)
+ XB3 c_msk = xt::zeros<bool>({n_srs, n_msk, n_tim});
+
+ for (std::size_t s = 0; s < n_srs; s++)
{
- xt::view(c_msk, m) =
- masks::generate_mask_from_conditions(
- m_cdt[0], xt::view(q_obs_, 0), q_prd_
- );
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ xt::view(c_msk, s, m) =
+ masks::generate_mask_from_conditions(
+ xt::view(m_cdt, m),
+ xt::view(q_obs_, 0),
+ xt::view(q_prd_, s, xt::newaxis())
+ );
+ }
}
return c_msk;
}
- // if neither t_msk nor m_cdt provided, generate dummy mask
else
{
- return XB2({xt::ones<bool>({std::size_t{1}, n_tim})});
+ return XB3({});
}
};
-
- auto msk = gen_msk();
+ const XB3 c_msk = gen_msk();
// apply streamflow transformation if requested
auto q_transform = [&](const XD2& q)
@@ -392,7 +397,11 @@ namespace evalhyd
}
// instantiate determinist evaluator
- determinist::Evaluator<XD2, XB2> evaluator(obs, prd, msk, exp);
+ determinist::Evaluator<XD2, XB3> evaluator(
+ obs, prd,
+ t_msk_.size() > 0 ? t_msk_: (m_cdt.size() > 0 ? c_msk : t_msk_),
+ exp
+ );
// retrieve or compute requested metrics
std::vector<xt::xarray<double>> r;
diff --git a/tests/test_determinist.cpp b/tests/test_determinist.cpp
index decdf97cef9b02b689d58be8f3ac2adc3d38dfe9..b7ee38e3364101e29d78b27393ac4632089ce7f4 100644
--- a/tests/test_determinist.cpp
+++ b/tests/test_determinist.cpp
@@ -152,9 +152,11 @@ TEST(DeterministTests, TestMasks)
std::tie(observed, predicted) = load_data_d();
// generate temporal subset by dropping 20 first time steps
- xt::xtensor<bool, 2> masks =
- xt::ones<bool>({std::size_t {1}, std::size_t {observed.size()}});
- xt::view(masks, 0, xt::range(0, 20)) = 0;
+ xt::xtensor<bool, 3> masks =
+ xt::ones<bool>({std::size_t {predicted.shape(0)},
+ std::size_t {1},
+ std::size_t {observed.size()}});
+ xt::view(masks, xt::all(), 0, xt::range(0, 20)) = 0;
// compute scores using masks to subset whole record
std::vector<xt::xarray<double>> metrics_masked =
@@ -183,7 +185,7 @@ TEST(DeterministTests, TestMaskingConditions)
std::tie(observed, predicted) = load_data_d();
// generate dummy empty masks required to access next optional argument
- xt::xtensor<bool, 2> masks;
+ xt::xtensor<bool, 3> masks;
// conditions on streamflow values _________________________________________
@@ -360,7 +362,7 @@ TEST(DeterministTests, TestBootstrap)
{}, // transform
{}, // exponent
{}, // epsilon
- xt::xtensor<bool, 2>({}), // t_msk
+ xt::xtensor<bool, 3>({}), // t_msk
{}, // m_cdt
bootstrap,
datetimes