add "sites" axis to temporal mask dimensions

different temporal subsets may be required of different sites, e.g.
if based on streamflow statistics, which are intrinsically site-specific

include/evalhyd/evalp.hpp

@@ -36,14 +36,14 @@ namespace evalhyd
     ///       Streamflow exceedance threshold(s).
     ///       shape: (thresholds,)
     ///
-    ///    t_msk: ``xt::xtensor<bool, 2>``, optional
+    ///    t_msk: ``xt::xtensor<bool, 3>``, 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). If not provided, no subset is
     ///       performed and only one set of metrics is returned corresponding
     ///       to the whole time series. If provided, as many sets of metrics are
     ///       returned as they are masks provided.
-    ///       shape: (subsets, time)
+    ///       shape: (sites, subsets, time)
     ///
     /// :Returns:
     ///
@@ -68,7 +68,7 @@ namespace evalhyd
     ///
     ///    .. code-block:: c++
     ///
-    ///       xt::xtensor<bool, 2> msk = {{ false, true, true, false, true }};
+    ///       xt::xtensor<bool, 3> msk = {{{ false, true, true, false, true }}};
     ///
     ///       evalhyd::evalp(obs, prd, {"BS"}, thr, msk);
     ///
@@ -82,7 +82,7 @@ namespace evalhyd
             const xt::xtensor<double, 4>& q_prd,
             const std::vector<std::string>& metrics,
             const xt::xtensor<double, 1>& q_thr = {},
-            const xt::xtensor<bool, 2>& t_msk = {}
+            const xt::xtensor<bool, 3>& t_msk = {}
     )
     {
         // check that the metrics to be computed are valid
@@ -140,15 +140,18 @@ namespace evalhyd
             r.emplace_back(xt::zeros<double>(dim[metric]));
 
         // compute variables one site at a time to minimise memory imprint
-        for (int s = 0; s < n_sit; s++)
-            // compute variables one lead time at a time to minimise memory imprint
-            for (int l = 0; l < n_ltm; l++)
+    for (int s = 0; s < n_sit; s++)
+        // compute variables one lead time at a time to minimise memory imprint
+        for (int l = 0; l < n_ltm; l++)
         {
             // instantiate probabilist evaluator
             const auto q_obs_v = xt::view(q_obs, s, xt::all());
             const auto q_prd_v = xt::view(q_prd, s, l, xt::all(), xt::all());
+            const auto t_msk_v = xt::view(t_msk, s, xt::all(), xt::all());
 
-            eh::probabilist::Evaluator evaluator(q_obs_v, q_prd_v, q_thr, t_msk);
+            eh::probabilist::Evaluator evaluator(
+                    q_obs_v, q_prd_v, q_thr, t_msk_v
+            );
 
             // pre-compute required elt
             for (const auto& element : req_elt)

tests/test_probabilist.cpp

@@ -120,9 +120,10 @@ TEST(ProbabilistTests, TestMasks)
     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;
+    xt::xtensor<double, 3> masks =
+            xt::ones<bool>({std::size_t {1}, std::size_t {1},
+                            std::size_t {observed.size()}});
+    xt::view(masks, 0, 0, xt::range(0, 20)) = 0;
 
     // compute scores using masks to subset whole record
     xt::xtensor<double, 1> thresholds = {690, 534, 445};