An error occurred while loading the file. Please try again.
-
Thibault Hallouin authored3eb041a7
// 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