From 0cedfa7c64f31c55a44ef2021d0d40d1017f9456 Mon Sep 17 00:00:00 2001 From: Thibault Hallouin <thibault.hallouin@inrae.fr> Date: Thu, 26 Jan 2023 16:30:57 +0100 Subject: [PATCH] add intervals-based metrics (CR, AW, AWN, AWI, WS, WSS) --- .../evalhyd/detail/probabilist/evaluator.hpp | 165 ++++- .../evalhyd/detail/probabilist/intervals.hpp | 622 ++++++++++++++++++ include/evalhyd/evalp.hpp | 49 +- 3 files changed, 833 insertions(+), 3 deletions(-) create mode 100644 include/evalhyd/detail/probabilist/intervals.hpp diff --git a/include/evalhyd/detail/probabilist/evaluator.hpp b/include/evalhyd/detail/probabilist/evaluator.hpp index c0c9601..4e5fb2c 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 0000000..dbc209b --- /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 59977c8..0aff6f6 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; -- GitLab