diff --git a/include/evalhyd/detail/determinist/efficiencies.hpp b/include/evalhyd/detail/determinist/efficiencies.hpp
index 034c84fdd30347cc390721aa57aaed48d430d029..93281ba80402a041cf10ae44dec048395a715dd9 100644
--- a/include/evalhyd/detail/determinist/efficiencies.hpp
+++ b/include/evalhyd/detail/determinist/efficiencies.hpp
@@ -83,6 +83,94 @@ namespace evalhyd
return r_pearson;
}
+ /// Compute the Spearman rank correlation coefficient.
+ ///
+ /// \param q_obs
+ /// Streamflow observations.
+ /// shape: (series, time)
+ /// \param q_prd
+ /// Streamflow predictions.
+ /// shape: (series, time)
+ /// \param t_msk
+ /// Temporal subsets of the whole record.
+ /// shape: (series, subsets, time)
+ /// \param b_exp
+ /// Boostrap samples.
+ /// shape: (samples, time slice)
+ /// \param n_srs
+ /// Number of prediction series.
+ /// \param n_msk
+ /// Number of temporal subsets.
+ /// \param n_exp
+ /// Number of bootstrap samples.
+ /// \return
+ /// Spearman rank correlation coefficients.
+ /// shape: (subsets, samples, series)
+ template <class XD2>
+ inline xt::xtensor<double, 3> calc_r_spearman(
+ const XD2& q_obs,
+ const XD2& q_prd,
+ const xt::xtensor<bool, 3>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_srs,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // calculate error in timing and dynamics $r_{spearman}$
+ // (Spearman's rank correlation coefficient)
+ xt::xtensor<double, 3> r_spearman =
+ xt::zeros<double>({n_msk, n_exp, n_srs});
+
+ for (std::size_t 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, xt::all(), m),
+ q_prd, NAN);
+ auto obs_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ q_obs, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t 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 prd_rank = xt::eval(xt::argsort(
+ xt::argsort(xt::eval(prd), {1}),
+ {1}
+ ));
+ auto obs_rank = xt::eval(xt::argsort(
+ xt::argsort(xt::eval(obs), {1}),
+ {1}
+ ));
+
+ auto mean_prd_rank =
+ xt::eval(xt::nanmean(prd_rank, {1}, xt::keep_dims));
+ auto mean_obs_rank =
+ xt::eval(xt::nanmean(obs_rank, {1}, xt::keep_dims));
+
+ auto prd_rank_err = xt::eval(prd_rank - mean_prd_rank);
+ auto obs_rank_err = xt::eval(obs_rank - mean_obs_rank);
+
+ auto r_num = xt::nansum(
+ prd_rank_err * obs_rank_err,
+ {1}
+ );
+
+ auto r_den = xt::sqrt(
+ xt::nansum(xt::square(prd_rank_err), {1})
+ * xt::nansum(xt::square(obs_rank_err), {1})
+ );
+
+ xt::view(r_spearman, m, e) = r_num / r_den;
+ }
+ }
+
+ return r_spearman;
+ }
+
/// Compute alpha.
///
/// \param q_obs
@@ -181,6 +269,85 @@ namespace evalhyd
return gamma;
}
+ /// Compute non-parametric alpha.
+ ///
+ /// \param q_obs
+ /// Streamflow observations.
+ /// shape: (series, time)
+ /// \param q_prd
+ /// Streamflow predictions.
+ /// shape: (series, time)
+ /// \param mean_obs
+ /// Mean observed streamflow.
+ /// shape: (subsets, samples, series, 1)
+ /// \param mean_prd
+ /// Mean predicted streamflow.
+ /// shape: (subsets, samples, series, 1)
+ /// \param t_msk
+ /// Temporal subsets of the whole record.
+ /// shape: (series, subsets, time)
+ /// \param b_exp
+ /// Boostrap samples.
+ /// shape: (samples, time slice)
+ /// \param n_srs
+ /// Number of prediction series.
+ /// \param n_msk
+ /// Number of temporal subsets.
+ /// \param n_exp
+ /// Number of bootstrap samples.
+ /// \return
+ /// Non-parametric alphas.
+ /// shape: (subsets, samples, series)
+ template <class XD2>
+ inline xt::xtensor<double, 3> calc_alpha_np(
+ const XD2& q_obs,
+ const XD2& q_prd,
+ const xt::xtensor<double, 4>& mean_obs,
+ const xt::xtensor<double, 4>& mean_prd,
+ const xt::xtensor<bool, 3>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_srs,
+ std::size_t n_msk,
+ std::size_t n_exp,
+ std::size_t n_tim
+ )
+ {
+ // calculate error in spread of flow $alpha$
+ xt::xtensor<double, 3> alpha_np =
+ xt::zeros<double>({n_msk, n_exp, n_srs});
+
+ for (std::size_t 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, xt::all(), m),
+ q_prd, NAN);
+ auto obs_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ q_obs, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t 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 prd_fdc = xt::sort(
+ xt::eval(prd / (n_tim * xt::view(mean_prd, m, e))),
+ {1}
+ );
+ auto obs_fdc = xt::sort(
+ xt::eval(obs / (n_tim * xt::view(mean_obs, m, e))),
+ {1}
+ );
+
+ xt::view(alpha_np, m, e) =
+ 1 - 0.5 * xt::nansum(xt::abs(prd_fdc - obs_fdc), {1});
+ }
+ }
+
+ return alpha_np;
+ }
+
/// Compute the bias.
///
/// \param q_obs
@@ -508,6 +675,106 @@ namespace evalhyd
return KGEPRIME_D;
}
+
+ /// Compute the non-parametric Kling-Gupta Efficiency (KGENP).
+ ///
+ /// \param r_spearman
+ /// Spearman rank correlation coefficients.
+ /// shape: (subsets, samples, series)
+ /// \param alpha_np
+ /// Non-parametric alphas.
+ /// shape: (subsets, samples, series)
+ /// \param bias
+ /// Biases.
+ /// shape: (subsets, samples, series)
+ /// \param n_srs
+ /// Number of prediction series.
+ /// \param n_msk
+ /// Number of temporal subsets.
+ /// \param n_exp
+ /// Number of bootstrap samples.
+ /// \return
+ /// Modified Kling-Gupta efficiencies.
+ /// shape: (series, subsets, samples)
+ inline xt::xtensor<double, 3> calc_KGENP(
+ const xt::xtensor<double, 3>& r_spearman,
+ const xt::xtensor<double, 3>& alpha_np,
+ const xt::xtensor<double, 3>& bias,
+ std::size_t n_srs,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ xt::xtensor<double, 3> KGENP =
+ xt::zeros<double>({n_srs, n_msk, n_exp});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // compute KGEPRIME
+ xt::view(KGENP, xt::all(), m, e) = 1 - xt::sqrt(
+ xt::square(xt::view(r_spearman, m, e) - 1)
+ + xt::square(xt::view(alpha_np, m, e) - 1)
+ + xt::square(xt::view(bias, m, e) - 1)
+ );
+ }
+ }
+
+ return KGENP;
+ }
+
+ /// Compute the non-parametric Kling-Gupta Efficiency
+ /// Decomposed (KGENP) into its three components that are the rank
+ /// correlation (r), the variability (non-parametric alpha), and
+ /// the bias (beta), in this order.
+ ///
+ /// \param r_spearman
+ /// Spearman correlation coefficients.
+ /// shape: (subsets, samples, series)
+ /// \param alpha_np
+ /// Non-parametric alphas.
+ /// shape: (subsets, samples, series)
+ /// \param bias
+ /// Biases.
+ /// shape: (subsets, samples, series)
+ /// \param n_srs
+ /// Number of prediction series.
+ /// \param n_msk
+ /// Number of temporal subsets.
+ /// \param n_exp
+ /// Number of bootstrap samples.
+ /// \return
+ /// Modified Kling-Gupta efficiencies.
+ /// shape: (series, subsets, samples)
+ inline xt::xtensor<double, 4> calc_KGENP_D(
+ const xt::xtensor<double, 3>& r_spearman,
+ const xt::xtensor<double, 3>& alpha_np,
+ const xt::xtensor<double, 3>& bias,
+ std::size_t n_srs,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ xt::xtensor<double, 4> KGENP_D =
+ xt::zeros<double>({n_srs, n_msk, n_exp, std::size_t {3}});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // put KGE components together
+ xt::view(KGENP_D, xt::all(), m, e, 0) =
+ xt::view(r_spearman, m, e);
+ xt::view(KGENP_D, xt::all(), m, e, 1) =
+ xt::view(alpha_np, m, e);
+ xt::view(KGENP_D, xt::all(), m, e, 2) =
+ xt::view(bias, m, e);
+ }
+ }
+
+ return KGENP_D;
+ }
}
}
}
diff --git a/include/evalhyd/detail/determinist/evaluator.hpp b/include/evalhyd/detail/determinist/evaluator.hpp
index 8276179cc3508b48ff73ce15cb86e88dadb07b2a..fb38f5bd6d4944289f385ecd28c071e6937427d0 100644
--- a/include/evalhyd/detail/determinist/evaluator.hpp
+++ b/include/evalhyd/detail/determinist/evaluator.hpp
@@ -47,8 +47,10 @@ namespace evalhyd
xtl::xoptional<xt::xtensor<double, 4>, bool> err_prd;
xtl::xoptional<xt::xtensor<double, 4>, bool> quad_err_prd;
xtl::xoptional<xt::xtensor<double, 3>, bool> r_pearson;
+ xtl::xoptional<xt::xtensor<double, 3>, bool> r_spearman;
xtl::xoptional<xt::xtensor<double, 3>, bool> alpha;
xtl::xoptional<xt::xtensor<double, 3>, bool> gamma;
+ xtl::xoptional<xt::xtensor<double, 3>, bool> alpha_np;
xtl::xoptional<xt::xtensor<double, 3>, bool> bias;
// members for evaluation metrics
@@ -61,6 +63,8 @@ namespace evalhyd
xtl::xoptional<xt::xtensor<double, 4>, bool> KGE_D;
xtl::xoptional<xt::xtensor<double, 3>, bool> KGEPRIME;
xtl::xoptional<xt::xtensor<double, 4>, bool> KGEPRIME_D;
+ xtl::xoptional<xt::xtensor<double, 3>, bool> KGENP;
+ xtl::xoptional<xt::xtensor<double, 4>, bool> KGENP_D;
// methods to compute elements
xt::xtensor<double, 4> get_mean_obs()
@@ -178,6 +182,18 @@ namespace evalhyd
return r_pearson.value();
};
+ xt::xtensor<double, 3> get_r_spearman()
+ {
+ if (!r_spearman.has_value())
+ {
+ r_spearman = elements::calc_r_spearman(
+ q_obs, q_prd, t_msk, b_exp,
+ n_srs, n_msk, n_exp
+ );
+ }
+ return r_spearman.value();
+ };
+
xt::xtensor<double, 3> get_alpha()
{
if (!alpha.has_value())
@@ -201,6 +217,18 @@ namespace evalhyd
return gamma.value();
};
+ xt::xtensor<double, 3> get_alpha_np()
+ {
+ if (!alpha_np.has_value())
+ {
+ alpha_np = elements::calc_alpha_np(
+ q_obs, q_prd, get_mean_obs(), get_mean_prd(),
+ t_msk, b_exp, n_srs, n_msk, n_exp, n_tim
+ );
+ }
+ return alpha_np.value();
+ };
+
xt::xtensor<double, 3> get_bias()
{
if (!bias.has_value())
@@ -351,6 +379,30 @@ namespace evalhyd
}
return KGEPRIME_D.value();
};
+
+ xt::xtensor<double, 3> get_KGENP()
+ {
+ if (!KGENP.has_value())
+ {
+ KGENP = metrics::calc_KGENP(
+ get_r_spearman(), get_alpha_np(), get_bias(),
+ n_srs, n_msk, n_exp
+ );
+ }
+ return KGENP.value();
+ };
+
+ xt::xtensor<double, 4> get_KGENP_D()
+ {
+ if (!KGENP_D.has_value())
+ {
+ KGENP_D = metrics::calc_KGENP_D(
+ get_r_spearman(), get_alpha_np(), get_bias(),
+ n_srs, n_msk, n_exp
+ );
+ }
+ return KGENP_D.value();
+ };
};
}
}
diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp
index ef858dfcfe97d058e7351155dd711f152a9f3d1b..76d2a8fcda8d8ca12e2f705fb6e9696cb7268123 100644
--- a/include/evalhyd/evald.hpp
+++ b/include/evalhyd/evald.hpp
@@ -209,7 +209,8 @@ namespace evalhyd
utils::check_metrics(
metrics,
{"MAE", "MARE", "MSE", "RMSE",
- "NSE", "KGE", "KGE_D", "KGEPRIME", "KGEPRIME_D"}
+ "NSE", "KGE", "KGE_D", "KGEPRIME", "KGEPRIME_D",
+ "KGENP", "KGENP_D"}
);
// check that optional parameters are valid
@@ -489,6 +490,18 @@ namespace evalhyd
uncertainty::summarise_d(evaluator.get_KGEPRIME_D(), summary)
);
}
+ else if ( metric == "KGENP" )
+ {
+ r.emplace_back(
+ uncertainty::summarise_d(evaluator.get_KGENP(), summary)
+ );
+ }
+ else if ( metric == "KGENP_D" )
+ {
+ r.emplace_back(
+ uncertainty::summarise_d(evaluator.get_KGENP_D(), summary)
+ );
+ }
}
return r;