Skip to content
GitLab
Commit f9c98efd authored by Thibault Hallouin's avatar Thibault Hallouin
Browse files

remove WSS and AWI 

resolves #7
parent 2915b9aa
main 11-check-the-computation-of-the-ranks-for-ties dev v0.1.2 v0.1.1
1 merge request!6release v0.1.1
Pipeline #47886 passed with stage
in 3 minutes and 43 seconds
Hide whitespace changes
Inline Side-by-side
Showing
with 10 additions and 287 deletions
+10 -287
changelog.rst
+ 6
− 0
View file @ f9c98efd
......@@ -5,6 +5,12 @@ latest
Yet to be versioned and released. Only available from *dev* branch until then.
.. rubric:: Scope changes
* remove `"WSS"` and `"AWI"` as probabilistic evaluation metric because of the
arbitrary nature of using the sample climatology as reference/benchmark
(`#7-CPP <https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-cpp/-/issues/7>`_)
.. rubric:: Bug fixes
* fix irrelevant rank check failure when passing arrays (and not tensors)
......
include/evalhyd/detail/probabilist/evaluator.hpp
+ 0
− 39
View file @ f9c98efd
......@@ -76,7 +76,6 @@ namespace evalhyd
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)
......@@ -124,9 +123,7 @@ namespace evalhyd
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;
// > Multi-variate
xtl::xoptional<xt::xtensor<double, 4>, bool> ES;
......@@ -367,19 +364,6 @@ namespace evalhyd
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()
{
......@@ -798,17 +782,6 @@ namespace evalhyd
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())
......@@ -821,18 +794,6 @@ namespace evalhyd
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(),
t_msk, b_exp, n_sit, n_ldt, n_itv, n_msk, n_exp
);
}
return WSS.value();
};
xt::xtensor<double, 4> get_ES()
{
if (!ES.has_value())
......
include/evalhyd/detail/probabilist/intervals.hpp
+ 0
− 230
View file @ f9c98efd
......@@ -123,119 +123,6 @@ namespace evalhyd
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
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++)
{
auto obs = xt::view(q_obs_masked_sampled,
s, l, xt::all());
auto obs_filtered =
xt::filter(obs, !xt::isnan(obs));
if (obs_filtered.size() > 0)
{
// lower bound
xt::view(clim_bnds, s, l, m, e, i, 0) =
xt::quantile(
obs_filtered,
{quantiles(i, 0)}
)();
// upper bound
xt::view(clim_bnds, s, l, m, e, i, 1) =
xt::quantile(
obs_filtered,
{quantiles(i, 1)}
)();
}
else
{
xt::view(clim_bnds, s, l, m, e, i, xt::all()) =
NAN;
}
}
}
}
}
}
return clim_bnds;
}
}
namespace intermediate
......@@ -487,34 +374,6 @@ namespace evalhyd
- std::numeric_limits<double>::infinity());
}
/// 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 xt::where(AW_clim > 0,
1 - (AW / AW_clim),
- std::numeric_limits<double>::infinity());
}
/// Compute the Winkler scores (WS), also known as interval score.
///
/// \param ws
......@@ -554,95 +413,6 @@ namespace evalhyd
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 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 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,
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" 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 l = 0; l < n_ldt; l++)
{
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(), l, xt::newaxis(),
m, e, xt::all(), xt::all(),
xt::newaxis())
);
// apply the mask
// (using NaN workaround until reducers work on masked_view)
auto ws_clim_masked = xt::where(
xt::view(t_msk, xt::all(), l, m, xt::newaxis(), xt::all()),
xt::view(ws_clim, xt::all(), l, xt::all(), xt::all()),
NAN
);
// apply the bootstrap sampling
auto ws_clim_masked_sampled =
xt::view(ws_clim_masked, xt::all(), xt::all(), b_exp[e]);
xt::view(WS_clim, xt::all(), l, m, e, xt::all()) =
xt::nanmean(ws_clim_masked_sampled, -1);
}
}
}
// compute the Winkler skill score
return xt::where(WS_clim > 0,
1 - (WS / WS_clim),
- std::numeric_limits<double>::infinity());
}
}
}
}
......
include/evalhyd/evalp.hpp
+ 1
− 13
View file @ f9c98efd
......@@ -240,7 +240,7 @@ namespace evalhyd
"QS", "CRPS_FROM_QS",
"POD", "POFD", "FAR", "CSI", "ROCSS",
"RANK_HIST", "DS", "AS",
"CR", "AW", "AWN", "AWI", "WS", "WSS",
"CR", "AW", "AWN", "WS",
"ES"}
);
......@@ -547,24 +547,12 @@ namespace evalhyd
uncertainty::summarise_p(evaluator.get_AWN(), summary)
);
}
else if ( metric == "AWI" )
{
r.emplace_back(
uncertainty::summarise_p(evaluator.get_AWI(), summary)
);
}
else if ( metric == "WS" )
{
r.emplace_back(
uncertainty::summarise_p(evaluator.get_WS(), summary)
);
}
else if ( metric == "WSS" )
{
r.emplace_back(
uncertainty::summarise_p(evaluator.get_WSS(), summary)
);
}
else if ( metric == "ES" )
{
r.emplace_back(
......
tests/expected/evalp/AWI.csv deleted 100644 → 0
+ 0
− 1
View file @ 2915b9aa
0.9821120161733,0.9880951944476
tests/expected/evalp/WSS.csv deleted 100644 → 0
+ 0
− 1
View file @ 2915b9aa
0.6621887740287,0.4360388849930
tests/test_probabilist.cpp
+ 3
− 3
View file @ f9c98efd
......@@ -35,7 +35,7 @@ std::vector<std::string> all_metrics_p = {
"QS", "CRPS_FROM_QS",
"POD", "POFD", "FAR", "CSI", "ROCSS",
"RANK_HIST", "DS", "AS",
"CR", "AW", "AWN", "AWI", "WS", "WSS",
"CR", "AW", "AWN", "WS",
"ES"
};
......@@ -263,7 +263,7 @@ TEST(ProbabilistTests, TestIntervals)
auto expected = load_expected_p();
// compute scores
std::vector<std::string> metrics = {"CR", "AW", "AWN", "AWI", "WS", "WSS"};
std::vector<std::string> metrics = {"CR", "AW", "AWN", "WS"};
std::vector<xt::xarray<double>> results =
evalhyd::evalp(
......@@ -271,7 +271,7 @@ TEST(ProbabilistTests, TestIntervals)
xt::eval(xt::view(observed, xt::newaxis(), xt::all())),
// shape: (sites [1], lead times [1], members [m], time [t])
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
{"CR", "AW", "AWN", "AWI", "WS", "WSS"},
{"CR", "AW", "AWN", "WS"},
xt::xtensor<double, 2>({}),
"", // events
{30., 80.} // c_lvl
......
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