From 1b82acb896df08217bf377811e26a70bcf93b376 Mon Sep 17 00:00:00 2001
From: Thibault Hallouin <thibault.hallouin@inrae.fr>
Date: Wed, 1 Feb 2023 18:27:06 +0100
Subject: [PATCH] add "mean+std" and "quantiles" as summaries for bootstrap
sampling
---
include/evalhyd/detail/uncertainty.hpp | 154 ++++++++++++++++++++-----
include/evalhyd/detail/utils.hpp | 3 +-
include/evalhyd/evald.hpp | 8 +-
include/evalhyd/evalp.hpp | 40 +++----
tests/test_determinist.cpp | 113 ++++++++++++++++++
tests/test_probabilist.cpp | 145 ++++++++++++++++++++++-
6 files changed, 408 insertions(+), 55 deletions(-)
diff --git a/include/evalhyd/detail/uncertainty.hpp b/include/evalhyd/detail/uncertainty.hpp
index de7d746..d3b086b 100644
--- a/include/evalhyd/detail/uncertainty.hpp
+++ b/include/evalhyd/detail/uncertainty.hpp
@@ -16,7 +16,6 @@
#include <xtensor/xtensor.hpp>
#include <xtensor/xadapt.hpp>
#include <xtensor/xrandom.hpp>
-#include <xtensor/xio.hpp>
#include "maths.hpp"
@@ -154,41 +153,70 @@ namespace evalhyd
return samples;
}
-
- inline auto summarise(const xt::xarray<double>& values, int summary)
+
+ inline auto summarise_d(const xt::xarray<double>& values, int summary)
{
- // TODO: wait for xt::quantile to be available
- // or implement it internally for n-dim expressions
+ // define axis along which samples are
+ std::size_t axis = 2;
+
+ // determine shape for output values
+ std::vector<std::size_t> shp;
+ std::size_t i = 0;
+ for (auto a : values.shape())
+ {
+ if (i != axis)
+ {
+ shp.push_back(a);
+ }
+ else
+ {
+ if (summary == 1)
+ {
+ shp.push_back(2);
+ }
+ else if (summary == 2)
+ {
+ shp.push_back(7);
+ }
+ }
+ i++;
+ }
+
// summary 2: series of quantiles across samples
if (summary == 2)
{
-// // adjust last axis size (from samples to number of statistics)
-// auto s = values.shape();
-// s.pop_back();
-// s.push_back(7);
-// xt::xarray<double> v = xt::zeros<double>(s);
-// // quantiles
-// xt::view(v, xt::all()) =
-// xt::quantile(values, {0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95}, -1);
-// return v;
- return values;
+ xt::xarray<double> v = xt::zeros<double>(shp);
+
+ // compute quantiles
+ auto quantiles = xt::quantile(
+ values,
+ {0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95},
+ axis
+ );
+
+ // transfer quantiles into correct axis
+ // (since xt::quantile puts the quantiles on the first axis)
+ for (std::size_t q = 0; q < 7; q++)
+ {
+ xt::view(v, xt::all(), xt::all(), q) =
+ xt::view(quantiles, q);
+ }
+
+ return v;
}
- // TODO: wait for xt::stddev to be fixed for rtensor
- // or implement it internally for n-dim expressions
// summary 1: mean and standard deviation across samples
else if (summary == 1)
{
-// // adjust last axis size (from samples to number of statistics)
-// auto s = values.shape();
-// s.pop_back();
-// s.push_back(2);
-// xt::xarray<double> v = xt::zeros<double>(s);
-// // mean
-// xt::strided_view(s, xt::ellipsis(), 0) = xt::mean(values, -1);
-// // standard deviation
-// xt::strided_view(s, xt::ellipsis(), 1) = xt::stddev(values, -1);
-// return v;
- return values;
+ xt::xarray<double> v = xt::zeros<double>(shp);
+
+ // compute mean
+ xt::view(v, xt::all(), xt::all(), 0) =
+ xt::mean(values, 2);
+ // compute standard deviation
+ xt::view(v, xt::all(), xt::all(), 1) =
+ xt::stddev(values, 2);
+
+ return v;
}
// summary 0: raw (keep all samples)
else
@@ -197,6 +225,76 @@ namespace evalhyd
}
}
+ inline auto summarise_p(const xt::xarray<double>& values, int summary)
+ {
+ // define axis along which samples are
+ std::size_t axis = 3;
+
+ // determine shape for output values
+ std::vector<std::size_t> shp;
+ std::size_t i = 0;
+ for (auto a : values.shape())
+ {
+ if (i != axis)
+ {
+ shp.push_back(a);
+ }
+ else
+ {
+ if (summary == 1)
+ {
+ shp.push_back(2);
+ }
+ else if (summary == 2)
+ {
+ shp.push_back(7);
+ }
+ }
+ i++;
+ }
+
+ // summary 2: series of quantiles across samples
+ if (summary == 2)
+ {
+ xt::xarray<double> v = xt::zeros<double>(shp);
+
+ // compute quantiles
+ auto quantiles = xt::quantile(
+ values,
+ {0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95},
+ axis
+ );
+
+ // transfer quantiles into correct axis
+ // (since xt::quantile puts the quantiles on the first axis)
+ for (std::size_t q = 0; q < 7; q++)
+ {
+ xt::view(v, xt::all(), xt::all(), xt::all(), q) =
+ xt::view(quantiles, q);
+ }
+
+ return v;
+ }
+ // summary 1: mean and standard deviation across samples
+ else if (summary == 1)
+ {
+ xt::xarray<double> v = xt::zeros<double>(shp);
+
+ // compute mean
+ xt::view(v, xt::all(), xt::all(), xt::all(), 0) =
+ xt::mean(values, axis);
+ // compute standard deviation
+ xt::view(v, xt::all(), xt::all(), xt::all(), 1) =
+ xt::stddev(values, axis);
+
+ return v;
+ }
+ // summary 0: raw (keep all samples)
+ else
+ {
+ return values;
+ }
+ }
}
}
diff --git a/include/evalhyd/detail/utils.hpp b/include/evalhyd/detail/utils.hpp
index 839b49d..94eae53 100644
--- a/include/evalhyd/detail/utils.hpp
+++ b/include/evalhyd/detail/utils.hpp
@@ -88,8 +88,7 @@ namespace evalhyd
);
}
auto summary = bootstrap.find("summary")->second;
- // TODO: change upper bound when mean+stddev and quantiles implemented
- if ((summary < 0) || (summary > 0))
+ if ((summary < 0) || (summary > 2))
{
throw std::runtime_error(
"invalid value for bootstrap summary"
diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp
index f6539c3..6826191 100644
--- a/include/evalhyd/evald.hpp
+++ b/include/evalhyd/evald.hpp
@@ -437,25 +437,25 @@ namespace evalhyd
if ( metric == "RMSE" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_RMSE(), summary)
+ uncertainty::summarise_d(evaluator.get_RMSE(), summary)
);
}
else if ( metric == "NSE" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_NSE(), summary)
+ uncertainty::summarise_d(evaluator.get_NSE(), summary)
);
}
else if ( metric == "KGE" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_KGE(), summary)
+ uncertainty::summarise_d(evaluator.get_KGE(), summary)
);
}
else if ( metric == "KGEPRIME" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_KGEPRIME(), summary)
+ uncertainty::summarise_d(evaluator.get_KGEPRIME(), summary)
);
}
}
diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp
index 8cf1fa7..80a14ce 100644
--- a/include/evalhyd/evalp.hpp
+++ b/include/evalhyd/evalp.hpp
@@ -400,121 +400,121 @@ namespace evalhyd
if ( metric == "BS" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_BS(), summary)
+ uncertainty::summarise_p(evaluator.get_BS(), summary)
);
}
else if ( metric == "BSS" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_BSS(), summary)
+ uncertainty::summarise_p(evaluator.get_BSS(), summary)
);
}
else if ( metric == "BS_CRD" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_BS_CRD(), summary)
+ uncertainty::summarise_p(evaluator.get_BS_CRD(), summary)
);
}
else if ( metric == "BS_LBD" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_BS_LBD(), summary)
+ uncertainty::summarise_p(evaluator.get_BS_LBD(), summary)
);
}
else if ( metric == "QS" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_QS(), summary)
+ uncertainty::summarise_p(evaluator.get_QS(), summary)
);
}
else if ( metric == "CRPS" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_CRPS(), summary)
+ uncertainty::summarise_p(evaluator.get_CRPS(), summary)
);
}
else if ( metric == "POD" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_POD(), summary)
+ uncertainty::summarise_p(evaluator.get_POD(), summary)
);
}
else if ( metric == "POFD" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_POFD(), summary)
+ uncertainty::summarise_p(evaluator.get_POFD(), summary)
);
}
else if ( metric == "FAR" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_FAR(), summary)
+ uncertainty::summarise_p(evaluator.get_FAR(), summary)
);
}
else if ( metric == "CSI" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_CSI(), summary)
+ uncertainty::summarise_p(evaluator.get_CSI(), summary)
);
}
else if ( metric == "ROCSS" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_ROCSS(), summary)
+ uncertainty::summarise_p(evaluator.get_ROCSS(), summary)
);
}
else if ( metric == "RANK_HIST" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_RANK_HIST(), summary)
+ uncertainty::summarise_p(evaluator.get_RANK_HIST(), summary)
);
}
else if ( metric == "DS" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_DS(), summary)
+ uncertainty::summarise_p(evaluator.get_DS(), summary)
);
}
else if ( metric == "AS" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_AS(), summary)
+ uncertainty::summarise_p(evaluator.get_AS(), summary)
);
}
else if ( metric == "CR" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_CR(), summary)
+ uncertainty::summarise_p(evaluator.get_CR(), summary)
);
}
else if ( metric == "AW" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_AW(), summary)
+ uncertainty::summarise_p(evaluator.get_AW(), summary)
);
}
else if ( metric == "AWN" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_AWN(), summary)
+ uncertainty::summarise_p(evaluator.get_AWN(), summary)
);
}
else if ( metric == "AWI" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_AWI(), summary)
+ uncertainty::summarise_p(evaluator.get_AWI(), summary)
);
}
else if ( metric == "WS" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_WS(), summary)
+ uncertainty::summarise_p(evaluator.get_WS(), summary)
);
}
else if ( metric == "WSS" )
{
r.emplace_back(
- uncertainty::summarise(evaluator.get_WSS(), summary)
+ uncertainty::summarise_p(evaluator.get_WSS(), summary)
);
}
}
diff --git a/tests/test_determinist.cpp b/tests/test_determinist.cpp
index 0ddcd18..c37e93f 100644
--- a/tests/test_determinist.cpp
+++ b/tests/test_determinist.cpp
@@ -12,6 +12,8 @@
#include <xtensor/xtensor.hpp>
#include <xtensor/xview.hpp>
#include <xtensor/xmanipulation.hpp>
+#include <xtensor/xmath.hpp>
+#include <xtensor/xsort.hpp>
#include <xtensor/xcsv.hpp>
#include "evalhyd/evald.hpp"
@@ -409,3 +411,114 @@ TEST(DeterministTests, TestBootstrap)
) << "Failure for (" << all_metrics_d[m] << ")";
}
}
+
+TEST(DeterministTests, TestBootstrapSummary)
+{
+ // read in data
+ std::ifstream ifs;
+
+ ifs.open(EVALHYD_DATA_DIR "/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 "/q_obs_1yr.csv");
+ xt::xtensor<double, 1> observed = xt::squeeze(xt::load_csv<double>(ifs, ',', 1));
+ ifs.close();
+
+ ifs.open(EVALHYD_DATA_DIR "/q_prd_1yr.csv");
+ xt::xtensor<double, 2> predicted = xt::load_csv<double>(ifs, ',', 1);
+ ifs.close();
+
+ // compute metrics via bootstrap with raw summary
+ std::unordered_map<std::string, int> bootstrap_0 =
+ {{"n_samples", 10}, {"len_sample", 3}, {"summary", 0}};
+
+ std::vector<xt::xarray<double>> metrics_raw =
+ evalhyd::evald(
+ xt::eval(xt::view(observed, xt::newaxis(), xt::all())),
+ predicted,
+ all_metrics_d,
+ {}, // transform
+ {}, // exponent
+ {}, // epsilon
+ xt::xtensor<bool, 3>({}), // 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::evald(
+ xt::eval(xt::view(observed, xt::newaxis(), xt::all())),
+ predicted,
+ all_metrics_d,
+ {}, // transform
+ {}, // exponent
+ {}, // epsilon
+ xt::xtensor<bool, 3>({}), // 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_d.size(); m++)
+ {
+ // mean
+ EXPECT_TRUE(
+ xt::all(xt::isclose(
+ xt::mean(metrics_raw[m], 2),
+ xt::view(metrics_mas[m], xt::all(), xt::all(), 0)
+ ))
+ ) << "Failure for (" << all_metrics_d[m] << ") on mean";
+ // standard deviation
+ EXPECT_TRUE(
+ xt::all(xt::isclose(
+ xt::stddev(metrics_raw[m], 2),
+ xt::view(metrics_mas[m], xt::all(), xt::all(), 1)
+ ))
+ ) << "Failure for (" << all_metrics_d[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::evald(
+ xt::eval(xt::view(observed, xt::newaxis(), xt::all())),
+ predicted,
+ all_metrics_d,
+ {}, // transform
+ {}, // exponent
+ {}, // epsilon
+ xt::xtensor<bool, 3>({}), // 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_d.size(); m++)
+ {
+ // 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}, 2),
+ xt::view(metrics_qtl[m], xt::all(), xt::all(), i)
+ ))
+ ) << "Failure for (" << all_metrics_d[m] << ") on quantile " << q;
+ i++;
+ }
+ }
+}
diff --git a/tests/test_probabilist.cpp b/tests/test_probabilist.cpp
index e0d8dd4..3196100 100644
--- a/tests/test_probabilist.cpp
+++ b/tests/test_probabilist.cpp
@@ -12,6 +12,7 @@
#include <xtl/xoptional.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xview.hpp>
+#include <xtensor/xmath.hpp>
#include <xtensor/xsort.hpp>
#include <xtensor/xmanipulation.hpp>
#include <xtensor/xcsv.hpp>
@@ -965,4 +966,146 @@ TEST(ProbabilistTests, TestBootstrap)
))
) << "Failure for (" << all_metrics_p[m] << ")";
}
-}
\ No newline at end of file
+}
+
+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 "/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 "/q_obs_1yr.csv");
+ xt::xtensor<double, 1> observed = xt::squeeze(xt::load_csv<double>(ifs, ',', 1));
+ ifs.close();
+
+ ifs.open(EVALHYD_DATA_DIR "/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())),
+ xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
+ all_metrics_p,
+ thresholds,
+ "high", // events
+ confidence_levels,
+ 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,
+ 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(
+ 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,
+ 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++;
+ }
+ }
+}
--
GitLab