diff --git a/include/evalhyd/detail/probabilist/brier.hpp b/include/evalhyd/detail/probabilist/brier.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..70293fc9c354f2d7db11ae122d68f33b8b2e1fdc
--- /dev/null
+++ b/include/evalhyd/detail/probabilist/brier.hpp
@@ -0,0 +1,639 @@
+#ifndef EVALHYD_PROBABILIST_BRIER_HPP
+#define EVALHYD_PROBABILIST_BRIER_HPP
+
+#include <xtensor/xtensor.hpp>
+#include <xtensor/xview.hpp>
+#include <xtensor/xmasked_view.hpp>
+
+
+namespace evalhyd
+{
+ namespace probabilist
+ {
+ namespace elements
+ {
+ // Determine observed realisation of threshold(s) exceedance.
+ //
+ // \param q_obs
+ // Streamflow observations.
+ // shape: (time,)
+ // \param q_thr
+ // Streamflow exceedance threshold(s).
+ // shape: (thresholds,)
+ // \return
+ // Event observed outcome.
+ // shape: (thresholds, time)
+ template<class V1D2>
+ inline xt::xtensor<double, 2> calc_o_k(
+ const V1D2& q_obs,
+ const V1D2& q_thr
+ )
+ {
+ // determine observed realisation of threshold(s) exceedance
+ return q_obs >= xt::view(q_thr, xt::all(), xt::newaxis());
+ }
+
+ // Determine mean observed realisation of threshold(s) exceedance.
+ //
+ // \param o_k
+ // Event observed outcome.
+ // shape: (thresholds, time)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, time)
+ // \param b_exp
+ // Boostrap samples.
+ // shape: (samples, time slice)
+ // \param n_thr
+ // Number of thresholds.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Mean event observed outcome.
+ // shape: (subsets, samples, thresholds)
+ inline xt::xtensor<double, 3> calc_bar_o(
+ const xt::xtensor<double, 2>& o_k,
+ const xt::xtensor<bool, 2>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_thr,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto o_k_masked = xt::where(
+ xt::view(t_msk, xt::all(), xt::newaxis(), xt::all()),
+ o_k, NAN
+ );
+
+ // compute variable one sample at a time
+ xt::xtensor<double, 3> bar_o = xt::zeros<double>({n_msk, n_exp, n_thr});
+
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // apply the bootstrap sampling
+ auto o_k_masked_sampled =
+ xt::view(o_k_masked, xt::all(), xt::all(), b_exp[e]);
+
+ // compute mean "climatology" relative frequency of the event
+ // $\bar{o} = \frac{1}{n} \sum_{k=1}^{n} o_k$
+ xt::view(bar_o, xt::all(), e, xt::all()) =
+ xt::nanmean(o_k_masked_sampled, -1);
+ }
+
+ return bar_o;
+ }
+
+ // Determine forecast probability of threshold(s) exceedance to occur.
+ //
+ // \param q_prd
+ // Streamflow predictions.
+ // shape: (members, time)
+ // \param q_thr
+ // Streamflow exceedance threshold(s).
+ // shape: (thresholds,)
+ // \param n_mbr
+ // Number of ensemble members.
+ // \return
+ // Event probability forecast.
+ // shape: (thresholds, time)
+ template<class V2D4, class V1D2>
+ inline xt::xtensor<double, 2> calc_y_k(
+ const V2D4& q_prd,
+ const V1D2& q_thr,
+ std::size_t n_mbr
+ )
+ {
+ // determine if members have exceeded threshold(s)
+ auto e_frc =
+ q_prd >= xt::view(q_thr, xt::all(),
+ xt::newaxis(), xt::newaxis());
+
+ // calculate how many members have exceeded threshold(s)
+ auto n_frc = xt::sum(e_frc, 1);
+
+ // determine probability of threshold(s) exceedance
+ // /!\ probability calculation dividing by n (the number of
+ // members), not n+1 (the number of ranks) like in other metrics
+ return xt::cast<double>(n_frc) / n_mbr;
+ }
+ }
+
+ namespace intermediate
+ {
+ // Compute the Brier score for each time step.
+ //
+ // \param o_k
+ // Observed event outcome.
+ // shape: (thresholds, time)
+ // \param y_k
+ // Event probability forecast.
+ // shape: (thresholds, time)
+ // \return
+ // Brier score for each threshold for each time step.
+ // shape: (thresholds, time)
+ inline xt::xtensor<double, 2> calc_bs(
+ const xt::xtensor<double, 2>& o_k,
+ const xt::xtensor<double, 2>& y_k
+ )
+ {
+ // return computed Brier score(s)
+ // $bs = (o_k - y_k)^2$
+ return xt::square(o_k - y_k);
+ }
+ }
+
+ namespace metrics
+ {
+ // Compute the Brier score (BS).
+ //
+ // \param bs
+ // Brier score for each threshold for each time step.
+ // shape: (thresholds, time)
+ // \param q_thr
+ // Streamflow exceedance threshold(s).
+ // shape: (thresholds,)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, time)
+ // \param b_exp
+ // Boostrap samples.
+ // shape: (samples, time slice)
+ // \param n_thr
+ // Number of thresholds.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Brier score for each subset and for each threshold.
+ // shape: (subsets, samples, thresholds)
+ template <class V1D2>
+ inline xt::xtensor<double, 3> calc_BS(
+ const xt::xtensor<double, 2>& bs,
+ const V1D2& q_thr,
+ const xt::xtensor<bool, 2>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_thr,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // initialise output variable
+ // shape: (subsets, thresholds)
+ xt::xtensor<double, 3> BS =
+ xt::zeros<double>({n_msk, n_exp, n_thr});
+
+ // compute variable one mask at a time to minimise memory imprint
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto bs_masked = xt::where(xt::row(t_msk, m), bs, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // apply the bootstrap sampling
+ auto bs_masked_sampled =
+ xt::view(bs_masked, xt::all(), b_exp[e]);
+
+ // calculate the mean over the time steps
+ // $BS = \frac{1}{n} \sum_{k=1}^{n} (o_k - y_k)^2$
+ xt::view(BS, m, e, xt::all()) =
+ xt::nanmean(bs_masked_sampled, -1);
+ }
+ }
+
+ // assign NaN where thresholds were not provided (i.e. set as NaN)
+ xt::masked_view(
+ BS,
+ xt::isnan(xt::view(q_thr, xt::newaxis(),
+ xt::newaxis(), xt::all()))
+ ) = NAN;
+
+ return BS;
+ }
+
+ // Compute the calibration-refinement decomposition of the Brier score
+ // into reliability, resolution, and uncertainty.
+ //
+ // BS = reliability - resolution + uncertainty
+ //
+ // \param q_thr
+ // Streamflow exceedance threshold(s).
+ // shape: (thresholds,)
+ // \param o_k
+ // Observed event outcome.
+ // shape: (thresholds, time)
+ // \param y_k
+ // Event probability forecast.
+ // shape: (thresholds, time)
+ // \param bar_o
+ // Mean event observed outcome.
+ // shape: (subsets, samples, thresholds)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, time)
+ // \param b_exp
+ // Boostrap samples.
+ // shape: (samples, time slice)
+ // \param n_thr
+ // Number of thresholds.
+ // \param n_mbr
+ // Number of ensemble members.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Brier score components (reliability, resolution, uncertainty)
+ // for each subset and for each threshold.
+ // shape: (subsets, samples, thresholds, components)
+ template <class V1D2>
+ inline xt::xtensor<double, 4> calc_BS_CRD(
+ const V1D2& q_thr,
+ const xt::xtensor<double, 2>& o_k,
+ const xt::xtensor<double, 2>& y_k,
+ const xt::xtensor<double, 3>& bar_o,
+ const xt::xtensor<bool, 2>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_thr,
+ std::size_t n_mbr,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // declare internal variables
+ // shape: (bins, thresholds, time)
+ xt::xtensor<double, 3> msk_bins;
+ // shape: (bins, thresholds)
+ xt::xtensor<double, 2> N_i, bar_o_i;
+ // shape: (bins,)
+ xt::xtensor<double, 1> y_i;
+
+ // compute range of forecast values $y_i$
+ y_i = xt::arange<double>(double(n_mbr + 1)) / n_mbr;
+
+ // initialise output variable
+ // shape: (subsets, thresholds, components)
+ xt::xtensor<double, 4> BS_CRD =
+ xt::zeros<double>({n_msk, n_exp, n_thr, std::size_t {3}});
+
+ // compute variable one mask at a time to minimise memory imprint
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto o_k_masked = xt::where(xt::row(t_msk, m), o_k, NAN);
+ auto y_k_masked = xt::where(xt::row(t_msk, m), y_k, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // apply the bootstrap sampling
+ auto o_k_masked_sampled =
+ xt::view(o_k_masked, xt::all(), b_exp[e]);
+ auto y_k_masked_sampled =
+ xt::view(y_k_masked, xt::all(), b_exp[e]);
+ auto t_msk_sampled =
+ xt::view(xt::row(t_msk, m), b_exp[e]);
+ auto bar_o_sampled =
+ xt::view(bar_o, xt::all(), e, xt::all());
+
+ // calculate length of subset
+ auto l = xt::sum(t_msk_sampled);
+
+ // 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(
+ // force evaluation to avoid segfault
+ xt::eval(y_k_masked_sampled),
+ xt::view(y_i, xt::all(), xt::newaxis(), xt::newaxis())
+ );
+
+ // compute number of forecasts in each forecast bin $N_i$
+ N_i = xt::nansum(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::nansum(
+ xt::where(
+ msk_bins,
+ xt::view(o_k_masked_sampled, xt::newaxis(),
+ xt::all(), xt::all()),
+ 0.
+ ),
+ -1
+ ) / N_i,
+ 0.
+ );
+
+ // calculate reliability =
+ // $\frac{1}{n} \sum_{i=1}^{I} N_i (y_i - \bar{o_i})^2$
+ xt::view(BS_CRD, m, e, xt::all(), 0) =
+ xt::nansum(
+ xt::square(
+ xt::view(y_i, xt::all(), xt::newaxis())
+ - bar_o_i
+ ) * N_i,
+ 0
+ ) / l;
+
+ // calculate resolution =
+ // $\frac{1}{n} \sum_{i=1}^{I} N_i (\bar{o_i} - \bar{o})^2$
+ xt::view(BS_CRD, m, e, xt::all(), 1) =
+ xt::nansum(
+ xt::square(
+ bar_o_i - xt::view(bar_o_sampled, m)
+ ) * N_i,
+ 0
+ ) / l;
+
+ // calculate uncertainty = $\bar{o} (1 - \bar{o})$
+ xt::view(BS_CRD, m, e, xt::all(), 2) =
+ xt::view(bar_o_sampled, m)
+ * (1 - xt::view(bar_o_sampled, m));
+ }
+ }
+
+ // assign NaN where thresholds were not provided (i.e. set as NaN)
+ xt::masked_view(
+ BS_CRD,
+ xt::isnan(xt::view(q_thr, xt::newaxis(), xt::newaxis(),
+ xt::all(), xt::newaxis()))
+ ) = NAN;
+
+ return BS_CRD;
+ }
+
+ // Compute the likelihood-base rate decomposition of the Brier score
+ // into type 2 bias, discrimination, and sharpness (a.k.a. refinement).
+ //
+ // BS = type 2 bias - discrimination + sharpness
+ //
+ // \param q_thr
+ // Streamflow exceedance threshold(s).
+ // shape: (thresholds,)
+ // \param o_k
+ // Observed event outcome.
+ // shape: (thresholds, time)
+ // \param y_k
+ // Event probability forecast.
+ // shape: (thresholds, time)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, time)
+ // \param b_exp
+ // Boostrap samples.
+ // shape: (samples, time slice)
+ // \param n_thr
+ // Number of thresholds.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Brier score components (type 2 bias, discrimination, sharpness)
+ // for each subset and for each threshold.
+ // shape: (subsets, samples, thresholds, components)
+ template <class V1D2>
+ inline xt::xtensor<double, 4> calc_BS_LBD(
+ const V1D2& q_thr,
+ const xt::xtensor<double, 2>& o_k,
+ const xt::xtensor<double, 2>& y_k,
+ const xt::xtensor<bool, 2>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_thr,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // 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;
+
+ // set the range of observed values $o_j$
+ o_j = {0., 1.};
+
+ // declare and initialise output variable
+ xt::xtensor<double, 4> BS_LBD =
+ xt::zeros<double>({n_msk, n_exp, n_thr, std::size_t {3}});
+
+ // compute variable one mask at a time to minimise memory imprint
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto o_k_masked = xt::where(xt::row(t_msk, m), o_k, NAN);
+ auto y_k_masked = xt::where(xt::row(t_msk, m), y_k, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // apply the bootstrap sampling
+ auto o_k_masked_sampled =
+ xt::view(o_k_masked, xt::all(), b_exp[e]);
+ auto y_k_masked_sampled =
+ xt::view(y_k_masked, xt::all(), b_exp[e]);
+ auto t_msk_sampled =
+ xt::view(t_msk, xt::all(), b_exp[e]);
+
+ // calculate length of subset
+ auto l = xt::sum(xt::row(t_msk_sampled, m));
+
+ // 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(
+ // force evaluation to avoid segfault
+ xt::eval(o_k_masked_sampled),
+ xt::view(o_j, xt::all(), xt::newaxis(), xt::newaxis())
+ );
+
+ // compute number of observations in each observation bin $M_j$
+ M_j = xt::nansum(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::nansum(
+ xt::where(
+ msk_bins,
+ xt::view(
+ y_k_masked_sampled,
+ 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::nanmean(y_k_masked_sampled, -1);
+
+ // calculate type 2 bias =
+ // $\frac{1}{n} \sum_{j=1}^{J} M_j (o_j - \bar{y_j})^2$
+ xt::view(BS_LBD, m, e, xt::all(), 0) =
+ xt::nansum(
+ xt::square(
+ xt::view(o_j, xt::all(), xt::newaxis())
+ - bar_y_j
+ ) * M_j,
+ 0
+ ) / l;
+
+ // calculate discrimination =
+ // $\frac{1}{n} \sum_{j=1}^{J} M_j (\bar{y_j} - \bar{y})^2$
+ xt::view(BS_LBD, m, e, xt::all(), 1) =
+ xt::nansum(
+ xt::square(
+ bar_y_j - bar_y
+ ) * M_j,
+ 0
+ ) / l;
+
+ // calculate sharpness =
+ // $\frac{1}{n} \sum_{k=1}^{n} (\bar{y_k} - \bar{y})^2$
+ xt::view(BS_LBD, m, e, xt::all(), 2) =
+ xt::nansum(
+ xt::square(
+ y_k_masked_sampled
+ - xt::view(bar_y, xt::all(), xt::newaxis())
+ ),
+ -1
+ ) / l;
+ }
+
+ }
+
+ // assign NaN where thresholds were not provided (i.e. set as NaN)
+ xt::masked_view(
+ BS_LBD,
+ xt::isnan(xt::view(q_thr, xt::newaxis(), xt::newaxis(),
+ xt::all(), xt::newaxis()))
+ ) = NAN;
+
+ return BS_LBD;
+ }
+
+ // Compute the Brier skill score (BSS).
+ //
+ // \param bs
+ // Brier score for each threshold for each time step.
+ // shape: (thresholds, time)
+ // \param q_thr
+ // Streamflow exceedance threshold(s).
+ // shape: (thresholds,)
+ // \param o_k
+ // Observed event outcome.
+ // shape: (thresholds, time)
+ // \param bar_o
+ // Mean event observed outcome.
+ // shape: (subsets, samples, thresholds)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, time)
+ // \param b_exp
+ // Boostrap samples.
+ // shape: (samples, time slice)
+ // \param n_thr
+ // Number of thresholds.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Brier skill score for each subset and for each threshold.
+ // shape: (subsets, samples, thresholds)
+ template <class V1D2>
+ inline xt::xtensor<double, 3> calc_BSS(
+ const xt::xtensor<double, 2>& bs,
+ const V1D2& q_thr,
+ const xt::xtensor<double, 2>& o_k,
+ const xt::xtensor<double, 3>& bar_o,
+ const xt::xtensor<bool, 2>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_thr,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // declare and initialise output variable
+ xt::xtensor<double, 3> BSS =
+ xt::zeros<double>({n_msk, n_exp, n_thr});
+
+ // compute variable one mask at a time to minimise memory imprint
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto o_k_masked = xt::where(xt::row(t_msk, m), o_k, NAN);
+ auto bs_masked = xt::where(xt::row(t_msk, m), bs, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // apply the bootstrap sampling
+ auto o_k_masked_sampled =
+ xt::view(o_k_masked, xt::all(), b_exp[e]);
+ auto bs_masked_sampled =
+ xt::view(bs_masked, xt::all(), b_exp[e]);
+ auto bar_o_sampled =
+ xt::view(bar_o, xt::all(), e, xt::all());
+
+ // calculate reference Brier score(s)
+ // $bs_{ref} = \frac{1}{n} \sum_{k=1}^{n} (o_k - \bar{o})^2$
+ xt::xtensor<double, 2> bs_ref =
+ xt::nanmean(
+ xt::square(
+ o_k_masked_sampled -
+ xt::view(
+ xt::view(bar_o_sampled, m),
+ xt::all(), xt::newaxis()
+ )
+ ),
+ -1, xt::keep_dims
+ );
+
+ // compute Brier skill score(s)
+ // $BSS = \frac{1}{n} \sum_{k=1}^{n} 1 - \frac{bs}{bs_{ref}}
+ xt::view(BSS, m, e, xt::all()) =
+ xt::nanmean(
+ xt::where(
+ xt::equal(bs_ref, 0),
+ 0, 1 - (bs_masked_sampled / bs_ref)
+ ),
+ -1
+ );
+ }
+ }
+
+ // assign NaN where thresholds were not provided (i.e. set as NaN)
+ xt::masked_view(
+ BSS,
+ xt::isnan(xt::view(q_thr, xt::newaxis(), xt::newaxis(),
+ xt::all()))
+ ) = NAN;
+
+ return BSS;
+ }
+ }
+ }
+}
+
+#endif //EVALHYD_PROBABILIST_BRIER_HPP
diff --git a/include/evalhyd/detail/probabilist/evaluator.hpp b/include/evalhyd/detail/probabilist/evaluator.hpp
index 1cce16193fffc677a526f394a77a9d9e86cdf943..7f0e2eb2bbf07da0ca9238c8a8b9192e1881890c 100644
--- a/include/evalhyd/detail/probabilist/evaluator.hpp
+++ b/include/evalhyd/detail/probabilist/evaluator.hpp
@@ -4,13 +4,12 @@
#include <stdexcept>
#include <vector>
+#include <xtl/xoptional.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xview.hpp>
-#include <xtensor/xmasked_view.hpp>
-#include <xtensor/xslice.hpp>
-#include <xtensor/xoperation.hpp>
-#include <xtensor/xmath.hpp>
-#include <xtensor/xsort.hpp>
+
+#include "brier.hpp"
+#include "quantiles.hpp"
namespace evalhyd
@@ -64,10 +63,104 @@ namespace evalhyd
std::size_t n_exp;
// members for computational elements
- xt::xtensor<double, 2> o_k;
- xt::xtensor<double, 3> bar_o;
- xt::xtensor<double, 2> y_k;
- xt::xtensor<double, 2> q_qnt;
+ xtl::xoptional<xt::xtensor<double, 2>, bool> o_k;
+ xtl::xoptional<xt::xtensor<double, 3>, bool> bar_o;
+ xtl::xoptional<xt::xtensor<double, 2>, bool> y_k;
+ xtl::xoptional<xt::xtensor<double, 2>, bool> q_qnt;
+
+ // members for intermediate evaluation metrics
+ // (i.e. before the reduction along the temporal axis)
+ xtl::xoptional<xt::xtensor<double, 2>, bool> bs;
+ xtl::xoptional<xt::xtensor<double, 2>, bool> qs;
+ xtl::xoptional<xt::xtensor<double, 2>, bool> crps;
+
+ // members for evaluation metrics
+ xtl::xoptional<xt::xtensor<double, 3>, bool> BS;
+ xtl::xoptional<xt::xtensor<double, 4>, bool> BS_CRD;
+ xtl::xoptional<xt::xtensor<double, 4>, bool> BS_LBD;
+ xtl::xoptional<xt::xtensor<double, 3>, bool> BSS;
+ xtl::xoptional<xt::xtensor<double, 3>, bool> QS;
+ xtl::xoptional<xt::xtensor<double, 2>, bool> CRPS;
+
+ // methods to compute elements
+ xt::xtensor<double, 2> get_o_k()
+ {
+ if (!o_k.has_value())
+ {
+ o_k = elements::calc_o_k<view1d_xtensor2d_double_type>(
+ q_obs, q_thr
+ );
+ }
+ return o_k.value();
+ };
+
+ xt::xtensor<double, 3> get_bar_o()
+ {
+ if (!bar_o.has_value())
+ {
+ bar_o = elements::calc_bar_o(
+ get_o_k(), t_msk, b_exp, n_thr, n_msk, n_exp
+ );
+ }
+ return bar_o.value();
+ };
+
+ xt::xtensor<double, 2> get_y_k()
+ {
+ if (!y_k.has_value())
+ {
+ y_k = elements::calc_y_k<view2d_xtensor4d_double_type,
+ view1d_xtensor2d_double_type>(
+ q_prd, q_thr, n_mbr
+ );
+ }
+ return y_k.value();
+ };
+
+ xt::xtensor<double, 2> get_q_qnt()
+ {
+ if (!q_qnt.has_value())
+ {
+ q_qnt = elements::calc_q_qnt<view2d_xtensor4d_double_type>(
+ q_prd
+ );
+ }
+ return q_qnt.value();
+ };
+
+ // methods to compute intermediate metrics
+ xt::xtensor<double, 2> get_bs()
+ {
+ if (!bs.has_value())
+ {
+ bs = intermediate::calc_bs(
+ get_o_k(), get_y_k()
+ );
+ }
+ return bs.value();
+ };
+
+ xt::xtensor<double, 2> get_qs()
+ {
+ if (!qs.has_value())
+ {
+ qs = intermediate::calc_qs<view1d_xtensor2d_double_type>(
+ q_obs, get_q_qnt(), n_mbr
+ );
+ }
+ return qs.value();
+ };;
+
+ xt::xtensor<double, 2> get_crps()
+ {
+ if (!crps.has_value())
+ {
+ crps = intermediate::calc_crps(
+ get_qs(), n_mbr
+ );
+ }
+ return crps.value();
+ };
public:
// constructor method
@@ -109,694 +202,76 @@ namespace evalhyd
xt::view(t_msk, xt::all(), xt::keep(msk_nan)) = false;
};
- // members for intermediate evaluation metrics
- // (i.e. before the reduction along the temporal axis)
- xt::xtensor<double, 2> bs;
- xt::xtensor<double, 2> qs;
- xt::xtensor<double, 2> crps;
-
- // members for evaluation metrics
- xt::xtensor<double, 3> BS;
- xt::xtensor<double, 4> BS_CRD;
- xt::xtensor<double, 4> BS_LBD;
- xt::xtensor<double, 3> BSS;
- xt::xtensor<double, 3> QS;
- xt::xtensor<double, 2> CRPS;
-
- // methods to compute derived data
- void calc_q_qnt();
-
- // methods to compute elements
- void calc_o_k();
- void calc_bar_o();
- void calc_y_k();
-
- // methods to compute intermediate metrics
- void calc_bs();
- void calc_qs();
- void calc_crps();
-
// methods to compute metrics
- void calc_BS();
- void calc_BS_CRD();
- void calc_BS_LBD();
- void calc_BSS();
- void calc_QS();
- void calc_CRPS();
- };
-
- // --------------------------------------------------------------------
- // ELEMENTS
- // --------------------------------------------------------------------
-
- // Determine observed realisation of threshold(s) exceedance.
- //
- // \require q_obs:
- // Streamflow observations.
- // shape: (time,)
- // \require q_thr:
- // Streamflow exceedance threshold(s).
- // shape: (thresholds,)
- // \assign o_k:
- // Event observed outcome.
- // shape: (thresholds, time)
- template <class D2, class D4, class B4>
- void Evaluator<D2, D4, B4>::calc_o_k()
- {
- // determine observed realisation of threshold(s) exceedance
- o_k = q_obs >= xt::view(q_thr, xt::all(), xt::newaxis());
- }
-
- // Determine mean observed realisation of threshold(s) exceedance.
- //
- // \require o_k:
- // Event observed outcome.
- // shape: (thresholds, time)
- // \require t_msk:
- // Temporal subsets of the whole record.
- // shape: (subsets, time)
- // \assign bar_o:
- // Mean event observed outcome.
- // shape: (subsets, samples, thresholds)
- template <class D2, class D4, class B4>
- void Evaluator<D2, D4, B4>::calc_bar_o()
- {
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto o_k_masked = xt::where(
- xt::view(t_msk, xt::all(), xt::newaxis(), xt::all()),
- o_k, NAN
- );
-
- // compute variable one sample at a time
- bar_o = xt::zeros<double>({n_msk, n_exp, n_thr});
-
- for (std::size_t e = 0; e < n_exp; e++)
- {
- // apply the bootstrap sampling
- auto o_k_masked_sampled =
- xt::view(o_k_masked, xt::all(), xt::all(), b_exp[e]);
-
- // compute mean "climatology" relative frequency of the event
- // $\bar{o} = \frac{1}{n} \sum_{k=1}^{n} o_k$
- xt::view(bar_o, xt::all(), e, xt::all()) =
- xt::nanmean(o_k_masked_sampled, -1);
- }
- }
-
- // Determine forecast probability of threshold(s) exceedance to occur.
- //
- // \require q_prd:
- // Streamflow predictions.
- // shape: (members, time)
- // \require q_thr:
- // Streamflow exceedance threshold(s).
- // shape: (thresholds,)
- // \assign y_k:
- // Event probability forecast.
- // shape: (thresholds, time)
- template <class D2, class D4, class B4>
- void Evaluator<D2, D4, B4>::calc_y_k()
- {
- // determine if members have exceeded threshold(s)
- auto e_frc =
- q_prd
- >= xt::view(q_thr, xt::all(), xt::newaxis(), xt::newaxis());
-
- // calculate how many members have exceeded threshold(s)
- auto n_frc = xt::sum(e_frc, 1);
-
- // determine probability of threshold(s) exceedance
- // /!\ probability calculation dividing by n (the number of
- // members), not n+1 (the number of ranks) like in other metrics
- y_k = xt::cast<double>(n_frc) / n_mbr;
- }
-
- // Compute the forecast quantiles from the ensemble members.
- //
- // \require q_prd:
- // Streamflow predictions.
- // shape: (members, time)
- // \assign q_qnt:
- // Streamflow forecast quantiles.
- // shape: (quantiles, time)
- template <class D2, class D4, class B4>
- void Evaluator<D2, D4, B4>::calc_q_qnt()
- {
- q_qnt = xt::sort(q_prd, 0);
- }
-
- // --------------------------------------------------------------------
- // BRIER
- // --------------------------------------------------------------------
-
- // Compute the Brier score for each time step.
- //
- // \require o_k:
- // Observed event outcome.
- // shape: (thresholds, time)
- // \require y_k:
- // Event probability forecast.
- // shape: (thresholds, time)
- // \assign bs:
- // Brier score for each threshold for each time step.
- // shape: (thresholds, time)
- template <class D2, class D4, class B4>
- void Evaluator<D2, D4, B4>::calc_bs()
- {
- // return computed Brier score(s)
- // $bs = (o_k - y_k)^2$
- bs = xt::square(o_k - y_k);
- }
-
- // Compute the Brier score (BS).
- //
- // \require t_msk:
- // Temporal subsets of the whole record.
- // shape: (subsets, time)
- // \require bs:
- // Brier score for each threshold for each time step.
- // shape: (thresholds, time)
- // \assign BS:
- // Brier score for each subset and for each threshold.
- // shape: (subsets, samples, thresholds)
- template <class D2, class D4, class B4>
- void Evaluator<D2, D4, B4>::calc_BS()
- {
- // initialise output variable
- // shape: (subsets, thresholds)
- BS = xt::zeros<double>({n_msk, n_exp, n_thr});
-
- // compute variable one mask at a time to minimise memory imprint
- for (std::size_t m = 0; m < n_msk; m++)
+ xt::xtensor<double, 3> get_BS()
{
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto bs_masked = xt::where(xt::row(t_msk, m), bs, NAN);
-
- // compute variable one sample at a time
- for (std::size_t e = 0; e < n_exp; e++)
+ if (!BS.has_value())
{
- // apply the bootstrap sampling
- auto bs_masked_sampled =
- xt::view(bs_masked, xt::all(), b_exp[e]);
-
- // calculate the mean over the time steps
- // $BS = \frac{1}{n} \sum_{k=1}^{n} (o_k - y_k)^2$
- xt::view(BS, m, e, xt::all()) =
- xt::nanmean(bs_masked_sampled, -1);
+ BS = metrics::calc_BS<view1d_xtensor2d_double_type>(
+ get_bs(), q_thr, t_msk, b_exp, n_thr, n_msk, n_exp
+ );
}
- }
-
- // assign NaN where thresholds were not provided (i.e. set as NaN)
- xt::masked_view(
- BS,
- xt::isnan(xt::view(q_thr, xt::newaxis(), xt::newaxis(), xt::all()))
- ) = NAN;
- }
-
- // Compute the calibration-refinement decomposition of the Brier score
- // into reliability, resolution, and uncertainty.
- //
- // BS = reliability - resolution + uncertainty
- //
- // \require q_thr:
- // Streamflow exceedance threshold(s).
- // shape: (thresholds,)
- // \require t_msk:
- // Temporal subsets of the whole record.
- // shape: (subsets, time)
- // \require o_k:
- // Observed event outcome.
- // shape: (thresholds, time)
- // \require y_k:
- // Event probability forecast.
- // shape: (thresholds, time)
- // \require bar_o:
- // Mean event observed outcome.
- // shape: (subsets, thresholds)
- // \assign BS_CRD:
- // Brier score components (reliability, resolution, uncertainty)
- // for each subset and for each threshold.
- // shape: (subsets, samples, thresholds, components)
- template <class D2, class D4, class B4>
- void Evaluator<D2, D4, B4>::calc_BS_CRD()
- {
- // declare internal variables
- // shape: (bins, thresholds, time)
- xt::xtensor<double, 3> msk_bins;
- // shape: (bins, thresholds)
- xt::xtensor<double, 2> N_i, bar_o_i;
- // shape: (bins,)
- xt::xtensor<double, 1> y_i;
-
- // compute range of forecast values $y_i$
- y_i = xt::arange<double>(double(n_mbr + 1)) / n_mbr;
-
- // initialise output variable
- // shape: (subsets, thresholds, components)
- BS_CRD = xt::zeros<double>({n_msk, n_exp, n_thr, std::size_t {3}});
-
- // compute variable one mask at a time to minimise memory imprint
- for (std::size_t m = 0; m < n_msk; m++)
- {
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto o_k_masked = xt::where(xt::row(t_msk, m), o_k, NAN);
- auto y_k_masked = xt::where(xt::row(t_msk, m), y_k, NAN);
+ return BS.value();
+ };
- // compute variable one sample at a time
- for (std::size_t e = 0; e < n_exp; e++)
+ xt::xtensor<double, 4> get_BS_CRD()
+ {
+ if (!BS_CRD.has_value())
{
- // apply the bootstrap sampling
- auto o_k_masked_sampled =
- xt::view(o_k_masked, xt::all(), b_exp[e]);
- auto y_k_masked_sampled =
- xt::view(y_k_masked, xt::all(), b_exp[e]);
- auto t_msk_sampled =
- xt::view(xt::row(t_msk, m), b_exp[e]);
- auto bar_o_sampled =
- xt::view(bar_o, xt::all(), e, xt::all());
-
- // calculate length of subset
- auto l = xt::sum(t_msk_sampled);
-
- // 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(
- // force evaluation to avoid segfault
- xt::eval(y_k_masked_sampled),
- xt::view(y_i, xt::all(), xt::newaxis(), xt::newaxis())
- );
-
- // compute number of forecasts in each forecast bin $N_i$
- N_i = xt::nansum(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::nansum(
- xt::where(
- msk_bins,
- xt::view(o_k_masked_sampled, xt::newaxis(),
- xt::all(), xt::all()),
- 0.
- ),
- -1
- ) / N_i,
- 0.
+ BS_CRD = metrics::calc_BS_CRD<view1d_xtensor2d_double_type>(
+ q_thr, get_o_k(), get_y_k(), get_bar_o(), t_msk,
+ b_exp, n_thr, n_mbr, n_msk, n_exp
);
-
- // calculate reliability =
- // $\frac{1}{n} \sum_{i=1}^{I} N_i (y_i - \bar{o_i})^2$
- xt::view(BS_CRD, m, e, xt::all(), 0) =
- xt::nansum(
- xt::square(
- xt::view(y_i, xt::all(), xt::newaxis())
- - bar_o_i
- ) * N_i,
- 0
- ) / l;
-
- // calculate resolution =
- // $\frac{1}{n} \sum_{i=1}^{I} N_i (\bar{o_i} - \bar{o})^2$
- xt::view(BS_CRD, m, e, xt::all(), 1) =
- xt::nansum(
- xt::square(
- bar_o_i - xt::view(bar_o_sampled, m)
- ) * N_i,
- 0
- ) / l;
-
- // calculate uncertainty = $\bar{o} (1 - \bar{o})$
- xt::view(BS_CRD, m, e, xt::all(), 2) =
- xt::view(bar_o_sampled, m)
- * (1 - xt::view(bar_o_sampled, m));
}
- }
-
- // assign NaN where thresholds were not provided (i.e. set as NaN)
- xt::masked_view(
- BS_CRD,
- xt::isnan(xt::view(q_thr, xt::newaxis(), xt::newaxis(),
- xt::all(), xt::newaxis()))
- ) = NAN;
- }
-
- // Compute the likelihood-base rate decomposition of the Brier score
- // into type 2 bias, discrimination, and sharpness (a.k.a. refinement).
- //
- // BS = type 2 bias - discrimination + sharpness
- //
- // \require t_msk:
- // Temporal subsets of the whole record.
- // shape: (subsets, time)
- // \require o_k:
- // Observed event outcome.
- // shape: (thresholds, time)
- // \require y_k:
- // Event probability forecast.
- // shape: (thresholds, time)
- // \return BS_LBD:
- // Brier score components (type 2 bias, discrimination, sharpness)
- // for each subset and for each threshold.
- // shape: (subsets, samples, thresholds, components)
- template <class D2, class D4, class B4>
- void Evaluator<D2, D4, B4>::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;
-
- // set the range of observed values $o_j$
- o_j = {0., 1.};
-
- // declare and initialise output variable
- // shape: (subsets, thresholds, components)
- BS_LBD = xt::zeros<double>({n_msk, n_exp, n_thr, std::size_t {3}});
-
- // compute variable one mask at a time to minimise memory imprint
- for (std::size_t m = 0; m < n_msk; m++)
- {
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto o_k_masked = xt::where(xt::row(t_msk, m), o_k, NAN);
- auto y_k_masked = xt::where(xt::row(t_msk, m), y_k, NAN);
+ return BS_CRD.value();
+ };
- // compute variable one sample at a time
- for (std::size_t e = 0; e < n_exp; e++)
+ xt::xtensor<double, 4> get_BS_LBD()
+ {
+ if (!BS_LBD.has_value())
{
- // apply the bootstrap sampling
- auto o_k_masked_sampled =
- xt::view(o_k_masked, xt::all(), b_exp[e]);
- auto y_k_masked_sampled =
- xt::view(y_k_masked, xt::all(), b_exp[e]);
- auto t_msk_sampled =
- xt::view(t_msk, xt::all(), b_exp[e]);
-
- // calculate length of subset
- auto l = xt::sum(xt::row(t_msk_sampled, m));
-
- // 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(
- // force evaluation to avoid segfault
- xt::eval(o_k_masked_sampled),
- xt::view(o_j, xt::all(), xt::newaxis(), xt::newaxis())
- );
-
- // compute number of observations in each observation bin $M_j$
- M_j = xt::nansum(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::nansum(
- xt::where(
- msk_bins,
- xt::view(
- y_k_masked_sampled,
- xt::newaxis(), xt::all(), xt::all()
- ),
- 0.
- ),
- -1
- ) / M_j,
- 0.
+ BS_LBD = metrics::calc_BS_LBD<view1d_xtensor2d_double_type>(
+ q_thr, get_o_k(), get_y_k(), t_msk,
+ b_exp, n_thr, n_msk, n_exp
);
-
- // compute mean "climatology" forecast probability
- // $\bar{y} = \frac{1}{n} \sum_{k=1}^{n} y_k$
- bar_y = xt::nanmean(y_k_masked_sampled, -1);
-
- // calculate type 2 bias =
- // $\frac{1}{n} \sum_{j=1}^{J} M_j (o_j - \bar{y_j})^2$
- xt::view(BS_LBD, m, e, xt::all(), 0) =
- xt::nansum(
- xt::square(
- xt::view(o_j, xt::all(), xt::newaxis())
- - bar_y_j
- ) * M_j,
- 0
- ) / l;
-
- // calculate discrimination =
- // $\frac{1}{n} \sum_{j=1}^{J} M_j (\bar{y_j} - \bar{y})^2$
- xt::view(BS_LBD, m, e, xt::all(), 1) =
- xt::nansum(
- xt::square(
- bar_y_j - bar_y
- ) * M_j,
- 0
- ) / l;
-
- // calculate sharpness =
- // $\frac{1}{n} \sum_{k=1}^{n} (\bar{y_k} - \bar{y})^2$
- xt::view(BS_LBD, m, e, xt::all(), 2) =
- xt::nansum(
- xt::square(
- y_k_masked_sampled
- - xt::view(bar_y, xt::all(), xt::newaxis())
- ),
- -1
- ) / l;
}
+ return BS_LBD.value();
+ };
- }
-
- // assign NaN where thresholds were not provided (i.e. set as NaN)
- xt::masked_view(
- BS_LBD,
- xt::isnan(xt::view(q_thr, xt::newaxis(), xt::newaxis(),
- xt::all(), xt::newaxis()))
- ) = NAN;
- }
-
- // Compute the Brier skill score (BSS).
- //
- // \require t_msk:
- // Temporal subsets of the whole record.
- // shape: (subsets, time)
- // \require o_k:
- // Observed event outcome.
- // shape: (thresholds, time)
- // \require bar_o:
- // Mean event observed outcome.
- // shape: (subsets, thresholds)
- // \require bs:
- // Brier scores for each time step.
- // shape: (thresholds, time)
- // \assign BSS:
- // Brier skill score for each subset and for each threshold.
- // shape: (subsets, samples, thresholds)
- template <class D2, class D4, class B4>
- void Evaluator<D2, D4, B4>::calc_BSS()
- {
- // declare and initialise output variable
- // shape: (subsets, thresholds)
- BSS = xt::zeros<double>({n_msk, n_exp, n_thr});
-
- // compute variable one mask at a time to minimise memory imprint
- for (std::size_t m = 0; m < n_msk; m++)
+ xt::xtensor<double, 3> get_BSS()
{
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto o_k_masked = xt::where(xt::row(t_msk, m), o_k, NAN);
- auto bs_masked = xt::where(xt::row(t_msk, m), bs, NAN);
-
- // compute variable one sample at a time
- for (std::size_t e = 0; e < n_exp; e++)
+ if (!BSS.has_value())
{
- // apply the bootstrap sampling
- auto o_k_masked_sampled =
- xt::view(o_k_masked, xt::all(), b_exp[e]);
- auto bs_masked_sampled =
- xt::view(bs_masked, xt::all(), b_exp[e]);
- auto bar_o_sampled =
- xt::view(bar_o, xt::all(), e, xt::all());
-
- // calculate reference Brier score(s)
- // $bs_{ref} = \frac{1}{n} \sum_{k=1}^{n} (o_k - \bar{o})^2$
- xt::xtensor<double, 2> bs_ref =
- xt::nanmean(
- xt::square(
- o_k_masked_sampled -
- xt::view(
- xt::view(bar_o_sampled, m),
- xt::all(), xt::newaxis()
- )
- ),
- -1, xt::keep_dims
- );
-
- // compute Brier skill score(s)
- // $BSS = \frac{1}{n} \sum_{k=1}^{n} 1 - \frac{bs}{bs_{ref}}
- xt::view(BSS, m, e, xt::all()) =
- xt::nanmean(
- xt::where(
- xt::equal(bs_ref, 0),
- 0, 1 - (bs_masked_sampled / bs_ref)
- ),
- -1
- );
+ BSS = metrics::calc_BSS<view1d_xtensor2d_double_type>(
+ get_bs(), q_thr, get_o_k(), get_bar_o(), t_msk,
+ b_exp, n_thr, n_msk, n_exp
+ );
}
- }
-
- // assign NaN where thresholds were not provided (i.e. set as NaN)
- xt::masked_view(
- BSS,
- xt::isnan(xt::view(q_thr, xt::newaxis(), xt::newaxis(),
- xt::all()))
- ) = NAN;
- }
-
- // --------------------------------------------------------------------
- // QUANTILES
- // --------------------------------------------------------------------
-
- // Compute the quantile scores for each time step.
- //
- // \require q_obs:
- // Streamflow observations.
- // shape: (time,)
- // \require q_qnt:
- // Streamflow quantiles.
- // shape: (quantiles, time)
- // \assign qs:
- // Quantile scores for each time step.
- // shape: (quantiles, time)
- template <class D2, class D4, class B4>
- void Evaluator<D2, D4, B4>::calc_qs()
- {
- // compute the quantile order $alpha$
- xt::xtensor<double, 1> alpha =
- xt::arange<double>(1., double(n_mbr + 1))
- / double(n_mbr + 1);
-
- // calculate the difference
- xt::xtensor<double, 2> diff = q_qnt - q_obs;
-
- // calculate the quantile scores
- qs = xt::where(
- diff > 0,
- 2 * (1 - xt::view(alpha, xt::all(), xt::newaxis())) * diff,
- - 2 * xt::view(alpha, xt::all(), xt::newaxis()) * diff
- );
- }
-
- // Compute the quantile score (QS).
- //
- // \require t_msk:
- // Temporal subsets of the whole record.
- // shape: (subsets, time)
- // \require qs:
- // Quantile scores for each time step.
- // shape: (quantiles, time)
- // \assign QS:
- // Quantile scores.
- // shape: (subsets, quantiles)
- template <class D2, class D4, class B4>
- void Evaluator<D2, D4, B4>::calc_QS()
- {
- // initialise output variable
- // shape: (subsets, quantiles)
- QS = xt::zeros<double>({n_msk, n_exp, n_mbr});
+ return BSS.value();
+ };
- // compute variable one mask at a time to minimise memory imprint
- for (std::size_t m = 0; m < n_msk; m++)
+ xt::xtensor<double, 3> get_QS()
{
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto qs_masked = xt::where(xt::row(t_msk, m), qs, NAN);
-
- // compute variable one sample at a time
- for (std::size_t e = 0; e < n_exp; e++)
+ if (!QS.has_value())
{
- // apply the bootstrap sampling
- auto qs_masked_sampled =
- xt::view(qs_masked, xt::all(), b_exp[e]);
-
- // calculate the mean over the time steps
- // $QS = \frac{1}{n} \sum_{k=1}^{n} qs$
- xt::view(QS, m, e, xt::all()) =
- xt::nanmean(qs_masked_sampled, -1);
+ QS = metrics::calc_QS(
+ get_qs(), t_msk, b_exp, n_mbr, n_msk, n_exp
+ );
}
- }
- }
-
- // Compute the continuous rank probability score(s) based
- // on quantile scores for each time step, and integrating using the
- // trapezoidal rule.
- //
- // /!\ The number of quantiles must be sufficiently large so that the
- // cumulative distribution is smooth enough for the numerical
- // integration to be accurate.
- //
- // \require qs:
- // Quantile scores for each time step.
- // shape: (quantiles, time)
- // \assign crps:
- // CRPS for each time step.
- // shape: (1, time)
- template <class D2, class D4, class B4>
- void Evaluator<D2, D4, B4>::calc_crps()
- {
- // integrate with trapezoidal rule
- crps = xt::view(
- // xt::trapz(y, dx=1/(n+1), axis=0)
- xt::trapz(qs, 1./(double(n_mbr) + 1.), 0),
- xt::newaxis(), xt::all()
- );
- }
-
- // Compute the continuous rank probability score (CRPS) based
- // on quantile scores.
- //
- // \require t_msk:
- // Temporal subsets of the whole record.
- // shape: (subsets, time)
- // \require crps:
- // CRPS for each time step.
- // shape: (1, time)
- // \assign CRPS:
- // CRPS.
- // shape: (subsets,)
- template <class D2, class D4, class B4>
- void Evaluator<D2, D4, B4>::calc_CRPS()
- {
- // initialise output variable
- // shape: (subsets,)
- CRPS = xt::zeros<double>({n_msk, n_exp});
+ return QS.value();
+ };
- // compute variable one mask at a time to minimise memory imprint
- for (std::size_t m = 0; m < n_msk; m++)
+ xt::xtensor<double, 2> get_CRPS()
{
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto crps_masked = xt::where(xt::row(t_msk, m), crps, NAN);
-
- // compute variable one sample at a time
- for (std::size_t e = 0; e < n_exp; e++)
+ if (!CRPS.has_value())
{
- // apply the bootstrap sampling
- auto crps_masked_sampled =
- xt::view(crps_masked, xt::all(), b_exp[e]);
-
- // calculate the mean over the time steps
- // $CRPS = \frac{1}{n} \sum_{k=1}^{n} crps$
- xt::view(CRPS, m, e) =
- xt::squeeze(xt::nanmean(crps_masked_sampled, -1));
+ CRPS = metrics::calc_CRPS(
+ get_crps(), t_msk, b_exp, n_msk, n_exp
+ );
}
- }
- }
+ return CRPS.value();
+ };
+ };
}
}
diff --git a/include/evalhyd/detail/probabilist/quantiles.hpp b/include/evalhyd/detail/probabilist/quantiles.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c05273d8ef8d408537484254ffa900e8a57fa353
--- /dev/null
+++ b/include/evalhyd/detail/probabilist/quantiles.hpp
@@ -0,0 +1,216 @@
+#ifndef EVALHYD_PROBABILIST_QUANTILES_HPP
+#define EVALHYD_PROBABILIST_QUANTILES_HPP
+
+#include <xtensor/xtensor.hpp>
+#include <xtensor/xview.hpp>
+#include <xtensor/xsort.hpp>
+
+
+namespace evalhyd
+{
+ namespace probabilist
+ {
+ namespace elements
+ {
+ // Compute the forecast quantiles from the ensemble members.
+ //
+ // \param q_prd
+ // Streamflow predictions.
+ // shape: (members, time)
+ // \return
+ // Streamflow forecast quantiles.
+ // shape: (quantiles, time)
+ template <class V2D4>
+ inline xt::xtensor<double, 2> calc_q_qnt(
+ const V2D4& q_prd
+ )
+ {
+ return xt::sort(q_prd, 0);
+ }
+ }
+
+ namespace intermediate
+ {
+ // Compute the quantile scores for each time step.
+ //
+ // \param q_obs
+ // Streamflow observations.
+ // shape: (time,)
+ // \param q_qnt
+ // Streamflow quantiles.
+ // shape: (quantiles, time)
+ // \return
+ // Quantile scores for each time step.
+ // shape: (quantiles, time)
+ template <class V1D2>
+ inline xt::xtensor<double, 2> calc_qs(
+ const V1D2 &q_obs,
+ const xt::xtensor<double, 2>& q_qnt,
+ std::size_t n_mbr
+ )
+ {
+ // compute the quantile order $alpha$
+ xt::xtensor<double, 1> alpha =
+ xt::arange<double>(1., double(n_mbr + 1))
+ / double(n_mbr + 1);
+
+ // calculate the difference
+ xt::xtensor<double, 2> diff = q_qnt - q_obs;
+
+ // calculate the quantile scores
+ xt::xtensor<double, 2> qs = xt::where(
+ diff > 0,
+ 2 * (1 - xt::view(alpha, xt::all(), xt::newaxis())) * diff,
+ - 2 * xt::view(alpha, xt::all(), xt::newaxis()) * diff
+ );
+
+ return qs;
+ }
+
+ // Compute the continuous rank probability score(s) based
+ // on quantile scores for each time step, and integrating using the
+ // trapezoidal rule.
+ //
+ // /!\ The number of quantiles must be sufficiently large so that the
+ // cumulative distribution is smooth enough for the numerical
+ // integration to be accurate.
+ //
+ // \param qs
+ // Quantile scores for each time step.
+ // shape: (quantiles, time)
+ // \return
+ // CRPS for each time step.
+ // shape: (1, time)
+ inline xt::xtensor<double, 2> calc_crps(
+ const xt::xtensor<double, 2>& qs,
+ std::size_t n_mbr
+ )
+ {
+ // integrate with trapezoidal rule
+ return xt::view(
+ // xt::trapz(y, dx=1/(n+1), axis=0)
+ xt::trapz(qs, 1./(double(n_mbr) + 1.), 0),
+ xt::newaxis(), xt::all()
+ );
+ }
+ }
+
+ namespace metrics
+ {
+ // Compute the quantile score (QS).
+ //
+ // \param qs
+ // Quantile scores for each time step.
+ // shape: (quantiles, time)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, time)
+ // \param b_exp
+ // Boostrap samples.
+ // shape: (samples, time slice)
+ // \param n_mbr
+ // Number of ensemble members.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Quantile scores.
+ // shape: (subsets, samples, quantiles)
+ inline xt::xtensor<double, 3> calc_QS(
+ const xt::xtensor<double, 2>& qs,
+ const xt::xtensor<bool, 2>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_mbr,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // initialise output variable
+ // shape: (subsets, quantiles)
+ xt::xtensor<double, 3> QS =
+ xt::zeros<double>({n_msk, n_exp, n_mbr});
+
+ // compute variable one mask at a time to minimise memory imprint
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto qs_masked = xt::where(xt::row(t_msk, m), qs, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // apply the bootstrap sampling
+ auto qs_masked_sampled =
+ xt::view(qs_masked, xt::all(), b_exp[e]);
+
+ // calculate the mean over the time steps
+ // $QS = \frac{1}{n} \sum_{k=1}^{n} qs$
+ xt::view(QS, m, e, xt::all()) =
+ xt::nanmean(qs_masked_sampled, -1);
+ }
+ }
+
+ return QS;
+ }
+
+ // Compute the continuous rank probability score (CRPS) based
+ // on quantile scores.
+ //
+ // \param crps
+ // CRPS for each time step.
+ // shape: (1, time)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, time)
+ // \param b_exp
+ // Boostrap samples.
+ // shape: (samples, time slice)
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // CRPS.
+ // shape: (subsets, samples)
+ inline xt::xtensor<double, 2> calc_CRPS(
+ const xt::xtensor<double, 2>& crps,
+ const xt::xtensor<bool, 2>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // initialise output variable
+ // shape: (subsets,)
+ xt::xtensor<double, 2> CRPS = xt::zeros<double>({n_msk, n_exp});
+
+ // compute variable one mask at a time to minimise memory imprint
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto crps_masked = xt::where(xt::row(t_msk, m), crps, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // apply the bootstrap sampling
+ auto crps_masked_sampled =
+ xt::view(crps_masked, xt::all(), b_exp[e]);
+
+ // calculate the mean over the time steps
+ // $CRPS = \frac{1}{n} \sum_{k=1}^{n} crps$
+ xt::view(CRPS, m, e) =
+ xt::squeeze(xt::nanmean(crps_masked_sampled, -1));
+ }
+ }
+
+ return CRPS;
+ }
+ }
+ }
+}
+
+#endif //EVALHYD_PROBABILIST_QUANTILES_HPP
diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp
index c41277ae5d25ce7a2dc876a3eefe80a58bef5425..350b923f70fb59a47cf2056e411e351ad8c69a7c 100644
--- a/include/evalhyd/evalp.hpp
+++ b/include/evalhyd/evalp.hpp
@@ -228,7 +228,6 @@ namespace evalhyd
}
}
-
// retrieve dimensions
std::size_t n_sit = q_prd_.shape(0);
std::size_t n_ltm = q_prd_.shape(1);
@@ -250,29 +249,6 @@ namespace evalhyd
dim["QS"] = {n_sit, n_ltm, n_msk, n_exp, n_mbr};
dim["CRPS"] = {n_sit, n_ltm, n_msk, n_exp};
- // declare maps for memoisation purposes
- std::unordered_map<std::string, std::vector<std::string>> elt;
- std::unordered_map<std::string, std::vector<std::string>> dep;
-
- // register potentially recurring computation elements across metrics
- elt["bs"] = {"o_k", "y_k"};
- elt["BSS"] = {"o_k", "bar_o"};
- elt["BS_CRD"] = {"o_k", "bar_o", "y_k"};
- elt["BS_LBD"] = {"o_k", "y_k"};
- elt["qs"] = {"q_qnt"};
-
- // register nested metrics (i.e. metric dependent on another metric)
- dep["BS"] = {"bs"};
- dep["BSS"] = {"bs"};
- dep["QS"] = {"qs"};
- dep["CRPS"] = {"qs", "crps"};
-
- // determine required elt/dep to be pre-computed
- std::vector<std::string> req_elt;
- std::vector<std::string> req_dep;
-
- utils::find_requirements(metrics, elt, dep, req_elt, req_dep);
-
// generate masks from conditions if provided
auto gen_msk = [&]()
{
@@ -353,44 +329,6 @@ namespace evalhyd
q_obs_v, q_prd_v, q_thr_v, t_msk_v, b_exp
);
- // pre-compute required elt
- for (const auto& element : req_elt)
- {
- if ( element == "o_k" )
- {
- evaluator.calc_o_k();
- }
- else if ( element == "bar_o" )
- {
- evaluator.calc_bar_o();
- }
- else if ( element == "y_k" )
- {
- evaluator.calc_y_k();
- }
- else if ( element == "q_qnt" )
- {
- evaluator.calc_q_qnt();
- }
- }
-
- // pre-compute required dep
- for (const auto& dependency : req_dep)
- {
- if ( dependency == "bs" )
- {
- evaluator.calc_bs();
- }
- else if ( dependency == "qs" )
- {
- evaluator.calc_qs();
- }
- else if ( dependency == "crps" )
- {
- evaluator.calc_crps();
- }
- }
-
// retrieve or compute requested metrics
for (std::size_t m = 0; m < metrics.size(); m++)
{
@@ -398,76 +336,45 @@ namespace evalhyd
if ( metric == "BS" )
{
- if (std::find(req_dep.begin(), req_dep.end(), metric)
- == req_dep.end())
- {
- evaluator.calc_BS();
- }
// (sites, lead times, subsets, samples, thresholds)
xt::view(r[m], s, l, xt::all(), xt::all(), xt::all()) =
- uncertainty::summarise(evaluator.BS, summary);
+ uncertainty::summarise(evaluator.get_BS(), summary);
}
else if ( metric == "BSS" )
{
- if (std::find(req_dep.begin(), req_dep.end(), metric)
- == req_dep.end())
- {
- evaluator.calc_BSS();
- }
// (sites, lead times, subsets, samples, thresholds)
xt::view(r[m], s, l, xt::all(), xt::all(), xt::all()) =
- uncertainty::summarise(evaluator.BSS, summary);
+ uncertainty::summarise(evaluator.get_BSS(), summary);
}
else if ( metric == "BS_CRD" )
{
- if (std::find(req_dep.begin(), req_dep.end(), metric)
- == req_dep.end())
- {
- evaluator.calc_BS_CRD();
- }
// (sites, lead times, subsets, samples, thresholds, components)
xt::view(r[m], s, l, xt::all(), xt::all(), xt::all(), xt::all()) =
- uncertainty::summarise(evaluator.BS_CRD, summary);
+ uncertainty::summarise(evaluator.get_BS_CRD(), summary);
}
else if ( metric == "BS_LBD" )
{
- if (std::find(req_dep.begin(), req_dep.end(), metric)
- == req_dep.end())
- {
- evaluator.calc_BS_LBD();
- }
-
// (sites, lead times, subsets, samples, thresholds, components)
xt::view(r[m], s, l, xt::all(), xt::all(), xt::all(), xt::all()) =
- uncertainty::summarise(evaluator.BS_LBD, summary);
+ uncertainty::summarise(evaluator.get_BS_LBD(), summary);
}
else if ( metric == "QS" )
{
- if (std::find(req_dep.begin(), req_dep.end(), metric)
- == req_dep.end())
- {
- evaluator.calc_QS();
- }
// (sites, lead times, subsets, samples, quantiles)
xt::view(r[m], s, l, xt::all(), xt::all(), xt::all()) =
- uncertainty::summarise(evaluator.QS, summary);
+ uncertainty::summarise(evaluator.get_QS(), summary);
}
else if ( metric == "CRPS" )
{
- if (std::find(req_dep.begin(), req_dep.end(), metric)
- == req_dep.end())
- {
- evaluator.calc_CRPS();
- }
// (sites, lead times, subsets, samples)
xt::view(r[m], s, l, xt::all(), xt::all()) =
- uncertainty::summarise(evaluator.CRPS, summary);
+ uncertainty::summarise(evaluator.get_CRPS(), summary);
}
}
}
}
return r;
- };
+ }
}
#endif //EVALHYD_EVALP_HPP