diff --git a/environment.yml b/environment.yml
index deaaf2428e75d27e30fe3e9fda2f490ad4f6ed81..b1f5e68243c28f695864e2d64abfde09012e617d 100644
--- a/environment.yml
+++ b/environment.yml
@@ -8,7 +8,7 @@ dependencies:
- make
# Host dependencies
- xtl
- - xtensor
+ - xtensor=0.25
# Test dependencies
- gtest
- gmock
diff --git a/include/evalhyd/detail/probabilist/evaluator.hpp b/include/evalhyd/detail/probabilist/evaluator.hpp
index 3f94fd6d85dcdd94fa5da2224bad5dc914a8397f..c01c54be9479593c98ced5f6a23c472f6c510fb4 100644
--- a/include/evalhyd/detail/probabilist/evaluator.hpp
+++ b/include/evalhyd/detail/probabilist/evaluator.hpp
@@ -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;
@@ -157,6 +158,11 @@ namespace evalhyd
}
}
+ auto get_q_lvl()
+ {
+ return _q_lvl;
+ }
+
bool is_high_flow_event()
{
if (_events.has_value())
@@ -336,7 +342,7 @@ namespace evalhyd
if (!itv_bnds.has_value())
{
itv_bnds = elements::calc_itv_bnds(
- q_prd, get_c_lvl(),
+ q_prd, get_c_lvl(), get_q_lvl(),
n_sit, n_ldt, n_itv, n_tim
);
}
@@ -393,7 +399,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, get_q_lvl()
);
}
return qs.value();
@@ -482,12 +488,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}
{
diff --git a/include/evalhyd/detail/probabilist/intervals.hpp b/include/evalhyd/detail/probabilist/intervals.hpp
index adbd487b592ab4903210072e9ae01836ba41e0ac..f52fd5a3f5a574ee0ad10b56ba94af6c0045d37c 100644
--- a/include/evalhyd/detail/probabilist/intervals.hpp
+++ b/include/evalhyd/detail/probabilist/intervals.hpp
@@ -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,41 @@ 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 / 100., 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, xt::all(), xt::all(), std::min(res[0][0], res[0][1]), xt::all());
+ xt::view(itv_bnds, xt::all(), xt::all(), i, 1, xt::all()) =
+ xt::view(q_prd, xt::all(), xt::all(), std::max(res[0][0], res[0][1]), xt::all());
+ }
+ }
+ }
+ 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;
diff --git a/include/evalhyd/detail/probabilist/quantiles.hpp b/include/evalhyd/detail/probabilist/quantiles.hpp
index 88bd5284b5c3983a2b562ce8e88e37bd92d8b95a..bb13751ea119c9a4538269c29d9b9abfb1e7d62f 100644
--- a/include/evalhyd/detail/probabilist/quantiles.hpp
+++ b/include/evalhyd/detail/probabilist/quantiles.hpp
@@ -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 =
diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp
index f7f4340f210a2027d49d30d6fcf0a8505c0e4110..0039677616398742211b14cb5ce05e63feea3e08 100644
--- a/include/evalhyd/evalp.hpp
+++ b/include/evalhyd/evalp.hpp
@@ -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
diff --git a/tests/expected/evalp/CR_QLVL.csv b/tests/expected/evalp/CR_QLVL.csv
new file mode 100644
index 0000000000000000000000000000000000000000..dfefd7e9ae92a09448178e081d471ade3fe9cc3f
--- /dev/null
+++ b/tests/expected/evalp/CR_QLVL.csv
@@ -0,0 +1 @@
+0.00643087,0.0514469
diff --git a/tests/test_probabilist.cpp b/tests/test_probabilist.cpp
index 87a0303a7387cc0fa14e515ccf4c10ac4ece21da..a4bd6b0c5973e5380009e835202bfad36d27bd7a 100644
--- a/tests/test_probabilist.cpp
+++ b/tests/test_probabilist.cpp
@@ -39,6 +39,17 @@ std::vector<std::string> all_metrics_p = {
"ES"
};
+std::vector<std::string> all_metrics_p_ = {
+ "BS", "BSS", "BS_CRD", "BS_LBD", "REL_DIAG", "CRPS_FROM_BS",
+ "CRPS_FROM_ECDF",
+ "QS", "CRPS_FROM_QS",
+ "CONT_TBL", "POD", "POFD", "FAR", "CSI", "ROCSS",
+ "RANK_HIST", "DS", "AS",
+ "CR", "AW", "AWN", "WS",
+ "ES",
+ "CR_QLVL"
+};
+
std::tuple<xt::xtensor<double, 1>, xt::xtensor<double, 2>> load_data_p()
{
// read in data
@@ -60,7 +71,7 @@ std::unordered_map<std::string, xt::xarray<double>> load_expected_p()
std::ifstream ifs;
std::unordered_map<std::string, xt::xarray<double>> expected;
- for (const auto& metric : all_metrics_p)
+ for (const auto& metric : all_metrics_p_)
{
ifs.open(EVALHYD_DATA_DIR "/expected/evalp/" + metric + ".csv");
expected[metric] = xt::view(
@@ -246,6 +257,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 +308,42 @@ 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<std::string> metrics_ = {"CR_QLVL"};
+
+ 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::keep(0, 15, 30, 50), xt::all())),
+ metrics,
+ xt::xtensor<double, 2>({}),
+ "", // events
+ {50., 80.}, // c_lvl
+ {10., 25., 75., 90.} // 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 +409,7 @@ TEST(ProbabilistTests, TestMasks)
thresholds,
"high",
confidence_levels,
+ {}, // q_lvl
// shape: (sites [1], lead times [1], subsets [1], time [t])
masks
);
@@ -431,6 +480,7 @@ TEST(ProbabilistTests, TestMaskingConditions)
thresholds,
"high",
confidence_levels,
+ {}, // q_lvl
masks,
q_conditions
);
@@ -488,6 +538,7 @@ TEST(ProbabilistTests, TestMaskingConditions)
thresholds,
"high",
confidence_levels,
+ {}, // q_lvl
masks,
q_conditions_
);
@@ -542,6 +593,7 @@ TEST(ProbabilistTests, TestMaskingConditions)
thresholds,
"high",
confidence_levels,
+ {}, // q_lvl
masks,
t_conditions
);
@@ -705,6 +757,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 +843,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 +862,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 +914,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 +991,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