diff --git a/include/evalhyd/detail/utils.hpp b/include/evalhyd/detail/utils.hpp
index 64934cab69a06e81fe19bbb923578244f2d06b51..a5bb9b3fe910951d93778c9577973ba9147fe059 100644
--- a/include/evalhyd/detail/utils.hpp
+++ b/include/evalhyd/detail/utils.hpp
@@ -52,6 +52,13 @@ namespace evalhyd
"number of samples missing for bootstrap"
);
}
+ auto n_samples = bootstrap.find("n_samples")->second;
+ if (n_samples < 1)
+ {
+ throw std::runtime_error(
+ "number of samples must be greater than zero"
+ );
+ }
// check len_sample
if (bootstrap.find("len_sample") == bootstrap.end())
{
@@ -59,6 +66,13 @@ namespace evalhyd
"length of sample missing for bootstrap"
);
}
+ auto len_sample = bootstrap.find("len_sample")->second;
+ if (len_sample < 1)
+ {
+ throw std::runtime_error(
+ "length of sample must be greater than zero"
+ );
+ }
// check summary
if (bootstrap.find("summary") == bootstrap.end())
{
@@ -66,9 +80,9 @@ namespace evalhyd
"summary missing for bootstrap"
);
}
- auto s = bootstrap.find("summary")->second;
+ auto summary = bootstrap.find("summary")->second;
// TODO: change upper bound when mean+stddev and quantiles implemented
- if ((s < 0) || (s > 0))
+ if ((summary < 0) || (summary > 0))
{
throw std::runtime_error(
"invalid value for bootstrap summary"
diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp
index 8b0041b6a9667f8d89532dc1ff4eb38d3006ae8a..88e0868ad057731c770ffe530a799082974ab658 100644
--- a/include/evalhyd/evald.hpp
+++ b/include/evalhyd/evald.hpp
@@ -4,6 +4,7 @@
#include <unordered_map>
#include <vector>
+#include <xtl/xoptional.hpp>
#include <xtensor/xexpression.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xarray.hpp>
@@ -58,18 +59,17 @@ namespace evalhyd
///
/// exponent: ``double``, optional
/// The value of the exponent n to use when the *transform* is the
- /// power function. If not provided (or set to default value 1),
- /// the streamflow observations and predictions remain untransformed.
+ /// power function. If not provided, the streamflow observations
+ /// and predictions remain untransformed.
///
/// epsilon: ``double``, optional
/// The value of the small constant ε to add to both the streamflow
/// observations and predictions prior to the calculation of the
/// *metrics* when the *transform* is the reciprocal function, the
/// natural logarithm, or the power function with a negative exponent
- /// (since none are defined for 0). If not provided (or set to default
- /// value -9), one hundredth of the mean of the streamflow
- /// observations is used as value for epsilon, as recommended by
- /// `Pushpalatha et al. (2012)
+ /// (since none are defined for 0). If not provided, one hundredth of
+ /// the mean of the streamflow observations is used as value for
+ /// epsilon, as recommended by `Pushpalatha et al. (2012)
/// <https://doi.org/10.1016/j.jhydrol.2011.11.055>`_.
///
/// t_msk: ``XB2``, optional
@@ -159,13 +159,16 @@ namespace evalhyd
const xt::xexpression<XD2>& q_obs,
const xt::xexpression<XD2>& q_prd,
const std::vector<std::string>& metrics,
- const std::string& transform = "none",
- double exponent = 1,
- double epsilon = -9,
+ xtl::xoptional<const std::string, bool> transform =
+ xtl::missing<const std::string>(),
+ xtl::xoptional<double, bool> exponent =
+ xtl::missing<double>(),
+ xtl::xoptional<double, bool> epsilon =
+ xtl::missing<double>(),
const xt::xexpression<XB2>& t_msk = XB2({}),
const xt::xtensor<std::array<char, 32>, 1>& m_cdt = {},
- const std::unordered_map<std::string, int>& bootstrap =
- {{"n_samples", -9}, {"len_sample", -9}, {"summary", 0}},
+ xtl::xoptional<const std::unordered_map<std::string, int>> bootstrap =
+ xtl::missing<const std::unordered_map<std::string, int>>(),
const std::vector<std::string>& dts = {}
)
{
@@ -196,7 +199,10 @@ namespace evalhyd
);
// check that optional parameters are valid
- utils::check_bootstrap(bootstrap);
+ if (bootstrap.has_value())
+ {
+ utils::check_bootstrap(bootstrap.value());
+ }
// check that data dimensions are compatible
// > time
@@ -277,68 +283,90 @@ namespace evalhyd
// apply streamflow transformation if requested
auto q_transform = [&](const XD2& q)
{
- if ( transform == "none" || (transform == "pow" && exponent == 1))
- {
- return q;
- }
- else if ( transform == "sqrt" )
- {
- return XD2(xt::sqrt(q));
- }
- else if ( transform == "inv" )
+ if (transform.has_value())
{
- if ( epsilon == -9 )
+ if ( transform.value() == "sqrt" )
{
- // determine an epsilon value to avoid zero divide
- epsilon = xt::mean(q_obs_)() * 0.01;
+ return XD2(xt::sqrt(q));
}
-
- return XD2(1. / (q + epsilon));
- }
- else if ( transform == "log" )
- {
- if ( epsilon == -9 )
+ else if ( transform.value() == "inv" )
{
- // determine an epsilon value to avoid log zero
- epsilon = xt::mean(q_obs_)() * 0.01;
- }
+ if ( !epsilon.has_value() )
+ {
+ // determine an epsilon value to avoid zero divide
+ epsilon = xt::mean(q_obs_)() * 0.01;
+ }
- return XD2(xt::log(q + epsilon));
- }
- else if ( transform == "pow" )
- {
- if ( exponent < 0 )
+ return XD2(1. / (q + epsilon.value()));
+ }
+ else if ( transform.value() == "log" )
{
- if ( epsilon == -9 )
+ if ( !epsilon.has_value() )
{
- // determine an epsilon value to avoid zero divide
+ // determine an epsilon value to avoid log zero
epsilon = xt::mean(q_obs_)() * 0.01;
}
- return XD2(xt::pow(q + epsilon, exponent));
+ return XD2(xt::log(q + epsilon.value()));
+ }
+ else if ( transform.value() == "pow" )
+ {
+ if ( exponent.has_value() )
+ {
+ if ( exponent.value() == 1)
+ {
+ return q;
+ }
+ else if ( exponent.value() < 0 )
+ {
+ if ( !epsilon.has_value() )
+ {
+ // determine an epsilon value to avoid zero divide
+ epsilon = xt::mean(q_obs_)() * 0.01;
+ }
+
+ return XD2(xt::pow(q + epsilon.value(),
+ exponent.value()));
+ }
+ else
+ {
+ return XD2(xt::pow(q, exponent.value()));
+ }
+ }
+ else
+ {
+ throw std::runtime_error(
+ "missing exponent for power transformation"
+ );
+ }
}
else
{
- return XD2(xt::pow(q, exponent));
+ throw std::runtime_error(
+ "invalid streamflow transformation: "
+ + transform.value()
+ );
}
}
else
{
- throw std::runtime_error(
- "invalid streamflow transformation: " + transform
- );
+ return q;
}
};
- const XD2 obs = q_transform(q_obs_);
- const XD2 prd = q_transform(q_prd_);
+ const XD2& obs = q_transform(q_obs_);
+ const XD2& prd = q_transform(q_prd_);
// generate bootstrap experiment if requested
std::vector<xt::xkeep_slice<int>> exp;
- auto n_samples = bootstrap.find("n_samples")->second;
- auto len_sample = bootstrap.find("len_sample")->second;
- if ((n_samples != -9) && (len_sample != -9))
+ int summary;
+
+ if (bootstrap.has_value())
{
+ auto n_samples = bootstrap.value().find("n_samples")->second;
+ auto len_sample = bootstrap.value().find("len_sample")->second;
+ summary = bootstrap.value().find("summary")->second;
+
if (dts.empty())
{
throw std::runtime_error(
@@ -354,6 +382,7 @@ namespace evalhyd
{
// if no bootstrap requested, generate one sample
// containing all the time indices once
+ summary = 0;
xt::xtensor<int, 1> all = xt::arange(n_tim);
exp.push_back(xt::keep(all));
}
@@ -364,8 +393,6 @@ namespace evalhyd
// retrieve or compute requested metrics
std::vector<xt::xarray<double>> r;
- auto summary = bootstrap.find("summary")->second;
-
for ( const auto& metric : metrics )
{
if ( metric == "RMSE" )
diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp
index a425a6692e23cfc9f1790511a7d51cd9a2730aee..02505cb8f2bb2d72dc112dc8c94e2961b3700b1a 100644
--- a/include/evalhyd/evalp.hpp
+++ b/include/evalhyd/evalp.hpp
@@ -4,6 +4,7 @@
#include <unordered_map>
#include <vector>
+#include <xtl/xoptional.hpp>
#include <xtensor/xexpression.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xarray.hpp>
@@ -144,8 +145,8 @@ namespace evalhyd
const xt::xexpression<XD2>& q_thr = XD2({}),
const xt::xexpression<XB4>& t_msk = XB4({}),
const xt::xtensor<std::array<char, 32>, 2>& m_cdt = {},
- const std::unordered_map<std::string, int>& bootstrap =
- {{"n_samples", -9}, {"len_sample", -9}, {"summary", 0}},
+ xtl::xoptional<const std::unordered_map<std::string, int>> bootstrap =
+ xtl::missing<const std::unordered_map<std::string, int>>(),
const std::vector<std::string>& dts = {}
)
{
@@ -184,7 +185,10 @@ namespace evalhyd
// check that optional parameters are given as arguments
utils::evalp::check_optionals(metrics, q_thr_);
- utils::check_bootstrap(bootstrap);
+ if (bootstrap.has_value())
+ {
+ utils::check_bootstrap(bootstrap.value());
+ }
// check that data dimensions are compatible
// > time
@@ -266,8 +270,8 @@ namespace evalhyd
std::size_t n_thr = q_thr_.shape(1);
std::size_t n_msk = t_msk_.size() > 0 ? t_msk_.shape(2) :
(m_cdt.size() > 0 ? m_cdt.shape(1) : 1);
- std::size_t n_exp = bootstrap.find("n_samples")->second == -9 ? 1:
- bootstrap.find("n_samples")->second;
+ std::size_t n_exp = !bootstrap.has_value() ? 1:
+ bootstrap.value().find("n_samples")->second;
// register metrics number of dimensions
std::unordered_map<std::string, std::vector<std::size_t>> dim;
@@ -308,10 +312,14 @@ namespace evalhyd
// generate bootstrap experiment if requested
std::vector<xt::xkeep_slice<int>> b_exp;
- auto n_samples = bootstrap.find("n_samples")->second;
- auto len_sample = bootstrap.find("len_sample")->second;
- if ((n_samples != -9) && (len_sample != -9))
+ int summary;
+
+ if (bootstrap.has_value())
{
+ auto n_samples = bootstrap.value().find("n_samples")->second;
+ auto len_sample = bootstrap.value().find("len_sample")->second;
+ summary = bootstrap.value().find("summary")->second;
+
if (dts.empty())
{
throw std::runtime_error(
@@ -325,6 +333,7 @@ namespace evalhyd
{
// if no bootstrap requested, generate one sample
// containing all the time indices once
+ summary = 0;
xt::xtensor<int, 1> all = xt::arange(n_tim);
b_exp.push_back(xt::keep(all));
}
@@ -336,8 +345,6 @@ namespace evalhyd
r.emplace_back(xt::zeros<double>(dim[metric]));
}
- auto summary = bootstrap.find("summary")->second;
-
// compute variables one site at a time to minimise memory imprint
for (std::size_t s = 0; s < n_sit; s++)
{