test_probabilist.cpp 30.80 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 <gtest/gtest.h>
#include <xtensor/xtensor.hpp>
#include <xtensor/xview.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",
        "QS", "CRPS",
        "POD", "POFD", "FAR", "CSI", "ROCSS"
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"
    // check results
    // Brier scores
    xt::xtensor<double, 5> bs =
            {{{{{0.10615136, 0.07395622, 0.08669186, NAN}}}}};
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
EXPECT_TRUE( xt::sum(xt::isclose(metrics[0], bs, 1e-05, 1e-08, true)) == xt::xscalar<double>(4) ); // Brier skill scores xt::xtensor<double, 5> bss = {{{{{0.5705594, 0.6661165, 0.5635126, NAN}}}}}; EXPECT_TRUE( xt::sum(xt::isclose(metrics[1], bss, 1e-05, 1e-08, true)) == xt::xscalar<double>(4) ); // 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::sum(xt::isclose(metrics[2], bs_crd, 1e-05, 1e-08, true)) == xt::xscalar<double>(12) ); // 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::sum(xt::isclose(metrics[3], bs_lbd, 1e-05, 1e-08, true)) == xt::xscalar<double>(12) ); } 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::allclose(metrics[0], qs)); // Continuous ranked probability scores