diff --git a/changelog.rst b/changelog.rst
index 2716fa599127e48f29b5ee884e0a2afab7c9e327..0ba34c1d6cc32c446cfbfd8ac9fd13111ac34235 100644
--- a/changelog.rst
+++ b/changelog.rst
@@ -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)
diff --git a/include/evalhyd/detail/probabilist/evaluator.hpp b/include/evalhyd/detail/probabilist/evaluator.hpp
index 6293ea29e546df8f9f3756608d213aa4e5773a9e..976145c58817f529c2786a246665c7eac55c5eb9 100644
--- a/include/evalhyd/detail/probabilist/evaluator.hpp
+++ b/include/evalhyd/detail/probabilist/evaluator.hpp
@@ -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())
diff --git a/include/evalhyd/detail/probabilist/intervals.hpp b/include/evalhyd/detail/probabilist/intervals.hpp
index c36a235838ead80eec395b4ccd342684217b5afc..adbd487b592ab4903210072e9ae01836ba41e0ac 100644
--- a/include/evalhyd/detail/probabilist/intervals.hpp
+++ b/include/evalhyd/detail/probabilist/intervals.hpp
@@ -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());
- }
}
}
}
diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp
index c79f47de421f3aa93084800d0fe8de581fdaa83d..c04412097f925f675cb5750037086c421dceda9e 100644
--- a/include/evalhyd/evalp.hpp
+++ b/include/evalhyd/evalp.hpp
@@ -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(
diff --git a/tests/expected/evalp/AWI.csv b/tests/expected/evalp/AWI.csv
deleted file mode 100644
index 8b3b7e7f9598b23911dc06f440d71572b99b48e5..0000000000000000000000000000000000000000
--- a/tests/expected/evalp/AWI.csv
+++ /dev/null
@@ -1 +0,0 @@
-0.9821120161733,0.9880951944476
diff --git a/tests/expected/evalp/WSS.csv b/tests/expected/evalp/WSS.csv
deleted file mode 100644
index fd929dabd8f58354be587455dcbd4402d6a9f83d..0000000000000000000000000000000000000000
--- a/tests/expected/evalp/WSS.csv
+++ /dev/null
@@ -1 +0,0 @@
-0.6621887740287,0.4360388849930
diff --git a/tests/test_probabilist.cpp b/tests/test_probabilist.cpp
index 65b16cdde5c4e0b30280bed9a96a8cfea6f36ba9..f4e53f71cc694abffe77c5d5287af2fa6124c4c9 100644
--- a/tests/test_probabilist.cpp
+++ b/tests/test_probabilist.cpp
@@ -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