implement temporal masking functionality for determinist evaluation

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;

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

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 =