Skip to content
GitLab
Select Git revision
  • dev
  • 11-check-the-computation-of-the-ranks-for-ties
  • main default
  • add-determinist-to-evalp
  • v0.1.2
  • v0.1.1
  • v0.1.0
  • v0.0.1
Find file BlameHistoryPermalink
  • Thibault Hallouin's avatar
    harmonise definitions for above/below threshold events · 4410a74f
    Thibault Hallouin authored 
    Low flow events were defined as the complement of high flow events,
meaning that Brier scores were symmetric. But it is more consistent to
include the threshold value in both the definitions of low flow and
high flow events since we either study one or the other, and rarely
both at the same time, especially not using the same threshold values.
Plus, the choice of including the threshold in high flow events rather
than low flow ones was really arbitrary.

Note, only the unit tests for contingency table-based metrics are
impacted because threshold-based metrics are using evants="high".
    4410a74f
test_probabilist.cpp 45.68 KiB
1 
// Copyright (c) 2023, INRAE.
2 
// Distributed under the terms of the GPL-3 Licence.
3 
// The full licence is in the file LICENCE, distributed with this software.
4
5 
#include <fstream>
6 
#include <vector>
7 
#include <tuple>
8 
#include <array>
9
10 
#include <gtest/gtest.h>
11
12 
#include <xtl/xoptional.hpp>
13 
#include <xtensor/xtensor.hpp>
14 
#include <xtensor/xview.hpp>
15 
#include <xtensor/xmath.hpp>
16 
#include <xtensor/xsort.hpp>
17 
#include <xtensor/xmanipulation.hpp>
18 
#include <xtensor/xcsv.hpp>
19
20 
#include "evalhyd/evalp.hpp"
21
22 
#ifndef EVALHYD_DATA_DIR
23 
#error "need to define data directory"
24 
#endif
25
26 
using namespace xt::placeholders;  // required for `_` to work
27
28
29 
std::vector<std::string> all_metrics_p = {
30 
        "BS", "BSS", "BS_CRD", "BS_LBD",
31 
        "QS", "CRPS",
32 
        "POD", "POFD", "FAR", "CSI", "ROCSS",
33 
        "RANK_HIST", "DS", "AS",
34 
        "CR", "AW", "AWN", "AWI", "WS", "WSS"
35 
};
36
37 
std::tuple<xt::xtensor<double, 1>, xt::xtensor<double, 2>> load_data_p()
38 
{
39 
    // read in data
40 
    std::ifstream ifs;
41 
    ifs.open(EVALHYD_DATA_DIR "/q_obs.csv");
42 
    xt::xtensor<double, 1> observed = xt::squeeze(xt::load_csv<int>(ifs));
43 
    ifs.close();
44
45 
    ifs.open(EVALHYD_DATA_DIR "/q_prd.csv");
46 
    xt::xtensor<double, 2> predicted = xt::load_csv<double>(ifs);
47 
    ifs.close();
48
49 
    return std::make_tuple(observed, predicted);
50 
}
51
52 
TEST(ProbabilistTests, TestBrier)
53 
{
54 
    // read in data
55 
    xt::xtensor<double, 1> observed;
56 
    xt::xtensor<double, 2> predicted;
57 
    std::tie(observed, predicted) = load_data_p();
58
59 
    // compute scores
60 
    xt::xtensor<double, 2> thresholds = {{690, 534, 445, NAN}};
61
62 
    std::vector<xt::xarray<double>> metrics =
63 
            evalhyd::evalp(
64 
                    // shape: (sites [1], time [t])
65 
                    xt::eval(xt::view(observed, xt::newaxis(), xt::all())),
66 
                    // shape: (sites [1], lead times [1], members [m], time [t])
67 
                    xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
68 
                    {"BS", "BSS", "BS_CRD", "BS_LBD"},
69 
                    thresholds,
70 
                    "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)));