diff --git a/include/evalhyd/detail/probabilist/evaluator.hpp b/include/evalhyd/detail/probabilist/evaluator.hpp
index c0c960158ada6581fb1d3f6f769b31f2fb8ffe87..4e5fb2c4a2b7fe107fe3b99aa80ebf629b45c41d 100644
--- a/include/evalhyd/detail/probabilist/evaluator.hpp
+++ b/include/evalhyd/detail/probabilist/evaluator.hpp
@@ -16,6 +16,7 @@
#include "quantiles.hpp"
#include "contingency.hpp"
#include "ranks.hpp"
+#include "intervals.hpp"
namespace evalhyd
@@ -31,6 +32,7 @@ namespace evalhyd
const XD4& q_prd;
// members for optional input data
const XD2& _q_thr;
+ const xt::xtensor<double, 1>& _c_lvl;
xtl::xoptional<const std::string, bool> _events;
XB4 t_msk;
const std::vector<xt::xkeep_slice<int>>& b_exp;
@@ -45,6 +47,7 @@ namespace evalhyd
std::size_t n_msk;
std::size_t n_mbr;
std::size_t n_thr;
+ std::size_t n_itv;
std::size_t n_exp;
// members for computational elements
@@ -64,6 +67,11 @@ namespace evalhyd
// > Ranks-based
xtl::xoptional<xt::xtensor<double, 3>, bool> r_k;
xtl::xoptional<xt::xtensor<double, 5>, bool> o_j;
+ // > Intervals-based
+ xtl::xoptional<xt::xtensor<double, 5>, bool> itv_bnds;
+ xtl::xoptional<xt::xtensor<double, 4>, bool> obs_in_itv;
+ xtl::xoptional<xt::xtensor<double, 4>, bool> itv_width;
+ xtl::xoptional<xt::xtensor<double, 6>, bool> clim_bnds;
// members for intermediate evaluation metrics
// (i.e. before the reduction along the temporal axis)
@@ -77,6 +85,8 @@ namespace evalhyd
xtl::xoptional<xt::xtensor<double, 5>, bool> pofd;
xtl::xoptional<xt::xtensor<double, 5>, bool> far;
xtl::xoptional<xt::xtensor<double, 5>, bool> csi;
+ // > Intervals-based
+ xtl::xoptional<xt::xtensor<double, 4>, bool> ws;
// members for evaluation metrics
// > Brier-based
@@ -97,6 +107,13 @@ namespace evalhyd
xtl::xoptional<xt::xtensor<double, 5>, bool> RANK_DIAG;
xtl::xoptional<xt::xtensor<double, 4>, bool> DS;
xtl::xoptional<xt::xtensor<double, 4>, bool> AS;
+ // > Intervals-based
+ xtl::xoptional<xt::xtensor<double, 5>, bool> CR;
+ xtl::xoptional<xt::xtensor<double, 5>, bool> AW;
+ xtl::xoptional<xt::xtensor<double, 5>, bool> AWN;
+ xtl::xoptional<xt::xtensor<double, 5>, bool> AWI;
+ xtl::xoptional<xt::xtensor<double, 5>, bool> WS;
+ xtl::xoptional<xt::xtensor<double, 5>, bool> WSS;
// methods to get optional parameters
auto get_q_thr()
@@ -113,6 +130,20 @@ namespace evalhyd
}
}
+ auto get_c_lvl()
+ {
+ if (_c_lvl.size() < 1)
+ {
+ throw std::runtime_error(
+ "interval-based metric requested, "
+ "but *c_lvl* not provided"
+ );
+ }
+ else{
+ return _c_lvl;
+ }
+ }
+
bool is_high_flow_event()
{
if (_events.has_value())
@@ -276,6 +307,53 @@ namespace evalhyd
return o_j.value();
};
+ xt::xtensor<double, 5> get_itv_bnds()
+ {
+ if (!itv_bnds.has_value())
+ {
+ itv_bnds = elements::calc_itv_bnds(
+ q_prd, get_c_lvl(),
+ n_sit, n_ldt, n_itv, n_tim
+ );
+ }
+ return itv_bnds.value();
+ };
+
+ xt::xtensor<double, 4> get_obs_in_itv()
+ {
+ if (!obs_in_itv.has_value())
+ {
+ obs_in_itv = elements::calc_obs_in_itv(
+ q_obs, get_itv_bnds()
+ );
+ }
+ return obs_in_itv.value();
+ };
+
+ xt::xtensor<double, 4> get_itv_width()
+ {
+ if (!itv_width.has_value())
+ {
+ itv_width = elements::calc_itv_width(
+ get_itv_bnds()
+ );
+ }
+ return itv_width.value();
+ };
+
+
+ xt::xtensor<double, 6> get_clim_bnds()
+ {
+ if (!clim_bnds.has_value())
+ {
+ clim_bnds = elements::calc_clim_bnds(
+ q_obs, get_c_lvl(), t_msk, b_exp,
+ n_sit, n_ldt, n_itv, n_msk, n_exp
+ );
+ }
+ return clim_bnds.value();
+ };
+
// methods to compute intermediate metrics
xt::xtensor<double, 4> get_bs()
{
@@ -354,17 +432,30 @@ namespace evalhyd
return csi.value();
};
+ xt::xtensor<double, 4> get_ws()
+ {
+ if (!ws.has_value())
+ {
+ ws = intermediate::calc_ws(
+ q_obs, get_c_lvl(), get_itv_bnds()
+ );
+ }
+ return ws.value();
+ };
+
public:
// constructor method
Evaluator(const XD2& obs,
const XD4& prd,
const XD2& thr,
+ const xt::xtensor<double, 1>& lvl,
xtl::xoptional<const std::string&, bool> events,
const XB4& msk,
const std::vector<xt::xkeep_slice<int>>& exp,
const long int seed) :
q_obs{obs}, q_prd{prd},
- _q_thr{thr}, _events{events}, t_msk(msk), b_exp(exp),
+ _q_thr{thr}, _c_lvl{lvl}, _events{events},
+ t_msk(msk), b_exp(exp),
random_seed{seed}
{
// initialise a mask if none provided
@@ -384,6 +475,7 @@ namespace evalhyd
n_tim = q_prd.shape(3);
n_msk = t_msk.shape(2);
n_thr = _q_thr.shape(1);
+ n_itv = _c_lvl.size();
n_exp = b_exp.size();
// drop time steps where observations and/or predictions are NaN
@@ -583,6 +675,77 @@ namespace evalhyd
}
return AS.value();
};
+
+ xt::xtensor<double, 5> get_CR()
+ {
+ if (!CR.has_value())
+ {
+ CR = metrics::calc_CR(
+ get_obs_in_itv(), t_msk, b_exp,
+ n_sit, n_ldt, n_itv, n_msk, n_exp
+ );
+ }
+ return CR.value();
+ };
+
+ xt::xtensor<double, 5> get_AW()
+ {
+ if (!AW.has_value())
+ {
+ AW = metrics::calc_AW(
+ get_itv_width(), t_msk, b_exp,
+ n_sit, n_ldt, n_itv, n_msk, n_exp
+ );
+ }
+ return AW.value();
+ };
+
+ xt::xtensor<double, 5> get_AWN()
+ {
+ if (!AWN.has_value())
+ {
+ AWN = metrics::calc_AWN(
+ q_obs, get_AW(), t_msk, b_exp,
+ n_sit, n_msk, n_exp
+ );
+ }
+ return AWN.value();
+ };
+
+ xt::xtensor<double, 5> get_AWI()
+ {
+ if (!AWI.has_value())
+ {
+ AWI = metrics::calc_AWI(
+ get_AW(), get_clim_bnds()
+ );
+ }
+ return AWI.value();
+ };
+
+ xt::xtensor<double, 5> get_WS()
+ {
+ if (!WS.has_value())
+ {
+ WS = metrics::calc_WS(
+ get_ws(), t_msk, b_exp,
+ n_sit, n_ldt, n_itv, n_msk, n_exp
+ );
+ }
+ return WS.value();
+ };
+
+ xt::xtensor<double, 5> get_WSS()
+ {
+ if (!WSS.has_value())
+ {
+ WSS = metrics::calc_WSS(
+ q_obs, get_c_lvl(), get_clim_bnds(), get_WS(),
+ n_sit, n_ldt, n_itv, n_msk, n_exp
+ );
+ }
+ return WSS.value();
+ };
};
}
}
diff --git a/include/evalhyd/detail/probabilist/intervals.hpp b/include/evalhyd/detail/probabilist/intervals.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..dbc209b66abbab50f6fa6dcc6b218ba7a519a423
--- /dev/null
+++ b/include/evalhyd/detail/probabilist/intervals.hpp
@@ -0,0 +1,622 @@
+// Copyright (c) 2023, INRAE.
+// Distributed under the terms of the GPL-3 Licence.
+// The full licence is in the file LICENCE, distributed with this software.
+
+#ifndef EVALHYD_PROBABILIST_INTERVALS_HPP
+#define EVALHYD_PROBABILIST_INTERVALS_HPP
+
+#include <xtensor/xtensor.hpp>
+#include <xtensor/xview.hpp>
+#include <xtensor/xsort.hpp>
+
+#include "../maths.hpp"
+
+
+namespace evalhyd
+{
+ namespace probabilist
+ {
+ namespace elements
+ {
+ /// Compute the bounds of the predictive intervals by computing
+ /// the quantiles of the predictive distribution corresponding
+ /// to the confidence intervals.
+ ///
+ /// \param q_prd
+ /// Streamflow predictions.
+ /// shape: (sites, lead times, members, time)
+ /// \param c_lvl
+ /// Confidence levels for the predictive intervals.
+ /// shape: (intervals,)
+ /// \param n_sit
+ /// Number of sites.
+ /// \param n_ldt
+ /// Number of lead times.
+ /// \param n_itv
+ /// Number of predictive intervals.
+ /// \param n_tim
+ /// Number of time steps.
+ /// \return
+ /// Bounds of the predictive intervals.
+ /// shape: (sites, lead times, intervals, bounds, time)
+ template <class XD4>
+ inline xt::xtensor<double, 5> calc_itv_bnds(
+ const XD4& q_prd,
+ const xt::xtensor<double, 1>& c_lvl,
+ std::size_t n_sit,
+ std::size_t n_ldt,
+ std::size_t n_itv,
+ std::size_t n_tim
+ )
+ {
+ xt::xtensor<double, 5> itv_bnds =
+ xt::zeros<double>({n_sit, n_ldt, n_itv, std::size_t {2}, n_tim});
+
+ // determine quantiles forming the predictive intervals
+ // from the confidence levels
+ xt::xtensor<double, 2> quantiles =
+ xt::zeros<double>({n_itv, std::size_t {2}});
+ xt::col(quantiles, 0) = 0.5 - c_lvl / 200.;
+ xt::col(quantiles, 1) = 0.5 + c_lvl / 200.;
+
+ // compute predictive interval bounds from quantiles
+
+ // TODO: replace with `xt::quantile` when available
+ for (std::size_t s = 0; s < n_sit; s++)
+ {
+ for (std::size_t l = 0; l < n_ldt; l++)
+ {
+ for (std::size_t i = 0; i < n_itv; i++)
+ {
+ for (std::size_t t = 0; t < n_tim; t++)
+ {
+ // lower bound
+ xt::view(itv_bnds, s, l, i, 0, t) =
+ evalhyd::maths::quantile(
+ xt::view(q_prd, s, l, xt::all(), t),
+ quantiles(i, 0)
+ );
+ // upper bound
+ xt::view(itv_bnds, s, l, i, 1, t) =
+ evalhyd::maths::quantile(
+ xt::view(q_prd, s, l, xt::all(), t),
+ quantiles(i, 1)
+ );
+ }
+ }
+ }
+ }
+
+ return itv_bnds;
+ }
+
+ /// Determine whether the observations are inside the predictive
+ /// intervals for each time step.
+ ///
+ /// \param q_obs
+ /// Streamflow predictions.
+ /// shape: (sites, time)
+ /// \param itv_bnds
+ /// Bounds of the predictive intervals.
+ /// shape: (sites, lead times, intervals, bounds, time)
+ /// \return
+ /// Boolean-like tensor evaluating to true where observations
+ /// are inside the predictive intervals.
+ /// shape: (sites, lead times, intervals, time)
+ template <class XD2>
+ inline xt::xtensor<double, 4> calc_obs_in_itv(
+ const XD2& q_obs,
+ const xt::xtensor<double, 5>& itv_bnds
+ )
+ {
+ // notations below follow Gneiting and Raftery (2007), sect 6.2
+ // https://doi.org/10.1198/016214506000001437
+
+ auto x = xt::view(q_obs, xt::all(), xt::newaxis(), xt::newaxis(), xt::all());
+ auto l = xt::view(itv_bnds, xt::all(), xt::all(), xt::all(), 0, xt::all());
+ auto u = xt::view(itv_bnds, xt::all(), xt::all(), xt::all(), 1, xt::all());
+
+ return ((x >= l) && (x <= u));
+ }
+
+ /// Compute the width of the predictive intervals for each time step.
+ ///
+ /// \param itv_bnds
+ /// Bounds of the predictive intervals.
+ /// shape: (sites, lead times, intervals, bounds, time)
+ /// \return
+ /// Interval scores for each time step.
+ /// shape: (sites, lead times, intervals, time)
+ inline xt::xtensor<double, 4> calc_itv_width(
+ const xt::xtensor<double, 5>& itv_bnds
+ )
+ {
+ // notations below follow Gneiting and Raftery (2007), sect 6.2
+ // https://doi.org/10.1198/016214506000001437
+
+ auto l = xt::view(itv_bnds, xt::all(), xt::all(), xt::all(), 0, xt::all());
+ auto u = xt::view(itv_bnds, xt::all(), xt::all(), xt::all(), 1, xt::all());
+
+ return (u - l);
+ }
+
+ /// Compute the bounds of the climatology corresponding to the
+ /// confidence levels.
+ ///
+ /// \param q_obs
+ /// Streamflow predictions.
+ /// shape: (sites, time)
+ /// \param c_lvl
+ /// Confidence levels for the predictive intervals.
+ /// shape: (intervals,)
+ /// \param t_msk
+ /// Temporal subsets of the whole record.
+ /// shape: (sites, lead times, subsets, time)
+ /// \param b_exp
+ /// Boostrap samples.
+ /// shape: (samples, time slice)
+ /// \param n_sit
+ /// Number of sites.
+ /// \param n_ldt
+ /// Number of lead times.
+ /// \param n_itv
+ /// Number of predictive intervals.
+ /// \param n_msk
+ /// Number of temporal subsets.
+ /// \param n_exp
+ /// Number of bootstrap samples.
+ /// \return
+ /// Climatology bounds.
+ /// shape: (sites, lead times, subsets, samples, intervals, bounds)
+ template <class XD2>
+ inline xt::xtensor<double, 6> calc_clim_bnds(
+ const XD2& q_obs,
+ const xt::xtensor<double, 1>& c_lvl,
+ const xt::xtensor<bool, 4>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_sit,
+ std::size_t n_ldt,
+ std::size_t n_itv,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // compute "climatology" average width
+ xt::xtensor<double, 6> clim_bnds =
+ xt::zeros<double>({n_sit, n_ldt, n_msk, n_exp,
+ n_itv, std::size_t {2}});
+
+ // determine quantiles forming the predictive intervals
+ // from the confidence levels
+ xt::xtensor<double, 2> quantiles =
+ xt::zeros<double>({n_itv, std::size_t {2}});
+ xt::col(quantiles, 0) = 0.5 - c_lvl / 200.;
+ xt::col(quantiles, 1) = 0.5 + c_lvl / 200.;
+
+ // 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 q_obs_masked = xt::where(
+ xt::view(t_msk, xt::all(), xt::all(), m, xt::all()),
+ xt::view(q_obs, xt::all(), xt::newaxis(), xt::all()),
+ NAN
+ );
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // apply the bootstrap sampling
+ auto q_obs_masked_sampled =
+ xt::view(q_obs_masked, xt::all(), xt::all(), b_exp[e]);
+
+ // compute "climatology" interval
+
+ // TODO: replace with `xt::quantile` when available
+ for (std::size_t s = 0; s < n_sit; s++)
+ {
+ for (std::size_t l = 0; l < n_ldt; l++)
+ {
+ for (std::size_t i = 0; i < n_itv; i++)
+ {
+ // lower bound
+ xt::view(clim_bnds, s, l, m, e, i, 0) =
+ evalhyd::maths::quantile(
+ xt::view(q_obs_masked_sampled,
+ s, l, xt::all()),
+ quantiles(i, 0)
+ );
+ // upper bound
+ xt::view(clim_bnds, s, l, m, e, i, 1) =
+ evalhyd::maths::quantile(
+ xt::view(q_obs_masked_sampled,
+ s, l, xt::all()),
+ quantiles(i, 1)
+ );
+
+ }
+ }
+ }
+ }
+ }
+
+ return clim_bnds;
+ }
+ }
+
+ namespace intermediate
+ {
+ /// Compute the Winkler score for each time step.
+ ///
+ /// \param q_obs
+ /// Streamflow predictions.
+ /// shape: (sites, time)
+ /// \param c_lvl
+ /// Confidence levels for the predictive intervals.
+ /// shape: (intervals,)
+ /// \param itv_bnds
+ /// Bounds of the predictive intervals.
+ /// shape: (sites, lead times, intervals, bounds, time)
+ /// \return
+ /// Interval scores for each time step.
+ /// shape: (sites, lead times, intervals, time)
+ template <class XD2>
+ inline xt::xtensor<double, 4> calc_ws(
+ const XD2& q_obs,
+ const xt::xtensor<double, 1>& c_lvl,
+ const xt::xtensor<double, 5>& itv_bnds
+ )
+ {
+ // notations below follow Gneiting and Raftery (2007), sect 6.2
+ // https://doi.org/10.1198/016214506000001437
+
+ auto x = xt::view(q_obs, xt::all(), xt::newaxis(), xt::newaxis(), xt::all());
+ auto alpha = 1 - xt::view(c_lvl, xt::all(), xt::newaxis()) / 100.;
+
+ // compute component corresponding to observations below interval
+ auto l = xt::view(itv_bnds, xt::all(), xt::all(), xt::all(), 0, xt::all());
+ // (l - x)ðŸ™{x < l}
+ auto ws_l = xt::where(x < l, l - x, 0.);
+
+ // compute component corresponding to observations above interval
+ auto u = xt::view(itv_bnds, xt::all(), xt::all(), xt::all(), 1, xt::all());
+ // (x - u)ðŸ™{x > u}
+ auto ws_u = xt::where(x > u, x - u, 0.);
+
+ // compute interval score
+ auto ws = (u - l) + 2. * (ws_l + ws_u) / alpha;
+
+ return ws;
+ }
+ }
+
+ namespace metrics
+ {
+ namespace detail {
+ inline xt::xtensor<double, 5> calc_METRIC_from_metric(
+ const xt::xtensor<double, 4>& metric,
+ const xt::xtensor<bool, 4>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_sit,
+ std::size_t n_ldt,
+ std::size_t n_itv,
+ std::size_t n_msk,
+ std::size_t n_exp
+ ) {
+ // initialise output variable
+ xt::xtensor<double, 5> METRIC =
+ xt::zeros<double>({n_sit, n_ldt, n_msk, n_exp, n_itv});
+
+ // 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 metric_masked = xt::where(
+ xt::view(t_msk, xt::all(), xt::all(), m,
+ xt::newaxis(), xt::all()),
+ metric,
+ NAN
+ );
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // apply the bootstrap sampling
+ auto metric_masked_sampled =
+ xt::view(metric_masked, xt::all(),
+ xt::all(), xt::all(), b_exp[e]);
+
+ // calculate the mean over the time steps
+ xt::view(METRIC, xt::all(), xt::all(), m, e, xt::all()) =
+ xt::nanmean(metric_masked_sampled, -1);
+ }
+ }
+
+ return METRIC;
+ }
+ }
+
+ /// Compute the Coverage Ratio (CR), i.e. the portion of
+ /// observations falling within the predictive intervals.
+ /// It is a measure of the reliability of the predictions.
+ ///
+ /// \param obs_in_itv
+ /// Boolean-like tensor evaluating to true where observations
+ /// are inside the predictive intervals.
+ /// shape: (sites, lead times, intervals, time)
+ /// \param t_msk
+ /// Temporal subsets of the whole record.
+ /// shape: (sites, lead times, subsets, time)
+ /// \param b_exp
+ /// Boostrap samples.
+ /// shape: (samples, time slice)
+ /// \param n_sit
+ /// Number of sites.
+ /// \param n_ldt
+ /// Number of lead times.
+ /// \param n_itv
+ /// Number of predictive intervals.
+ /// \param n_msk
+ /// Number of temporal subsets.
+ /// \param n_exp
+ /// Number of bootstrap samples.
+ /// \return
+ /// Coverage ratios.
+ /// shape: (sites, lead times, subsets, samples, intervals)
+ inline xt::xtensor<double, 5> calc_CR(
+ const xt::xtensor<double, 4>& obs_in_itv,
+ const xt::xtensor<bool, 4>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_sit,
+ std::size_t n_ldt,
+ std::size_t n_itv,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ return detail::calc_METRIC_from_metric(
+ obs_in_itv, t_msk, b_exp,
+ n_sit, n_ldt, n_itv, n_msk, n_exp
+ );
+ }
+
+ /// Compute the Average Width (AW) of the predictive intervals.
+ /// It is a measure of the sharpness of the predictions.
+ ///
+ /// \param itv_width
+ /// Widths of predictive intervals for each time step.
+ /// shape: (sites, lead times, intervals, time)
+ /// \param t_msk
+ /// Temporal subsets of the whole record.
+ /// shape: (sites, lead times, subsets, time)
+ /// \param b_exp
+ /// Boostrap samples.
+ /// shape: (samples, time slice)
+ /// \param n_sit
+ /// Number of sites.
+ /// \param n_ldt
+ /// Number of lead times.
+ /// \param n_itv
+ /// Number of predictive intervals.
+ /// \param n_msk
+ /// Number of temporal subsets.
+ /// \param n_exp
+ /// Number of bootstrap samples.
+ /// \return
+ /// Average widths.
+ /// shape: (sites, lead times, subsets, samples, intervals)
+ inline xt::xtensor<double, 5> calc_AW(
+ const xt::xtensor<double, 4>& itv_width,
+ const xt::xtensor<bool, 4>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_sit,
+ std::size_t n_ldt,
+ std::size_t n_itv,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ return detail::calc_METRIC_from_metric(
+ itv_width, t_msk, b_exp, n_sit, n_ldt, n_itv, n_msk, n_exp
+ );
+ }
+
+ /// Compute the Average Width Normalised (AWN).
+ ///
+ /// \param q_obs
+ /// Streamflow predictions.
+ /// shape: (sites, time)
+ /// \param AW
+ /// Average widths.
+ /// shape: (sites, lead times, subsets, samples, intervals)
+ /// \param t_msk
+ /// Temporal subsets of the whole record.
+ /// shape: (sites, lead times, subsets, time)
+ /// \param b_exp
+ /// Boostrap samples.
+ /// shape: (samples, time slice)
+ /// \param n_sit
+ /// Number of sites.
+ /// \param n_msk
+ /// Number of temporal subsets.
+ /// \param n_exp
+ /// Number of bootstrap samples.
+ /// \return
+ /// Average widths normalised with mean observations.
+ /// shape: (sites, lead times, subsets, samples, intervals)
+ template <class XD2>
+ inline xt::xtensor<double, 5> calc_AWN(
+ const XD2& q_obs,
+ const xt::xtensor<double, 5>& AW,
+ const xt::xtensor<bool, 4>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_sit,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // compute "climatology" average width
+ xt::xtensor<double, 5> mean_obs =
+ xt::zeros<double>({n_sit, std::size_t {1}, n_msk,
+ n_exp, std::size_t {1}});
+
+ // 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 q_obs_masked = xt::where(
+ xt::view(t_msk, xt::all(), xt::all(), m, xt::all()),
+ xt::view(q_obs, xt::all(), xt::newaxis(), xt::all()),
+ NAN
+ );
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // apply the bootstrap sampling
+ auto q_obs_masked_sampled =
+ xt::view(q_obs_masked, xt::all(), xt::all(), b_exp[e]);
+
+ // compute mean observation
+ xt::view(mean_obs, xt::all(), xt::all(), m, e, 0) =
+ xt::nanmean(q_obs_masked_sampled, -1);
+ }
+ }
+
+ return (AW / mean_obs);
+ }
+
+ /// Compute the Average Width Index (AWI).
+ ///
+ /// \param AW
+ /// Average widths.
+ /// shape: (sites, lead times, subsets, samples, intervals)
+ /// \param clim_bnds
+ /// Climatology bounds.
+ /// shape: (sites, lead times, subsets, samples, intervals, bounds)
+ /// \return
+ /// Average width indices.
+ /// shape: (sites, lead times, subsets, samples, intervals)
+ inline xt::xtensor<double, 5> calc_AWI(
+ const xt::xtensor<double, 5>& AW,
+ const xt::xtensor<double, 6>& clim_bnds
+ )
+ {
+ // compute "climatology" average width
+ auto AW_clim =
+ xt::view(clim_bnds, xt::all(), xt::all(), xt::all(),
+ xt::all(), xt::all(), 1)
+ - xt::view(clim_bnds, xt::all(), xt::all(), xt::all(),
+ xt::all(), xt::all(), 0);
+
+ return 1 - (AW / AW_clim);
+ }
+
+ /// Compute the Winkler scores (WS), also known as interval score.
+ ///
+ /// \param ws
+ /// Winkler scores for each time step.
+ /// shape: (sites, lead times, intervals, time)
+ /// \param t_msk
+ /// Temporal subsets of the whole record.
+ /// shape: (sites, lead times, subsets, time)
+ /// \param b_exp
+ /// Boostrap samples.
+ /// shape: (samples, time slice)
+ /// \param n_sit
+ /// Number of sites.
+ /// \param n_ldt
+ /// Number of lead times.
+ /// \param n_itv
+ /// Number of predictive intervals.
+ /// \param n_msk
+ /// Number of temporal subsets.
+ /// \param n_exp
+ /// Number of bootstrap samples.
+ /// \return
+ /// Winkler scores.
+ /// shape: (sites, lead times, subsets, samples, intervals)
+ inline xt::xtensor<double, 5> calc_WS(
+ const xt::xtensor<double, 4>& ws,
+ const xt::xtensor<bool, 4>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_sit,
+ std::size_t n_ldt,
+ std::size_t n_itv,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ return detail::calc_METRIC_from_metric(
+ ws, t_msk, b_exp, n_sit, n_ldt, n_itv, n_msk, n_exp
+ );
+ }
+
+ /// Compute the Winkler skill scores (WSS).
+ ///
+ /// \param q_obs
+ /// Streamflow predictions.
+ /// shape: (sites, time)
+ /// \param c_lvl
+ /// Confidence levels for the predictive intervals.
+ /// shape: (intervals,)
+ /// \param clim_bnds
+ /// Climatology bounds.
+ /// shape: (sites, lead times, subsets, samples, intervals, bounds)
+ /// \param WS
+ /// Winkler scores.
+ /// shape: (sites, lead times, subsets, samples, intervals)
+ /// \param n_sit
+ /// Number of sites.
+ /// \param n_ldt
+ /// Number of lead times.
+ /// \param n_itv
+ /// Number of predictive intervals.
+ /// \param n_msk
+ /// Number of temporal subsets.
+ /// \param n_exp
+ /// Number of bootstrap samples.
+ /// \return
+ /// Winkler skill scores.
+ /// shape: (sites, lead times, subsets, samples, intervals)
+ template <class XD2>
+ inline xt::xtensor<double, 5> calc_WSS(
+ const XD2& q_obs,
+ const xt::xtensor<double, 1>& c_lvl,
+ const xt::xtensor<double, 6>& clim_bnds,
+ const xt::xtensor<double, 5>& WS,
+ std::size_t n_sit,
+ std::size_t n_ldt,
+ std::size_t n_itv,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // compute "climatology" Winkler score
+ xt::xtensor<double, 5> WS_clim =
+ xt::zeros<double>({n_sit, n_ldt, n_msk, n_exp, n_itv});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ auto ws_clim = intermediate::calc_ws(
+ q_obs, c_lvl,
+ xt::view(clim_bnds, xt::all(), xt::all(), m, e,
+ xt::all(), xt::all(), xt::newaxis())
+ );
+
+ xt::view(WS_clim, xt::all(), xt::all(), m, e, xt::all()) =
+ xt::nanmean(ws_clim, -1);
+ }
+ }
+
+ return 1 - (WS / WS_clim);
+ }
+ }
+ }
+}
+
+#endif //EVALHYD_PROBABILIST_INTERVALS_HPP
diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp
index 59977c80086e021b1714fff620ad7e84928c56aa..0aff6f6392a1254a455d0f6203d78119801ece49 100644
--- a/include/evalhyd/evalp.hpp
+++ b/include/evalhyd/evalp.hpp
@@ -12,6 +12,7 @@
#include <xtensor/xexpression.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xarray.hpp>
+#include <xtensor/xadapt.hpp>
#include "detail/utils.hpp"
#include "detail/masks.hpp"
@@ -73,6 +74,9 @@ namespace evalhyd
/// occurring when streamflow goes below threshold). It must be
/// provided if *q_thr* is provided.
///
+ /// c_lvl: ``std::vector<double>``, optional
+ /// Confidence interval(s).
+ ///
/// t_msk: ``XB4``, optional
/// Mask(s) used to generate temporal subsets of the whole streamflow
/// time series (where True/False is used for the time steps to
@@ -162,6 +166,7 @@ namespace evalhyd
const xt::xexpression<XD2>& q_thr = XD2({}),
xtl::xoptional<const std::string, bool> events =
xtl::missing<const std::string>(),
+ const std::vector<double>& c_lvl = {},
const xt::xexpression<XB4>& t_msk = XB4({}),
const xt::xtensor<std::array<char, 32>, 2>& m_cdt = {},
xtl::xoptional<const std::unordered_map<std::string, int>, bool> bootstrap =
@@ -198,13 +203,17 @@ namespace evalhyd
const XB4& t_msk_ = t_msk.derived_cast();
+ // adapt vector to tensor
+ const xt::xtensor<double, 1> c_lvl_ = xt::adapt(c_lvl);
+
// check that the metrics to be computed are valid
utils::check_metrics(
metrics,
{"BS", "BSS", "BS_CRD", "BS_LBD",
"QS", "CRPS",
"POD", "POFD", "FAR", "CSI", "ROCSS",
- "RANK_DIAG", "DS", "AS"}
+ "RANK_DIAG", "DS", "AS",
+ "CR", "AW", "AWN", "AWI", "WS", "WSS"}
);
// check optional parameters
@@ -371,7 +380,7 @@ namespace evalhyd
// instantiate determinist evaluator
probabilist::Evaluator<XD2, XD4, XB4> evaluator(
- q_obs_, q_prd_, q_thr_, events,
+ q_obs_, q_prd_, q_thr_, c_lvl_, events,
t_msk_.size() > 0 ? t_msk_: (m_cdt.size() > 0 ? c_msk : t_msk_),
b_exp,
random_seed
@@ -466,6 +475,42 @@ namespace evalhyd
uncertainty::summarise(evaluator.get_AS(), summary)
);
}
+ else if ( metric == "CR" )
+ {
+ r.emplace_back(
+ uncertainty::summarise(evaluator.get_CR(), summary)
+ );
+ }
+ else if ( metric == "AW" )
+ {
+ r.emplace_back(
+ uncertainty::summarise(evaluator.get_AW(), summary)
+ );
+ }
+ else if ( metric == "AWN" )
+ {
+ r.emplace_back(
+ uncertainty::summarise(evaluator.get_AWN(), summary)
+ );
+ }
+ else if ( metric == "AWI" )
+ {
+ r.emplace_back(
+ uncertainty::summarise(evaluator.get_AWI(), summary)
+ );
+ }
+ else if ( metric == "WS" )
+ {
+ r.emplace_back(
+ uncertainty::summarise(evaluator.get_WS(), summary)
+ );
+ }
+ else if ( metric == "WSS" )
+ {
+ r.emplace_back(
+ uncertainty::summarise(evaluator.get_WSS(), summary)
+ );
+ }
}
return r;