test_probabilist.cpp

#include <fstream>
#include <vector>
#include <array>
#include <gtest/gtest.h>
#include <xtensor/xtensor.hpp>
#include <xtensor/xmanipulation.hpp>
#include <xtensor/xcsv.hpp>

#include "evalhyd/evalp.hpp"

using namespace xt::placeholders;  // required for `_` to work

TEST(ProbabilistTests, TestBrier)
{
    // read in data
    std::ifstream ifs;
    ifs.open("./data/q_obs.csv");
    xt::xtensor<double, 1> observed = xt::squeeze(xt::load_csv<int>(ifs));
    ifs.close();

    ifs.open("./data/q_prd.csv");
    xt::xtensor<double, 2> predicted = xt::load_csv<double>(ifs);
    ifs.close();

    // 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::view(observed, xt::newaxis(), xt::all()),
                    // shape: (sites [1], lead times [1], members [m], time [t])
                    xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all()),
                    {"BS", "BSS", "BS_CRD", "BS_LBD"},
                    thresholds
            );

    // check results
    // Brier scores
    xt::xtensor<double, 4> bs =
            {{{{0.10615136, 0.07395622, 0.08669186, NAN}}}};
    EXPECT_TRUE(
            xt::sum(xt::isclose(metrics[0], bs, 1e-05, 1e-08, true))
            == xt::xscalar<double>(4)
    );

    // Brier skill scores
    xt::xtensor<double, 4> 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, 5> 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, 5> bs_lbd =
            {{{{{0.012159881, 0.1506234, 0.2446149},
                {0.008031746, 0.1473869, 0.2133114},
                {0.017191279, 0.1048221, 0.1743227},