An error occurred while loading the file. Please try again.
-
Thibault Hallouin authored24fcb9c4
#ifndef EVALHYD_DETERMINIST_EVALUATOR_HPP
#define EVALHYD_DETERMINIST_EVALUATOR_HPP
#include <vector>
#include <xtensor/xexpression.hpp>
#include <xtensor/xtensor.hpp>
#include "../maths.hpp"
namespace evalhyd
{
namespace determinist
{
class Evaluator
{
private:
// members for input data
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_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;
std::vector<std::size_t> mean_dims;
std::vector<std::size_t> final_dims;
// members for computational elements
xt::xtensor<double, 4> mean_obs;
xt::xtensor<double, 4> mean_prd;
xt::xtensor<double, 2> quad_err;
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
Evaluator(const xt::xtensor<double, 2>& obs,
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}, b_exp{exp}
{
// 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, 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);
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
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, 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();
void calc_mean_prd();
void calc_quad_err();
void calc_quad_obs();
void calc_quad_prd();
void calc_r_pearson();
void calc_alpha();
void calc_bias();
// methods to compute metrics
void calc_RMSE();
void calc_NSE();
void calc_KGE();
void calc_KGEPRIME();
};
// Compute the mean of the observations.
//
// \require q_obs:
// Streamflow observations.
// shape: (series, time)
// \assign mean_obs:
// Mean observed streamflow.
// shape: (subsets, samples, series, 1)
void Evaluator::calc_mean_obs()
{
mean_obs = xt::zeros<double>(mean_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);
}
}
}
// Compute the mean of the predictions.
//
// \require q_prd:
// Streamflow predictions.
// shape: (series, time)
// \assign mean_prd:
// Mean predicted streamflow.
// shape: (subsets, samples, series, 1)
void Evaluator::calc_mean_prd()
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
{
mean_prd = xt::zeros<double>(mean_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);
}
}
}
// Compute the quadratic error between observations and predictions.
//
// \require q_obs:
// Streamflow observations.
// shape: (series, time)
// \require q_prd:
// Streamflow predictions.
// shape: (series, time)
// \assign quad_err:
// Quadratic errors between observations and predictions.
// shape: (series, time)
void Evaluator::calc_quad_err()
{
quad_err = xt::square(q_obs - q_prd);
}
// Compute the quadratic error between observations and mean observation.
//
// \require q_obs:
// Streamflow observations.
// shape: (series, time)
// \require mean_obs:
// Mean observed streamflow.
// shape: (samples, series, 1)
// \assign quad_obs:
// Quadratic errors between observations and mean observation.
// shape: (subsets, samples, series, time)
void Evaluator::calc_quad_obs()
{
quad_obs = xt::zeros<double>(inter_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);
for (int e = 0; e < n_exp; e++) {
xt::view(quad_obs, m, e) = xt::square(
obs_masked - xt::view(mean_obs, m, e)
);
}
}
}
// Compute the quadratic error between predictions and mean prediction.
//
// \require q_prd:
// Streamflow predictions.
// shape: (series, time)
// \require mean_prd:
// Mean predicted streamflow.
// shape: (samples, series, 1)
// \assign quad_prd:
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
// Quadratic errors between predictions and mean prediction.
// shape: (subsets, samples, series, time)
void Evaluator::calc_quad_prd()
{
quad_prd = xt::zeros<double>(inter_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);
for (int e = 0; e < n_exp; e++) {
xt::view(quad_prd, m, e) = xt::square(
prd_masked - xt::view(mean_prd, m, e)
);
}
}
}
// Compute the Pearson correlation coefficient.
//
// \require q_obs:
// Streamflow observations.
// shape: (series, time)
// \require q_prd:
// Streamflow predictions.
// shape: (series, time)
// \require mean_obs:
// Mean observed streamflow.
// shape: (samples, series, 1)
// \require mean_prd:
// Mean predicted streamflow.
// shape: (samples, series, 1)
// \require quad_obs:
// Quadratic errors between observations and mean observation.
// shape: (samples, series, time)
// \require quad_prd:
// Quadratic errors between predictions and mean prediction.
// shape: (samples, series, time)
// \assign r_pearson:
// Pearson correlation coefficients.
// 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 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;
}
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
}
}
// Compute alpha.
//
// \require q_obs:
// Streamflow observations.
// shape: (series, time)
// \require q_prd:
// Streamflow predictions.
// shape: (series, time)
// \require mean_obs:
// Mean observed streamflow.
// shape: (samples, series, 1)
// \require mean_prd:
// Mean predicted streamflow.
// shape: (samples, series, 1)
// \assign alpha:
// Alphas, ratios of standard deviations.
// shape: (subsets, samples, series)
void Evaluator::calc_alpha()
{
// calculate error in spread of flow $alpha$
alpha = xt::zeros<double>(inner_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);
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) =
evalhyd::maths::nanstd(prd, xt::view(mean_prd, m, e))
/ evalhyd::maths::nanstd(obs, xt::view(mean_obs, m, e));
}
}
}
// Compute the bias.
//
// \require q_obs:
// Streamflow observations.
// shape: (series, time)
// \require q_prd:
// Streamflow predictions.
// shape: (series, time)
// \assign bias:
// Biases.
// shape: (subsets, samples, series)
void Evaluator::calc_bias()
{
// calculate $bias$
bias = xt::zeros<double>(inner_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);
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);
}
}
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
}
// Compute the root-mean-square error (RMSE).
//
// \require quad_err:
// Quadratic errors between observations and predictions.
// shape: (series, time)
// \assign RMSE:
// Root-mean-square errors.
// shape: (series, subsets, samples)
void Evaluator::calc_RMSE()
{
// compute RMSE
RMSE = xt::zeros<double>(final_dims);
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));
}
}
}
// Compute the Nash-Sutcliffe Efficiency (NSE).
//
// \require quad_err:
// Quadratic errors between observations and predictions.
// shape: (series, time)
// \require quad_obs:
// Quadratic errors between observations and mean observation.
// shape: (samples, series, time)
// \assign NSE:
// Nash-Sutcliffe efficiencies.
// shape: (series, subsets, samples)
void Evaluator::calc_NSE()
{
NSE = xt::zeros<double>(final_dims);
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);
}
}
}
// Compute the Kling-Gupta Efficiency (KGE).
//
// \require r_pearson:
// Pearson correlation coefficients.
// shape: (samples, series)
// \require alpha:
// Alphas, ratios of standard deviations.
// shape: (samples, series)
421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485
// \require bias:
// Biases.
// shape: (samples, series)
// \assign KGE:
// Kling-Gupta efficiencies.
// shape: (series, subsets, samples)
void Evaluator::calc_KGE()
{
KGE = xt::zeros<double>(final_dims);
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)
);
}
}
}
// Compute the modified Kling-Gupta Efficiency (KGEPRIME).
//
// \require mean_obs:
// Mean observed streamflow.
// shape: (samples, series, 1)
// \require mean_prd:
// Mean predicted streamflow.
// shape: (samples, series, 1)
// \require r_pearson:
// Pearson correlation coefficients.
// shape: (samples, series)
// \require alpha:
// Alphas, ratios of standard deviations.
// shape: (samples, series)
// \require bias:
// Biases.
// shape: (samples, series)
// \assign KGEPRIME:
// Modified Kling-Gupta efficiencies.
// shape: (series, subsets, samples)
void Evaluator::calc_KGEPRIME()
{
KGEPRIME = xt::zeros<double>(final_dims);
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)
);
}
}
}
}
}
#endif //EVALHYD_DETERMINIST_EVALUATOR_HPP