diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp
index 3e37e432927a005aec0288bcf71ea15d26fe5151..56926350ad3f9c87c7f8a69eaf0cdbcc20c4e481 100644
--- a/include/evalhyd/evald.hpp
+++ b/include/evalhyd/evald.hpp
@@ -28,17 +28,15 @@ namespace evalhyd
/// Streamflow observations. Time steps with missing observations
/// must be assigned `NAN` values. Those time steps will be ignored
/// both in the observations and the predictions before the
- /// *metrics* are computed. Observations and predictions must be
- /// broadcastable.
- /// shape: (1 or series, time)
+ /// *metrics* are computed.
+ /// shape: (1, time)
///
/// q_prd: ``xt::xtensor<double, 2>``
/// Streamflow predictions. Time steps with missing predictions
/// must be assigned `NAN` values. Those time steps will be ignored
/// both in the observations and the predictions before the
- /// *metrics* are computed. Observations and predictions must feature
- /// be broadcastable.
- /// shape: (1 or series, time)
+ /// *metrics* are computed.
+ /// shape: (series, time)
///
/// metrics: ``std::vector<std::string>``
/// The sequence of evaluation metrics to be computed.
@@ -85,11 +83,11 @@ namespace evalhyd
/// the subset). If provided, masks must feature the same number of
/// dimensions as observations and predictions, and it must
/// broadcastable with both of them.
- /// shape: (1 or series, time)
+ /// shape: (subsets, time)
///
/// .. seealso:: :doc:`../../functionalities/temporal-masking`
///
- /// m_cdt: ``xt::xtensor<std::array<char, 32>, 2>``, optional
+ /// m_cdt: ``xt::xtensor<std::array<char, 32>, 1>``, optional
/// Masking conditions to use to generate temporal subsets. Each
/// condition consists in a string and can be specified on observed
/// streamflow values/statistics (mean, median, quantile), or on time
@@ -97,7 +95,7 @@ namespace evalhyd
/// precedence. If not provided and neither is *t_msk*, no subset is
/// performed. If provided, there must be as many conditions as there
/// are time series of observations.
- /// shape: (1 or series, 1)
+ /// shape: (subsets,)
///
/// .. seealso:: :doc:`../../functionalities/conditional-masking`
///
@@ -125,7 +123,7 @@ namespace evalhyd
/// ``std::vector<xt::xarray<double>>``
/// The sequence of evaluation metrics computed
/// in the same order as given in *metrics*.
- /// shape: (metrics,)<(1 or series, samples)>
+ /// shape: (metrics,)<(series, subsets, samples)>
///
/// :Examples:
///
@@ -168,7 +166,7 @@ namespace evalhyd
const double exponent = 1,
double epsilon = -9,
const xt::xtensor<bool, 2>& t_msk = {},
- const xt::xtensor<std::array<char, 32>, 2>& m_cdt = {},
+ const xt::xtensor<std::array<char, 32>, 1>& m_cdt = {},
const std::unordered_map<std::string, int>& bootstrap =
{{"n_samples", -9}, {"len_sample", -9}},
const std::vector<std::string>& dts = {}
@@ -183,6 +181,30 @@ namespace evalhyd
// check that optional parameters are valid
eh::utils::check_bootstrap(bootstrap);
+ // check that data dimensions are compatible
+ // > time
+ if (q_obs.shape(1) != q_prd.shape(1))
+ throw std::runtime_error(
+ "observations and predictions feature different "
+ "temporal lengths"
+ );
+ if (t_msk.size() > 0)
+ if (q_obs.shape(1) != t_msk.shape(1))
+ throw std::runtime_error(
+ "observations and masks feature different "
+ "temporal lengths"
+ );
+ // > series
+ if (q_obs.shape(0) != 1)
+ throw std::runtime_error(
+ "observations contain more than one time series"
+ );
+
+ // retrieve dimensions
+ std::size_t n_tim = q_obs.shape(1);
+ std::size_t n_msk = t_msk.size() > 0 ? t_msk.shape(0) :
+ (m_cdt.size() > 0 ? m_cdt.shape(0) : 1);
+
// initialise a mask if none provided
// (corresponding to no temporal subset)
auto gen_msk = [&]() {
@@ -192,95 +214,22 @@ namespace evalhyd
// else if m_cdt provided, use them to generate t_msk
else if (m_cdt.size() > 0)
{
- xt::xtensor<bool, 2> c_msk = xt::zeros<bool>(q_obs.shape());
-
- // flatten arrays to bypass n-dim considerations
- // (possible because shapes are constrained to be the same)
- if (m_cdt.shape(m_cdt.dimension() - 1) != 1)
- throw std::runtime_error("length of last axis in masking "
- "conditions must be equal to one");
- for (int a = 0; a < m_cdt.dimension() - 1; a++)
- if (q_obs.shape(a) != m_cdt.shape(a))
- throw std::runtime_error("masking conditions and observations "
- "feature incompatible shapes");
-
- auto f_cdt = xt::flatten(m_cdt);
- auto f_msk = xt::flatten(c_msk);
- auto f_obs = xt::flatten(q_obs);
-
- // determine length of temporal axis
- auto nt = q_obs.shape(q_obs.dimension() - 1);
+ xt::xtensor<bool, 2> c_msk = xt::zeros<bool>({n_msk, n_tim});
- // generate mask gradually, one condition at a time
- for (int i = 0; i < m_cdt.size(); i++)
- {
- xt::view(f_msk, xt::range(i*nt, (i+1)*nt)) =
- evalhyd::masks::generate_mask_from_conditions(
- xt::view(f_cdt, i),
- xt::view(f_obs, xt::range(i*nt, (i+1)*nt))
+ for (int m = 0; m < n_msk; m++)
+ xt::view(c_msk, m) =
+ eh::masks::generate_mask_from_conditions(
+ m_cdt[0], xt::view(q_obs, 0), q_prd
);
- }
return c_msk;
}
// if neither t_msk nor m_cdt provided, generate dummy mask
else
- return xt::xtensor<bool, 2>{xt::ones<bool>(q_obs.shape())};
+ return xt::xtensor<bool, 2>{xt::ones<bool>({std::size_t{1}, n_tim})};
};
- const xt::xtensor<bool, 2> msk = gen_msk();
-
- // check that observations, predictions, and masks dimensions are compatible
- if (q_obs.dimension() != q_prd.dimension())
- {
- throw std::runtime_error(
- "observations and predictions feature "
- "different numbers of dimensions"
- );
- }
- if (q_obs.dimension() != msk.dimension())
- {
- throw std::runtime_error(
- "observations and masks feature "
- "different numbers of dimensions"
- );
- }
- if (!dts.empty())
- if (q_obs.shape(q_obs.dimension() - 1) != dts.size())
- throw std::runtime_error(
- "observations and datetimes feature different "
- "temporal lengths"
- );
-
- auto obs_shp = q_obs.shape();
- auto prd_shp = q_prd.shape();
- auto msk_shp = msk.shape();
-
- // check observations, predictions, and masks are broadcastable
- if (obs_shp != prd_shp)
- for (int a = 0; a < q_obs.dimension(); a++)
- if (obs_shp[a] != prd_shp[a])
- if ((obs_shp[a] != 1) and (prd_shp[a] != 1))
- throw std::runtime_error(
- "observations and predictions are not "
- "broadcastable"
- );
-
- if (obs_shp != msk_shp)
- for (int a = 0; a < q_obs.dimension(); a++)
- if (obs_shp[a] != msk_shp[a])
- if ((obs_shp[a] != 1) and (msk_shp[a] != 1))
- throw std::runtime_error(
- "observations and masks are not broadcastable"
- );
-
- if (prd_shp != msk_shp)
- for (int a = 0; a < q_prd.dimension(); a++)
- if (prd_shp[a] != msk_shp[a])
- if ((prd_shp[a] != 1) and (msk_shp[a] != 1))
- throw std::runtime_error(
- "predictions and masks are not broadcastable"
- );
+ auto msk = gen_msk();
// apply streamflow transformation if requested
auto q_transform = [&](const xt::xtensor<double, 2>& q)
@@ -349,7 +298,7 @@ namespace evalhyd
{
// if no bootstrap requested, generate one sample
// containing all the time indices once
- xt::xtensor<int, 1> all = xt::arange(obs.shape(obs.dimension() - 1));
+ xt::xtensor<int, 1> all = xt::arange(n_tim);
exp.push_back(xt::keep(all));
}
diff --git a/src/determinist/evaluator.hpp b/src/determinist/evaluator.hpp
index 0f09078732323942465291f74f4e2ab69b7cdf62..bec074d92f4ddc70558382dff26241ff33006d35 100644
--- a/src/determinist/evaluator.hpp
+++ b/src/determinist/evaluator.hpp
@@ -27,36 +27,17 @@ namespace evalhyd
{
class Evaluator
{
- // define type for function to be used to "mask" (i.e. yield NaNs) the
- // time steps for which the observations/predictions are missing data
- typedef xt::xfunction<
- xt::detail::conditional_ternary,
- xt::xfunction<
- xt::detail::bitwise_or,
- xt::xfunction<
- xt::math::isnan_fun,
- const xt::xtensor<double, 2>&
- >,
- xt::xfunction<
- xt::detail::logical_not,
- const xt::xtensor<bool, 2>&
- >
- >,
- xt::xscalar<float>,
- const xt::xtensor<double, 2>&
- > nan_function;
-
private:
// members for input data
- const xt::xtensor<double, 2>& _q_obs;
- const xt::xtensor<double, 2>& _q_prd;
- nan_function q_obs;
- nan_function q_prd;
+ const xt::xtensor<double, 2>& q_obs;
+ const xt::xtensor<double, 2>& q_prd;
+ xt::xtensor<bool, 3> t_msk;
const std::vector<xt::xkeep_slice<int>>& b_exp;
// members for dimensions
- std::size_t n_srs;
std::size_t n_tim;
+ std::size_t n_msk;
+ std::size_t n_srs;
std::size_t n_exp;
std::vector<std::size_t> inner_dims;
std::vector<std::size_t> inter_dims;
@@ -64,14 +45,14 @@ namespace evalhyd
std::vector<std::size_t> final_dims;
// members for computational elements
- xt::xtensor<double, 3> mean_obs;
- xt::xtensor<double, 3> mean_prd;
+ xt::xtensor<double, 4> mean_obs;
+ xt::xtensor<double, 4> mean_prd;
xt::xtensor<double, 2> quad_err;
- xt::xtensor<double, 3> quad_obs;
- xt::xtensor<double, 3> quad_prd;
- xt::xtensor<double, 2> r_pearson;
- xt::xtensor<double, 2> alpha;
- xt::xtensor<double, 2> bias;
+ xt::xtensor<double, 4> quad_obs;
+ xt::xtensor<double, 4> quad_prd;
+ xt::xtensor<double, 3> r_pearson;
+ xt::xtensor<double, 3> alpha;
+ xt::xtensor<double, 3> bias;
public:
// constructor method
@@ -79,45 +60,44 @@ namespace evalhyd
const xt::xtensor<double, 2>& prd,
const xt::xtensor<bool, 2>& msk,
const std::vector<xt::xkeep_slice<int>>& exp) :
- _q_obs{obs},
- _q_prd{prd},
- // "mask" (i.e. use NaN) observations where predictions are
- // missing (i.e. NaN)
- q_obs{xt::where(xt::isnan(prd) | !msk, NAN, obs)},
- // "mask" (i.e. use NaN) predictions where observations are
- // missing (i.e. NaN)
- q_prd{xt::where(xt::isnan(obs) | !msk, NAN, prd)},
- b_exp{exp}
+ q_obs{obs}, q_prd{prd}, b_exp{exp}
{
- if (q_obs.shape(0) > q_prd.shape(0))
- n_srs = q_obs.shape(0);
- else
- n_srs = q_prd.shape(0);
-
- if (q_obs.shape(1) > q_prd.shape(1))
- n_tim = q_obs.shape(1);
- else
- n_tim = q_prd.shape(1);
-
+ // determine size for recurring dimensions
+ n_tim = q_prd.shape(1);
+ n_srs = q_prd.shape(0);
+ n_msk = msk.shape(0);
n_exp = b_exp.size();
// determine dimensions for inner components, intermediate
// elements, for time average elements, and for final metrics:
- // -> inner components shape: (samples, series)
- inner_dims = {n_exp, n_srs};
- // -> intermediate elements shape: (samples, series, time)
- inter_dims = {n_exp, n_srs, n_tim};
- // -> time average elements shape: (samples, series, 1)
- mean_dims = {n_exp, n_srs, 1};
- // -> final metrics shape: (series, samples)
- final_dims = {n_srs, n_exp};
+ // -> inner components shape: (samples, subsets, series)
+ inner_dims = {n_msk, n_exp, n_srs};
+ // -> intermediate elements shape: (samples, subsets, series, time)
+ inter_dims = {n_msk, n_exp, n_srs, n_tim};
+ // -> time average elements shape: (samples, subsets, series, 1)
+ mean_dims = {n_msk, n_exp, n_srs, 1};
+ // -> final metrics shape: (series, subsets, samples)
+ final_dims = {n_srs, n_msk, n_exp};
+
+ // drop time steps where observations or predictions are NaN
+ auto obs_nan = xt::isnan(q_obs);
+ auto prd_nan = xt::isnan(q_prd);
+
+ t_msk = xt::ones<bool>({n_msk, n_srs, n_tim});
+ xt::view(t_msk, xt::all()) =
+ xt::view(msk, xt::all(), xt::newaxis(), xt::all());
+ for (int m = 0; m < n_msk; m++) {
+ xt::view(t_msk, m) =
+ xt::where(obs_nan | prd_nan,
+ false, xt::view(t_msk, m));
+ }
};
// members for evaluation metrics
- xt::xtensor<double, 2> RMSE;
- xt::xtensor<double, 2> NSE;
- xt::xtensor<double, 2> KGE;
- xt::xtensor<double, 2> KGEPRIME;
+ xt::xtensor<double, 3> RMSE;
+ xt::xtensor<double, 3> NSE;
+ xt::xtensor<double, 3> KGE;
+ xt::xtensor<double, 3> KGEPRIME;
// methods to compute elements
void calc_mean_obs();
@@ -143,13 +123,22 @@ namespace evalhyd
// shape: (series, time)
// \assign mean_obs:
// Mean observed streamflow.
- // shape: (samples, series, 1)
+ // shape: (subsets, samples, series, 1)
void Evaluator::calc_mean_obs()
{
mean_obs = xt::zeros<double>(mean_dims);
- for (int e = 0; e < n_exp; e++) {
- auto obs = xt::view(q_obs, xt::all(), b_exp[e]);
- xt::view(mean_obs, e) = xt::nanmean(obs, -1, xt::keep_dims);
+ for (int m = 0; m < n_msk; m++) {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+
+ // compute variable one sample at a time
+ for (int e = 0; e < n_exp; e++) {
+ // apply the bootstrap sampling
+ auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
+ xt::view(mean_obs, m, e) =
+ xt::nanmean(obs, -1, xt::keep_dims);
+ }
}
}
@@ -160,13 +149,22 @@ namespace evalhyd
// shape: (series, time)
// \assign mean_prd:
// Mean predicted streamflow.
- // shape: (samples, series, 1)
+ // shape: (subsets, samples, series, 1)
void Evaluator::calc_mean_prd()
{
mean_prd = xt::zeros<double>(mean_dims);
- for (int e = 0; e < n_exp; e++) {
- auto prd = xt::view(q_prd, xt::all(), b_exp[e]);
- xt::view(mean_prd, e) = xt::nanmean(prd, -1, xt::keep_dims);
+ for (int m = 0; m < n_msk; m++) {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
+
+ // compute variable one sample at a time
+ for (int e = 0; e < n_exp; e++) {
+ // apply the bootstrap sampling
+ auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
+ xt::view(mean_prd, m, e) =
+ xt::nanmean(prd, -1, xt::keep_dims);
+ }
}
}
@@ -197,14 +195,20 @@ namespace evalhyd
// shape: (samples, series, 1)
// \assign quad_obs:
// Quadratic errors between observations and mean observation.
- // shape: (samples, series, time)
+ // shape: (subsets, samples, series, time)
void Evaluator::calc_quad_obs()
{
quad_obs = xt::zeros<double>(inter_dims);
- for (int e = 0; e < n_exp; e++) {
- xt::view(quad_obs, e) = xt::square(
- q_obs - xt::view(mean_obs, e)
- );
+ for (int m = 0; m < n_msk; m++) {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+
+ for (int e = 0; e < n_exp; e++) {
+ xt::view(quad_obs, m, e) = xt::square(
+ obs_masked - xt::view(mean_obs, m, e)
+ );
+ }
}
}
@@ -218,14 +222,20 @@ namespace evalhyd
// shape: (samples, series, 1)
// \assign quad_prd:
// Quadratic errors between predictions and mean prediction.
- // shape: (samples, series, time)
+ // shape: (subsets, samples, series, time)
void Evaluator::calc_quad_prd()
{
quad_prd = xt::zeros<double>(inter_dims);
- for (int e = 0; e < n_exp; e++) {
- xt::view(quad_prd, e) = xt::square(
- q_prd - xt::view(mean_prd, e)
- );
+ for (int m = 0; m < n_msk; m++) {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
+
+ for (int e = 0; e < n_exp; e++) {
+ xt::view(quad_prd, m, e) = xt::square(
+ prd_masked - xt::view(mean_prd, m, e)
+ );
+ }
}
}
@@ -251,28 +261,36 @@ namespace evalhyd
// shape: (samples, series, time)
// \assign r_pearson:
// Pearson correlation coefficients.
- // shape: (samples, series)
+ // shape: (subsets, samples, series)
void Evaluator::calc_r_pearson()
{
// calculate error in timing and dynamics $r_{pearson}$
// (Pearson's correlation coefficient)
r_pearson = xt::zeros<double>(inner_dims);
- for (int e = 0; e < n_exp; e++) {
- auto prd = xt::view(q_prd, xt::all(), b_exp[e]);
- auto obs = xt::view(q_obs, xt::all(), b_exp[e]);
- auto r_num = xt::nansum(
- (prd - xt::view(mean_prd, e))
- * (obs - xt::view(mean_obs, e)),
- -1
- );
-
- auto prd2 = xt::view(quad_prd, e, xt::all(), b_exp[e]);
- auto obs2 = xt::view(quad_obs, e, xt::all(), b_exp[e]);
- auto r_den = xt::sqrt(
- xt::nansum(prd2, -1) * xt::nansum(obs2, -1)
- );
-
- xt::view(r_pearson, e) = r_num / r_den;
+ for (int m = 0; m < n_msk; m++) {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
+ auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+
+ // compute variable one sample at a time
+ for (int e = 0; e < n_exp; e++) {
+ auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
+ auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
+ auto r_num = xt::nansum(
+ (prd - xt::view(mean_prd, m, e))
+ * (obs - xt::view(mean_obs, m, e)),
+ -1
+ );
+
+ auto prd2 = xt::view(quad_prd, m, e, xt::all(), b_exp[e]);
+ auto obs2 = xt::view(quad_obs, m, e, xt::all(), b_exp[e]);
+ auto r_den = xt::sqrt(
+ xt::nansum(prd2, -1) * xt::nansum(obs2, -1)
+ );
+
+ xt::view(r_pearson, m, e) = r_num / r_den;
+ }
}
}
@@ -292,17 +310,25 @@ namespace evalhyd
// shape: (samples, series, 1)
// \assign alpha:
// Alphas, ratios of standard deviations.
- // shape: (samples, series)
+ // shape: (subsets, samples, series)
void Evaluator::calc_alpha()
{
// calculate error in spread of flow $alpha$
alpha = xt::zeros<double>(inner_dims);
- for (int e = 0; e < n_exp; e++) {
- auto prd = xt::view(q_prd, xt::all(), b_exp[e]);
- auto obs = xt::view(q_obs, xt::all(), b_exp[e]);
- xt::view(alpha, e) =
- detail::std_dev(prd, xt::view(mean_prd, e))
- / detail::std_dev(obs, xt::view(mean_obs, e));
+ for (int m = 0; m < n_msk; m++) {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
+ auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+
+ // compute variable one sample at a time
+ for (int e = 0; e < n_exp; e++) {
+ auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
+ auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
+ xt::view(alpha, m, e) =
+ detail::std_dev(prd, xt::view(mean_prd, m, e))
+ / detail::std_dev(obs, xt::view(mean_obs, m, e));
+ }
}
}
@@ -316,16 +342,24 @@ namespace evalhyd
// shape: (series, time)
// \assign bias:
// Biases.
- // shape: (samples, series)
+ // shape: (subsets, samples, series)
void Evaluator::calc_bias()
{
// calculate $bias$
bias = xt::zeros<double>(inner_dims);
- for (int e = 0; e < n_exp; e++) {
- auto prd = xt::view(q_prd, xt::all(), b_exp[e]);
- auto obs = xt::view(q_obs, xt::all(), b_exp[e]);
- xt::view(bias, e) =
- xt::nansum(prd, -1) / xt::nansum(obs, -1);
+ for (int m = 0; m < n_msk; m++) {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
+ auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+
+ // compute variable one sample at a time
+ for (int e = 0; e < n_exp; e++) {
+ auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
+ auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
+ xt::view(bias, m, e) =
+ xt::nansum(prd, -1) / xt::nansum(obs, -1);
+ }
}
}
@@ -336,14 +370,21 @@ namespace evalhyd
// shape: (series, time)
// \assign RMSE:
// Root-mean-square errors.
- // shape: (series, samples)
+ // shape: (series, subsets, samples)
void Evaluator::calc_RMSE()
{
// compute RMSE
RMSE = xt::zeros<double>(final_dims);
- for (int e = 0; e < n_exp; e++) {
- auto err2 = xt::view(quad_err, xt::all(), b_exp[e]);
- xt::view(RMSE, xt::all(), e) = xt::sqrt(xt::nanmean(err2, -1));
+ for (int m = 0; m < n_msk; m++) {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto quad_err_masked = xt::where(xt::view(t_msk, m), quad_err, NAN);
+
+ // compute variable one sample at a time
+ for (int e = 0; e < n_exp; e++) {
+ auto err2 = xt::view(quad_err_masked, xt::all(), b_exp[e]);
+ xt::view(RMSE, xt::all(), m, e) = xt::sqrt(xt::nanmean(err2, -1));
+ }
}
}
@@ -357,21 +398,28 @@ namespace evalhyd
// shape: (samples, series, time)
// \assign NSE:
// Nash-Sutcliffe efficiencies.
- // shape: (series, samples)
+ // shape: (series, subsets, samples)
void Evaluator::calc_NSE()
{
NSE = xt::zeros<double>(final_dims);
- for (int e = 0; e < n_exp; e++) {
- // compute squared errors operands
- auto err2 = xt::view(quad_err, xt::all(), b_exp[e]);
- xt::xtensor<double, 1> f_num =
- xt::nansum(err2, -1);
- auto obs2 = xt::view(quad_obs, e, xt::all(), b_exp[e]);
- xt::xtensor<double, 1> f_den =
- xt::nansum(obs2, -1);
-
- // compute NSE
- xt::view(NSE, xt::all(), e) = 1 - (f_num / f_den);
+ for (int m = 0; m < n_msk; m++) {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto quad_err_masked = xt::where(xt::view(t_msk, m), quad_err, NAN);
+
+ // compute variable one sample at a time
+ for (int e = 0; e < n_exp; e++) {
+ // compute squared errors operands
+ auto err2 = xt::view(quad_err_masked, xt::all(), b_exp[e]);
+ xt::xtensor<double, 1> f_num =
+ xt::nansum(err2, -1);
+ auto obs2 = xt::view(quad_obs, m, e, xt::all(), b_exp[e]);
+ xt::xtensor<double, 1> f_den =
+ xt::nansum(obs2, -1);
+
+ // compute NSE
+ xt::view(NSE, xt::all(), m, e) = 1 - (f_num / f_den);
+ }
}
}
@@ -388,17 +436,19 @@ namespace evalhyd
// shape: (samples, series)
// \assign KGE:
// Kling-Gupta efficiencies.
- // shape: (series, samples)
+ // shape: (series, subsets, samples)
void Evaluator::calc_KGE()
{
KGE = xt::zeros<double>(final_dims);
- for (int e = 0; e < n_exp; e++) {
- // compute KGE
- xt::view(KGE, xt::all(), e) = 1 - xt::sqrt(
- xt::square(xt::view(r_pearson, e) - 1)
- + xt::square(xt::view(alpha, e) - 1)
- + xt::square(xt::view(bias, e) - 1)
- );
+ for (int m = 0; m < n_msk; m++) {
+ for (int e = 0; e < n_exp; e++) {
+ // compute KGE
+ xt::view(KGE, xt::all(), m, e) = 1 - xt::sqrt(
+ xt::square(xt::view(r_pearson, m, e) - 1)
+ + xt::square(xt::view(alpha, m, e) - 1)
+ + xt::square(xt::view(bias, m, e) - 1)
+ );
+ }
}
}
@@ -421,22 +471,24 @@ namespace evalhyd
// shape: (samples, series)
// \assign KGEPRIME:
// Modified Kling-Gupta efficiencies.
- // shape: (series, samples)
+ // shape: (series, subsets, samples)
void Evaluator::calc_KGEPRIME()
{
KGEPRIME = xt::zeros<double>(final_dims);
- for (int e = 0; e < n_exp; e++) {
- // calculate error in spread of flow $gamma$
- auto gamma = xt::view(alpha, e)
- * (xt::view(mean_obs, e, xt::all(), 0)
- / xt::view(mean_prd, e, xt::all(), 0));
-
- // compute KGEPRIME
- xt::view(KGEPRIME, xt::all(), e) = 1 - xt::sqrt(
- xt::square(xt::view(r_pearson, e) - 1)
- + xt::square(gamma - 1)
- + xt::square(xt::view(bias, e) - 1)
- );
+ for (int m = 0; m < n_msk; m++) {
+ for (int e = 0; e < n_exp; e++) {
+ // calculate error in spread of flow $gamma$
+ auto gamma = xt::view(alpha, m, e)
+ * (xt::view(mean_obs, m, e, xt::all(), 0)
+ / xt::view(mean_prd, m, e, xt::all(), 0));
+
+ // compute KGEPRIME
+ xt::view(KGEPRIME, xt::all(), m, e) = 1 - xt::sqrt(
+ xt::square(xt::view(r_pearson, m, e) - 1)
+ + xt::square(gamma - 1)
+ + xt::square(xt::view(bias, m, e) - 1)
+ );
+ }
}
}
}
diff --git a/tests/test_determinist.cpp b/tests/test_determinist.cpp
index 76dbfdbca912c1e8378088415448bcfc2bacb23b..eb5ef2dfcdc99f0f341619720df1bcdd669225f7 100644
--- a/tests/test_determinist.cpp
+++ b/tests/test_determinist.cpp
@@ -43,36 +43,36 @@ TEST(DeterministTests, TestMetrics)
);
// check results on all metrics
- xt::xtensor<double, 2> rmse =
- {{777.034272},
- {776.878479},
- {777.800217},
- {778.151082},
- {778.61487}};
+ xt::xtensor<double, 3> rmse =
+ {{{777.034272}},
+ {{776.878479}},
+ {{777.800217}},
+ {{778.151082}},
+ {{778.61487 }}};
EXPECT_TRUE(xt::allclose(metrics[0], rmse));
- xt::xtensor<double, 2> nse =
- {{0.718912},
- {0.719025},
- {0.718358},
- {0.718104},
- {0.717767}};
+ xt::xtensor<double, 3> nse =
+ {{{0.718912}},
+ {{0.719025}},
+ {{0.718358}},
+ {{0.718104}},
+ {{0.717767}}};
EXPECT_TRUE(xt::allclose(metrics[1], nse));
- xt::xtensor<double, 2> kge =
- {{0.748088},
- {0.746106},
- {0.744111},
- {0.743011},
- {0.741768}};
+ xt::xtensor<double, 3> kge =
+ {{{0.748088}},
+ {{0.746106}},
+ {{0.744111}},
+ {{0.743011}},
+ {{0.741768}}};
EXPECT_TRUE(xt::allclose(metrics[2], kge));
- xt::xtensor<double, 2> kgeprime =
- {{0.813141},
- {0.812775},
- {0.812032},
- {0.811787},
- {0.811387}};
+ xt::xtensor<double, 3> kgeprime =
+ {{{0.813141}},
+ {{0.812775}},
+ {{0.812032}},
+ {{0.811787}},
+ {{0.811387}}};
EXPECT_TRUE(xt::allclose(metrics[3], kgeprime));
}
@@ -87,45 +87,45 @@ TEST(DeterministTests, TestTransform)
std::vector<xt::xarray<double>> metrics =
evalhyd::evald(observed, predicted, {"NSE"}, "sqrt");
- xt::xtensor<double, 2> nse_sqrt =
- {{0.882817},
- {0.883023},
- {0.883019},
- {0.883029},
- {0.882972}};
+ xt::xtensor<double, 3> nse_sqrt =
+ {{{0.882817}},
+ {{0.883023}},
+ {{0.883019}},
+ {{0.883029}},
+ {{0.882972}}};
EXPECT_TRUE(xt::allclose(metrics[0], nse_sqrt));
// compute and check results on inverted streamflow series
metrics = evalhyd::evald(observed, predicted, {"NSE"}, "inv");
- xt::xtensor<double, 2> nse_inv =
- {{0.737323},
- {0.737404},
- {0.737429},
- {0.737546},
- {0.737595}};
+ xt::xtensor<double, 3> nse_inv =
+ {{{0.737323}},
+ {{0.737404}},
+ {{0.737429}},
+ {{0.737546}},
+ {{0.737595}}};
EXPECT_TRUE(xt::allclose(metrics[0], nse_inv));
// compute and check results on square-rooted streamflow series
metrics = evalhyd::evald(observed, predicted, {"NSE"}, "log");
- xt::xtensor<double, 2> nse_log =
- {{0.893344},
- {0.893523},
- {0.893585},
- {0.893758},
- {0.893793}};
+ xt::xtensor<double, 3> nse_log =
+ {{{0.893344}},
+ {{0.893523}},
+ {{0.893585}},
+ {{0.893758}},
+ {{0.893793}}};
EXPECT_TRUE(xt::allclose(metrics[0], nse_log));
// compute and check results on power-transformed streamflow series
metrics = evalhyd::evald(observed, predicted, {"NSE"}, "pow", 0.2);
- xt::xtensor<double, 2> nse_pow =
- {{0.899207},
- {0.899395},
- {0.899451},
- {0.899578},
- {0.899588}};
+ xt::xtensor<double, 3> nse_pow =
+ {{{0.899207}},
+ {{0.899395}},
+ {{0.899451}},
+ {{0.899578}},
+ {{0.899588}}};
EXPECT_TRUE(xt::allclose(metrics[0], nse_pow));
}
@@ -180,8 +180,8 @@ TEST(DeterministTests, TestMaskingConditions)
// conditions on streamflow values _________________________________________
// compute scores using masking conditions on streamflow to subset whole record
- xt::xtensor<std::array<char, 32>, 2> q_conditions ={
- {{"q_obs{<2000,>3000}"}}
+ xt::xtensor<std::array<char, 32>, 1> q_conditions ={
+ {"q_obs{<2000,>3000}"}
};
std::vector<xt::xarray<double>> metrics_q_conditioned =
@@ -211,8 +211,8 @@ TEST(DeterministTests, TestMaskingConditions)
// conditions on streamflow statistics _____________________________________
// compute scores using masking conditions on streamflow to subset whole record
- xt::xtensor<std::array<char, 32>, 2> q_conditions_ ={
- {{"q_obs{>=mean}"}}
+ xt::xtensor<std::array<char, 32>, 1> q_conditions_ ={
+ {"q_obs{>=mean}"}
};
double mean = xt::mean(observed, {1})();
@@ -244,8 +244,8 @@ TEST(DeterministTests, TestMaskingConditions)
// conditions on temporal indices __________________________________________
// compute scores using masking conditions on time indices to subset whole record
- xt::xtensor<std::array<char, 32>, 2> t_conditions = {
- {{"t{0,1,2,3,4,5:97,97,98,99}"}}
+ xt::xtensor<std::array<char, 32>, 1> t_conditions = {
+ {"t{0,1,2,3,4,5:97,97,98,99}"}
};
std::vector<xt::xarray<double>> metrics_t_conditioned =