drop support for customisable array orientation

Since arrays are row-major, and that the most costly reduction is along the time dimension, it makes sense to always have it as the last axis, so having a customised array orientation would not make this the case all the time. Dropping such support will also help with the documentation to state which axis corresponds to what.

include/evalhyd/deterministic.hpp

@@ -10,29 +10,25 @@ namespace evalhyd
     /// Compute the Nash-Sutcliffe Efficiency (NSE).
     ///
     /// If multi-dimensional arrays are provided, the arrays of simulations
-    /// and observations must feature the same number of dimensions, and
-    /// they must be broadcastable. The optional parameter *axis* may be used
-    /// to indicate which axis corresponds to the temporal dimension on which
-    /// to perform the reductions. The temporal dimension should not be one
-    /// to be broadcast.
+    /// and observations must feature the same number of dimensions, they must
+    /// be broadcastable, and their temporal dimensions must be along their
+    /// last axis.
     ///
     /// \tparam A:
     ///     The type of data structures (e.g. `xt::xarray`, `xt::xtensor`).
     /// \param [in] sim:
     ///     Array of streamflow simulations.
+    ///     shape: (..., time)
     /// \param [in] obs:
     ///     Array of streamflow observations.
-    /// \param [in] axis (optional):
-    ///     Array axis along which the temporal dimension is for both *sim*
-    ///     arrays are multi-dimensional. If not provided, set to default
-    ///     value 0.
+    ///     shape: (1, time)
     /// \return
     ///     Array of computed efficiencies.
+    ///     shape: (...)
     template <class A>
     inline A nse(
             const xt::xexpression<A>& sim,
-            const xt::xexpression<A>& obs,
-            int axis = 0
+            const xt::xexpression<A>& obs
     )
     {
         const A& q_sim = sim.derived_cast();
@@ -46,15 +42,12 @@ namespace evalhyd
             );
         }
 
-        // turn axis into a vector
-        std::vector<int> axes = { axis };
-
         // compute average observed flow
-        A q_avg = xt::mean(q_obs, axes, xt::keep_dims);
+        A q_avg = xt::mean(q_obs, -1, xt::keep_dims);
 
         // compute squared errors operands
-        A f_num = xt::sum(xt::square(q_obs - q_sim), axes, xt::keep_dims);
-        A f_den = xt::sum(xt::square(q_obs - q_avg), axes, xt::keep_dims);
+        A f_num = xt::sum(xt::square(q_obs - q_sim), -1, xt::keep_dims);
+        A f_den = xt::sum(xt::square(q_obs - q_avg), -1, xt::keep_dims);
 
         // return computed NSE
         return 1 - (f_num / f_den);

include/evalhyd/probabilistic.hpp

