evaluator.hpp 18.42 KiB
#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