An error occurred while loading the file. Please try again.
-
Thibault Hallouin authoredin `xt::stddev``, when passing a scalar, it is passed to *ddof* instead of *axes*, so to compute the standard deviation on a given axis, the axis must be passed in curly braces
64eca8fa
// Copyright (c) 2023, INRAE.// Distributed under the terms of the GPL-3 Licence.// The full licence is in the file LICENCE, distributed with this software.#include <fstream>#include <vector>#include <tuple>#include <array>#include <gtest/gtest.h>#include <xtl/xoptional.hpp>#include <xtensor/xtensor.hpp>#include <xtensor/xview.hpp>#include <xtensor/xmath.hpp>#include <xtensor/xsort.hpp>#include <xtensor/xmanipulation.hpp>#include <xtensor/xcsv.hpp>#include "evalhyd/evalp.hpp"#ifndef EVALHYD_DATA_DIR#error "need to define data directory"#endifusing namespace xt::placeholders; // required for `_` to workstd::vector<std::string> all_metrics_p = { "BS", "BSS", "BS_CRD", "BS_LBD", "QS", "CRPS", "POD", "POFD", "FAR", "CSI", "ROCSS", "RANK_HIST", "DS", "AS", "CR", "AW", "AWN", "AWI", "WS", "WSS"};std::tuple<xt::xtensor<double, 1>, xt::xtensor<double, 2>> load_data_p(){ // read in data std::ifstream ifs; ifs.open(EVALHYD_DATA_DIR "/q_obs.csv"); xt::xtensor<double, 1> observed = xt::squeeze(xt::load_csv<int>(ifs)); ifs.close(); ifs.open(EVALHYD_DATA_DIR "/q_prd.csv"); xt::xtensor<double, 2> predicted = xt::load_csv<double>(ifs); ifs.close(); return std::make_tuple(observed, predicted);}TEST(ProbabilistTests, TestBrier){ // read in data xt::xtensor<double, 1> observed; xt::xtensor<double, 2> predicted; std::tie(observed, predicted) = load_data_p(); // compute scores xt::xtensor<double, 2> thresholds = {{690, 534, 445, NAN}}; std::vector<xt::xarray<double>> metrics = 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::all(), xt::all())), {"BS", "BSS", "BS_CRD", "BS_LBD"}, thresholds, "high"7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
);
// check results
// Brier scores
xt::xtensor<double, 5> bs =
{{{{{0.10615136, 0.07395622, 0.08669186, NAN}}}}};
EXPECT_TRUE(
xt::all(xt::isclose(metrics[0], bs, 1e-05, 1e-08, true))
);
// Brier skill scores
xt::xtensor<double, 5> bss =
{{{{{0.5705594, 0.6661165, 0.5635126, NAN}}}}};
EXPECT_TRUE(
xt::all(xt::isclose(metrics[1], bss, 1e-05, 1e-08, true))
);
// Brier calibration-refinement decompositions
xt::xtensor<double, 6> bs_crd =
{{{{{{0.011411758, 0.1524456, 0.2471852},
{0.005532413, 0.1530793, 0.2215031},
{0.010139431, 0.1220601, 0.1986125},
{NAN, NAN, NAN}}}}}};
EXPECT_TRUE(
xt::all(xt::isclose(metrics[2], bs_crd, 1e-05, 1e-08, true))
);
// Brier likelihood-base rate decompositions
xt::xtensor<double, 6> bs_lbd =
{{{{{{0.012159881, 0.1506234, 0.2446149},
{0.008031746, 0.1473869, 0.2133114},
{0.017191279, 0.1048221, 0.1743227},
{NAN, NAN, NAN}}}}}};
EXPECT_TRUE(
xt::all(xt::isclose(metrics[3], bs_lbd, 1e-05, 1e-08, true))
);
}
TEST(ProbabilistTests, TestQuantiles)
{
// read in data
xt::xtensor<double, 1> observed;
xt::xtensor<double, 2> predicted;
std::tie(observed, predicted) = load_data_p();
// compute scores
std::vector<xt::xarray<double>> metrics =
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::all(), xt::all())),
{"QS", "CRPS"}
);
// check results
// Quantile scores
xt::xtensor<double, 5> 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::all(xt::isclose(metrics[0], qs)));
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
// Continuous ranked probability scores
xt::xtensor<double, 4> crps =
{{{{252.956919}}}};
EXPECT_TRUE(xt::all(xt::isclose(metrics[1], crps)));
}
TEST(ProbabilistTests, TestContingency)
{
// read in data
xt::xtensor<double, 1> observed;
xt::xtensor<double, 2> predicted;
std::tie(observed, predicted) = load_data_p();
// compute scores
xt::xtensor<double, 2> thresholds = {{690, 534, 445, NAN}};
std::vector<xt::xarray<double>> metrics =
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::all(), xt::all())),
{"POD", "POFD", "FAR", "CSI", "ROCSS"},
thresholds,
"low"
);
// check results
// POD
xt::xtensor<double, 6> pod =
{{{{{{ 1. , 1. , 1. , NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.863309, 0.854369, 0.752941, NAN},
{ 0.848921, 0.854369, 0.752941, NAN},
{ 0.848921, 0.854369, 0.752941, NAN},
{ 0.848921, 0.84466 , 0.752941, NAN}}}}}};
EXPECT_TRUE(
xt::all(xt::isclose(metrics[0], pod, 1e-05, 1e-08, true))
);
// POFD
xt::xtensor<double, 6> pofd =
{{{{{{ 1. , 1. , 1. , NAN},
{ 0.087209, 0.038462, 0.026549, NAN},
{ 0.087209, 0.038462, 0.026549, NAN},
{ 0.087209, 0.038462, 0.026549, NAN},
{ 0.087209, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
{ 0.081395, 0.038462, 0.026549, NAN},
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
{ 0.081395, 0.038462, 0.022124, NAN}}}}}};
EXPECT_TRUE(
xt::all(xt::isclose(metrics[1], pofd, 1e-04, 1e-07, true))
);
// FAR
xt::xtensor<double, 6> far =
{{{{{{ 0.553055, 0.66881 , 0.726688, NAN},
{ 0.111111, 0.083333, 0.085714, NAN},
{ 0.111111, 0.083333, 0.085714, NAN},
{ 0.111111, 0.083333, 0.085714, NAN},
{ 0.111111, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.104478, 0.083333, 0.085714, NAN},
{ 0.106061, 0.083333, 0.085714, NAN},
{ 0.106061, 0.083333, 0.085714, NAN},
{ 0.106061, 0.084211, 0.072464, NAN}}}}}};
EXPECT_TRUE(
xt::all(xt::isclose(metrics[2], far, 1e-05, 1e-08, true))
);
// CSI
xt::xtensor<double, 6> csi =
{{{{{{ 0.446945, 0.33119 , 0.273312, NAN},
{ 0.779221, 0.792793, 0.703297, NAN},
{ 0.779221, 0.792793, 0.703297, NAN},
{ 0.779221, 0.792793, 0.703297, NAN},
{ 0.779221, 0.792793, 0.703297, NAN},
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.784314, 0.792793, 0.703297, NAN},
{ 0.771242, 0.792793, 0.703297, NAN},
{ 0.771242, 0.792793, 0.703297, NAN},
{ 0.771242, 0.783784, 0.711111, NAN}}}}}}
;
EXPECT_TRUE(
xt::all(xt::isclose(metrics[3], csi, 1e-05, 1e-08, true))
);
// ROC skill scores
xt::xtensor<double, 5> rocss =
{{{{{ 0.71085 , 0.783047, 0.713066, NAN}}}}};
EXPECT_TRUE(
xt::all(xt::isclose(metrics[4], rocss, 1e-05, 1e-08, true))
);
}
TEST(ProbabilistTests, TestRanks)
{
// read in data
xt::xtensor<double, 1> observed;
xt::xtensor<double, 2> predicted;
std::tie(observed, predicted) = load_data_p();
// compute scores
std::vector<xt::xarray<double>> metrics =
evalhyd::evalp(
421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490
// 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::all(), xt::all())),
{"RANK_HIST", "DS", "AS"},
xt::xtensor<double, 2>({}),
"high", // events
{}, // c_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
{}, // dts
7 // seed
);
// check results
// Rank histogram
xt::xtensor<double, 5> rank_hist;
#if EVALHYD_TESTING_OS == WINDOWS
rank_hist = {{{{{ 0.607717, 0. , 0. , 0. , 0. , 0.003215,
0. , 0. , 0. , 0. , 0. , 0. ,
0.003215, 0. , 0.003215, 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. , 0.003215,
0. , 0. , 0. , 0. , 0. , 0.006431,
0. , 0. , 0.003215, 0.006431, 0. , 0. ,
0. , 0.003215, 0. , 0. , 0.003215, 0.003215,
0.003215, 0. , 0.006431, 0.344051}}}}};
#elif EVALHYD_TESTING_OS == MACOS
rank_hist = {{{{{ 0.607717, 0. , 0. , 0. , 0. , 0.003215,
0. , 0. , 0. , 0. , 0. , 0. ,
0.003215, 0. , 0.003215, 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. , 0.003215,
0. , 0. , 0. , 0. , 0. , 0.006431,
0. , 0. , 0.003215, 0.006431, 0. , 0. ,
0. , 0.003215, 0. , 0. , 0.003215, 0.003215,
0.003215, 0. , 0.006431, 0.344051}}}}};
#elif EVALHYD_TESTING_OS == LINUX
rank_hist = {{{{{ 0.607717, 0. , 0. , 0. , 0. , 0.003215,
0. , 0. , 0. , 0. , 0. , 0. ,
0.003215, 0. , 0.003215, 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. , 0.003215,
0. , 0. , 0. , 0. , 0. , 0.006431,
0. , 0. , 0.003215, 0.006431, 0. , 0. ,
0. , 0.003215, 0. , 0. , 0.003215, 0.003215,
0.003215, 0. , 0.006431, 0.344051}}}}};
#endif
EXPECT_TRUE(
xt::all(xt::isclose(metrics[0], rank_hist, 1e-04, 1e-06, true))
);
// Delta scores
xt::xtensor<double, 4> ds;
#if EVALHYD_TESTING_OS == WINDOWS
ds = {{{{ 148.790164}}}};
#elif EVALHYD_TESTING_OS == MACOS
ds = {{{{ 148.790164}}}};
#elif EVALHYD_TESTING_OS == LINUX
ds = {{{{ 148.790164}}}};
#endif
EXPECT_TRUE(
xt::all(xt::isclose(metrics[1], ds, 1e-04, 1e-07, true))
);
// Alpha scores
xt::xtensor<double, 4> as;
#if EVALHYD_TESTING_OS == WINDOWS
491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560
as = {{{{ 0.491481}}}};
#elif EVALHYD_TESTING_OS == MACOS
as = {{{{ 0.491481}}}};
#elif EVALHYD_TESTING_OS == LINUX
as = {{{{ 0.491481}}}};
#endif
EXPECT_TRUE(
xt::all(xt::isclose(metrics[2], as, 1e-04, 1e-07, true))
);
}
TEST(ProbabilistTests, TestIntervals)
{
// read in data
xt::xtensor<double, 1> observed;
xt::xtensor<double, 2> predicted;
std::tie(observed, predicted) = load_data_p();
// compute scores
std::vector<xt::xarray<double>> metrics =
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::all(), xt::all())),
{"CR", "AW", "AWN", "AWI", "WS", "WSS"},
xt::xtensor<double, 2>({}),
"", // events
{30., 80.} // c_lvl
);
// check results
// coverage ratios
xt::xtensor<double, 5> cr = {{{{{ 0.006431, 0.03537 }}}}};
EXPECT_TRUE(xt::all(xt::isclose(metrics[0], cr, 1e-05, 1e-06, true)));
// average widths
xt::xtensor<double, 5> aw = {{{{{ 9.27492, 31.321543}}}}};
EXPECT_TRUE(xt::all(xt::isclose(metrics[1], aw, 1e-05, 1e-06, true)));
// average widths normalised
xt::xtensor<double, 5> awn = {{{{{ 0.007383, 0.024931}}}}};
EXPECT_TRUE(xt::all(xt::isclose(metrics[2], awn, 1e-05, 1e-06, true)));
// average widths indices
xt::xtensor<double, 5> awi = {{{{{ 0.982112, 0.988095}}}}};
EXPECT_TRUE(xt::all(xt::isclose(metrics[3], awi, 1e-05, 1e-06, true)));
// Winkler scores
xt::xtensor<double, 5> ws = {{{{{ 764.447175, 2578.138264}}}}};
EXPECT_TRUE(xt::all(xt::isclose(metrics[4], ws, 1e-05, 1e-06, true)));
// Winkler skill scores
xt::xtensor<double, 5> wss = {{{{{ 0.662189, 0.436039}}}}};
EXPECT_TRUE(xt::all(xt::isclose(metrics[5], wss, 1e-05, 1e-06, true)));
}
TEST(ProbabilistTests, TestMasks)
{
// read in data
xt::xtensor<double, 1> observed;
xt::xtensor<double, 2> predicted;
std::tie(observed, predicted) = load_data_p();
// generate temporal subset by dropping 20 first time steps
xt::xtensor<bool, 4> masks =
xt::ones<bool>({std::size_t {1}, std::size_t {1}, std::size_t {1},
std::size_t {observed.size()}});
xt::view(masks, 0, xt::all(), 0, xt::range(0, 20)) = 0;
561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630
// compute scores using masks to subset whole record
xt::xtensor<double, 2> thresholds = {{690, 534, 445}};
std::vector<double> confidence_levels = {30., 80.};
std::vector<xt::xarray<double>> metrics_masked =
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::all(), xt::all())),
all_metrics_p,
thresholds,
"high",
confidence_levels,
// shape: (sites [1], lead times [1], subsets [1], time [t])
masks
);
// compute scores on pre-computed subset of whole record
std::vector<xt::xarray<double>> metrics_subset =
evalhyd::evalp(
// shape: (sites [1], time [t-20])
xt::eval(xt::view(observed, xt::newaxis(), xt::range(20, _))),
// shape: (sites [1], lead times [1], members [m], time [t-20])
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::range(20, _))),
all_metrics_p,
thresholds,
"high",
confidence_levels
);
// check results are identical
for (std::size_t m = 0; m < all_metrics_p.size(); m++)
{
// ---------------------------------------------------------------------
// /!\ skip ranks-based metrics because it contains a random process
// for which setting the seed will not work because the time series
// lengths are different between "masked" and "subset", which
// results in different tensor shapes, and hence in different
// random ranks for ties
if ((all_metrics_p[m] == "RANK_HIST")
|| (all_metrics_p[m] == "DS")
|| (all_metrics_p[m] == "AS"))
{
continue;
}
// ---------------------------------------------------------------------
EXPECT_TRUE(xt::all(xt::isclose(metrics_masked[m], metrics_subset[m], 1e-04, 1e-07, true)))
<< "Failure for (" << all_metrics_p[m] << ")";
}
}
TEST(ProbabilistTests, TestMaskingConditions)
{
xt::xtensor<double, 2> thresholds = {{690, 534, 445}};
std::vector<double> confidence_levels = {30., 80.};
// read in data
xt::xtensor<double, 1> observed_;
xt::xtensor<double, 2> predicted;
std::tie(observed_, predicted) = load_data_p();
// turn observed into 2D view (to simplify syntax later on)
auto observed = xt::view(observed_, xt::newaxis(), xt::all());
// generate dummy empty masks required to access next optional argument
xt::xtensor<bool, 4> masks;
// conditions on streamflow values _________________________________________
631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700
// compute scores using masking conditions on streamflow to subset whole record
xt::xtensor<std::array<char, 32>, 2> q_conditions = {
{std::array<char, 32> {"q_obs{<2000,>3000}"}}
};
std::vector<xt::xarray<double>> metrics_q_conditioned =
evalhyd::evalp(
xt::eval(observed),
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
all_metrics_p,
thresholds,
"high",
confidence_levels,
masks,
q_conditions
);
// compute scores using "NaN-ed" time indices where conditions on streamflow met
std::vector<xt::xarray<double>> metrics_q_preconditioned =
evalhyd::evalp(
xt::eval(xt::where((observed < 2000) | (observed > 3000), observed, NAN)),
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
all_metrics_p,
thresholds,
"high",
confidence_levels
);
// check results are identical
for (std::size_t m = 0; m < all_metrics_p.size(); m++)
{
// ---------------------------------------------------------------------
// /!\ skip ranks-based metrics because it contains a random process
// for which setting the seed will not work because the time series
// lengths are different between "conditioned" and "preconditioned",
// which results in different tensor shapes, and hence in different
// random ranks for ties
if ((all_metrics_p[m] == "RANK_HIST")
|| (all_metrics_p[m] == "DS")
|| (all_metrics_p[m] == "AS"))
{
continue;
}
// ---------------------------------------------------------------------
EXPECT_TRUE(
xt::all(xt::isclose(metrics_q_conditioned[m],
metrics_q_preconditioned[m],
1e-05, 1e-08, true))
) << "Failure for (" << all_metrics_p[m] << ")";
}
// conditions on streamflow statistics _____________________________________
// compute scores using masking conditions on streamflow to subset whole record
xt::xtensor<std::array<char, 32>, 2> q_conditions_ = {
{std::array<char, 32> {"q_prd_mean{>=median}"}}
};
auto q_prd_mean = xt::mean(predicted, {0}, xt::keep_dims);
double median = xt::median(q_prd_mean);
std::vector<xt::xarray<double>> metrics_q_conditioned_ =
evalhyd::evalp(
xt::eval(observed),
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
all_metrics_p,
thresholds,
"high",
701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770
confidence_levels,
masks,
q_conditions_
);
// compute scores using "NaN-ed" time indices where conditions on streamflow met
std::vector<xt::xarray<double>> metrics_q_preconditioned_ =
evalhyd::evalp(
xt::eval(xt::where(q_prd_mean >= median, observed, NAN)),
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
all_metrics_p,
thresholds,
"high",
confidence_levels
);
// check results are identical
for (std::size_t m = 0; m < all_metrics_p.size(); m++)
{
// ---------------------------------------------------------------------
// /!\ skip ranks-based metrics because it contains a random process
// for which setting the seed will not work because the time series
// lengths are different between "conditioned" and "preconditioned",
// which results in different tensor shapes, and hence in different
// random ranks for ties
if ((all_metrics_p[m] == "RANK_HIST")
|| (all_metrics_p[m] == "DS")
|| (all_metrics_p[m] == "AS"))
{
continue;
}
// ---------------------------------------------------------------------
EXPECT_TRUE(
xt::all(xt::isclose(metrics_q_conditioned_[m],
metrics_q_preconditioned_[m],
1e-05, 1e-08, true))
) << "Failure for (" << all_metrics_p[m] << ")";
}
// conditions on temporal indices __________________________________________
// compute scores using masking conditions on time indices to subset whole record
xt::xtensor<std::array<char, 32>, 2> t_conditions = {
{std::array<char, 32> {"t{0,1,2,3,4,5:97,97,98,99}"}}
};
std::vector<xt::xarray<double>> metrics_t_conditioned =
evalhyd::evalp(
xt::eval(observed),
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
all_metrics_p,
thresholds,
"high",
confidence_levels,
masks,
t_conditions
);
// compute scores on already subset time series
std::vector<xt::xarray<double>> metrics_t_subset =
evalhyd::evalp(
xt::eval(xt::view(observed_, xt::newaxis(), xt::range(0, 100))),
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::range(0, 100))),
all_metrics_p,
thresholds,
"high",
confidence_levels
);
771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840
// check results are identical
for (std::size_t m = 0; m < all_metrics_p.size(); m++)
{
// ---------------------------------------------------------------------
// /!\ skip ranks-based metrics because it contains a random process
// for which setting the seed will not work because the time series
// lengths are different between "conditioned" and "subset",
// which results in different tensor shapes, and hence in different
// random ranks for ties
if ((all_metrics_p[m] == "RANK_HIST")
|| (all_metrics_p[m] == "DS")
|| (all_metrics_p[m] == "AS"))
{
continue;
}
// ---------------------------------------------------------------------
EXPECT_TRUE(
xt::all(xt::isclose(metrics_t_conditioned[m],
metrics_t_subset[m],
1e-05, 1e-08, true))
) << "Failure for (" << all_metrics_p[m] << ")";
}
}
TEST(ProbabilistTests, TestMissingData)
{
xt::xtensor<double, 2> thresholds = {{ 4., 5. }};
std::vector<double> confidence_levels = {30., 80.};
// compute metrics on series with NaN
xt::xtensor<double, 4> forecast_nan {{
// leadtime 1
{{ 5.3, 4.2, 5.7, 2.3, NAN },
{ 4.3, 4.2, 4.7, 4.3, NAN },
{ 5.3, 5.2, 5.7, 2.3, NAN }},
// leadtime 2
{{ NAN, 4.2, 5.7, 2.3, 3.1 },
{ NAN, 4.2, 4.7, 4.3, 3.3 },
{ NAN, 5.2, 5.7, 2.3, 3.9 }}
}};
xt::xtensor<double, 2> observed_nan
{{ 4.7, 4.3, NAN, 2.7, 4.1 }};
std::vector<xt::xarray<double>> metrics_nan =
evalhyd::evalp(
observed_nan,
forecast_nan,
all_metrics_p,
thresholds,
"high",
confidence_levels
);
// compute metrics on manually subset series (one leadtime at a time)
xt::xtensor<double, 4> forecast_pp1 {{
// leadtime 1
{{ 5.3, 4.2, 2.3 },
{ 4.3, 4.2, 4.3 },
{ 5.3, 5.2, 2.3 }},
}};
xt::xtensor<double, 2> observed_pp1
{{ 4.7, 4.3, 2.7 }};
std::vector<xt::xarray<double>> metrics_pp1 =
evalhyd::evalp(
observed_pp1,
forecast_pp1,
841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910
all_metrics_p,
thresholds,
"high",
confidence_levels
);
xt::xtensor<double, 4> forecast_pp2 {{
// leadtime 2
{{ 4.2, 2.3, 3.1 },
{ 4.2, 4.3, 3.3 },
{ 5.2, 2.3, 3.9 }}
}};
xt::xtensor<double, 2> observed_pp2
{{ 4.3, 2.7, 4.1 }};
std::vector<xt::xarray<double>> metrics_pp2 =
evalhyd::evalp(
observed_pp2,
forecast_pp2,
all_metrics_p,
thresholds,
"high",
confidence_levels
);
// check that numerical results are identical
for (std::size_t m = 0; m < all_metrics_p.size(); m++)
{
// for leadtime 1
EXPECT_TRUE(
xt::all(xt::isclose(xt::view(metrics_nan[m], xt::all(), 0),
xt::view(metrics_pp1[m], xt::all(), 0),
1e-05, 1e-08, true))
) << "Failure for (" << all_metrics_p[m] << ", " << "leadtime 1)";
// for leadtime 2
EXPECT_TRUE(
xt::all(xt::isclose(xt::view(metrics_nan[m], xt::all(), 1),
xt::view(metrics_pp2[m], xt::all(), 0),
1e-05, 1e-08, true))
) << "Failure for (" << all_metrics_p[m] << ", " << "leadtime 2)";
}
}
TEST(ProbabilistTests, TestBootstrap)
{
xt::xtensor<double, 2> thresholds = {{ 33.87, 55.67 }};
std::vector<double> confidence_levels = {30., 80.};
// read in data
std::ifstream ifs;
ifs.open(EVALHYD_DATA_DIR "/q_obs_1yr.csv");
xt::xtensor<std::string, 1> x_dts = xt::squeeze(xt::load_csv<std::string>(ifs, ',', 0, 1));
ifs.close();
std::vector<std::string> datetimes (x_dts.begin(), x_dts.end());
ifs.open(EVALHYD_DATA_DIR "/q_obs_1yr.csv");
xt::xtensor<double, 1> observed = xt::squeeze(xt::load_csv<double>(ifs, ',', 1));
ifs.close();
ifs.open(EVALHYD_DATA_DIR "/q_prd_1yr.csv");
xt::xtensor<double, 2> predicted = xt::load_csv<double>(ifs, ',', 1);
ifs.close();
// compute metrics via bootstrap
std::unordered_map<std::string, int> bootstrap =
{{"n_samples", 10}, {"len_sample", 3}, {"summary", 0}};
911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980
std::vector<xt::xarray<double>> metrics_bts =
evalhyd::evalp(
xt::eval(xt::view(observed, xt::newaxis(), xt::all())),
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
all_metrics_p,
thresholds,
"high", // events
confidence_levels,
xt::xtensor<bool, 4>({}), // t_msk
xt::xtensor<std::array<char, 32>, 2>({}), // m_cdt
bootstrap,
datetimes
);
// compute metrics by repeating year of data 3 times
// (since there is only one year of data, and that the bootstrap works on
// one-year blocks, it can only select that given year to form samples,
// and the length of the sample corresponds to how many times this year
// is repeated in the sample, so that repeating the input data this many
// times should result in the same numerical results)
xt::xtensor<double, 1> observed_x3 =
xt::concatenate(xt::xtuple(observed, observed, observed), 0);
xt::xtensor<double, 2> predicted_x3 =
xt::concatenate(xt::xtuple(predicted, predicted, predicted), 1);
std::vector<xt::xarray<double>> metrics_rep =
evalhyd::evalp(
xt::eval(xt::view(observed_x3, xt::newaxis(), xt::all())),
xt::eval(xt::view(predicted_x3, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
all_metrics_p,
thresholds,
"high",
confidence_levels
);
// check results are identical
for (std::size_t m = 0; m < all_metrics_p.size(); m++)
{
// ---------------------------------------------------------------------
// /!\ skip ranks-based metrics because it contains a random process
// for which setting the seed will not work because the time series
// lengths are different between "bts" and "rep", which
// results in different tensor shapes, and hence in different
// random ranks for ties
if ((all_metrics_p[m] == "RANK_HIST")
|| (all_metrics_p[m] == "DS")
|| (all_metrics_p[m] == "AS"))
{
continue;
}
// ---------------------------------------------------------------------
EXPECT_TRUE(
xt::all(xt::isclose(
metrics_bts[m], metrics_rep[m]
))
) << "Failure for (" << all_metrics_p[m] << ")";
}
}
TEST(ProbabilistTests, TestBootstrapSummary)
{
xt::xtensor<double, 2> thresholds = {{ 33.87, 55.67 }};
std::vector<double> confidence_levels = {30., 80.};
// read in data
std::ifstream ifs;
ifs.open(EVALHYD_DATA_DIR "/q_obs_1yr.csv");
xt::xtensor<std::string, 1> x_dts = xt::squeeze(xt::load_csv<std::string>(ifs, ',', 0, 1));
981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050
ifs.close();
std::vector<std::string> datetimes (x_dts.begin(), x_dts.end());
ifs.open(EVALHYD_DATA_DIR "/q_obs_1yr.csv");
xt::xtensor<double, 1> observed = xt::squeeze(xt::load_csv<double>(ifs, ',', 1));
ifs.close();
ifs.open(EVALHYD_DATA_DIR "/q_prd_1yr.csv");
xt::xtensor<double, 2> predicted = xt::load_csv<double>(ifs, ',', 1);
ifs.close();
// compute metrics via bootstrap
std::unordered_map<std::string, int> bootstrap_0 =
{{"n_samples", 10}, {"len_sample", 3}, {"summary", 0}};
std::vector<xt::xarray<double>> metrics_raw =
evalhyd::evalp(
xt::eval(xt::view(observed, xt::newaxis(), xt::all())),
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
all_metrics_p,
thresholds,
"high", // events
confidence_levels,
xt::xtensor<bool, 4>({}), // t_msk
xt::xtensor<std::array<char, 32>, 2>({}), // m_cdt
bootstrap_0,
datetimes
);
// compute metrics via bootstrap with mean and standard deviation summary
std::unordered_map<std::string, int> bootstrap_1 =
{{"n_samples", 10}, {"len_sample", 3}, {"summary", 1}};
std::vector<xt::xarray<double>> metrics_mas =
evalhyd::evalp(
xt::eval(xt::view(observed, xt::newaxis(), xt::all())),
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
all_metrics_p,
thresholds,
"high", // events
confidence_levels,
xt::xtensor<bool, 4>({}), // t_msk
xt::xtensor<std::array<char, 32>, 2>({}), // m_cdt
bootstrap_1,
datetimes
);
// check results are identical
for (std::size_t m = 0; m < all_metrics_p.size(); m++)
{
// ---------------------------------------------------------------------
// /!\ skip ranks-based metrics because it contains a random process
// for which setting the seed will not work because the time series
// lengths are different between "bts" and "rep", which
// results in different tensor shapes, and hence in different
// random ranks for ties
if ((all_metrics_p[m] == "RANK_HIST")
|| (all_metrics_p[m] == "DS")
|| (all_metrics_p[m] == "AS"))
{
continue;
}
// ---------------------------------------------------------------------
// mean
EXPECT_TRUE(
xt::all(xt::isclose(
xt::mean(metrics_raw[m], {3}),
xt::view(metrics_mas[m], xt::all(), xt::all(), xt::all(), 0)
))
10511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112
) << "Failure for (" << all_metrics_p[m] << ") on mean";
// standard deviation
EXPECT_TRUE(
xt::all(xt::isclose(
xt::stddev(metrics_raw[m], {3}),
xt::view(metrics_mas[m], xt::all(), xt::all(), xt::all(), 1)
))
) << "Failure for (" << all_metrics_p[m] << ") on standard deviation";
}
// compute metrics via bootstrap with quantiles summary
std::unordered_map<std::string, int> bootstrap_2 =
{{"n_samples", 10}, {"len_sample", 3}, {"summary", 2}};
std::vector<xt::xarray<double>> metrics_qtl =
evalhyd::evalp(
xt::eval(xt::view(observed, xt::newaxis(), xt::all())),
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
all_metrics_p,
thresholds,
"high", // events
confidence_levels,
xt::xtensor<bool, 4>({}), // t_msk
xt::xtensor<std::array<char, 32>, 2>({}), // m_cdt
bootstrap_2,
datetimes
);
// check results are identical
for (std::size_t m = 0; m < all_metrics_p.size(); m++)
{
// ---------------------------------------------------------------------
// /!\ skip ranks-based metrics because it contains a random process
// for which setting the seed will not work because the time series
// lengths are different between "bts" and "rep", which
// results in different tensor shapes, and hence in different
// random ranks for ties
if ((all_metrics_p[m] == "RANK_HIST")
|| (all_metrics_p[m] == "DS")
|| (all_metrics_p[m] == "AS"))
{
continue;
}
// ---------------------------------------------------------------------
// quantiles
std::vector<double> quantiles = {0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95};
std::size_t i = 0;
for (auto q : quantiles)
{
EXPECT_TRUE(
xt::all(xt::isclose(
xt::quantile(metrics_raw[m], {q}, 3),
xt::view(metrics_qtl[m], xt::all(), xt::all(), xt::all(), i)
))
) << "Failure for (" << all_metrics_p[m] << ") on quantile " << q;
i++;
}
}
}