test_probabilist.cpp 36.56 KiB
// 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 <string>
#include <unordered_map>
#include <gtest/gtest.h>
#include <xtl/xoptional.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xarray.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"
#endif
using namespace xt::placeholders;  // required for `_` to work
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"
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
    std::ifstream ifs;
    ifs.open(EVALHYD_DATA_DIR "/data/q_obs.csv");
    xt::xtensor<double, 1> observed = xt::squeeze(xt::load_csv<double>(ifs));
    ifs.close();
    ifs.open(EVALHYD_DATA_DIR "/data/q_prd.csv");
    xt::xtensor<double, 2> predicted = xt::load_csv<double>(ifs);
    ifs.close();
    return std::make_tuple(observed, predicted);
std::unordered_map<std::string, xt::xarray<double>> load_expected_p()
    // read in expected results
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
std::ifstream ifs; std::unordered_map<std::string, xt::xarray<double>> expected; for (const auto& metric : all_metrics_p_) { ifs.open(EVALHYD_DATA_DIR "/expected/evalp/" + metric + ".csv"); expected[metric] = xt::view( xt::squeeze(xt::load_csv<double>(ifs)), xt::newaxis(), xt::newaxis(), xt::newaxis(), xt::newaxis(), xt::all() ); ifs.close(); } return expected; } TEST(ProbabilistTests, TestBrier) { // 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 xt::xtensor<double, 2> thresholds = {{690, 534, 445, NAN}}; std::vector<std::string> metrics = {"BS", "BSS", "BS_CRD", "BS_LBD", "REL_DIAG", "CRPS_FROM_BS"}; 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::all(), xt::all())), metrics, thresholds, "high" ); // check results for (std::size_t m = 0; m < metrics.size(); m++) { if ( metrics[m] == "REL_DIAG" ) { // /!\ stacked-up thresholds in CSV file because 7D metric, // so need to resize array expected[metrics[m]].resize( {std::size_t {1}, std::size_t {1}, std::size_t {1}, std::size_t {1}, thresholds.shape(1), predicted.shape(0) + 1, std::size_t {3}} ); } EXPECT_TRUE(xt::all(xt::isclose( results[m], expected[metrics[m]], 1e-05, 1e-08, true ))) << "Failure for (" << metrics[m] << ")"; } } TEST(ProbabilistTests, TestCDF) { // 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();
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
// compute scores std::vector<std::string> metrics = {"CRPS_FROM_ECDF"}; 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::all(), xt::all())), metrics ); // 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, TestQuantiles) { // 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 = {"QS", "CRPS_FROM_QS"}; 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::all(), xt::all())), metrics ); // 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, TestContingency) { // 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 xt::xtensor<double, 2> thresholds = {{690, 534, 445, NAN}}; std::vector<std::string> metrics = {"CONT_TBL", "POD", "POFD", "FAR", "CSI", "ROCSS"}; std::vector<xt::xarray<double>> results = evalhyd::evalp( // shape: (sites [1], time [t])
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
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())), metrics, thresholds, "low" ); // check results for (std::size_t m = 0; m < metrics.size(); m++) { if (metrics[m] == "CONT_TBL") { // /!\ stacked-up thresholds and cells in CSV file because 7D metric, // so need to resize array accordingly expected[metrics[m]].resize( {std::size_t {1}, std::size_t {1}, std::size_t {1}, std::size_t {1}, predicted.shape(0) + 1, thresholds.shape(1), std::size_t {4}} ); } EXPECT_TRUE(xt::all(xt::isclose( results[m], expected[metrics[m]], 1e-05, 1e-08, true ))) << "Failure for (" << metrics[m] << ")"; } } TEST(ProbabilistTests, TestRanks) { // 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(); std::vector<std::string> metrics = {"RANK_HIST", "DS", "AS"}; // compute scores 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::all(), xt::all())), {"RANK_HIST", "DS", "AS"}, 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 {}, // dts 7 // seed ); // 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, TestIntervals) { // read in data xt::xtensor<double, 1> observed;
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
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", "AW", "AWN", "WS"}; 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::all(), xt::all())), {"CR", "AW", "AWN", "WS"}, xt::xtensor<double, 2>({}), "", // events {30., 80.} // c_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, 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 xt::xtensor<double, 1> observed;
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
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 = {"ES"}; xt::xtensor<double, 2> obs = xt::repeat( xt::view(observed, xt::newaxis(), xt::all()), 5, 0 ); xt::xtensor<double, 4> prd = xt::repeat( xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all()), 5, 0 ); std::vector<xt::xarray<double>> results = evalhyd::evalp( // shape: (sites [5], time [t]) obs, // shape: (sites [5], lead times [1], members [m], time [t]) prd, metrics ); // 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, 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; // 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, {}, // q_lvl // 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])
421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490
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 _________________________________________ // 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, {}, // q_lvl 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(
491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560
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", confidence_levels, {}, // q_lvl 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++) { // ---------------------------------------------------------------------
561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630
// /!\ 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, {}, // q_lvl 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 ); // 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],
631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700
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, 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,
701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770
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 "/data/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 "/data/q_obs_1yr.csv"); xt::xtensor<double, 1> observed = xt::squeeze(xt::load_csv<double>(ifs, ',', 1)); ifs.close(); ifs.open(EVALHYD_DATA_DIR "/data/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}}; 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, {}, // q_lvl 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
771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840
// 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 "/data/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 "/data/q_obs_1yr.csv"); xt::xtensor<double, 1> observed = xt::squeeze(xt::load_csv<double>(ifs, ',', 1)); ifs.close(); ifs.open(EVALHYD_DATA_DIR "/data/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())),
841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())), all_metrics_p, thresholds, "high", // events confidence_levels, {}, // q_lvl 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, {}, // q_lvl 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) )) ) << "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(
911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980
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, {}, // q_lvl 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++; } } } TEST(ProbabilistTests, TestCompleteness) { std::vector<std::string> diags = {"completeness"}; // compute metrics on series with NaN xt::xtensor<double, 4> prd = {{ // leadtime 1 {{ 5.3, NAN, 5.7, 2.3, 3.3, NAN }, { 4.3, NAN, 4.7, 4.3, 3.4, NAN }, { 5.3, NAN, 5.7, 2.3, 3.8, NAN }}, // leadtime 2 {{ NAN, 4.2, 5.7, 2.3, 3.1, 4.1 }, { NAN, 4.2, 4.7, 4.3, 3.3, 2.8 }, { NAN, 5.2, 5.7, 2.3, 3.9, 3.5 }} }}; xt::xtensor<double, 2> obs = {{ 4.7, 4.3, NAN, 2.7, 4.1, 5.0 }}; xt::xtensor<bool, 4> msk = {{ // leadtime 1 {{ true, true, true, false, true, true }, { true, true, true, true, true, true }},
9819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018
// leadtime 2 {{ true, true, true, true, true, false }, { true, true, true, true, true, true }}, }}; std::vector<xt::xarray<double>> results = evalhyd::evalp( obs, prd, std::vector<std::string> {}, // metrics 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 {}, // dts xtl::missing<const int>(), // seed diags ); // check that numerical results are identical xt::xtensor<double, 4> expected = {{ // leadtime 1 {{ 2. }, { 3. }}, // leadtime 2 {{ 3. }, { 4. }}, }}; EXPECT_TRUE( xt::all(xt::isclose(results[0], expected, 1e-05, 1e-08, true)) ); }