move away from fully generic API for evald

The use of `xexpression` and template parameter was showing its limits
as the data type could not be different between the different n-dim
inputs (e.g. masks had to be as double rather than bool). Since there is
no real need to be able to use `evald` on `xarray`, the genericity of
its API is reduced so that it now expects inputs to be tensor only. The
number of dimensions for the tensors becomes a template parameter.

include/evalhyd/evald.hpp

@@ -136,38 +136,34 @@ namespace evalhyd
     ///       evalhyd::evald(obs, prd, {"NSE"}, "none", 1, -9, msk);
     ///
     /// \endrst
-    template <class C>
+    template <std::size_t N>
     std::vector<xt::xarray<double>> evald(
-            const xt::xexpression<C>& q_obs,
-            const xt::xexpression<C>& q_prd,
+            const xt::xtensor<double, N>& q_obs,
+            const xt::xtensor<double, N>& q_prd,
             const std::vector<std::string>& metrics,
             const std::string& transform = "none",
             const double exponent = 1,
             double epsilon = -9,
-            const xt::xexpression<C>& t_msk = C()
+            const xt::xtensor<bool, N>& t_msk = {}
     )
     {
-        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)
-        C msk;
-        if (t_msk_.size() < 1)
-            msk = xt::ones<bool>(q_obs_.shape());
+        xt::xtensor<bool, N> msk;
+        if (t_msk.size() < 1)
+            msk = xt::ones<bool>(q_obs.shape());
         else
             msk = std::move(t_msk.derived_cast());
 
         // check that observations, predictions, and masks dimensions are compatible
-        if (q_obs_.dimension() != q_prd_.dimension())
+        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())
+        if (q_obs.dimension() != msk.dimension())
         {
             throw std::runtime_error(
                     "observations and masks feature different numbers "
@@ -175,65 +171,67 @@ namespace evalhyd
             );
         }
 
-        auto obs_shp = q_obs_.shape();
-        auto prd_shp = q_prd_.shape();
+        auto obs_shp = q_obs.shape();
+        auto prd_shp = q_prd.shape();
         auto msk_shp = msk.shape();
 
+        // check observations, predictions, and masks are broadcastable
         if (obs_shp != prd_shp)
-        {
-            // check observations, predictions, and masks are broadcastable
-            for (int a = 0; a < obs_shp.size(); a++)
-            {
+            for (int a = 0; a < q_obs.dimension(); a++)
                 if (obs_shp[a] != prd_shp[a])
                     if ((obs_shp[a] != 1) and (prd_shp[a] != 1))
                         throw std::runtime_error(
                                 "observations and predictions are not "
                                 "broadcastable"
                         );
+
+        if (obs_shp != msk_shp)
+            for (int a = 0; a < q_obs.dimension(); a++)
                 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 != msk_shp)
+            for (int a = 0; a < q_prd.dimension(); a++)
                 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"
                         );
-            }
-        }
 
         // apply streamflow transformation if requested
-        C obs;
-        C prd;
+        xt::xtensor<double, N> obs;
+        xt::xtensor<double, N> prd;
 
         if ( transform == "none" || (transform == "pow" && exponent == 1))
         {
-            obs = std::move(q_obs_);
-            prd = std::move(q_prd_);
+            obs = std::move(q_obs);
+            prd = std::move(q_prd);
         }
         else if ( transform == "sqrt" )
         {
-            obs = xt::sqrt(q_obs_);
-            prd = xt::sqrt(q_prd_);
+            obs = xt::sqrt(q_obs);
+            prd = xt::sqrt(q_prd);
         }
         else if ( transform == "inv" )
         {
             if ( epsilon == -9 )
                 // determine an epsilon value to avoid zero divide
-                epsilon = xt::mean(q_obs_)() * 0.01;
+                epsilon = xt::mean(q_obs)() * 0.01;
 
-            obs = 1. / (q_obs_ + epsilon);
-            prd = 1. / (q_prd_ + epsilon);
+            obs = 1. / (q_obs + epsilon);
+            prd = 1. / (q_prd + epsilon);
         }
         else if ( transform == "log" )
         {
             if ( epsilon == -9 )
                 // determine an epsilon value to avoid log zero
-                epsilon = xt::mean(q_obs_)() * 0.01;
+                epsilon = xt::mean(q_obs)() * 0.01;
 
-            obs = xt::log(q_obs_ + epsilon);
-            prd = xt::log(q_prd_ + epsilon);
+            obs = xt::log(q_obs + epsilon);
+            prd = xt::log(q_prd + epsilon);
         }
         else if ( transform == "pow" )
         {
@@ -241,15 +239,15 @@ namespace evalhyd
             {
                 if ( epsilon == -9 )
                     // determine an epsilon value to avoid zero divide
-                    epsilon = xt::mean(q_obs_)() * 0.01;
+                    epsilon = xt::mean(q_obs)() * 0.01;
 
-                obs = xt::pow(q_obs_ + epsilon, exponent);
-                prd = xt::pow(q_prd_ + epsilon, exponent);
+                obs = xt::pow(q_obs + epsilon, exponent);
+                prd = xt::pow(q_prd + epsilon, exponent);
             }
             else
             {
-                obs = xt::pow(q_obs_, exponent);
-                prd = xt::pow(q_prd_, exponent);
+                obs = xt::pow(q_obs, exponent);
+                prd = xt::pow(q_prd, exponent);
             }
         }
         else
