diff --git a/include/evalhyd/determinist.hpp b/include/evalhyd/determinist.hpp
index 205c668e569c7aa90aaed63cdc55b2cf7d84988d..7bb8068c212b1c8ca47cfd2544e4ae22a2176493 100644
--- a/include/evalhyd/determinist.hpp
+++ b/include/evalhyd/determinist.hpp
@@ -6,7 +6,7 @@
#include <xtensor/xexpression.hpp>
#include "utils.hpp"
-#include "deterministic/nse.hpp"
+#include "deterministic/evaluator.hpp"
namespace eh = evalhyd;
@@ -36,11 +36,11 @@ namespace evalhyd
const A& obs = q_obs.derived_cast();
const A& sim = q_sim.derived_cast();
+ eh::determinist::Evaluator<A> evaluator(obs, sim);
+
// declare maps for memoisation purposes
std::unordered_map<std::string, std::vector<std::string>> elt;
- std::unordered_map<std::string, A> mem_elt;
std::unordered_map<std::string, std::vector<std::string>> dep;
- std::unordered_map<std::string, A> mem_dep;
// register potentially recurring computation elt across metrics
// TODO
@@ -71,15 +71,11 @@ namespace evalhyd
for ( const auto& metric : metrics )
{
- // if exists, retrieve pre-computed metric
- if ( mem_dep.find(metric) != mem_dep.end() )
- {
- r.emplace_back(mem_dep[metric]);
- }
- // else, compute relevant metric
- else if ( metric == "nse" )
+ if ( metric == "nse" )
{
- r.emplace_back(eh::nse<A>(obs, sim));
+ if ( req_dep.find(metric) == req_dep.end() )
+ evaluator.calc_nse();
+ r.emplace_back(evaluator.nse);
}
}
diff --git a/include/evalhyd/deterministic/evaluator.hpp b/include/evalhyd/deterministic/evaluator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a0ae42dde690421a7991ad0da0ffb09d0c7b7ce0
--- /dev/null
+++ b/include/evalhyd/deterministic/evaluator.hpp
@@ -0,0 +1,72 @@
+#ifndef EVALHYD_DETERMINIST_EVALUATOR_HPP
+#define EVALHYD_DETERMINIST_EVALUATOR_HPP
+
+#include <xtensor/xtensor.hpp>
+
+namespace evalhyd
+{
+ namespace determinist
+ {
+ template <class A>
+ class Evaluator
+ {
+ private:
+ // members for input data
+ const A& q_obs;
+ const A& q_sim;
+ // members for computational elements
+ // TODO
+
+ public:
+ // constructor method
+ Evaluator(const A& obs, const A& sim) : q_obs{obs}, q_sim{sim} {
+ // check that data dimensions are compatible
+ if (q_sim.dimension() != q_obs.dimension())
+ {
+ throw std::runtime_error(
+ "number of dimensions of 'sim' and 'obs' "
+ "must be identical"
+ );
+ }
+ };
+ // members for evaluation metrics
+ A nse;
+ // methods to compute elements
+ // TODO
+ // methods to compute metrics
+ void calc_nse();
+ };
+
+ // 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, they must
+ // be broadcastable, and their temporal dimensions must be along their
+ // last axis.
+ //
+ // \require q_obs:
+ // Array of streamflow observations.
+ // shape: (1, time)
+ // \require q_sim:
+ // Array of streamflow simulations.
+ // shape: (..., time)
+ // \assign nse
+ // Array of computed Nash-Sutcliffe efficiencies.
+ // shape: (...)
+ template <class A>
+ void Evaluator<A>::calc_nse()
+ {
+ // compute average observed flow
+ 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), -1, xt::keep_dims);
+ A f_den = xt::sum(xt::square(q_obs - q_avg), -1, xt::keep_dims);
+
+ // return computed NSE
+ nse = 1 - (f_num / f_den);
+ }
+ }
+}
+
+#endif //EVALHYD_DETERMINIST_EVALUATOR_HPP
diff --git a/include/evalhyd/deterministic/nse.hpp b/include/evalhyd/deterministic/nse.hpp
deleted file mode 100644
index d35d7fe5b4428065ea29c6bfad0bd9d089dde155..0000000000000000000000000000000000000000
--- a/include/evalhyd/deterministic/nse.hpp
+++ /dev/null
@@ -1,54 +0,0 @@
-#ifndef EVALHYD_NSE_HPP
-#define EVALHYD_NSE_HPP
-
-#include <vector>
-#include <xtensor/xexpression.hpp>
-#include <xtensor/xmath.hpp>
-
-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, 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] obs:
- /// Array of streamflow observations.
- /// shape: (1, time)
- /// \param [in] sim:
- /// Array of streamflow simulations.
- /// shape: (..., time)
- /// \return
- /// Array of computed efficiencies.
- /// shape: (...)
- template <class A>
- inline A nse(
- const A& obs,
- const A& sim
- )
- {
- // check that data dimensions are compatible
- if (sim.dimension() != obs.dimension())
- {
- throw std::runtime_error(
- "number of dimensions of 'sim' and 'obs' must be identical"
- );
- }
-
- // compute average observed flow
- A avg = xt::mean(obs, -1, xt::keep_dims);
-
- // compute squared errors operands
- A f_num = xt::sum(xt::square(obs - sim), -1, xt::keep_dims);
- A f_den = xt::sum(xt::square(obs - avg), -1, xt::keep_dims);
-
- // return computed NSE
- return 1 - (f_num / f_den);
- }
-}
-
-#endif //EVALHYD_NSE_HPP
diff --git a/include/evalhyd/probabilist.hpp b/include/evalhyd/probabilist.hpp
index 2743b920b90cec9f2e371e326757d0899a28daec..00a305b69846b75272cef5af9687ff4ad5dc6905 100644
--- a/include/evalhyd/probabilist.hpp
+++ b/include/evalhyd/probabilist.hpp
@@ -1,19 +1,21 @@
#ifndef EVALHYD_PROBABILIST_HPP
#define EVALHYD_PROBABILIST_HPP
+#include <utility>
#include <unordered_map>
#include <unordered_set>
#include <xtensor/xtensor.hpp>
+#include <xtensor/xview.hpp>
#include "utils.hpp"
-#include "probabilistic/elements.hpp"
-#include "probabilistic/brier.hpp"
+#include "probabilistic/evaluator.h"
namespace eh = evalhyd;
namespace evalhyd
{
- namespace probabilist {
+ namespace probabilist
+ {
/// Function allowing the evaluation of streamflow forecasts using a
/// range of relevant metrics.
///
@@ -37,11 +39,11 @@ namespace evalhyd
const xt::xtensor<double, 1>& q_thr
)
{
+ eh::probabilist::Evaluator evaluator(q_obs, q_frc, q_thr);
+
// declare maps for memoisation purposes
std::unordered_map<std::string, std::vector<std::string>> elt;
- std::unordered_map<std::string, xt::xtensor<double, 2>> mem_elt;
std::unordered_map<std::string, std::vector<std::string>> dep;
- std::unordered_map<std::string, xt::xtensor<double, 1>> mem_dep;
// register potentially recurring computation elt across metrics
elt["bs"] = {"o_k", "y_k"};
@@ -62,18 +64,16 @@ namespace evalhyd
for (const auto& element : req_elt)
{
if ( element == "o_k" )
- mem_elt["o_k"] =
- eh::elements::event_observed_outcome(q_obs, q_thr);
+ evaluator.calc_o_k();
else if ( element == "y_k" )
- mem_elt["y_k"] =
- eh::elements::event_probability_forecast(q_frc, q_thr);
+ evaluator.calc_y_k();
}
// pre-compute required dep
for (const auto& dependency : req_dep)
{
if ( dependency == "bs" )
- mem_dep["bs"] = eh::bs(mem_elt);
+ evaluator.calc_bs();
}
// retrieve or compute requested metrics
@@ -81,33 +81,33 @@ namespace evalhyd
for (const auto& metric : metrics)
{
- // if exists, retrieve pre-computed metric
- if (mem_dep.find(metric) != mem_dep.end())
- {
- r.emplace_back(mem_dep[metric]);
- }
- // else, compute relevant metric
- else if ( metric == "bs" )
+ if ( metric == "bs" )
{
- r.emplace_back(eh::bs(mem_elt));
+ if (req_dep.find(metric) == req_dep.end())
+ evaluator.calc_bs();
+ r.emplace_back(evaluator.bs);
}
else if ( metric == "bss" )
{
- r.emplace_back(eh::bss(mem_elt, mem_dep));
+ if (req_dep.find(metric) == req_dep.end())
+ evaluator.calc_bss();
+ r.emplace_back(evaluator.bss);
}
else if ( metric == "bs_crd" )
{
- auto c = eh::bs_crd(mem_elt, q_frc.shape(0));
- r.emplace_back(xt::row(c, 0));
- r.emplace_back(xt::row(c, 1));
- r.emplace_back(xt::row(c, 2));
+ if (req_dep.find(metric) == req_dep.end())
+ evaluator.calc_bs_crd();
+ r.emplace_back(xt::row(evaluator.bs_crd, 0));
+ r.emplace_back(xt::row(evaluator.bs_crd, 1));
+ r.emplace_back(xt::row(evaluator.bs_crd, 2));
}
else if ( metric == "bs_lbd" )
{
- auto c = eh::bs_lbd(mem_elt);
- r.emplace_back(xt::row(c, 0));
- r.emplace_back(xt::row(c, 1));
- r.emplace_back(xt::row(c, 2));
+ if (req_dep.find(metric) == req_dep.end())
+ evaluator.calc_bs_lbd();
+ r.emplace_back(xt::row(evaluator.bs_lbd, 0));
+ r.emplace_back(xt::row(evaluator.bs_lbd, 1));
+ r.emplace_back(xt::row(evaluator.bs_lbd, 2));
}
}
diff --git a/include/evalhyd/probabilistic/brier.hpp b/include/evalhyd/probabilistic/brier.hpp
deleted file mode 100644
index 2164d1bdbb596380092b4510669977136624dd75..0000000000000000000000000000000000000000
--- a/include/evalhyd/probabilistic/brier.hpp
+++ /dev/null
@@ -1,289 +0,0 @@
-#ifndef EVALHYD_BRIER_HPP
-#define EVALHYD_BRIER_HPP
-
-#include <unordered_map>
-#include <xtensor/xtensor.hpp>
-#include <xtensor/xmath.hpp>
-#include <xtensor/xoperation.hpp>
-
-#include "elements.hpp"
-
-namespace eh = evalhyd;
-
-namespace evalhyd
-{
- // NOTE --------------------------------------------------------------------
- // All equations in metrics below are following notations from
- // "Wilks, D. S. (2011). Statistical methods in the atmospheric sciences.
- // Amsterdam; Boston: Elsevier Academic Press. ISBN: 9780123850225".
- // In particular, pp. 302-303, 332-333.
- // -------------------------------------------------------------------------
-
- /// Compute the Brier score (BS).
- ///
- /// \param [in] elt:
- /// Map of pre-computed elements, including required elements:
- /// o_k:
- /// Observed event outcome.
- /// shape: (thresholds, time)
- /// y_k:
- /// Event probability forecast.
- /// shape: (thresholds, time)
- /// \return
- /// 1D array of Brier Score for each threshold.
- /// shape: (thresholds,)
- xt::xtensor<double, 1> bs(
- std::unordered_map<std::string, xt::xtensor<double, 2>>& elt
- )
- {
- // return computed Brier score(s)
- // $BS = \frac{1}{n} \sum_{k=1}^{n} (o_k - y_k)^2$
- return xt::mean(xt::square(elt["o_k"] - elt["y_k"]), -1);
- }
-
- /// Compute the calibration-refinement decomposition of the Brier score
- /// into reliability, resolution, and uncertainty.
- ///
- /// BS = reliability - resolution + uncertainty
- ///
- /// \param [in] elt:
- /// Map of pre-computed elements, including required elements:
- /// o_k:
- /// Observed event outcome.
- /// shape: (thresholds, time)
- /// y_k:
- /// Event probability forecast.
- /// shape: (thresholds, time)
- /// \param [in] n_mbr:
- /// Number of forecast members.
- /// \return
- /// 2D array of Brier score components (reliability, resolution,
- /// uncertainty) for each threshold.
- /// shape: (components, thresholds)
- xt::xtensor<double, 2> bs_crd(
- std::unordered_map<std::string, xt::xtensor<double, 2>>& elt,
- std::size_t n_mbr
- )
- {
- // declare internal variables
- // shape: (bins, thresholds, time)
- xt::xtensor<double, 3> msk_bins;
- // shape: (thresholds,)
- xt::xtensor<double, 1> bar_o;
- // shape: (bins, thresholds)
- xt::xtensor<double, 2> N_i, bar_o_i;
- // shape: (bins,)
- xt::xtensor<double, 1> y_i;
-
- // define some dimensions
- std::size_t n = elt["o_k"].shape(1);
- std::size_t n_thr = elt["o_k"].shape(0);
- std::size_t n_cmp = 3;
-
- // declare and initialise output variable
- // shape: (components, thresholds)
- xt::xtensor<double, 2> bs_components = xt::zeros<double>({n_cmp, n_thr});
-
- // compute range of forecast values $y_i$
- y_i = xt::arange<double>(double (n_mbr + 1)) / n_mbr;
-
- // compute mask to subsample time steps belonging to same forecast bin
- // (where bins are defined as the range of forecast values)
- msk_bins = xt::equal(
- elt["y_k"],
- xt::view(y_i, xt::all(), xt::newaxis(), xt::newaxis())
- );
-
- // compute number of forecasts in each forecast bin $N_i$
- N_i = xt::sum(msk_bins, -1);
-
- // compute subsample relative frequency
- // $\bar{o_i} = \frac{1}{N_i} \sum_{k \in N_i} o_k$
- bar_o_i = xt::where(
- N_i > 0,
- xt::sum(
- xt::where(
- msk_bins,
- xt::view(elt["o_k"], xt::newaxis(), xt::all(), xt::all()),
- 0.
- ),
- -1
- ) / N_i,
- 0.
- );
-
- // compute mean "climatology" relative frequency of the event
- // $\bar{o} = \frac{1}{n} \sum_{k=1}^{n} o_k$
- bar_o = xt::mean(elt["o_k"], -1);
-
- // calculate reliability =
- // $\frac{1}{n} \sum_{i=1}^{I} N_i (y_i - \bar{o_i})^2$
- xt::row(bs_components, 0) =
- xt::sum(
- xt::square(
- xt::view(y_i, xt::all(), xt::newaxis())
- - bar_o_i
- ) * N_i,
- 0
- ) / n;
-
- // calculate resolution =
- // $\frac{1}{n} \sum_{i=1}^{I} N_i (\bar{o_i} - \bar{o})^2$
- xt::row(bs_components, 1) =
- xt::sum(
- xt::square(
- bar_o_i - bar_o
- ) * N_i,
- 0
- ) / n;
-
- // calculate uncertainty = $\bar{o} (1 - \bar{o})$
- xt::row(bs_components, 2) = bar_o * (1 - bar_o);
-
- return bs_components;
- }
-
- /// Compute the likelihood-based decomposition of the Brier score
- /// into type 2 bias, discrimination, and sharpness (a.k.a. refinement).
- ///
- /// BS = type 2 bias - discrimination + sharpness
- ///
- /// \param [in] elt:
- /// Map of pre-computed elements, including required elements:
- /// o_k:
- /// Observed event outcome.
- /// shape: (thresholds, time)
- /// y_k:
- /// Event probability forecast.
- /// shape: (thresholds, time)
- /// \return
- /// 2D array of Brier score components (type 2 bias, discrimination,
- /// sharpness) for each threshold.
- /// shape: (components, thresholds)
- xt::xtensor<double, 2> bs_lbd(
- std::unordered_map<std::string, xt::xtensor<double, 2>>& elt
- )
- {
- // declare internal variables
- // shape: (bins, thresholds, time)
- xt::xtensor<double, 3> msk_bins;
- // shape: (thresholds,)
- xt::xtensor<double, 1> bar_y;
- // shape: (bins, thresholds)
- xt::xtensor<double, 2> M_j, bar_y_j;
- // shape: (bins,)
- xt::xtensor<double, 1> o_j;
-
- // define some dimensions
- std::size_t n = elt["o_k"].shape(1);
- std::size_t n_thr = elt["o_k"].shape(0);
- std::size_t n_cmp = 3;
-
- // declare and initialise output variable
- // shape: (components, thresholds)
- xt::xtensor<double, 2> bs_components = xt::zeros<double>({n_cmp, n_thr});
-
- // set the range of observed values $o_j$
- o_j = {0., 1.};
-
- // compute mask to subsample time steps belonging to same observation bin
- // (where bins are defined as the range of forecast values)
- msk_bins = xt::equal(
- elt["o_k"],
- xt::view(o_j, xt::all(), xt::newaxis(), xt::newaxis())
- );
-
- // compute number of observations in each observation bin $M_j$
- M_j = xt::sum(msk_bins, -1);
-
- // compute subsample relative frequency
- // $\bar{y_j} = \frac{1}{M_j} \sum_{k \in M_j} y_k$
- bar_y_j = xt::where(
- M_j > 0,
- xt::sum(
- xt::where(
- msk_bins,
- xt::view(elt["y_k"], xt::newaxis(), xt::all(), xt::all()),
- 0.
- ),
- -1
- ) / M_j,
- 0.
- );
-
- // compute mean "climatology" forecast probability
- // $\bar{y} = \frac{1}{n} \sum_{k=1}^{n} y_k$
- bar_y = xt::mean(elt["y_k"], -1);
-
- // calculate type 2 bias =
- // $\frac{1}{n} \sum_{j=1}^{J} M_j (o_j - \bar{y_j})^2$
- xt::row(bs_components, 0) =
- xt::sum(
- xt::square(
- xt::view(o_j, xt::all(), xt::newaxis())
- - bar_y_j
- ) * M_j,
- 0
- ) / n;
-
- // calculate discrimination =
- // $\frac{1}{n} \sum_{j=1}^{J} M_j (\bar{y_j} - \bar{y})^2$
- xt::row(bs_components, 1) =
- xt::sum(
- xt::square(
- bar_y_j - bar_y
- ) * M_j,
- 0
- ) / n;
-
- // calculate sharpness =
- // $\frac{1}{n} \sum_{k=1}^{n} (\bar{y_k} - \bar{y})^2$
- xt::row(bs_components, 2) =
- xt::sum(
- xt::square(
- elt["y_k"] -
- xt::view(bar_y, xt::all(), xt::newaxis())
- ),
- -1
- ) / n;
-
- return bs_components;
- }
-
- /// Compute the Brier skill score (BSS).
- ///
- /// \param [in] elt:
- /// Map of pre-computed elements, including required elements:
- /// o_k:
- /// Observed event outcome.
- /// shape: (thresholds, time)
- /// \param [in] dep:
- /// Map of pre-computed metrics, including required metrics:
- /// bs:
- /// Brier score(s).
- /// shape: (thresholds,)
- /// \return
- /// 1D array of Brier skill score for each threshold.
- /// shape: (thresholds,)
- xt::xtensor<double, 1> bss(
- std::unordered_map<std::string, xt::xtensor<double, 2>>& elt,
- std::unordered_map<std::string, xt::xtensor<double, 1>>& dep
- )
- {
- // calculate mean observed event outcome
- // $\bar{o_k} = \frac{1}{n} \sum_{k=1}^{n} o_k$
- xt::xtensor<double, 2> bar_o_k = xt::mean(elt["o_k"], -1, xt::keep_dims);
-
- // calculate reference Brier score(s)
- // $BS_{ref} = \frac{1}{n} \sum_{k=1}^{n} (o_k - \bar{o_k})^2$
- xt::xtensor<double, 1> bs_ref = xt::mean(
- xt::square(elt["o_k"] - bar_o_k), -1
- );
-
- // return computed Brier skill score(s)
- // $BSS = 1 - \frac{BS}{BS_{ref}}
- return 1 - (dep["bs"] / bs_ref);
- }
-}
-
-#endif //EVALHYD_BRIER_HPP
diff --git a/include/evalhyd/probabilistic/elements.hpp b/include/evalhyd/probabilistic/elements.hpp
deleted file mode 100644
index 9253f9018e51244f3319f402ec0294f99b5bb21c..0000000000000000000000000000000000000000
--- a/include/evalhyd/probabilistic/elements.hpp
+++ /dev/null
@@ -1,99 +0,0 @@
-#ifndef EVALHYD_ELEMENTS_HPP
-#define EVALHYD_ELEMENTS_HPP
-
-#include <xtensor/xtensor.hpp>
-#include "xtensor/xview.hpp"
-#include <xtensor/xmath.hpp>
-
-
-namespace eh = evalhyd;
-
-namespace evalhyd
-{
- namespace elements
- {
- namespace detail
- {
- // determine whether flows `q` are greater than given threshold(s) `thr`
- // q shape: (1+, time)
- // thr shape: (thresholds,)
- // returned shape: (thresholds, 1+, time)
- xt::xtensor<double, 3> is_above_threshold(
- const xt::xtensor<double, 2>& q,
- const xt::xtensor<double, 1>& thr
- )
- {
- // return a boolean-like array
- return q >= xt::view(thr, xt::all(), xt::newaxis(), xt::newaxis());
- }
-
- // determine whether flows `q` are strictly lower than given threshold(s) `thr`
- // q shape: (1+, time)
- // thr shape: (thresholds,)
- // returned shape: (thresholds, 1+, time)
- xt::xtensor<double, 3> is_below_threshold(
- const xt::xtensor<double, 2>& q,
- const xt::xtensor<double, 1>& thr
- )
- {
- // return a boolean-like array
- return q < xt::view(thr, xt::all(), xt::newaxis(), xt::newaxis());
- }
- }
-
- /// Determine whether exceedance event has occurred.
- ///
- /// \param [in] q_obs:
- /// 2D array of streamflow observations.
- /// shape: (1, time)
- /// \param [in] q_thr:
- /// 1D array of streamflow exceedance threshold(s).
- /// shape: (thresholds,)
- /// \return
- /// 2D array of event observed outcome
- /// shape: (thresholds, time)
- xt::xtensor<double, 2> event_observed_outcome(
- const xt::xtensor<double, 2>& q_obs,
- const xt::xtensor<double, 1>& q_thr
- )
- {
- // determine observed realisation of threshold(s) exceedance
- return xt::squeeze(detail::is_above_threshold(q_obs, q_thr));
- }
-
- /// Determine forecast probability of event to occur.
- ///
- /// \param [in] q_frc:
- /// 2D array of streamflow forecasts.
- /// shape: (members, time)
- /// \param [in] q_thr:
- /// 1D array of streamflow exceedance threshold(s).
- /// shape: (thresholds,)
- /// \return
- /// 2D array of event probability forecast
- /// shape: (thresholds, time)
- xt::xtensor<double, 2> event_probability_forecast(
- const xt::xtensor<double, 2>& q_frc,
- const xt::xtensor<double, 1>& q_thr
- )
- {
- // determine if members have exceeded threshold(s)
- xt::xtensor<double, 3> e_frc =
- detail::is_above_threshold(q_frc, q_thr);
-
- // calculate how many members have exceeded threshold(s)
- 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 in other metrics
- std::size_t n_mbr = q_frc.shape(0);
- xt::xtensor<double, 2> p_frc = n_frc / n_mbr;
-
- return p_frc;
- }
- }
-}
-
-#endif //EVALHYD_ELEMENTS_HPP
diff --git a/include/evalhyd/probabilistic/evaluator.h b/include/evalhyd/probabilistic/evaluator.h
new file mode 100644
index 0000000000000000000000000000000000000000..f2d8b72ec802a8003b7607a9a28cac6bb70cc622
--- /dev/null
+++ b/include/evalhyd/probabilistic/evaluator.h
@@ -0,0 +1,53 @@
+#ifndef EVALHYD_PROBABILIST_EVALUATOR_H
+#define EVALHYD_PROBABILIST_EVALUATOR_H
+
+#include <xtensor/xtensor.hpp>
+
+namespace evalhyd
+{
+ namespace probabilist
+ {
+ class Evaluator
+ {
+ private:
+ // members for input data
+ const xt::xtensor<double, 2>& q_obs;
+ const xt::xtensor<double, 2>& q_frc;
+ const xt::xtensor<double, 1>& q_thr;
+ // members for computational elements
+ xt::xtensor<double, 2> o_k;
+ xt::xtensor<double, 2> y_k;
+
+ public:
+ // constructor method
+ Evaluator(const xt::xtensor<double, 2>& obs,
+ const xt::xtensor<double, 2>& frc,
+ const xt::xtensor<double, 1>& thr) :
+ q_obs{obs}, q_frc{frc}, q_thr{thr} {};
+ // members for evaluation metrics
+ xt::xtensor<double, 1> bs;
+ xt::xtensor<double, 2> bs_crd;
+ xt::xtensor<double, 2> bs_lbd;
+ xt::xtensor<double, 1> bss;
+ // methods to compute elements
+ void calc_o_k();
+ void calc_y_k();
+ // methods to compute metrics
+ void calc_bs();
+ void calc_bs_crd();
+ void calc_bs_lbd();
+ void calc_bss();
+ };
+
+ xt::xtensor<double, 3> is_above_threshold(
+ const xt::xtensor<double, 2>&, const xt::xtensor<double, 1>&
+ );
+
+ xt::xtensor<double, 3> is_below_threshold(
+ const xt::xtensor<double, 2>&, const xt::xtensor<double, 1>&
+ );
+
+ }
+}
+
+#endif //EVALHYD_PROBABILIST_EVALUATOR_H
diff --git a/include/evalhyd/probabilistic/evaluator_brier.cpp b/include/evalhyd/probabilistic/evaluator_brier.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..251cef469bd5994445ab4417c78c67b7c1e7228e
--- /dev/null
+++ b/include/evalhyd/probabilistic/evaluator_brier.cpp
@@ -0,0 +1,266 @@
+#include <unordered_map>
+#include <xtensor/xmath.hpp>
+#include <xtensor/xview.hpp>
+#include <xtensor/xoperation.hpp>
+
+#include "evaluator.h"
+
+namespace eh = evalhyd;
+
+// NOTE ------------------------------------------------------------------------
+// All equations in metrics below are following notations from
+// "Wilks, D. S. (2011). Statistical methods in the atmospheric sciences.
+// Amsterdam; Boston: Elsevier Academic Press. ISBN: 9780123850225".
+// In particular, pp. 302-303, 332-333.
+// -----------------------------------------------------------------------------
+
+namespace evalhyd
+{
+ namespace probabilist
+ {
+ // Compute the Brier score (BS).
+ //
+ // \require o_k:
+ // Observed event outcome.
+ // shape: (thresholds, time)
+ // \require y_k:
+ // Event probability forecast.
+ // shape: (thresholds, time)
+ // \assign bs:
+ // 1D array of Brier Score for each threshold.
+ // shape: (thresholds,)
+ void Evaluator::calc_bs()
+ {
+ // return computed Brier score(s)
+ // $BS = \frac{1}{n} \sum_{k=1}^{n} (o_k - y_k)^2$
+ bs = xt::mean(xt::square(o_k - y_k), -1);
+ }
+
+ // Compute the calibration-refinement decomposition of the Brier score
+ // into reliability, resolution, and uncertainty.
+ //
+ // BS = reliability - resolution + uncertainty
+ //
+ // \require q_thr:
+ // 1D array of streamflow exceedance threshold(s).
+ // shape: (thresholds,)
+ // \require o_k:
+ // Observed event outcome.
+ // shape: (thresholds, time)
+ // \require y_k:
+ // Event probability forecast.
+ // shape: (thresholds, time)
+ // \assign bs_crd:
+ // 2D array of Brier score components (reliability, resolution,
+ // uncertainty) for each threshold.
+ // shape: (components, thresholds)
+ void Evaluator::calc_bs_crd()
+ {
+ // declare internal variables
+ // shape: (bins, thresholds, time)
+ xt::xtensor<double, 3> msk_bins;
+ // shape: (thresholds,)
+ xt::xtensor<double, 1> bar_o;
+ // shape: (bins, thresholds)
+ xt::xtensor<double, 2> N_i, bar_o_i;
+ // shape: (bins,)
+ xt::xtensor<double, 1> y_i;
+
+ // define some dimensions
+ std::size_t n = o_k.shape(1);
+ std::size_t n_thr = o_k.shape(0);
+ std::size_t n_mbr = q_frc.shape(0);
+ std::size_t n_cmp = 3;
+
+ // initialise output variable
+ // shape: (components, thresholds)
+ bs_crd = xt::zeros<double>({n_cmp, n_thr});
+
+ // compute range of forecast values $y_i$
+ y_i = xt::arange<double>(double(n_mbr + 1)) / n_mbr;
+
+ // compute mask to subsample time steps belonging to same forecast bin
+ // (where bins are defined as the range of forecast values)
+ msk_bins = xt::equal(
+ y_k, xt::view(y_i, xt::all(), xt::newaxis(), xt::newaxis())
+ );
+
+ // compute number of forecasts in each forecast bin $N_i$
+ N_i = xt::sum(msk_bins, -1);
+
+ // compute subsample relative frequency
+ // $\bar{o_i} = \frac{1}{N_i} \sum_{k \in N_i} o_k$
+ bar_o_i = xt::where(
+ N_i > 0,
+ xt::sum(
+ xt::where(
+ msk_bins,
+ xt::view(o_k, xt::newaxis(),
+ xt::all(), xt::all()),
+ 0.
+ ),
+ -1
+ ) / N_i,
+ 0.
+ );
+
+ // compute mean "climatology" relative frequency of the event
+ // $\bar{o} = \frac{1}{n} \sum_{k=1}^{n} o_k$
+ bar_o = xt::mean(o_k, -1);
+
+ // calculate reliability =
+ // $\frac{1}{n} \sum_{i=1}^{I} N_i (y_i - \bar{o_i})^2$
+ xt::row(bs_crd, 0) =
+ xt::sum(
+ xt::square(
+ xt::view(y_i, xt::all(), xt::newaxis())
+ - bar_o_i
+ ) * N_i,
+ 0
+ ) / n;
+
+ // calculate resolution =
+ // $\frac{1}{n} \sum_{i=1}^{I} N_i (\bar{o_i} - \bar{o})^2$
+ xt::row(bs_crd, 1) =
+ xt::sum(
+ xt::square(
+ bar_o_i - bar_o
+ ) * N_i,
+ 0
+ ) / n;
+
+ // calculate uncertainty = $\bar{o} (1 - \bar{o})$
+ xt::row(bs_crd, 2) = bar_o * (1 - bar_o);
+ }
+
+ // Compute the likelihood-based decomposition of the Brier score
+ // into type 2 bias, discrimination, and sharpness (a.k.a. refinement).
+ //
+ // BS = type 2 bias - discrimination + sharpness
+ //
+ // \require o_k:
+ // Observed event outcome.
+ // shape: (thresholds, time)
+ // \require y_k:
+ // Event probability forecast.
+ // shape: (thresholds, time)
+ // \return bs_lbd:
+ // 2D array of Brier score components (type 2 bias, discrimination,
+ // sharpness) for each threshold.
+ // shape: (components, thresholds)
+ void Evaluator::calc_bs_lbd()
+ {
+ // declare internal variables
+ // shape: (bins, thresholds, time)
+ xt::xtensor<double, 3> msk_bins;
+ // shape: (thresholds,)
+ xt::xtensor<double, 1> bar_y;
+ // shape: (bins, thresholds)
+ xt::xtensor<double, 2> M_j, bar_y_j;
+ // shape: (bins,)
+ xt::xtensor<double, 1> o_j;
+
+ // define some dimensions
+ std::size_t n = o_k.shape(1);
+ std::size_t n_thr = o_k.shape(0);
+ std::size_t n_cmp = 3;
+
+ // declare and initialise output variable
+ // shape: (components, thresholds)
+ bs_lbd = xt::zeros<double>({n_cmp, n_thr});
+
+ // set the range of observed values $o_j$
+ o_j = {0., 1.};
+
+ // compute mask to subsample time steps belonging to same observation bin
+ // (where bins are defined as the range of forecast values)
+ msk_bins = xt::equal(
+ o_k, xt::view(o_j, xt::all(), xt::newaxis(), xt::newaxis())
+ );
+
+ // compute number of observations in each observation bin $M_j$
+ M_j = xt::sum(msk_bins, -1);
+
+ // compute subsample relative frequency
+ // $\bar{y_j} = \frac{1}{M_j} \sum_{k \in M_j} y_k$
+ bar_y_j = xt::where(
+ M_j > 0,
+ xt::sum(
+ xt::where(
+ msk_bins,
+ xt::view(y_k, xt::newaxis(),
+ xt::all(), xt::all()),
+ 0.
+ ),
+ -1
+ ) / M_j,
+ 0.
+ );
+
+ // compute mean "climatology" forecast probability
+ // $\bar{y} = \frac{1}{n} \sum_{k=1}^{n} y_k$
+ bar_y = xt::mean(y_k, -1);
+
+ // calculate type 2 bias =
+ // $\frac{1}{n} \sum_{j=1}^{J} M_j (o_j - \bar{y_j})^2$
+ xt::row(bs_lbd, 0) =
+ xt::sum(
+ xt::square(
+ xt::view(o_j, xt::all(), xt::newaxis())
+ - bar_y_j
+ ) * M_j,
+ 0
+ ) / n;
+
+ // calculate discrimination =
+ // $\frac{1}{n} \sum_{j=1}^{J} M_j (\bar{y_j} - \bar{y})^2$
+ xt::row(bs_lbd, 1) =
+ xt::sum(
+ xt::square(
+ bar_y_j - bar_y
+ ) * M_j,
+ 0
+ ) / n;
+
+ // calculate sharpness =
+ // $\frac{1}{n} \sum_{k=1}^{n} (\bar{y_k} - \bar{y})^2$
+ xt::row(bs_lbd, 2) =
+ xt::sum(
+ xt::square(
+ y_k -
+ xt::view(bar_y, xt::all(), xt::newaxis())
+ ),
+ -1
+ ) / n;
+ }
+
+ // Compute the Brier skill score (BSS).
+ //
+ // \require o_k:
+ // Observed event outcome.
+ // shape: (thresholds, time)
+ // \require bs:
+ // Brier score(s).
+ // shape: (thresholds,)
+ // \assign bss:
+ // 1D array of Brier skill score for each threshold.
+ // shape: (thresholds,)
+ void Evaluator::calc_bss()
+ {
+ // calculate mean observed event outcome
+ // $\bar{o_k} = \frac{1}{n} \sum_{k=1}^{n} o_k$
+ xt::xtensor<double, 2> bar_o_k = xt::mean(o_k, -1,
+ xt::keep_dims);
+
+ // calculate reference Brier score(s)
+ // $BS_{ref} = \frac{1}{n} \sum_{k=1}^{n} (o_k - \bar{o_k})^2$
+ xt::xtensor<double, 1> bs_ref = xt::mean(
+ xt::square(o_k - bar_o_k), -1
+ );
+
+ // return computed Brier skill score(s)
+ // $BSS = 1 - \frac{BS}{BS_{ref}}
+ bss = 1 - (bs / bs_ref);
+ }
+ }
+}
diff --git a/include/evalhyd/probabilistic/evaluator_elements.cpp b/include/evalhyd/probabilistic/evaluator_elements.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7d36c5d390675f58e438eb5b66360add17c6c959
--- /dev/null
+++ b/include/evalhyd/probabilistic/evaluator_elements.cpp
@@ -0,0 +1,56 @@
+#include <xtensor/xmath.hpp>
+
+#include "evaluator.h"
+
+namespace evalhyd
+{
+ namespace probabilist
+ {
+ // Determine observed realisation of threshold(s) exceedance.
+ //
+ // \require q_obs:
+ // 2D array of streamflow observations.
+ // shape: (1, time)
+ // \require q_thr:
+ // 1D array of streamflow exceedance threshold(s).
+ // shape: (thresholds,)
+ // \assign o_k:
+ // 2D array of event observed outcome.
+ // shape: (thresholds, time)
+ void Evaluator::calc_o_k()
+ {
+ // determine observed realisation of threshold(s) exceedance
+ o_k = xt::squeeze(is_above_threshold(q_obs, q_thr));
+ }
+
+ // Determine forecast probability of threshold(s) exceedance to occur.
+ //
+ // \require q_frc:
+ // 2D array of streamflow forecasts.
+ // shape: (members, time)
+ // \require q_thr:
+ // 1D array of streamflow exceedance threshold(s).
+ // shape: (thresholds,)
+ // \assign y_k
+ // 2D array of event probability forecast.
+ // shape: (thresholds, time)
+ void Evaluator::calc_y_k()
+ {
+ // determine if members have exceeded threshold(s)
+ xt::xtensor<double, 3> e_frc =
+ is_above_threshold(q_frc, q_thr);
+
+ // calculate how many members have exceeded threshold(s)
+ 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 in other metrics
+ std::size_t n_mbr = q_frc.shape(0);
+
+ // determine probability of threshold(s) exceedance
+ y_k = n_frc / n_mbr;
+ }
+ }
+}
diff --git a/include/evalhyd/probabilistic/evaluator_utils.cpp b/include/evalhyd/probabilistic/evaluator_utils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5999e1fe6ad25a135c7d487fdaad687da1d11d72
--- /dev/null
+++ b/include/evalhyd/probabilistic/evaluator_utils.cpp
@@ -0,0 +1,44 @@
+#include <xtensor/xtensor.hpp>
+#include <xtensor/xview.hpp>
+
+namespace evalhyd
+{
+ namespace probabilist
+ {
+ // Determine whether flows are greater than given threshold(s).
+ //
+ // \param q:
+ // Array of streamflow data.
+ // shape: (1+, time)
+ // \param thr:
+ // Array of streamflow thresholds.
+ // shape: (thresholds,)
+ // \return
+ // Array of boolean-like threshold(s) exceedance.
+ // shape: (thresholds, 1+, time)
+ xt::xtensor<double, 3> is_above_threshold(
+ const xt::xtensor<double, 2> &q,
+ const xt::xtensor<double, 1> &thr
+ ) {
+ return q >= xt::view(thr, xt::all(), xt::newaxis(), xt::newaxis());
+ }
+
+ // Determine whether flows are strictly lower than given threshold(s)
+ //
+ // \param q:
+ // Array of streamflow data.
+ // shape: (1+, time)
+ // \param thr:
+ // Array of streamflow thresholds.
+ // shape: (thresholds,)
+ // \return
+ // Array of boolean-like threshold(s) exceedance.
+ // shape: (thresholds, 1+, time)
+ xt::xtensor<double, 3> is_below_threshold(
+ const xt::xtensor<double, 2> &q,
+ const xt::xtensor<double, 1> &thr
+ ) {
+ return q < xt::view(thr, xt::all(), xt::newaxis(), xt::newaxis());
+ }
+ }
+}