@@ -14,19 +14,10 @@ namespace evalhyd
 {
     namespace detail
     {
-        /// check whether the axis is compatible with 2D arrays
-        void check_axis(int axis)
-        {
-            // check that axis is a valid index for a 2D array
-            if ((axis != 0) && (axis != 1))
-            {
-                throw std::out_of_range(
-                        "array axis must be 0 or 1 for 2D arrays"
-                );
-            }
-        }
-
-        /// determine whether exceedance event occurred
+        // determine whether exceedance event has occurred
+        //     q_obs shape: (1, time)
+        //     q_thr shape: (thresholds,)
+        //     returned shape: (thresholds, time)
         xt::xtensor<double, 2> event_observed_outcome(
                 const xt::xtensor<double, 2>& q_obs,
                 const xt::xtensor<double, 1>& q_thr,
@@ -43,11 +34,13 @@ namespace evalhyd
             }
         }
 
-        /// determine forecast probability of event to occur
+        // determine forecast probability of event to occur
+        //     q_obs shape: (members, time)
+        //     q_thr shape: (thresholds,)
+        //     returned shape: (thresholds, members, time)
         xt::xtensor<double, 2> event_probability_forecast(
                 const xt::xtensor<double, 2>& q_frc,
                 const xt::xtensor<double, 1>& q_thr,
-                int axis,
                 bool above
         )
         {
@@ -57,14 +50,13 @@ namespace evalhyd
             xt::xtensor<double, 3> e_frc = f(q_frc, q_thr);
 
             // calculate how many members have exceeded threshold(s)
-            int axis_mbr = (axis == 0) ? 1 : 0;
-            xt::xtensor<double, 2> n_frc = xt::sum(e_frc, axis_mbr);
+            xt::xtensor<double, 2> n_frc = xt::sum(e_frc, 1);
 
             // determine correspondence between number of exceeding members
             // and forecast probabilities of exceedance
             // /!\ probability calculation dividing by n (the number of
             //     members), not n+1 (the number of ranks) like other metrics
-            std::size_t n_mbr = q_frc.shape(axis_mbr);
+            std::size_t n_mbr = q_frc.shape(0);
             xt::xtensor<double, 2> p_frc = n_frc / n_mbr;
 
             return p_frc;
@@ -75,14 +67,13 @@ namespace evalhyd
     ///
     /// \param [in] q_obs:
     ///     2D array of streamflow observations.
+    ///     shape: (1, time)
     /// \param [in] q_frc:
     ///     2D array of streamflow forecasts.
+    ///     shape: (members, time)
     /// \param [in] q_thr:
     ///     1D array of streamflow exceedance threshold(s).
-    /// \param [in] axis (optional):
-    ///     Array axis along which the temporal dimension is for both
-    ///     *frc* and *obs* (either 0 or 1). If not provided, set to
-    ///     default value 0.
+    ///     shape: (thresholds,)
     /// \param [in] above (optional):
     ///     Boolean to indicate whether exceedance corresponds to being greater
     ///     (above) the thresholds or being lower than the thresholds (below).
@@ -90,39 +81,36 @@ namespace evalhyd
     ///     default value True.
     /// \return
     ///     1D array of Brier Score for each threshold.
+    ///     shape: (thresholds,)
     xt::xtensor<double, 1> bs(
             const xt::xtensor<double, 2>& q_obs,
             const xt::xtensor<double, 2>& q_frc,
             const xt::xtensor<double, 1>& q_thr,
-            int axis = 0,
             bool above = true
     )
     {
-        detail::check_axis(axis);
-
         // get observed event and forecast probability of event
         xt::xtensor<double, 2> p_obs =
                 detail::event_observed_outcome(q_obs, q_thr, above);
 
         xt::xtensor<double, 2> p_frc =
-                detail::event_probability_forecast(q_frc, q_thr, axis, above);
+                detail::event_probability_forecast(q_frc, q_thr, above);
 
         // return computed BS
-        return xt::mean(xt::square(p_obs - p_frc), std::vector<int> { axis });
+        return xt::mean(xt::square(p_obs - p_frc), -1);
     }
 
     /// Compute the Brier Skill Score (BSS).
     ///
     /// \param [in] q_obs:
     ///     2D array of streamflow observations.
+    ///     shape: (1, time)
     /// \param [in] q_frc:
     ///     2D array of streamflow forecasts.
+    ///     shape: (members, time)
     /// \param [in] q_thr:
     ///     1D array of streamflow exceedance threshold(s).
-    /// \param [in] axis (optional):
-    ///     Array axis along which the temporal dimension is for both
-    ///     *frc* and *obs* (either 0 or 1). If not provided, set to
-    ///     default value 0.
+    ///     shape: (thresholds,)
     /// \param [in] above (optional):
     ///     Boolean to indicate whether exceedance corresponds to being greater
     ///     (above) the thresholds or being lower than the thresholds (below).
@@ -130,11 +118,11 @@ namespace evalhyd
     ///     default value True.
     /// \return
     ///     1D array of Brier Skill Score for each threshold.
+    ///     shape: (thresholds,)
     xt::xtensor<double, 1> bss(
             const xt::xtensor<double, 2>& q_obs,
             const xt::xtensor<double, 2>& q_frc,
             const xt::xtensor<double, 1>& q_thr,
-            int axis = 0,
             bool above = true
     )
     {
@@ -143,16 +131,15 @@ namespace evalhyd
                 detail::event_observed_outcome(q_obs, q_thr, above);
 
         // calculate mean observed realisation of threshold(s) exceedance
-        xt::xtensor<double, 1> p_obs_mean =
-                xt::mean(p_obs, std::vector<int> { axis });
+        xt::xtensor<double, 2> p_obs_mean = xt::mean(p_obs, -1, xt::keep_dims);
 
         // calculate reference Brier Score
         xt::xtensor<double, 1> bs_ref = xt::mean(
-                xt::square(p_obs - p_obs_mean), std::vector<int> { axis }
+                xt::square(p_obs - p_obs_mean), -1
         );
 
         // return computed BSS
-        return 1 - (bs(q_obs, q_frc, q_thr, axis) / bs_ref);
+        return 1 - (bs(q_obs, q_frc, q_thr) / bs_ref);
     }
 }
 

include/evalhyd/utils.hpp

@@ -22,25 +22,21 @@ namespace evalhyd
         /// the threshold is exceeded, and zeros otherwise.
         ///
         /// \param [in] q:
-        ///     2D array of streamflow data
+        ///     2D array of streamflow data.
+        ///     shape: (1+, time)
         /// \param [in] thr:
-        ///     1D array of streamflow threshold(s)
+        ///     1D array of streamflow threshold(s).
+        ///     shape: (thresholds,)
         /// \return
-        ///     3D array of ones and zeros
+        ///     3D array of ones and zeros.
+        ///     shape: (thresholds, 1+, time)
         xt::xtensor<double, 3> is_above_threshold(
                 const xt::xtensor<double, 2>& q,
                 const xt::xtensor<double, 1>& thr
         )
         {
-            // add one dimension to the streamflow data
-            auto new_shape = {q.shape(0), q.shape(1), thr.size()};
-            auto q2 = xt::broadcast(
-                    xt::view(q, xt::all(), xt::all(), xt::newaxis()),
-                    new_shape
-            );
-
             // return a boolean-like array
-            return q2 >= thr;
+            return q >= xt::view(thr, xt::all(), xt::newaxis(), xt::newaxis());
         }
 
         /// Determine whether flows are strictly lower than given threshold(s)
@@ -54,25 +50,21 @@ namespace evalhyd
         /// the threshold is strictly not exceeded, and zeros otherwise.
         ///
         /// \param [in] q:
-        ///     2D array of streamflow data
+        ///     2D array of streamflow data.
+        ///     shape: (1+, time)
         /// \param [in] thr:
-        ///     1D array of streamflow threshold(s)
+        ///     1D array of streamflow threshold(s).
+        ///     shape: (thresholds,)
         /// \return
-        ///     3D array of ones and zeros
+        ///     3D array of ones and zeros.
+        ///     shape: (thresholds, 1+, time)
         xt::xtensor<double, 3> is_below_threshold(
                 const xt::xtensor<double, 2>& q,
                 const xt::xtensor<double, 1>& thr
         )
         {
-            // add one dimension to the streamflow data
-            auto new_shape = {q.shape(0), q.shape(1), thr.size()};
-            auto q2 = xt::broadcast(
-                    xt::view(q, xt::all(), xt::all(), xt::newaxis()),
-                    new_shape
-            );
-
             // return a boolean-like array
-            return q2 < thr;
+            return q < xt::view(thr, xt::all(), xt::newaxis(), xt::newaxis());
         }
     }
 }