@@ -266,7 +264,7 @@ namespace evalhyd
         );
 
         // instantiate determinist evaluator
-        eh::determinist::Evaluator<C> evaluator(obs, prd, msk);
+        eh::determinist::Evaluator<N> evaluator(obs, prd, msk);
 
         // declare maps for memoisation purposes
         std::unordered_map<std::string, std::vector<std::string>> elt;

src/determinist/evaluator.hpp

@@ -23,7 +23,7 @@ namespace evalhyd
 
     namespace determinist
     {
-        template <class C>
+        template <std::size_t N>
         class Evaluator
         {
         // define type for function to be used to "mask" (i.e. yield NaNs) the
@@ -34,63 +34,53 @@ namespace evalhyd
                 xt::detail::bitwise_or,
                 xt::xfunction<
                     xt::math::isnan_fun,
-                    const C&
+                    const xt::xtensor<double, N>&
                 >,
                 xt::xfunction<
                     xt::detail::logical_not,
-                    const C&
+                    const xt::xtensor<bool, N>&
                 >
             >,
             xt::xscalar<float>,
-            const C&
+            const xt::xtensor<double, N>&
         > nan_function;
 
         private:
             // members for input data
-            const C& _q_obs;
-            const C& _q_prd;
+            const xt::xtensor<double, N>& _q_obs;
+            const xt::xtensor<double, N>& _q_prd;
             nan_function q_obs;
             nan_function q_prd;
 
             // members for computational elements
-            C mean_obs;
-            C mean_prd;
-            C quad_err;
-            C quad_obs;
-            C quad_prd;
-            C r_pearson;
-            C bias;
+            xt::xtensor<double, N> mean_obs;
+            xt::xtensor<double, N> mean_prd;
+            xt::xtensor<double, N> quad_err;
+            xt::xtensor<double, N> quad_obs;
+            xt::xtensor<double, N> quad_prd;
+            xt::xtensor<double, N> r_pearson;
+            xt::xtensor<double, N> bias;
 
         public:
             // constructor method
-            Evaluator(const xt::xexpression<C>& obs,
-                      const xt::xexpression<C>& prd,
-                      const xt::xexpression<C>& msk) :
-                    _q_obs{obs.derived_cast()},
-                    _q_prd{prd.derived_cast()},
+            Evaluator(const xt::xtensor<double, N>& obs,
+                      const xt::xtensor<double, N>& prd,
+                      const xt::xtensor<bool, N>& msk) :
+                    _q_obs{obs},
+                    _q_prd{prd},
                     // "mask" (i.e. use NaN) observations where predictions are
                     // missing (i.e. NaN)
-                    q_obs{
-                        xt::where(
-                            xt::isnan(prd.derived_cast()) | !msk.derived_cast(),
-                            NAN, obs.derived_cast()
-                        )
-                    },
+                    q_obs{xt::where(xt::isnan(prd) | !msk, NAN, obs)},
                     // "mask" (i.e. use NaN) predictions where observations are
                     // missing (i.e. NaN)
-                    q_prd{
-                        xt::where(
-                            xt::isnan(obs.derived_cast()) | !msk.derived_cast(),
-                            NAN, prd.derived_cast()
-                        )
-                    }
+                    q_prd{xt::where(xt::isnan(obs) | !msk, NAN, prd)}
             {};
 
             // members for evaluation metrics
-            C RMSE;
-            C NSE;
-            C KGE;
-            C KGEPRIME;
+            xt::xtensor<double, N> RMSE;
+            xt::xtensor<double, N> NSE;
+            xt::xtensor<double, N> KGE;
+            xt::xtensor<double, N> KGEPRIME;
 
             // methods to compute elements
             void calc_mean_obs();
@@ -116,8 +106,8 @@ namespace evalhyd
         // \assign mean_obs:
         //     Mean observed streamflow.
         //     shape: ({... ,} 1)
-        template <class C>
-        void Evaluator<C>::calc_mean_obs()
+        template <std::size_t N>
+        void Evaluator<N>::calc_mean_obs()
         {
             mean_obs = xt::nanmean(q_obs, -1, xt::keep_dims);
         }
@@ -130,8 +120,8 @@ namespace evalhyd
         // \assign mean_prd:
         //     Mean predicted streamflow.
         //     shape: ({... ,} 1)
-        template <class C>
-        void Evaluator<C>::calc_mean_prd()
+        template <std::size_t N>
+        void Evaluator<N>::calc_mean_prd()
         {
             mean_prd = xt::nanmean(q_prd, -1, xt::keep_dims);
         }
@@ -147,8 +137,8 @@ namespace evalhyd
         // \assign quad_err:
         //     Quadratic errors between observations and predictions.
         //     shape: ({... ,} time)
-        template <class C>
-        void Evaluator<C>::calc_quad_err()
+        template <std::size_t N>
+        void Evaluator<N>::calc_quad_err()
         {
             quad_err = xt::square(q_obs - q_prd);
         }
@@ -164,8 +154,8 @@ namespace evalhyd
         // \assign quad_obs:
         //     Quadratic errors between observations and mean observation.
         //     shape: ({... ,} time)
-        template <class C>
-        void Evaluator<C>::calc_quad_obs()
+        template <std::size_t N>
+        void Evaluator<N>::calc_quad_obs()
         {
             quad_obs = xt::square(q_obs - mean_obs);
         }
@@ -181,8 +171,8 @@ namespace evalhyd
         // \assign quad_prd:
         //     Quadratic errors between predictions and mean prediction.
         //     shape: ({... ,} time)
-        template <class C>
-        void Evaluator<C>::calc_quad_prd()
+        template <std::size_t N>
+        void Evaluator<N>::calc_quad_prd()
         {
             quad_prd = xt::square(q_prd - mean_prd);
         }
@@ -210,8 +200,8 @@ namespace evalhyd
         // \assign r_pearson:
         //     Pearson correlation coefficients.
         //     shape: ({... ,} time)
-        template <class C>
-        void Evaluator<C>::calc_r_pearson()
+        template <std::size_t N>
+        void Evaluator<N>::calc_r_pearson()
         {
             // calculate error in timing and dynamics $r_{pearson}$
             // (Pearson's correlation coefficient)
@@ -238,8 +228,8 @@ namespace evalhyd
         // \assign bias:
         //     Biases.
         //     shape: ({... ,} 1)
-        template <class C>
-        void Evaluator<C>::calc_bias()
+        template <std::size_t N>
+        void Evaluator<N>::calc_bias()
         {
             // calculate $bias$
             bias = xt::nansum(q_prd, -1, xt::keep_dims)
@@ -254,8 +244,8 @@ namespace evalhyd
         // \assign RMSE:
         //     Root-mean-square errors.
         //     shape: ({... ,} 1)
-        template <class C>
-        void Evaluator<C>::calc_RMSE()
+        template <std::size_t N>
+        void Evaluator<N>::calc_RMSE()
         {
            // compute RMSE
            RMSE = xt::sqrt(xt::nanmean(quad_err, -1, xt::keep_dims));
@@ -272,12 +262,12 @@ namespace evalhyd
         // \assign NSE:
         //     Nash-Sutcliffe efficiencies.
         //     shape: ({... ,} 1)
-        template <class C>
-        void Evaluator<C>::calc_NSE()
+        template <std::size_t N>
+        void Evaluator<N>::calc_NSE()
         {
             // compute squared errors operands
-            C f_num = xt::nansum(quad_err, -1, xt::keep_dims);
-            C f_den = xt::nansum(quad_obs, -1, xt::keep_dims);
+            xt::xtensor<double, N> f_num = xt::nansum(quad_err, -1, xt::keep_dims);
+            xt::xtensor<double, N> f_den = xt::nansum(quad_obs, -1, xt::keep_dims);
 
             // compute NSE
             NSE = 1 - (f_num / f_den);
@@ -306,8 +296,8 @@ namespace evalhyd
         // \assign KGE:
         //     Kling-Gupta efficiencies.
         //     shape: ({... ,} 1)
-        template <class C>
-        void Evaluator<C>::calc_KGE()
+        template <std::size_t N>
+        void Evaluator<N>::calc_KGE()
         {
             // calculate error in spread of flow $alpha$
             auto alpha = detail::std_dev(q_prd, mean_prd)
@@ -344,8 +334,8 @@ namespace evalhyd
         // \assign KGEPRIME:
         //     Modified Kling-Gupta efficiencies.
         //     shape: ({... ,} 1)
-        template <class C>
-        void Evaluator<C>::calc_KGEPRIME()
+        template <std::size_t N>
+        void Evaluator<N>::calc_KGEPRIME()
         {
             // calculate error in spread of flow $gamma$
             auto gamma = (detail::std_dev(q_prd, mean_prd) / mean_prd)

tests/test_determinist.cpp

@@ -26,7 +26,7 @@ TEST(DeterministTests, TestMetrics2D)
 
     // compute scores (with 2D tensors)
     std::vector<xt::xarray<double>> metrics =
-            evalhyd::evald<xt::xtensor<double, 2>>(
+            evalhyd::evald<2>(
                     observed, predicted, {"RMSE", "NSE", "KGE", "KGEPRIME"}
             );
 
@@ -80,7 +80,7 @@ TEST(DeterministTests, TestMetrics1D)
 
     // compute scores (with 1D tensors)
     std::vector<xt::xarray<double>> metrics =
-            evalhyd::evald<xt::xtensor<double, 1>>(
+            evalhyd::evald<1>(
                     observed, predicted, {"RMSE", "NSE", "KGE", "KGEPRIME"}
             );
 
@@ -114,9 +114,7 @@ TEST(DeterministTests, TestTransform)
 
     // compute and check results on square-rooted streamflow series
     std::vector<xt::xarray<double>> metrics =
-            evalhyd::evald<xt::xtensor<double, 2>>(
-                    observed, predicted, {"NSE"}, "sqrt"
-            );
+            evalhyd::evald<2>(observed, predicted, {"NSE"}, "sqrt");
 
     xt::xtensor<double, 2> nse_sqrt =
             {{0.882817},
@@ -127,10 +125,7 @@ TEST(DeterministTests, TestTransform)
     EXPECT_TRUE(xt::allclose(metrics[0], nse_sqrt));
 
     // compute and check results on inverted streamflow series
-    metrics =
-            evalhyd::evald<xt::xtensor<double, 2>>(
-                    observed, predicted, {"NSE"}, "inv"
-            );
+    metrics = evalhyd::evald<2>(observed, predicted, {"NSE"}, "inv");
 
     xt::xtensor<double, 2> nse_inv =
             {{0.737323},
@@ -141,10 +136,7 @@ TEST(DeterministTests, TestTransform)
     EXPECT_TRUE(xt::allclose(metrics[0], nse_inv));
 
     // compute and check results on square-rooted streamflow series
-    metrics =
-            evalhyd::evald<xt::xtensor<double, 2>>(
-                    observed, predicted, {"NSE"}, "log"
-            );
+    metrics = evalhyd::evald<2>(observed, predicted, {"NSE"}, "log");
 
     xt::xtensor<double, 2> nse_log =
             {{0.893344},
@@ -155,10 +147,7 @@ TEST(DeterministTests, TestTransform)
     EXPECT_TRUE(xt::allclose(metrics[0], nse_log));
 
     // compute and check results on power-transformed streamflow series
-    metrics =
-    evalhyd::evald<xt::xtensor<double, 2>>(
-            observed, predicted, {"NSE"}, "pow", 0.2
-    );
+    metrics = evalhyd::evald<2>(observed, predicted, {"NSE"}, "pow", 0.2);
 
     xt::xtensor<double, 2> nse_pow =
             {{0.899207},
@@ -192,14 +181,14 @@ TEST(DeterministTests, TestMasks)
             {"RMSE", "NSE", "KGE", "KGEPRIME"};
 
     std::vector<xt::xarray<double>> metrics_masked =
-            evalhyd::evald(observed, predicted, metrics, "none", 1, -9, masks);
+            evalhyd::evald<2>(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);
+            evalhyd::evald<2>(obs, prd, metrics);
 
     // check results are identical
     for (int m = 0; m < metrics.size(); m++)
@@ -237,9 +226,7 @@ TEST(DeterministTests, TestMissingData)
 
     // compute metrics with observations containing NaN values
     std::vector<xt::xarray<double>> metrics_nan =
-            evalhyd::evald<xt::xtensor<double, 2>>(
-                    observed, predicted, metrics
-            );
+            evalhyd::evald<2>(observed, predicted, metrics);
 
     for (int m = 0; m < metrics.size(); m++) {
         for (int p = 0; p < predicted.shape(0); p++) {
@@ -251,9 +238,7 @@ TEST(DeterministTests, TestMissingData)
                     xt::view(predicted, p, xt::range(20+(3*(p+1)), _));
 
             std::vector<xt::xarray<double>> metrics_sbs =
-                    evalhyd::evald<xt::xtensor<double, 1>>(
-                            obs, prd, {metrics[m]}
-                    );
+                    evalhyd::evald<1>(obs, prd, {metrics[m]});
 
             // compare to check results are the same
             EXPECT_TRUE(