Commit 92de7a26 authored by Thibault Hallouin's avatar Thibault Hallouin
Browse files

add unittest on QS/CRPS

No related merge requests found
Pipeline #37099 passed with stages
in 9 seconds
Showing with 59 additions and 8 deletions
+59 -8
...@@ -20,18 +20,16 @@ TEST(ProbabilistTests, TestBrier) { ...@@ -20,18 +20,16 @@ TEST(ProbabilistTests, TestBrier) {
xt::xtensor<double, 2> predicted = xt::load_csv<double>(ifs); xt::xtensor<double, 2> predicted = xt::load_csv<double>(ifs);
ifs.close(); ifs.close();
// add "artificial" site dimension (size 1) to observations
auto obs = xt::view(observed, xt::newaxis(), xt::all());
// add "artificial" site dimension (size 1) and lead time dimension (size 1)
// to predictions
auto prd = xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all());
// compute scores // compute scores
xt::xtensor<double, 1> thresholds = {690, 534, 445}; xt::xtensor<double, 1> thresholds = {690, 534, 445};
std::vector<xt::xarray<double>> metrics = std::vector<xt::xarray<double>> metrics =
evalhyd::evalp( evalhyd::evalp(
obs, prd, {"BS", "BSS", "BS_CRD", "BS_LBD"}, // shape: (sites [1], time [t])
xt::view(observed, xt::newaxis(), xt::all()),
// shape: (sites [1], lead times [1], members [m], time [t])
xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all()),
{"BS", "BSS", "BS_CRD", "BS_LBD"},
thresholds thresholds
); );
...@@ -61,6 +59,53 @@ TEST(ProbabilistTests, TestBrier) { ...@@ -61,6 +59,53 @@ TEST(ProbabilistTests, TestBrier) {
EXPECT_TRUE(xt::allclose(metrics[3], bs_lbd)); EXPECT_TRUE(xt::allclose(metrics[3], bs_lbd));
} }
TEST(ProbabilistTests, TestQuantiles) {
// read in data
std::ifstream ifs;
ifs.open("./data/q_obs.csv");
xt::xtensor<double, 1> observed = xt::squeeze(xt::load_csv<int>(ifs));
ifs.close();
ifs.open("./data/q_prd.csv");
xt::xtensor<double, 2> predicted = xt::load_csv<double>(ifs);
ifs.close();
// compute scores
xt::xtensor<double, 1> thresholds = {690, 534, 445};
std::vector<xt::xarray<double>> metrics =
evalhyd::evalp(
// shape: (sites [1], time [t])
xt::view(observed, xt::newaxis(), xt::all()),
// shape: (sites [1], lead times [1], members [m], time [t])
xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all()),
{"QS", "CRPS"},
thresholds
);
// check results
// Quantile scores
xt::xtensor<double, 4> qs =
{{{{ 345.91578 , 345.069256, 343.129359, 340.709869, 338.281598,
335.973535, 333.555157, 330.332426, 327.333539, 324.325996,
321.190082, 318.175117, 315.122186, 311.97205 , 308.644942,
305.612169, 302.169552, 298.445956, 294.974648, 291.273807,
287.724586, 284.101905, 280.235592, 276.21865 , 272.501484,
268.652733, 264.740168, 260.8558 , 256.90329 , 252.926292,
248.931239, 244.986396, 240.662998, 236.328964, 232.089785,
227.387089, 222.976008, 218.699975, 214.099678, 209.67252 ,
205.189587, 200.395746, 195.2372 , 190.080139, 185.384244,
180.617858, 174.58323 , 169.154093, 163.110932, 156.274796,
147.575315}}}};
EXPECT_TRUE(xt::allclose(metrics[0], qs));
// Continuous ranked probability scores
xt::xtensor<double, 3> crps =
{{{ 257.412129}}};
EXPECT_TRUE(xt::allclose(metrics[1], crps));
}
TEST(ProbabilistTests, TestMasks) { TEST(ProbabilistTests, TestMasks) {
// read in data // read in data
std::ifstream ifs; std::ifstream ifs;
...@@ -79,21 +124,27 @@ TEST(ProbabilistTests, TestMasks) { ...@@ -79,21 +124,27 @@ TEST(ProbabilistTests, TestMasks) {
// compute scores using masks to subset whole record // compute scores using masks to subset whole record
xt::xtensor<double, 1> thresholds = {690, 534, 445}; xt::xtensor<double, 1> thresholds = {690, 534, 445};
std::vector<std::string> metrics = {"BS", "BSS", "BS_CRD", "BS_LBD"}; std::vector<std::string> metrics =
{"BS", "BSS", "BS_CRD", "BS_LBD", "QS", "CRPS"};
std::vector<xt::xarray<double>> metrics_masked = std::vector<xt::xarray<double>> metrics_masked =
evalhyd::evalp( evalhyd::evalp(
// shape: (sites [1], time [t])
xt::view(observed, xt::newaxis(), xt::all()), xt::view(observed, xt::newaxis(), xt::all()),
// shape: (sites [1], lead times [1], members [m], time [t])
xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all()), xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all()),
metrics, metrics,
thresholds, thresholds,
// shape: (subsets [1], time [t])
masks masks
); );
// compute scores on pre-computed subset of whole record // compute scores on pre-computed subset of whole record
std::vector<xt::xarray<double>> metrics_subset = std::vector<xt::xarray<double>> metrics_subset =
evalhyd::evalp( evalhyd::evalp(
// shape: (sites [1], time [t-20])
xt::view(observed, xt::newaxis(), xt::range(20, _)), xt::view(observed, xt::newaxis(), xt::range(20, _)),
// shape: (sites [1], lead times [1], members [m], time [t-20])
xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::range(20, _)), xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::range(20, _)),
metrics, metrics,
thresholds thresholds
......
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