-
Thibault Hallouin authorede785de1d
#ifndef EVALHYD_DETERMINIST_EVALUATOR_HPP#define EVALHYD_DETERMINIST_EVALUATOR_HPP#include <xtensor/xexpression.hpp>#include <xtensor/xtensor.hpp>namespace evalhyd{ namespace determinist { template <class A> class Evaluator { private: // members for input data const A& q_obs; const A& q_prd; // members for computational elements A mean_obs; A mean_prd; A quad_err; A quad_obs; A quad_prd; A r_pearson; A bias; public: // constructor method Evaluator(const xt::xexpression<A>& obs, const xt::xexpression<A>& prd) : q_obs{obs.derived_cast()}, q_prd{prd.derived_cast()} { // check that data dimensions are compatible if (q_prd.dimension() != q_obs.dimension()) { throw std::runtime_error( "number of dimensions of 'sim' and 'obs' " "must be identical" ); } }; // members for evaluation metrics A RMSE; A NSE; A KGE; A 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_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: ({1, ...}, time) // \assign mean_obs: // Mean observed streamflow. // shape: ({1, ...}, 1)7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
template <class A>
void Evaluator<A>::calc_mean_obs()
{
mean_obs = xt::mean(q_obs, -1, xt::keep_dims);
}
// Compute the mean of the predictions.
//
// \require q_prd:
// Streamflow predictions.
// shape: ({1+, ...}, time)
// \assign mean_prd:
// Mean predicted streamflow.
// shape: ({1+, ...}, 1)
template <class A>
void Evaluator<A>::calc_mean_prd()
{
mean_prd = xt::mean(q_prd, -1, xt::keep_dims);
}
// Compute the quadratic error between observations and predictions.
//
// \require q_obs:
// Streamflow observations.
// shape: ({1, ...}, time)
// \require q_prd:
// Streamflow predictions.
// shape: ({1+, ...}, time)
// \assign quad_err:
// Quadratic errors between observations and predictions.
// shape: ({1+, ...}, time)
template <class A>
void Evaluator<A>::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: ({1, ...}, time)
// \require mean_obs:
// Mean observed streamflow.
// shape: ({1, ...}, 1)
// \assign quad_obs:
// Quadratic errors between observations and mean observation.
// shape: ({1, ...}, time)
template <class A>
void Evaluator<A>::calc_quad_obs()
{
quad_obs = xt::square(q_obs - mean_obs);
}
// Compute the quadratic error between predictions and mean prediction.
//
// \require q_prd:
// Streamflow predictions.
// shape: ({1+, ...}, time)
// \require mean_prd:
// Mean predicted streamflow.
// shape: ({1+, ...}, 1)
// \assign quad_prd:
// Quadratic errors between predictions and mean prediction.
// shape: ({1+, ...}, time)
template <class A>
void Evaluator<A>::calc_quad_prd()
{
quad_prd = xt::square(q_prd - mean_prd);
}
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
// Compute the Pearson correlation coefficient.
//
// \require q_obs:
// Streamflow observations.
// shape: ({1, ...}, time)
// \require mean_obs:
// Mean observed streamflow.
// shape: ({1, ...}, 1)
// \require q_prd:
// Streamflow predictions.
// shape: ({1+, ...}, time)
// \require mean_prd:
// Mean predicted streamflow.
// shape: ({1+, ...}, 1)
// \require quad_obs:
// Quadratic errors between observations and mean observation.
// shape: ({1, ...}, time)
// \require quad_prd:
// Quadratic errors between predictions and mean prediction.
// shape: ({1+, ...}, time)
// \assign r_pearson:
// Pearson correlation coefficients.
// shape: ({1+, ...}, time)
template <class A>
void Evaluator<A>::calc_r_pearson()
{
// calculate error in timing and dynamics $r_{pearson}$
// (Pearson's correlation coefficient)
auto r_num = xt::sum(
(q_prd - mean_prd) * (q_obs - mean_obs), -1, xt::keep_dims
);
auto r_den = xt::sqrt(
xt::sum(quad_prd, -1, xt::keep_dims)
* xt::sum(quad_obs, -1, xt::keep_dims)
);
r_pearson = r_num / r_den;
}
// Compute the bias.
//
// \require q_obs:
// Streamflow observations.
// shape: ({1, ...}, time)
// \require q_prd:
// Streamflow predictions.
// shape: ({1+, ...}, time)
// \assign bias:
// Biases.
// shape: ({1+, ...}, 1)
template <class A>
void Evaluator<A>::calc_bias()
{
// calculate $bias$
bias = xt::sum(q_prd, -1, xt::keep_dims)
/ xt::sum(q_obs, -1, xt::keep_dims);
}
// Compute the root-mean-square error (RMSE).
//
// \require quad_err:
// Quadratic errors between observations and predictions.
// shape: ({1+, ...}, time)
// \assign RMSE:
// Root-mean-square errors.
// shape: ({1+, ...}, 1)
template <class A>
void Evaluator<A>::calc_RMSE()
{
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
// compute RMSE
RMSE = xt::sqrt(xt::mean(quad_err, -1, xt::keep_dims));
}
// Compute the Nash-Sutcliffe Efficiency (NSE).
//
// If multi-dimensional arrays are provided, the arrays of predictions
// and observations must feature the same number of dimensions, they
// must be broadcastable, and their temporal dimensions must be along
// their last axis.
//
// \require quad_err:
// Quadratic errors between observations and predictions.
// shape: ({1+, ...}, time)
// \require quad_obs:
// Quadratic errors between observations and mean observation.
// shape: ({1, ...}, time)
// \assign NSE:
// Nash-Sutcliffe efficiencies.
// shape: ({1+, ...}, 1)
template <class A>
void Evaluator<A>::calc_NSE()
{
// compute squared errors operands
A f_num = xt::sum(quad_err, -1, xt::keep_dims);
A f_den = xt::sum(quad_obs, -1, xt::keep_dims);
// compute NSE
NSE = 1 - (f_num / f_den);
}
// Compute the Kling-Gupta Efficiency (KGE).
//
// If multi-dimensional arrays are provided, the arrays of predictions
// and observations must feature the same number of dimensions, they
// must be broadcastable, and their temporal dimensions must be along
// their last axis.
//
// \require q_obs:
// Streamflow observations.
// shape: ({1, ...}, time)
// \require q_prd:
// Streamflow predictions.
// shape: ({1+, ...}, time)
// \require r_pearson:
// Pearson correlation coefficients.
// shape: ({1+, ...}, time)
// \require bias:
// Biases.
// shape: ({1+, ...}, 1)
// \assign KGE:
// Kling-Gupta efficiencies.
// shape: ({1+, ...}, 1)
template <class A>
void Evaluator<A>::calc_KGE()
{
// calculate error in spread of flow $alpha$
auto alpha =
xt::stddev(q_prd, {q_prd.dimension() - 1}, xt::keep_dims)
/ xt::stddev(q_obs, {q_obs.dimension() - 1}, xt::keep_dims);
// compute KGE
KGE = 1 - xt::sqrt(
xt::square(r_pearson - 1)
+ xt::square(alpha - 1)
+ xt::square(bias - 1)
);
}
// Compute the modified Kling-Gupta Efficiency (KGEPRIME).
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327
//
// If multi-dimensional arrays are provided, the arrays of predictions
// and observations must feature the same number of dimensions, they
// must be broadcastable, and their temporal dimensions must be along
// their last axis.
//
// \require q_obs:
// Streamflow observations.
// shape: ({1, ...}, time)
// \require mean_obs:
// Mean observed streamflow.
// shape: ({1, ...}, 1)
// \require q_prd:
// Streamflow predictions.
// shape: ({1+, ...}, time)
// \require mean_prd:
// Mean predicted streamflow.
// shape: ({1+, ...}, 1)
// \require r_pearson:
// Pearson correlation coefficients.
// shape: ({1+, ...}, time)
// \require bias:
// Biases.
// shape: ({1+, ...}, 1)
// \assign KGEPRIME:
// Modified Kling-Gupta efficiencies.
// shape: ({1+, ...}, 1)
template <class A>
void Evaluator<A>::calc_KGEPRIME()
{
// calculate error in spread of flow $gamma$
auto gamma =
(xt::stddev(q_prd, {q_prd.dimension() - 1}, xt::keep_dims) / mean_prd)
/ (xt::stddev(q_obs, {q_obs.dimension() - 1}, xt::keep_dims) / mean_obs);
// compute KGE
KGEPRIME = 1 - xt::sqrt(
xt::square(r_pearson - 1)
+ xt::square(gamma - 1)
+ xt::square(bias - 1)
);
}
}
}
#endif //EVALHYD_DETERMINIST_EVALUATOR_HPP