evaluator.hpp 4.31 KiB
#ifndef EVALHYD_PROBABILIST_EVALUATOR_H
#define EVALHYD_PROBABILIST_EVALUATOR_H
#include <stdexcept>
#include <vector>
#include <xtensor/xtensor.hpp>
#include <xtensor/xview.hpp>
#include <xtensor/xslice.hpp>
using view1d_xtensor2d_double_type = decltype(
    xt::view(
            std::declval<const xt::xtensor<double, 2>&>(),
            std::declval<int>(),
            xt::all()
using view2d_xtensor4d_double_type = decltype(
    xt::view(
            std::declval<const xt::xtensor<double, 4>&>(),
            std::declval<int>(),
            std::declval<int>(),
            xt::all(),
            xt::all()
using view2d_xtensor4d_bool_type = decltype(
    xt::view(
            std::declval<const xt::xtensor<bool, 4>&>(),
            std::declval<int>(),
            std::declval<int>(),
            xt::all(),
            xt::all()
namespace evalhyd
    namespace probabilist
        class Evaluator
        private:
            // members for input data
            const view1d_xtensor2d_double_type& q_obs;
            const view2d_xtensor4d_double_type& q_prd;
            const view1d_xtensor2d_double_type& q_thr;
            xt::xtensor<bool, 2> t_msk;
            const std::vector<xt::xkeep_slice<int>>& b_exp;
            // members for dimensions
            std::size_t n;
            std::size_t n_msk;
            std::size_t n_mbr;
            std::size_t n_thr;
            std::size_t n_exp;
            // members for computational elements
            xt::xtensor<double, 2> o_k;
            xt::xtensor<double, 3> bar_o;
            xt::xtensor<double, 2> y_k;
            xt::xtensor<double, 2> q_qnt;
        public:
            // constructor method
            Evaluator(const view1d_xtensor2d_double_type& obs,
                      const view2d_xtensor4d_double_type& prd,
                      const view1d_xtensor2d_double_type& thr,
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
const view2d_xtensor4d_bool_type& msk, const std::vector<xt::xkeep_slice<int>>& exp) : q_obs{obs}, q_prd{prd}, q_thr{thr}, t_msk(msk), b_exp(exp) { // initialise a mask if none provided // (corresponding to no temporal subset) if (msk.size() < 1) t_msk = xt::ones<bool>({std::size_t {1}, q_obs.size()}); // determine size for recurring dimensions n = q_obs.size(); n_msk = t_msk.shape(0); n_mbr = q_prd.shape(0); n_thr = q_thr.size(); n_exp = b_exp.size(); // drop time steps where observations and/or predictions are NaN auto obs_nan = xt::isnan(q_obs); auto prd_nan = xt::isnan(q_prd); auto sum_nan = xt::eval(xt::sum(prd_nan, -1)); if (xt::amin(sum_nan) != xt::amax(sum_nan)) throw std::runtime_error( "predictions across members feature non-equal lengths" ); auto msk_nan = xt::where(obs_nan | xt::row(prd_nan, 0))[0]; xt::view(t_msk, xt::all(), xt::keep(msk_nan)) = false; }; // members for intermediate evaluation metrics // (i.e. before the reduction along the temporal axis) xt::xtensor<double, 2> bs; xt::xtensor<double, 2> qs; xt::xtensor<double, 2> crps; // members for evaluation metrics xt::xtensor<double, 3> BS; xt::xtensor<double, 4> BS_CRD; xt::xtensor<double, 4> BS_LBD; xt::xtensor<double, 3> BSS; xt::xtensor<double, 3> QS; xt::xtensor<double, 2> CRPS; // methods to compute derived data void calc_q_qnt(); // methods to compute elements void calc_o_k(); void calc_bar_o(); void calc_y_k(); // methods to compute intermediate metrics void calc_bs(); void calc_qs(); void calc_crps(); // methods to compute metrics void calc_BS(); void calc_BS_CRD(); void calc_BS_LBD(); void calc_BSS(); void calc_QS(); void calc_CRPS(); }; } } #endif //EVALHYD_PROBABILIST_EVALUATOR_H
141