diff --git a/include/evalhyd/detail/probabilist/brier.hpp b/include/evalhyd/detail/probabilist/brier.hpp
index b8d80eb927751567ff8c156d117890d804266c07..1da4fdd01dc3982591d6fba9bbf18cb33d478d98 100644
--- a/include/evalhyd/detail/probabilist/brier.hpp
+++ b/include/evalhyd/detail/probabilist/brier.hpp
@@ -37,11 +37,21 @@ namespace evalhyd
template<class XV1D2>
inline xt::xtensor<double, 2> calc_o_k(
const XV1D2& q_obs,
- const XV1D2& q_thr
+ const XV1D2& q_thr,
+ bool high_flow_events
)
{
- // determine observed realisation of threshold(s) exceedance
- return q_obs >= xt::view(q_thr, xt::all(), xt::newaxis());
+ if (high_flow_events)
+ {
+ // observations above threshold(s)
+ return q_obs >= xt::view(q_thr, xt::all(), xt::newaxis());
+ }
+ else
+ {
+ // observations below threshold(s)
+ return q_obs < xt::view(q_thr, xt::all(), xt::newaxis());
+ }
+
}
// Determine mean observed realisation of threshold(s) exceedance.
@@ -112,15 +122,28 @@ namespace evalhyd
template<class XV2D4, class XV1D2>
inline xt::xtensor<double, 2> calc_sum_f_k(
const XV2D4& q_prd,
- const XV1D2& q_thr
+ const XV1D2& q_thr,
+ bool high_flow_events
)
{
- // determine if members have exceeded threshold(s)
- auto f_k = q_prd >=
- xt::view(q_thr, xt::all(), xt::newaxis(), xt::newaxis());
+ if (high_flow_events)
+ {
+ // determine if members are above threshold(s)
+ auto f_k = q_prd >=
+ xt::view(q_thr, xt::all(), xt::newaxis(), xt::newaxis());
+
+ // calculate how many members are above threshold(s)
+ return xt::sum(f_k, 1);
+ }
+ else
+ {
+ // determine if members are below threshold(s)
+ auto f_k = q_prd <
+ xt::view(q_thr, xt::all(), xt::newaxis(), xt::newaxis());
- // calculate how many members have exceeded threshold(s)
- return xt::sum(f_k, 1);
+ // calculate how many members are below threshold(s)
+ return xt::sum(f_k, 1);
+ }
}
// Determine forecast probability of threshold(s) exceedance to occur.
diff --git a/include/evalhyd/detail/probabilist/evaluator.hpp b/include/evalhyd/detail/probabilist/evaluator.hpp
index d4d4c1537be9dd892298d0dff1d0048f92c5df02..90c65d037838a86f9fda98c842dbed90c88155c4 100644
--- a/include/evalhyd/detail/probabilist/evaluator.hpp
+++ b/include/evalhyd/detail/probabilist/evaluator.hpp
@@ -57,6 +57,7 @@ namespace evalhyd
const view1d_xtensor2d_double_type& q_obs;
const view2d_xtensor4d_double_type& q_prd;
const view1d_xtensor2d_double_type& q_thr;
+ const bool high_flow_events;
xt::xtensor<bool, 2> t_msk;
const std::vector<xt::xkeep_slice<int>>& b_exp;
@@ -112,7 +113,7 @@ namespace evalhyd
if (!o_k.has_value())
{
o_k = elements::calc_o_k(
- q_obs, q_thr
+ q_obs, q_thr, high_flow_events
);
}
return o_k.value();
@@ -134,7 +135,7 @@ namespace evalhyd
if (!sum_f_k.has_value())
{
sum_f_k = elements::calc_sum_f_k(
- q_prd, q_thr
+ q_prd, q_thr, high_flow_events
);
}
return sum_f_k.value();
@@ -300,9 +301,12 @@ namespace evalhyd
Evaluator(const view1d_xtensor2d_double_type& obs,
const view2d_xtensor4d_double_type& prd,
const view1d_xtensor2d_double_type& thr,
+ const bool high_flow_events,
const view2d_xtensor4d_bool_type& msk,
const std::vector<xt::xkeep_slice<int>>& exp) :
- q_obs{obs}, q_prd{prd}, q_thr{thr}, t_msk(msk), b_exp(exp)
+ q_obs{obs}, q_prd{prd}, q_thr{thr},
+ high_flow_events{high_flow_events},
+ t_msk(msk), b_exp(exp)
{
// initialise a mask if none provided
// (corresponding to no temporal subset)
diff --git a/include/evalhyd/detail/utils.hpp b/include/evalhyd/detail/utils.hpp
index 8d966610bdff6fd4e935b33decb3e5c173d26a12..1c55572163fa0200ff91771a0746f17fe356ef6e 100644
--- a/include/evalhyd/detail/utils.hpp
+++ b/include/evalhyd/detail/utils.hpp
@@ -9,20 +9,23 @@
#include <unordered_set>
#include <vector>
#include <stdexcept>
+
+#include <xtl/xoptional.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xrandom.hpp>
+
namespace evalhyd
{
namespace utils
{
- /// Procedure to check that all elements in the list of metrics are
- /// valid metrics.
- ///
- /// \param [in] requested_metrics:
- /// Vector of strings for the metric(s) to be computed.
- /// \param [in] valid_metrics:
- /// Vector of strings for the metric(s) to can be computed.
+ // Procedure to check that all elements in the list of metrics are
+ // valid metrics.
+ //
+ // \param requested_metrics
+ // Vector of strings for the metric(s) to be computed.
+ // \param valid_metrics
+ // Vector of strings for the metric(s) to can be computed.
inline void check_metrics (
const std::vector<std::string>& requested_metrics,
const std::vector<std::string>& valid_metrics
@@ -40,11 +43,11 @@ namespace evalhyd
}
}
- /// Procedure to check that all elements for a bootstrap experiment
- /// are provided and valid.
- ///
- /// \param [in] bootstrap:
- /// Map of parameters for the bootstrap experiment.
+ // Procedure to check that all elements for a bootstrap experiment
+ // are provided and valid.
+ //
+ // \param bootstrap
+ // Map of parameters for the bootstrap experiment.
inline void check_bootstrap (
const std::unordered_map<std::string, int>& bootstrap
)
@@ -105,27 +108,34 @@ namespace evalhyd
namespace evalp
{
- /// Procedure to check that optional parameters are provided
- /// as arguments when required metrics need them.
- ///
- /// \param [in] metrics:
- /// Vector of strings for the metric(s) to be computed.
- /// \param [in] thresholds:
- /// Array of thresholds for metrics based on exceedance events.
- inline void check_optionals (
+ // Procedure to check that optional parameters are provided
+ // as arguments when required metrics need them.
+ //
+ // \param metrics
+ // Vector of strings for the metric(s) to be computed.
+ // \param thresholds
+ // Array of thresholds for metrics based on exceedance events.
+ // \param events
+ // Kind of streamflow exceedance events.
+ // \return
+ // Whether high flow events are considered.
+ inline bool check_optionals (
const std::vector<std::string>& metrics,
- const xt::xtensor<double, 2>& thresholds
+ const xt::xtensor<double, 2>& thresholds,
+ xtl::xoptional<const std::string, bool> events
)
{
std::vector<std::string>threshold_metrics =
{"BS", "BS_CRD", "BS_LBD", "BSS",
"POD", "POFD", "FAR", "CSI", "ROCSS"};
+ bool thresholds_required_and_provided = false;
for (const auto& metric : metrics)
{
if (std::find(threshold_metrics.begin(), threshold_metrics.end(),
metric) != threshold_metrics.end())
{
+ // check thresholds
if (thresholds.size() < 1)
{
throw std::runtime_error(
@@ -133,8 +143,45 @@ namespace evalhyd
"compute " + metric
);
}
+ else
+ {
+ thresholds_required_and_provided = true;
+ break;
+ }
}
}
+
+ // check events
+ if (thresholds_required_and_provided)
+ {
+ if (events.has_value())
+ {
+ if (events.value() == "high")
+ {
+ return true;
+ }
+ else if (events.value() == "low")
+ {
+ return false;
+ }
+ else
+ {
+ throw std::runtime_error(
+ "invalid value for streamflow *events*"
+ );
+ }
+ }
+ else
+ {
+ throw std::runtime_error(
+ "*q_thr* provided but *events* is missing"
+ );
+ }
+ }
+ else
+ {
+ return true;
+ }
}
}
}
diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp
index e56e5b14ef009da1d67cb0e158949f51e3a57d83..5582e1f1090404751a054dc4cf92980b74c9c409 100644
--- a/include/evalhyd/evalp.hpp
+++ b/include/evalhyd/evalp.hpp
@@ -60,9 +60,19 @@ namespace evalhyd
/// The sequence of evaluation metrics to be computed.
///
/// q_thr: ``XD2``, optional
- /// Streamflow exceedance threshold(s).
+ /// Streamflow exceedance threshold(s). If provided, *events* must
+ /// also be provided.
/// shape: (sites, thresholds)
///
+ /// events: ``std::string``, optional
+ /// The type of streamflow events to consider for threshold
+ /// exceedance-based metrics. It can either be set as "high" when
+ /// flooding conditions/high flow events are evaluated (i.e. event
+ /// occurring when streamflow goes above threshold) or as "low" when
+ /// drought conditions/low flow events are evaluated (i.e. event
+ /// occurring when streamflow goes below threshold). It must be
+ /// provided if *q_thr* is provided.
+ ///
/// t_msk: ``XB4``, optional
/// Mask(s) used to generate temporal subsets of the whole streamflow
/// time series (where True/False is used for the time steps to
@@ -147,6 +157,8 @@ namespace evalhyd
const xt::xexpression<XD4>& q_prd,
const std::vector<std::string>& metrics,
const xt::xexpression<XD2>& q_thr = XD2({}),
+ xtl::xoptional<const std::string, bool> events =
+ xtl::missing<const std::string>(),
const xt::xexpression<XB4>& t_msk = XB4({}),
const xt::xtensor<std::array<char, 32>, 2>& m_cdt = {},
xtl::xoptional<const std::unordered_map<std::string, int>, bool> bootstrap =
@@ -190,7 +202,9 @@ namespace evalhyd
);
// check that optional parameters are given as arguments
- utils::evalp::check_optionals(metrics, q_thr_);
+ bool high_flow_events =
+ utils::evalp::check_optionals(metrics, q_thr_, events);
+
if (bootstrap.has_value())
{
utils::check_bootstrap(bootstrap.value());
@@ -374,7 +388,9 @@ namespace evalhyd
xt::view(t_msk_, s, l, xt::all(), xt::all()));
probabilist::Evaluator<XD2, XD4, XB4> evaluator(
- q_obs_v, q_prd_v, q_thr_v, t_msk_v, b_exp
+ q_obs_v, q_prd_v, q_thr_v,
+ high_flow_events,
+ t_msk_v, b_exp
);
// retrieve or compute requested metrics
diff --git a/tests/test_probabilist.cpp b/tests/test_probabilist.cpp
index d406c29b2da95f60c7bf6a23332a8bf6d2a83dd4..ee863a7ab165ebb7e07a924649917bbd26a32ecf 100644
--- a/tests/test_probabilist.cpp
+++ b/tests/test_probabilist.cpp
@@ -53,7 +53,8 @@ TEST(ProbabilistTests, TestBrier)
// 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
+ thresholds,
+ "high"
);
// check results
@@ -96,6 +97,78 @@ TEST(ProbabilistTests, TestBrier)
);
}
+TEST(ProbabilistTests, TestContingency)
+{
+ // 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())),
+ {"POD", "POFD", "FAR", "CSI", "ROCSS"},
+ thresholds,
+ "high"
+ );
+
+ std::cout << "POD" << std::endl;
+ std::cout << metrics[0] << std::endl;
+ std::cout << "POFD" << std::endl;
+ std::cout << metrics[1] << std::endl;
+ std::cout << "FAR" << std::endl;
+ std::cout << metrics[2] << std::endl;
+ std::cout << "CSI" << std::endl;
+ std::cout << metrics[3] << std::endl;
+ std::cout << "ROCSS" << std::endl;
+ std::cout << metrics[4] << std::endl;
+
+// // check results
+// // POD
+// xt::xtensor<double, 6> pod =
+// {{{{{{0.10615136, 0.07395622, 0.08669186, NAN}}}}}};
+// EXPECT_TRUE(
+// xt::sum(xt::isclose(metrics[0], pod, 1e-05, 1e-08, true))
+// == xt::xscalar<double>(4)
+// );
+//
+// // POFD
+// xt::xtensor<double, 6> pofd =
+// {{{{{{0.5705594, 0.6661165, 0.5635126, NAN}}}}}};
+// EXPECT_TRUE(
+// xt::sum(xt::isclose(metrics[1], pofd, 1e-05, 1e-08, true))
+// == xt::xscalar<double>(4)
+// );
+//
+// // FAR
+// xt::xtensor<double, 6> far =
+// {{{{{{{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], far, 1e-05, 1e-08, true))
+// == xt::xscalar<double>(12)
+// );
+//
+// // CSI
+// xt::xtensor<double, 6> csi =
+// {{{{{{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], csi, 1e-05, 1e-08, true))
+// == xt::xscalar<double>(12)
+// );
+}
+
TEST(ProbabilistTests, TestQuantiles)
{
// read in data
@@ -161,6 +234,7 @@ TEST(ProbabilistTests, TestMasks)
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
metrics,
thresholds,
+ "high",
// shape: (sites [1], lead times [1], subsets [1], time [t])
masks
);
@@ -173,7 +247,8 @@ TEST(ProbabilistTests, TestMasks)
// shape: (sites [1], lead times [1], members [m], time [t-20])
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::range(20, _))),
metrics,
- thresholds
+ thresholds,
+ "high"
);
// check results are identical
@@ -214,6 +289,7 @@ TEST(ProbabilistTests, TestMaskingConditions)
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
metrics,
thresholds,
+ "high",
masks,
q_conditions
);
@@ -224,7 +300,8 @@ TEST(ProbabilistTests, TestMaskingConditions)
xt::eval(xt::where((observed < 2000) | (observed > 3000), observed, NAN)),
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
metrics,
- thresholds
+ thresholds,
+ "high"
);
// check results are identical
@@ -253,6 +330,7 @@ TEST(ProbabilistTests, TestMaskingConditions)
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
metrics,
thresholds,
+ "high",
masks,
q_conditions_
);
@@ -263,7 +341,8 @@ TEST(ProbabilistTests, TestMaskingConditions)
xt::eval(xt::where(q_prd_mean >= median, observed, NAN)),
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
metrics,
- thresholds
+ thresholds,
+ "high"
);
// check results are identical
@@ -289,6 +368,7 @@ TEST(ProbabilistTests, TestMaskingConditions)
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
metrics,
thresholds,
+ "high",
masks,
t_conditions
);
@@ -299,7 +379,8 @@ TEST(ProbabilistTests, TestMaskingConditions)
xt::eval(xt::view(observed_, xt::newaxis(), xt::range(0, 100))),
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::range(0, 100))),
metrics,
- thresholds
+ thresholds,
+ "high"
);
// check results are identical
@@ -340,7 +421,8 @@ TEST(ProbabilistTests, TestMissingData)
observed_nan,
forecast_nan,
metrics,
- thresholds
+ thresholds,
+ "high"
);
// compute metrics on manually subset series (one leadtime at a time)
@@ -359,7 +441,8 @@ TEST(ProbabilistTests, TestMissingData)
observed_pp1,
forecast_pp1,
metrics,
- thresholds
+ thresholds,
+ "high"
);
xt::xtensor<double, 4> forecast_pp2 {{
@@ -377,7 +460,8 @@ TEST(ProbabilistTests, TestMissingData)
observed_pp2,
forecast_pp2,
metrics,
- thresholds
+ thresholds,
+ "high"
);
// check that numerical results are identical
@@ -434,6 +518,7 @@ TEST(ProbabilistTests, TestBootstrap)
xt::eval(xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
metrics,
thresholds,
+ "high", // events
xt::xtensor<bool, 4>({}), // t_msk
{}, // m_cdt
bootstrap,
@@ -456,7 +541,8 @@ TEST(ProbabilistTests, TestBootstrap)
xt::eval(xt::view(observed_x3, xt::newaxis(), xt::all())),
xt::eval(xt::view(predicted_x3, xt::newaxis(), xt::newaxis(), xt::all(), xt::all())),
metrics,
- thresholds
+ thresholds,
+ "high"
);
// check results are identical