Commit 07ceab5a authored by Thibault Hallouin's avatar Thibault Hallouin
Browse files

add non-parametric KGE metric

1 merge request!3release v0.1.0
Pipeline #44369 passed with stage
in 3 minutes and 15 seconds
Showing with 333 additions and 1 deletion
+333 -1
......@@ -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;
}
}
}
}
......
......@@ -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();
};
};
}
}
......
......@@ -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;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment