implement Brier Score and Brier Skill Score

include/evalhyd/probabilistic.hpp

+#ifndef EVALHYD_PROBABILISTIC_HPP
+#define EVALHYD_PROBABILISTIC_HPP
+
+#include <xtensor/xexpression.hpp>
+#include <xtensor/xmath.hpp>
+#include <xtensor/xindex_view.hpp>
+#include <xtensor/xio.hpp>
+
+#include "utils.hpp"
+
+namespace eh = evalhyd;
+
+namespace evalhyd
+{
+    /// Compute the Brier Score (BS).
+    ///
+    /// \param [in] q_obs:
+    ///     2D array of streamflow observations.
+    /// \param [in] q_frc:
+    ///     2D array of streamflow forecasts.
+    /// \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.
+    /// \return
+    ///     1D array of Brier Score for each threshold.
+    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
+    )
+    {
+        // 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 observed realisation of threshold(s) exceedance
+        xt::xtensor<double, 2> p_obs =
+                xt::squeeze(eh::utils::is_above_threshold(q_obs, q_thr));
+
+        // determine forecast probability of threshold(s) exceedance
+        // > determine if members have exceeded threshold(s)
+        xt::xtensor<double, 3> e_frc =
+                eh::utils::is_above_threshold(q_frc, q_thr);
+
+        // > calculate how many members have exceeded threshold(s)
+        int axis_mbr = (axis == 0) ? 1 : 0;
+        xt::xtensor<int, 2> n_frc = xt::sum(e_frc, axis_mbr);
+
+        // > determine correspondence between number of exceeding members
+        //   and forecast probabilities of exceedance assuming members are
+        //   equiprobable
+        std::size_t n_mbr = q_frc.shape(axis_mbr);
+        xt::xtensor<double, 1> p_mbr =
+                xt::arange<double>(double (n_mbr + 1)) / (n_mbr + 1);
+
+        // > replace number exceeding members with corresponding probability
+        xt::xtensor<double, 2> p_frc = xt::zeros<double>(n_frc.shape());
+        for (int i = 0; i < q_thr.size(); i++)
+        {
+            xt::col(p_frc, i) = xt::index_view(p_mbr, xt::col(n_frc, i));
+        }
+
+        // return computed BS
+        return xt::mean(xt::square(p_obs - p_frc), std::vector<int> { axis });
+    }
+
+    /// Compute the Brier Skill Score (BSS).
+    ///
+    /// \param [in] q_obs:
+    ///     2D array of streamflow observations.
+    /// \param [in] q_frc:
+    ///     2D array of streamflow forecasts.
+    /// \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.
+    /// \return
+    ///     1D array of Brier Skill Score for each threshold.
+    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
+    )
+    {
+        // determine observed realisation of threshold(s) exceedance
+        xt::xtensor<double, 2> p_obs =
+                xt::squeeze(eh::utils::is_above_threshold(q_obs, q_thr));
+
+        // calculate mean observed realisation of threshold(s) exceedance
+        xt::xtensor<double, 1> p_obs_mean =
+                xt::mean(p_obs, std::vector<int> { axis });
+
+        // calculate reference Brier Score
+        xt::xtensor<double, 1> bs_ref = xt::mean(
+                xt::square(p_obs - p_obs_mean), std::vector<int> { axis }
+        );
+
+        // return computed BSS
+        return 1 - (bs(q_obs, q_frc, q_thr, axis) / bs_ref);
+    }
+
+}
+
+#endif //EVALHYD_PROBABILISTIC_HPP

include/evalhyd/utils.hpp

+#ifndef EVALHYD_UTILS_HPP
+#define EVALHYD_UTILS_HPP
+
+#include <xtensor/xtensor.hpp>
+#include <xtensor/xbroadcast.hpp>
+#include <xtensor/xadapt.hpp>
+#include <xtensor/xview.hpp>
+#include <xtensor/xmath.hpp>
+
+namespace evalhyd
+{
+    namespace utils
+    {
+        /// Determine whether flows are strictly greater than given threshold(s)
+        ///
+        /// Streamflow data is always expected as a 2D array (even if only for
+        /// one time series). Threshold(s) given is (are) always expected as a
+        /// 1D array (even if only one threshold is given). A 3D array is
+        /// returned whose first two dimensions are of the same sizes as the
+        /// streamflow data and whose third dimension is of size equal to the
+        /// number of thresholds given. The returned array contains ones where
+        /// the threshold is strictly exceeded, and zeros otherwise.
+        ///
+        /// \param [in] q:
+        ///     2D array of streamflow data
+        /// \param [in] thr:
+        ///     1D array of streamflow threshold(s)
+        /// \return
+        ///     3D array of ones and zeros
+        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;
+        }
+
+        /// Determine whether flows are strictly lower than given threshold(s)
+        ///
+        /// Streamflow data is always expected as a 2D array (even if only for
+        /// one time series). Threshold(s) given is (are) always expected as a
+        /// 1D array (even if only one threshold is given). A 3D array is
+        /// returned whose first two dimensions are of the same sizes as the
+        /// streamflow data and whose third dimension is of size equal to the
+        /// number of thresholds given. The returned array contains ones where
+        /// the threshold is strictly not exceeded, and zeros otherwise.
+        ///
+        /// \param [in] q:
+        ///     2D array of streamflow data
+        /// \param [in] thr:
+        ///     1D array of streamflow threshold(s)
+        /// \return
+        ///     3D array of ones and zeros
+        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;
+        }
+    }
+}
+
+#endif //EVALHYD_UTILS_HPP