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: