diff --git a/include/evalhyd/deterministic.hpp b/include/evalhyd/deterministic.hpp
index 4691796bae9fa956f876114edd2e65a0de1755a0..a7f62771908f413877afa5c45fbf4ad99d567977 100644
--- a/include/evalhyd/deterministic.hpp
+++ b/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);
diff --git a/include/evalhyd/probabilistic.hpp b/include/evalhyd/probabilistic.hpp
index 350c82c05615dfd9662b29b9225f0eefd291c726..881bfeb10517290257e79f29d6824f971837b07c 100644
--- a/include/evalhyd/probabilistic.hpp
+++ b/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);
}
}
diff --git a/include/evalhyd/utils.hpp b/include/evalhyd/utils.hpp
index aad268c2f05d13cd018379cf3523a6d8cf158efc..3a83c6550218cc7d86ad49d5528a09b2bfadd095 100644
--- a/include/evalhyd/utils.hpp
+++ b/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());
}
}
}