An error occurred while loading the file. Please try again.
-
Thibault Hallouin authoredfixes #5
85ecd375
// Copyright (c) 2023, INRAE.// Distributed under the terms of the GPL-3 Licence.// The full licence is in the file LICENCE, distributed with this software.#ifndef EVALHYD_DETERMINIST_EFFICIENCIES_HPP#define EVALHYD_DETERMINIST_EFFICIENCIES_HPP#include <xtensor/xtensor.hpp>#include <xtensor/xview.hpp>#include <xtensor/xoperation.hpp>#include "../maths.hpp"namespace evalhyd{ namespace determinist { namespace elements { /// Compute the Pearson correlation coefficient. /// /// \param err_obs /// Errors between observations and mean observation. /// shape: (subsets, samples, series, time) /// \param err_prd /// Errors between predictions and mean prediction. /// shape: (subsets, samples, series, time) /// \param quad_err_obs /// Quadratic errors between observations and mean observation. /// shape: (subsets, samples, series, time) /// \param quad_err_prd /// Quadratic errors between predictions and mean prediction. /// shape: (subsets, samples, 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 /// Pearson correlation coefficients. /// shape: (subsets, samples, series) inline xt::xtensor<double, 3> calc_r_pearson( const xt::xtensor<double, 4>& err_obs, const xt::xtensor<double, 4>& err_prd, const xt::xtensor<double, 4>& quad_err_obs, const xt::xtensor<double, 4>& quad_err_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_{pearson}$ // (Pearson's correlation coefficient) xt::xtensor<double, 3> r_pearson = xt::zeros<double>({n_msk, n_exp, n_srs}); for (std::size_t m = 0; m < n_msk; m++) { // compute variable one sample at a time for (std::size_t e = 0; e < n_exp; e++)7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
{
auto prd = xt::view(err_prd, m, e, xt::all(), b_exp[e]);
auto obs = xt::view(err_obs, m, e, xt::all(), b_exp[e]);
auto r_num = xt::nansum(prd * obs, -1);
auto prd2 = xt::view(quad_err_prd, m, e, xt::all(), b_exp[e]);
auto obs2 = xt::view(quad_err_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;
}
}
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++)
{
// compute one series at a time because xt::sort does not
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
// consistently put NaN values at the end/beginning, so
// need to eliminate them before the sorting
for (std::size_t s = 0; s < n_srs; s++)
{
auto prd = xt::view(prd_masked, s, b_exp[e]);
auto obs = xt::view(obs_masked, s, b_exp[e]);
auto prd_filtered =
xt::filter(prd, !xt::isnan(prd));
auto obs_filtered =
xt::filter(obs, !xt::isnan(obs));
auto prd_sort = xt::argsort(
xt::eval(prd_filtered), {0},
xt::sorting_method::stable
);
auto obs_sort = xt::argsort(
xt::eval(obs_filtered), {0},
xt::sorting_method::stable
);
auto prd_rank = xt::eval(xt::argsort(prd_sort));
auto obs_rank = xt::eval(xt::argsort(obs_sort));
auto mean_rank = (prd_rank.size() - 1) / 2.;
auto prd_rank_err = xt::eval(prd_rank - mean_rank);
auto obs_rank_err = xt::eval(obs_rank - mean_rank);
auto r_num = xt::nansum(prd_rank_err * obs_rank_err);
auto r_den = xt::sqrt(
xt::nansum(xt::square(prd_rank_err))
* xt::nansum(xt::square(obs_rank_err))
);
xt::view(r_spearman, m, e, s) = r_num / r_den;
}
}
}
return r_spearman;
}
/// Compute 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.
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
/// \return
/// Alphas, ratios of standard deviations.
/// shape: (subsets, samples, series)
template <class XD2>
inline xt::xtensor<double, 3> calc_alpha(
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
)
{
// calculate error in spread of flow $alpha$
xt::xtensor<double, 3> alpha =
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]);
xt::view(alpha, m, e) =
maths::nanstd(prd, xt::view(mean_prd, m, e))
/ maths::nanstd(obs, xt::view(mean_obs, m, e));
}
}
return alpha;
}
/// Compute gamma.
///
/// \param mean_obs
/// Mean observed streamflow.
/// shape: (subsets, samples, series, 1)
/// \param mean_prd
/// Mean predicted streamflow.
/// shape: (subsets, samples, series, 1)
/// \param alpha
/// Alphas, ratios of standard deviations.
/// shape: (subsets, samples, series)
/// \return
/// Gammas, ratios of standard deviations normalised by
/// their means.
/// shape: (subsets, samples, series)
inline xt::xtensor<double, 3> calc_gamma(
const xt::xtensor<double, 4>& mean_obs,
const xt::xtensor<double, 4>& mean_prd,
const xt::xtensor<double, 3>& alpha
)
{
// calculate normalised error in spread of flow $gamma$
xt::xtensor<double, 3> gamma =
alpha * (xt::view(mean_obs, xt::all(), xt::all(), xt::all(), 0)
/ xt::view(mean_prd, xt::all(), xt::all(), xt::all(), 0));
return gamma;
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
}
/// 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
)
{
// 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++)
{
// compute one series at a time because xt::sort does not
// consistently put NaN values at the end/beginning, so
// need to eliminate them before the sorting
for (std::size_t s = 0; s < n_srs; s++)
{
auto prd = xt::view(prd_masked, s, b_exp[e]);
auto obs = xt::view(obs_masked, s, b_exp[e]);
auto prd_filtered =
xt::filter(prd, !xt::isnan(prd));
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
auto obs_filtered =
xt::filter(obs, !xt::isnan(obs));
auto prd_fdc = xt::sort(
xt::eval(prd_filtered
/ (prd_filtered.size()
* xt::view(mean_prd, m, e, s)))
);
auto obs_fdc = xt::sort(
xt::eval(obs_filtered
/ (obs_filtered.size()
* xt::view(mean_obs, m, e, s)))
);
xt::view(alpha_np, m, e, s) =
1 - 0.5 * xt::nansum(xt::abs(prd_fdc - obs_fdc));
}
}
}
return alpha_np;
}
/// Compute the bias.
///
/// \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
/// Biases.
/// shape: (subsets, samples, series)
template <class XD2>
inline xt::xtensor<double, 3> calc_bias(
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 $bias$
xt::xtensor<double, 3> bias =
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);
421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490
// 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]);
xt::view(bias, m, e) =
xt::nansum(prd, -1) / xt::nansum(obs, -1);
}
}
return bias;
}
}
namespace metrics
{
/// Compute the Nash-Sutcliffe Efficiency (NSE).
///
/// \param quad_err
/// Quadratic errors between observations and predictions.
/// shape: (series, time)
/// \param quad_err_obs
/// Quadratic errors between observations and mean observation.
/// shape: (subsets, samples, 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
/// Nash-Sutcliffe efficiencies.
/// shape: (series, subsets, samples)
inline xt::xtensor<double, 3> calc_NSE(
const xt::xtensor<double, 2>& quad_err,
const xt::xtensor<double, 4>& quad_err_obs,
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
)
{
xt::xtensor<double, 3> NSE =
xt::zeros<double>({n_srs, n_msk, n_exp});
for (std::size_t 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, xt::all(), m),
quad_err, NAN);
// compute variable one sample at a time
for (std::size_t 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_err_obs, m, e, xt::all(), b_exp[e]);
xt::xtensor<double, 1> f_den =
xt::nansum(obs2, -1);
491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560
// compute NSE
xt::view(NSE, xt::all(), m, e) = 1 - (f_num / f_den);
}
}
return NSE;
}
/// Compute the Kling-Gupta Efficiency (KGE).
///
/// \param r_pearson
/// Pearson correlation coefficients.
/// shape: (subsets, samples, series)
/// \param alpha
/// Alphas, ratios of standard deviations.
/// 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
/// Kling-Gupta efficiencies.
/// shape: (series, subsets, samples)
inline xt::xtensor<double, 3> calc_KGE(
const xt::xtensor<double, 3>& r_pearson,
const xt::xtensor<double, 3>& alpha,
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> KGE =
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 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)
);
}
}
return KGE;
}
/// Compute the Kling-Gupta Efficiency Decomposed (KGE_D) into
/// its three components that are the linear correlation (r),
/// the variability (alpha), and the bias (beta), in this order.
///
/// \param r_pearson
/// Pearson correlation coefficients.
/// shape: (subsets, samples, series)
/// \param alpha
/// Alphas, ratios of standard deviations.
/// shape: (subsets, samples, series)
/// \param bias
/// Biases.
/// shape: (subsets, samples, series)
/// \param n_srs
561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630
/// Number of prediction series.
/// \param n_msk
/// Number of temporal subsets.
/// \param n_exp
/// Number of bootstrap samples.
/// \return
/// KGE components (r, alpha, beta) for each subset
/// and for each threshold.
/// shape: (series, subsets, samples, 3)
inline xt::xtensor<double, 4> calc_KGE_D(
const xt::xtensor<double, 3>& r_pearson,
const xt::xtensor<double, 3>& alpha,
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> KGE_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(KGE_D, xt::all(), m, e, 0) =
xt::view(r_pearson, m, e);
xt::view(KGE_D, xt::all(), m, e, 1) =
xt::view(alpha, m, e);
xt::view(KGE_D, xt::all(), m, e, 2) =
xt::view(bias, m, e);
}
}
return KGE_D;
}
/// Compute the modified Kling-Gupta Efficiency (KGEPRIME).
///
/// \param r_pearson
/// Pearson correlation coefficients.
/// shape: (subsets, samples, series)
/// \param gamma
/// Gammas, ratios of standard deviations normalised by
/// their means.
/// 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_KGEPRIME(
const xt::xtensor<double, 3>& r_pearson,
const xt::xtensor<double, 3>& gamma,
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> KGEPRIME =
xt::zeros<double>({n_srs, n_msk, n_exp});
631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700
for (std::size_t m = 0; m < n_msk; m++)
{
for (std::size_t e = 0; e < n_exp; e++)
{
// compute KGEPRIME
xt::view(KGEPRIME, xt::all(), m, e) = 1 - xt::sqrt(
xt::square(xt::view(r_pearson, m, e) - 1)
+ xt::square(xt::view(gamma, m, e) - 1)
+ xt::square(xt::view(bias, m, e) - 1)
);
}
}
return KGEPRIME;
}
/// Compute the modified Kling-Gupta Efficiency Decomposed
/// (KGEPRIME_D) into its three components that are the linear
/// correlation (r), the variability (gamma), and the bias (beta),
/// in this order.
///
/// \param r_pearson
/// Pearson correlation coefficients.
/// shape: (subsets, samples, series)
/// \param gamma
/// Gammas, ratios of standard deviations normalised by
/// their means.
/// 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 KGE components (r, gamma, beta) for each subset
/// and for each threshold.
/// shape: (series, subsets, samples, 3)
inline xt::xtensor<double, 4> calc_KGEPRIME_D(
const xt::xtensor<double, 3>& r_pearson,
const xt::xtensor<double, 3>& gamma,
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> KGEPRIME_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(KGEPRIME_D, xt::all(), m, e, 0) =
xt::view(r_pearson, m, e);
xt::view(KGEPRIME_D, xt::all(), m, e, 1) =
xt::view(gamma, m, e);
xt::view(KGEPRIME_D, xt::all(), m, e, 2) =
xt::view(bias, m, e);
}
}
return KGEPRIME_D;
}
701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770
/// 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
771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806
/// 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;
}
}
}
}
#endif //EVALHYD_DETERMINIST_EFFICIENCIES_HPP