Commit 0cedfa7c authored by Thibault Hallouin's avatar Thibault Hallouin
Browse files

add intervals-based metrics (CR, AW, AWN, AWI, WS, WSS)

1 merge request!3release v0.1.0
Pipeline #43737 failed with stage
in 1 minute and 50 seconds
Showing with 833 additions and 3 deletions
+833 -3
......@@ -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();
};
};
}
}
......
// 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
......@@ -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;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment