diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp
index 236a8768f1f34d3462bc84d11815779ca7ac130f..9dd828578453084b4b1d3ec07d525844ef953095 100644
--- a/include/evalhyd/evalp.hpp
+++ b/include/evalhyd/evalp.hpp
@@ -23,16 +23,16 @@ namespace evalhyd
///
/// q_obs: :cpp:`xt::xtensor<double, 2>`
/// Streamflow observations.
- /// shape: (1, time)
+ /// shape: (sites, time)
///
- /// q_prd: :cpp:`xt::xtensor<double, 2>`
+ /// q_prd: :cpp:`xt::xtensor<double, 4>`
/// Streamflow predictions.
- /// shape: (members, time)
+ /// shape: (sites, lead times, members, time)
///
/// metrics: :cpp:`std::vector<std::string>`
/// The sequence of evaluation metrics to be computed.
///
- /// q_thr: :cpp:`xt::xtensor<double, 1>`, optional
+ /// q_thr: :cpp:`xt::xtensor<double, 2>`, optional
/// Streamflow exceedance threshold(s).
/// shape: (thresholds,)
///
@@ -43,11 +43,11 @@ namespace evalhyd
/// performed and only one set of metrics is returned corresponding
/// to the whole time series. If provided, as many sets of metrics are
/// returned as they are masks provided.
- /// shape: (1+, time)
+ /// shape: (subsets, time)
///
/// :Returns:
///
- /// :cpp:`std::vector<xt::xarray<double>>`
+ /// :cpp:`std::vector<xt::xtensor<double, 5>>`
/// The sequence of evaluation metrics computed
/// in the same order as given in *metrics*.
///
@@ -69,7 +69,7 @@ namespace evalhyd
/// \endrst
std::vector<xt::xarray<double>> evalp(
const xt::xtensor<double, 2>& q_obs,
- const xt::xtensor<double, 2>& q_prd,
+ const xt::xtensor<double, 4>& q_prd,
const std::vector<std::string>& metrics,
const xt::xtensor<double, 1>& q_thr = {},
const xt::xtensor<bool, 2>& t_msk = {}
@@ -84,8 +84,22 @@ namespace evalhyd
// check that optional parameters are given as arguments
eh::utils::check_optionals(metrics, q_thr);
- // instantiate probabilist evaluator
- eh::probabilist::Evaluator evaluator(q_obs, q_prd, q_thr, t_msk);
+ // retrieve dimensions
+ std::size_t n_sit = q_prd.shape(0);
+ std::size_t n_ltm = q_prd.shape(1);
+ std::size_t n_mbr = q_prd.shape(2);
+ std::size_t n_thr = q_thr.size();
+ std::size_t n_msk = t_msk.size() < 1 ? 1 : t_msk.shape(0);
+
+ // register metrics number of dimensions
+ std::unordered_map<std::string, std::vector<std::size_t>> dim;
+
+ dim["BS"] = {n_sit, n_ltm, n_msk, n_thr};
+ dim["BSS"] = {n_sit, n_ltm, n_msk, n_thr};
+ dim["BS_CRD"] = {n_sit, n_ltm, n_msk, n_thr, 3};
+ dim["BS_LBD"] = {n_sit, n_ltm, n_msk, n_thr, 3};
+ dim["QS"] = {n_sit, n_ltm, n_msk, n_mbr};
+ dim["CRPS"] = {n_sit, n_ltm, n_msk};
// declare maps for memoisation purposes
std::unordered_map<std::string, std::vector<std::string>> elt;
@@ -110,76 +124,105 @@ namespace evalhyd
eh::utils::find_requirements(metrics, elt, dep, req_elt, req_dep);
- // pre-compute required elt
- for (const auto& element : req_elt)
- {
- if ( element == "o_k" )
- evaluator.calc_o_k();
- else if ( element == "bar_o" )
- evaluator.calc_bar_o();
- else if ( element == "y_k" )
- evaluator.calc_y_k();
- else if ( element == "q_qnt" )
- evaluator.calc_q_qnt();
- }
+ // initialise data structure for outputs
+ std::vector<xt::xarray<double>> r;
+ for (const auto& metric : metrics)
+ r.emplace_back(xt::zeros<double>(dim[metric]));
- // pre-compute required dep
- for (const auto& dependency : req_dep)
+ // compute variables one site at a time to minimise memory imprint
+ for (int s = 0; s < n_sit; s++)
+ // compute variables one lead time at a time to minimise memory imprint
+ for (int l = 0; l < n_ltm; l++)
{
- if ( dependency == "bs" )
- evaluator.calc_bs();
- else if ( dependency == "qs" )
- evaluator.calc_qs();
- else if ( dependency == "crps" )
- evaluator.calc_crps();
- }
+ // instantiate probabilist evaluator
+ const auto q_obs_v = xt::view(q_obs, s, xt::all());
+ const auto q_prd_v = xt::view(q_prd, s, l, xt::all(), xt::all());
- // retrieve or compute requested metrics
- std::vector<xt::xarray<double>> r;
+ eh::probabilist::Evaluator evaluator(q_obs_v, q_prd_v, q_thr, t_msk);
- for (const auto& metric : metrics)
- {
- if ( metric == "BS" )
+ // pre-compute required elt
+ for (const auto& element : req_elt)
{
- if (std::find(req_dep.begin(), req_dep.end(), metric)
- == req_dep.end())
- evaluator.calc_BS();
- r.emplace_back(xt::squeeze(evaluator.BS));
+ if ( element == "o_k" )
+ evaluator.calc_o_k();
+ else if ( element == "bar_o" )
+ evaluator.calc_bar_o();
+ else if ( element == "y_k" )
+ evaluator.calc_y_k();
+ else if ( element == "q_qnt" )
+ evaluator.calc_q_qnt();
}
- else if ( metric == "BSS" )
+
+ // pre-compute required dep
+ for (const auto& dependency : req_dep)
{
- if (std::find(req_dep.begin(), req_dep.end(), metric)
- == req_dep.end())
- evaluator.calc_BSS();
- r.emplace_back(xt::squeeze(evaluator.BSS));
+ if ( dependency == "bs" )
+ evaluator.calc_bs();
+ else if ( dependency == "qs" )
+ evaluator.calc_qs();
+ else if ( dependency == "crps" )
+ evaluator.calc_crps();
}
- else if ( metric == "BS_CRD" )
+
+ // retrieve or compute requested metrics
+ for (int m = 0; m < metrics.size(); m++)
{
- if (std::find(req_dep.begin(), req_dep.end(), metric)
+ const auto& metric = metrics[m];
+
+ if ( metric == "BS" )
+ {
+ if (std::find(req_dep.begin(), req_dep.end(), metric)
== req_dep.end())
- evaluator.calc_BS_CRD();
- r.emplace_back(xt::squeeze(evaluator.BS_CRD));
- }
- else if ( metric == "BS_LBD" )
- {
- if (std::find(req_dep.begin(), req_dep.end(), metric)
+ evaluator.calc_BS();
+ // (sites, lead times, subsets, thresholds)
+ xt::view(r[m], s, l, xt::all(), xt::all()) =
+ evaluator.BS;
+ }
+ else if ( metric == "BSS" )
+ {
+ if (std::find(req_dep.begin(), req_dep.end(), metric)
== req_dep.end())
- evaluator.calc_BS_LBD();
- r.emplace_back(xt::squeeze(evaluator.BS_LBD));
- }
- else if ( metric == "QS" )
- {
- if (std::find(req_dep.begin(), req_dep.end(), metric)
+ evaluator.calc_BSS();
+ // (sites, lead times, subsets, thresholds)
+ xt::view(r[m], s, l, xt::all(), xt::all()) =
+ evaluator.BSS;
+ }
+ else if ( metric == "BS_CRD" )
+ {
+ if (std::find(req_dep.begin(), req_dep.end(), metric)
== req_dep.end())
- evaluator.calc_QS();
- r.emplace_back(xt::squeeze(evaluator.QS));
- }
- else if ( metric == "CRPS" )
- {
- if (std::find(req_dep.begin(), req_dep.end(), metric)
+ evaluator.calc_BS_CRD();
+ // (sites, lead times, subsets, thresholds, components)
+ xt::view(r[m], s, l, xt::all(), xt::all(), xt::all()) =
+ evaluator.BS_CRD;
+ }
+ else if ( metric == "BS_LBD" )
+ {
+ if (std::find(req_dep.begin(), req_dep.end(), metric)
+ == req_dep.end())
+ evaluator.calc_BS_LBD();
+ // (sites, lead times, subsets, thresholds, components)
+ xt::view(r[m], s, l, xt::all(), xt::all(), xt::all()) =
+ evaluator.BS_LBD;
+ }
+ else if ( metric == "QS" )
+ {
+ if (std::find(req_dep.begin(), req_dep.end(), metric)
+ == req_dep.end())
+ evaluator.calc_QS();
+ // (sites, lead times, subsets, quantiles)
+ xt::view(r[m], s, l, xt::all(), xt::all()) =
+ evaluator.QS;
+ }
+ else if ( metric == "CRPS" )
+ {
+ if (std::find(req_dep.begin(), req_dep.end(), metric)
== req_dep.end())
- evaluator.calc_CRPS();
- r.emplace_back(xt::squeeze(evaluator.CRPS));
+ evaluator.calc_CRPS();
+ // (sites, lead times, subsets)
+ xt::view(r[m], s, l, xt::all()) =
+ evaluator.CRPS;
+ }
}
}
diff --git a/src/probabilist/evaluator.h b/src/probabilist/evaluator.h
index 26f0223430efc75759aa13eaa3461a4044e86fa8..0a622cef15e98c24b4dd38e4433974830ec60656 100644
--- a/src/probabilist/evaluator.h
+++ b/src/probabilist/evaluator.h
@@ -2,6 +2,25 @@
#define EVALHYD_PROBABILIST_EVALUATOR_H
#include <xtensor/xtensor.hpp>
+#include <xtensor/xview.hpp>
+
+using view1d_xtensor2d_type = decltype(
+ xt::view(
+ std::declval<const xt::xtensor<double, 2>&>(),
+ std::declval<int>(),
+ xt::all()
+ )
+);
+
+using view2d_xtensor4d_type = decltype(
+ xt::view(
+ std::declval<const xt::xtensor<double, 4>&>(),
+ std::declval<int>(),
+ std::declval<int>(),
+ xt::all(),
+ xt::all()
+ )
+);
namespace evalhyd
{
@@ -11,8 +30,8 @@ namespace evalhyd
{
private:
// members for input data
- const xt::xtensor<double, 2>& q_obs;
- const xt::xtensor<double, 2>& q_prd;
+ const view1d_xtensor2d_type& q_obs;
+ const view2d_xtensor4d_type& q_prd;
const xt::xtensor<double, 1>& q_thr;
xt::xtensor<bool, 2> t_msk;
@@ -30,8 +49,8 @@ namespace evalhyd
public:
// constructor method
- Evaluator(const xt::xtensor<double, 2>& obs,
- const xt::xtensor<double, 2>& prd,
+ Evaluator(const view1d_xtensor2d_type& obs,
+ const view2d_xtensor4d_type& prd,
const xt::xtensor<double, 1>& thr,
const xt::xtensor<bool, 2>& msk) :
q_obs{obs}, q_prd{prd}, q_thr{thr}, t_msk(msk)
@@ -39,13 +58,13 @@ namespace evalhyd
// initialise a mask if none provided
// (corresponding to no temporal subset)
if (msk.size() < 1)
- t_msk = xt::ones<bool>(obs.shape());
+ t_msk = xt::ones<bool>({std::size_t {1}, q_obs.size()});
// determine size for recurring dimensions
- n = obs.shape(1);
+ n = q_obs.size();
n_msk = t_msk.shape(0);
- n_mbr = prd.shape(0);
- n_thr = thr.size();
+ n_mbr = q_prd.shape(0);
+ n_thr = q_thr.size();
};
// members for intermediate evaluation metrics
diff --git a/src/probabilist/evaluator_brier.cpp b/src/probabilist/evaluator_brier.cpp
index dcc8b7e17a6d63d20f9e4a7e8c47cde04b7db6cb..c250f17d870e4d722be8fb57e322e78cfa175ec4 100644
--- a/src/probabilist/evaluator_brier.cpp
+++ b/src/probabilist/evaluator_brier.cpp
@@ -26,7 +26,7 @@ namespace evalhyd
// Event probability forecast.
// shape: (thresholds, time)
// \assign bs:
- // 2D array of Brier score for each threshold for each time step.
+ // Brier score for each threshold for each time step.
// shape: (thresholds, time)
void Evaluator::calc_bs()
{
diff --git a/src/probabilist/evaluator_elements.cpp b/src/probabilist/evaluator_elements.cpp
index 5f354e30634b33d1e7687adad95446693b953312..2df49be8cf19c1e496b089ac40df26108425b6f7 100644
--- a/src/probabilist/evaluator_elements.cpp
+++ b/src/probabilist/evaluator_elements.cpp
@@ -12,7 +12,7 @@ namespace evalhyd
//
// \require q_obs:
// Streamflow observations.
- // shape: (1, time)
+ // shape: (time,)
// \require q_thr:
// Streamflow exceedance threshold(s).
// shape: (thresholds,)
@@ -22,7 +22,7 @@ namespace evalhyd
void Evaluator::calc_o_k()
{
// determine observed realisation of threshold(s) exceedance
- o_k = xt::squeeze(is_above_threshold(q_obs, q_thr));
+ o_k = q_obs >= xt::view(q_thr, xt::all(), xt::newaxis());
}
// Determine mean observed realisation of threshold(s) exceedance.
@@ -65,7 +65,8 @@ namespace evalhyd
{
// determine if members have exceeded threshold(s)
xt::xtensor<double, 3> e_frc =
- is_above_threshold(q_prd, q_thr);
+ q_prd
+ >= xt::view(q_thr, xt::all(), xt::newaxis(), xt::newaxis());
// calculate how many members have exceeded threshold(s)
xt::xtensor<double, 2> n_frc = xt::sum(e_frc, 1);
diff --git a/src/probabilist/evaluator_quantiles.cpp b/src/probabilist/evaluator_quantiles.cpp
index 726d14d8f0a0a8ac0bcd2f36c485bb302e1146af..a92c76a5dd4d60dd252ed0fbd59e4e973bd2cb6d 100644
--- a/src/probabilist/evaluator_quantiles.cpp
+++ b/src/probabilist/evaluator_quantiles.cpp
@@ -15,7 +15,7 @@ namespace evalhyd
//
// \require q_obs:
// Streamflow observations.
- // shape: (1, time)
+ // shape: (time,)
// \require q_qnt:
// Streamflow quantiles.
// shape: (quantiles, time)
@@ -50,11 +50,11 @@ namespace evalhyd
// shape: (quantiles, time)
// \assign QS:
// Quantile scores.
- // shape: (subsets, quantiles, 1)
+ // shape: (subsets, quantiles)
void Evaluator::calc_QS()
{
// initialise output variable
- // shape: (subsets, quantiles, 1)
+ // shape: (subsets, quantiles)
QS = xt::zeros<double>({n_msk, n_mbr});
// compute variable one mask at a time to minimise memory imprint
@@ -103,14 +103,14 @@ namespace evalhyd
// shape: (subsets, time)
// \require crps:
// CRPS for each time step.
- // shape: (quantiles, time)
+ // shape: (1, time)
// \assign CRPS:
// CRPS.
- // shape: (subsets, 1, 1)
+ // shape: (subsets,)
void Evaluator::calc_CRPS()
{
// initialise output variable
- // shape: (subsets, thresholds, 1)
+ // shape: (subsets,)
CRPS = xt::zeros<double>({n_msk});
// compute variable one mask at a time to minimise memory imprint
@@ -123,7 +123,7 @@ namespace evalhyd
// calculate the mean over the time steps
// $CRPS = \frac{1}{n} \sum_{k=1}^{n} crps$
xt::view(CRPS, m) =
- xt::nanmean(crps_masked, -1);
+ xt::squeeze(xt::nanmean(crps_masked, -1));
}
}
}
diff --git a/src/probabilist/evaluator_utils.cpp b/src/probabilist/evaluator_utils.cpp
deleted file mode 100644
index 923cc3aef2ce64446d56188104e4e14861a9b74e..0000000000000000000000000000000000000000
--- a/src/probabilist/evaluator_utils.cpp
+++ /dev/null
@@ -1,46 +0,0 @@
-#include <xtensor/xtensor.hpp>
-#include <xtensor/xview.hpp>
-
-namespace evalhyd
-{
- namespace probabilist
- {
- // Determine whether flows are greater than given threshold(s).
- //
- // \param q:
- // Streamflow data.
- // shape: (1+, time)
- // \param thr:
- // Streamflow thresholds.
- // shape: (thresholds,)
- // \return
- // Boolean-like threshold(s) exceedance.
- // shape: (thresholds, 1+, time)
- xt::xtensor<double, 3> is_above_threshold(
- const xt::xtensor<double, 2> &q,
- const xt::xtensor<double, 1> &thr
- )
- {
- return q >= xt::view(thr, xt::all(), xt::newaxis(), xt::newaxis());
- }
-
- // Determine whether flows are strictly lower than given threshold(s)
- //
- // \param q:
- // Streamflow data.
- // shape: (1+, time)
- // \param thr:
- // Streamflow thresholds.
- // shape: (thresholds,)
- // \return
- // Boolean-like threshold(s) exceedance.
- // shape: (thresholds, 1+, time)
- xt::xtensor<double, 3> is_below_threshold(
- const xt::xtensor<double, 2> &q,
- const xt::xtensor<double, 1> &thr
- )
- {
- return q < xt::view(thr, xt::all(), xt::newaxis(), xt::newaxis());
- }
- }
-}
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index d60239f47d28c0cb28cd1b6e85e65adfb87fdaf0..a616f8fedd137ee672156b1dd1b86de0fde7b5e6 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -37,6 +37,5 @@ add_executable(
../src/probabilist/evaluator_brier.cpp
../src/probabilist/evaluator_quantiles.cpp
../src/probabilist/evaluator_elements.cpp
- ../src/probabilist/evaluator_utils.cpp
)
target_link_libraries(evalhyd_tests gtest gtest_main)
diff --git a/tests/test_probabilist.cpp b/tests/test_probabilist.cpp
index 880993fdd5ebbc77efab7b734870fb0f6d80cfa3..ae694f1fb39ac3e1270bda5e4bd560d964416532 100644
--- a/tests/test_probabilist.cpp
+++ b/tests/test_probabilist.cpp
@@ -11,44 +11,50 @@ TEST(ProbabilistTests, TestBrier) {
// read in data
std::ifstream ifs;
ifs.open("./data/q_obs.csv");
- xt::xtensor<double, 2> observed = xt::load_csv<int>(ifs);
+ 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();
+ // add "artificial" site dimension (size 1) to observations
+ auto obs = xt::view(observed, xt::newaxis(), xt::all());
+ // add "artificial" site dimension (size 1) and lead time dimension (size 1)
+ // to predictions
+ auto prd = xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all());
+
// compute scores
xt::xtensor<double, 1> thresholds = {690, 534, 445};
std::vector<xt::xarray<double>> metrics =
evalhyd::evalp(
- observed, predicted, {"BS", "BSS", "BS_CRD", "BS_LBD"},
+ obs, prd, {"BS", "BSS", "BS_CRD", "BS_LBD"},
thresholds
);
// check results
// Brier scores
- xt::xtensor<double, 1> bs =
- {0.10615136, 0.07395622, 0.08669186};
+ xt::xtensor<double, 4> bs =
+ {{{{0.10615136, 0.07395622, 0.08669186}}}};
EXPECT_TRUE(xt::allclose(metrics[0], bs));
// Brier skill scores
- xt::xtensor<double, 1> bss =
- {0.5705594, 0.6661165, 0.5635126};
+ xt::xtensor<double, 4> bss =
+ {{{{0.5705594, 0.6661165, 0.5635126}}}};
EXPECT_TRUE(xt::allclose(metrics[1], bss));
// Brier calibration-refinement decompositions
- xt::xtensor<double, 2> bs_crd =
- {{0.011411758, 0.1524456, 0.2471852},
- {0.005532413, 0.1530793, 0.2215031},
- {0.010139431, 0.1220601, 0.1986125}};
+ xt::xtensor<double, 5> bs_crd =
+ {{{{{0.011411758, 0.1524456, 0.2471852},
+ {0.005532413, 0.1530793, 0.2215031},
+ {0.010139431, 0.1220601, 0.1986125}}}}};
EXPECT_TRUE(xt::allclose(metrics[2], bs_crd));
// Brier likelihood-base rate decompositions
- xt::xtensor<double, 2> bs_lbd =
- {{0.012159881, 0.1506234, 0.2446149},
- {0.008031746, 0.1473869, 0.2133114},
- {0.017191279, 0.1048221, 0.1743227}};
+ xt::xtensor<double, 5> bs_lbd =
+ {{{{{0.012159881, 0.1506234, 0.2446149},
+ {0.008031746, 0.1473869, 0.2133114},
+ {0.017191279, 0.1048221, 0.1743227}}}}};
EXPECT_TRUE(xt::allclose(metrics[3], bs_lbd));
}