Skip to content
GitLab

Resolve "evaluation of quantiles prediction"

Merged François Bourgin requested to merge 12-evaluation-of-quantiles-prediction into dev
Compare
and
6 files
+ 94
− 10
include/evalhyd/detail/probabilist/evaluator.hpp
+ 5
− 3
@@ -36,6 +36,7 @@ namespace evalhyd
// members for optional input data
const XD2& _q_thr;
const xt::xtensor<double, 1>& _c_lvl;
const xt::xtensor<double, 1>& _q_lvl;
xtl::xoptional<const std::string, bool> _events;
xt::xtensor<bool, 4> t_msk;
const std::vector<xt::xkeep_slice<int>>& b_exp;
@@ -336,7 +337,7 @@ namespace evalhyd
if (!itv_bnds.has_value())
{
itv_bnds = elements::calc_itv_bnds(
q_prd, get_c_lvl(),
q_prd, get_c_lvl(), _q_lvl,
n_sit, n_ldt, n_itv, n_tim
);
}
@@ -393,7 +394,7 @@ namespace evalhyd
if (!qs.has_value())
{
qs = intermediate::calc_qs(
q_obs, get_q_qnt(), n_mbr
q_obs, get_q_qnt(), n_mbr, _q_lvl
);
}
return qs.value();
@@ -482,12 +483,13 @@ namespace evalhyd
const XD4& prd,
const XD2& thr,
const xt::xtensor<double, 1>& lvl,
const xt::xtensor<double, 1>& qlvl,
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}, _c_lvl{lvl}, _events{events},
_q_thr{thr}, _c_lvl{lvl}, _q_lvl{qlvl}, _events{events},
t_msk(msk), b_exp(exp),
random_seed{seed}
{
include/evalhyd/detail/probabilist/intervals.hpp
+ 28
− 2
@@ -44,6 +44,7 @@ namespace evalhyd
inline xt::xtensor<double, 5> calc_itv_bnds(
const XD4& q_prd,
const xt::xtensor<double, 1>& c_lvl,
const xt::xtensor<double, 1>& q_lvl,
std::size_t n_sit,
std::size_t n_ldt,
std::size_t n_itv,
@@ -60,15 +61,40 @@ namespace evalhyd
xt::col(quantiles, 0) = 0.5 - c_lvl / 200.;
xt::col(quantiles, 1) = 0.5 + c_lvl / 200.;
// compute predictive interval bounds from quantiles
for (std::size_t i = 0; i < n_itv; i++)
// get or compute predictive interval bounds from quantiles
if (q_lvl.size() > 0)
{
for (std::size_t i = 0; i < n_itv; i++)
{
auto a = xt::broadcast(xt::view(quantiles, i), std::vector<std::size_t>({q_lvl.size(), 2}));
auto b = xt::broadcast(q_lvl, std::vector<std::size_t>({2, q_lvl.size()}));
auto res = xt::where(xt::equal(a, xt::transpose(b)));
if (res.size() != 2)
{
throw std::runtime_error(
"interval-based metric requested, "
"but *c_lvl* not matching *q_lvl"
);
} else
{
xt::view(itv_bnds, xt::all(), xt::all(), i, 0, xt::all()) =
xt::view(q_prd, std::min(res[0][1], res[1][1]));
xt::view(itv_bnds, xt::all(), xt::all(), i, 1, xt::all()) =
xt::view(q_prd, std::max(res[0][1], res[1][1]));
}
}
}
else
{
for (std::size_t i = 0; i < n_itv; i++)
{
auto q = xt::quantile(q_prd, xt::view(quantiles, i), 2);
xt::view(itv_bnds, xt::all(), xt::all(), i, 0, xt::all()) =
xt::view(q, 0);
xt::view(itv_bnds, xt::all(), xt::all(), i, 1, xt::all()) =
xt::view(q, 1);
}
}
return itv_bnds;
include/evalhyd/detail/probabilist/quantiles.hpp
+ 12
− 3
@@ -60,13 +60,22 @@ namespace evalhyd
inline xt::xtensor<double, 4> calc_qs(
const XD2 &q_obs,
const xt::xtensor<double, 4>& q_qnt,
std::size_t n_mbr
std::size_t n_mbr,
const xt::xtensor<double, 1>& q_lvl
)
{
// compute the quantile order $alpha$
xt::xtensor<double, 1> alpha =
xt::xtensor<double, 1> alpha;
// get or compute the quantile order $alpha$
if (q_lvl.size() > 0)
{
alpha = xt::sort(q_lvl);
}
else
{
alpha =
xt::arange<double>(1., double(n_mbr + 1))
/ double(n_mbr + 1);
}
// calculate the difference
xt::xtensor<double, 4> diff =
include/evalhyd/evalp.hpp
+ 3
− 1
@@ -181,6 +181,7 @@ namespace evalhyd
xtl::xoptional<const std::string, bool> events =
xtl::missing<const std::string>(),
const std::vector<double>& c_lvl = {},
const std::vector<double>& q_lvl = {},
const xt::xexpression<XB4>& t_msk = XB4({}),
const xt::xexpression<XS2>& m_cdt = XS2({}),
xtl::xoptional<const std::unordered_map<std::string, int>, bool> bootstrap =
@@ -233,6 +234,7 @@ namespace evalhyd
// adapt vector to tensor
const xt::xtensor<double, 1> c_lvl_ = xt::adapt(c_lvl);
const xt::xtensor<double, 1> q_lvl_ = xt::adapt(q_lvl);
// check that the metrics/diagnostics to be computed are valid
utils::check_metrics(
@@ -418,7 +420,7 @@ namespace evalhyd
// instantiate determinist evaluator
probabilist::Evaluator<XD2, XD4, XB4> evaluator(
q_obs_, q_prd_, q_thr_, c_lvl_, events,
q_obs_, q_prd_, q_thr_, c_lvl_, q_lvl_, events,
t_msk_.size() > 0 ? t_msk_: (m_cdt_.size() > 0 ? c_msk : t_msk_),
b_exp,
random_seed
tests/test_probabilist.cpp
+ 45
− 0
@@ -246,6 +246,7 @@ TEST(ProbabilistTests, TestRanks)
xt::xtensor<double, 2>({}),
"high", // events
{}, // c_lvl
{}, // q_lvl
xt::xtensor<bool, 4>({}), // t_msk
xt::xtensor<std::array<char, 32>, 2>({}), // m_cdt
xtl::missing<const std::unordered_map<std::string, int>>(), // bootstrap
@@ -296,6 +297,41 @@ TEST(ProbabilistTests, TestIntervals)
}
}
TEST(ProbabilistTests, TestIntervalsQLVL)
{
// read in data
xt::xtensor<double, 1> observed;
xt::xtensor<double, 2> predicted;
std::tie(observed, predicted) = load_data_p();
// read in expected results
auto expected = load_expected_p();
// compute scores
std::vector<std::string> metrics = {"CR"};
std::vector<xt::xarray<double>> results =
evalhyd::evalp(
// shape: (sites [1], time [t])
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::range(0, 4), xt::all())),
metrics,
xt::xtensor<double, 2>({}),
"", // events
{50., 80.}, // c_lvl
{0.1, 0.25, 0.75, 0.9} // q_lvl
);
// check results
for (std::size_t m = 0; m < metrics.size(); m++)
{
EXPECT_TRUE(xt::all(xt::isclose(
results[m], expected[metrics[m]], 1e-05, 1e-08, true
))) << "Failure for (" << metrics[m] << ")";
}
}
TEST(ProbabilistTests, TestMultiVariate)
{
// read in data
@@ -361,6 +397,7 @@ TEST(ProbabilistTests, TestMasks)
thresholds,
"high",
confidence_levels,
{}, // q_lvl
// shape: (sites [1], lead times [1], subsets [1], time [t])
masks
);
@@ -431,6 +468,7 @@ TEST(ProbabilistTests, TestMaskingConditions)
thresholds,
"high",
confidence_levels,
{}, // q_lvl
masks,
q_conditions
);
@@ -488,6 +526,7 @@ TEST(ProbabilistTests, TestMaskingConditions)
thresholds,
"high",
confidence_levels,
{}, // q_lvl
masks,
q_conditions_
);
@@ -542,6 +581,7 @@ TEST(ProbabilistTests, TestMaskingConditions)
thresholds,
"high",
confidence_levels,
{}, // q_lvl
masks,
t_conditions
);
@@ -705,6 +745,7 @@ TEST(ProbabilistTests, TestBootstrap)
thresholds,
"high", // events
confidence_levels,
{}, // q_lvl
xt::xtensor<bool, 4>({}), // t_msk
xt::xtensor<std::array<char, 32>, 2>({}), // m_cdt
bootstrap,
@@ -790,6 +831,7 @@ TEST(ProbabilistTests, TestBootstrapSummary)
thresholds,
"high", // events
confidence_levels,
{}, // q_lvl
xt::xtensor<bool, 4>({}), // t_msk
xt::xtensor<std::array<char, 32>, 2>({}), // m_cdt
bootstrap_0,
@@ -808,6 +850,7 @@ TEST(ProbabilistTests, TestBootstrapSummary)
thresholds,
"high", // events
confidence_levels,
{}, // q_lvl
xt::xtensor<bool, 4>({}), // t_msk
xt::xtensor<std::array<char, 32>, 2>({}), // m_cdt
bootstrap_1,
@@ -859,6 +902,7 @@ TEST(ProbabilistTests, TestBootstrapSummary)
thresholds,
"high", // events
confidence_levels,
{}, // q_lvl
xt::xtensor<bool, 4>({}), // t_msk
xt::xtensor<std::array<char, 32>, 2>({}), // m_cdt
bootstrap_2,
@@ -935,6 +979,7 @@ TEST(ProbabilistTests, TestCompleteness)
xt::xtensor<double, 2>({}), // thresholds
xtl::missing<const std::string>(), // events
{},
{}, // q_lvl
msk, // t_msk
xt::xtensor<std::array<char, 32>, 2>({}), // m_cdt
xtl::missing<const std::unordered_map<std::string, int>>(), // bootstrap
environment.yml
+ 1
− 1
@@ -8,7 +8,7 @@ dependencies:
- make
# Host dependencies
- xtl
- xtensor
- xtensor=0.25
# Test dependencies
- gtest
- gmock