evaluator_elements.cpp 3.51 KiB
#include <xtensor/xmath.hpp>
#include <xtensor/xview.hpp>
#include <xtensor/xsort.hpp>
#include "probabilist/evaluator.hpp"
namespace evalhyd
    namespace probabilist
        // Determine observed realisation of threshold(s) exceedance.
        // \require q_obs:
        //     Streamflow observations.
        //     shape: (time,)
        // \require q_thr:
        //     Streamflow exceedance threshold(s).
        //     shape: (thresholds,)
        // \assign o_k:
        //     Event observed outcome.
        //     shape: (thresholds, time)
        void Evaluator::calc_o_k()
            // determine observed realisation of threshold(s) exceedance
            o_k = q_obs >= xt::view(q_thr, xt::all(), xt::newaxis());
        // Determine mean observed realisation of threshold(s) exceedance.
        // \require o_k:
        //     Event observed outcome.
        //     shape: (thresholds, time)
        // \require t_msk:
        //     Temporal subsets of the whole record.
        //     shape: (subsets, time)
        // \assign bar_o:
        //     Mean event observed outcome.
        //     shape: (subsets, samples, thresholds)
        void Evaluator::calc_bar_o()
            // apply the mask
            // (using NaN workaround until reducers work on masked_view)
            auto o_k_masked = xt::where(
                    xt::view(t_msk, xt::all(), xt::newaxis(), xt::all()),
                    o_k, NAN
            // compute variable one sample at a time
            bar_o = xt::zeros<double>({n_msk, n_exp, n_thr});
            for (int e = 0; e < n_exp; e++)
                // apply the bootstrap sampling
                auto o_k_masked_sampled =
                        xt::view(o_k_masked, xt::all(), xt::all(), b_exp[e]);
                // compute mean "climatology" relative frequency of the event
                // $\bar{o} = \frac{1}{n} \sum_{k=1}^{n} o_k$
                xt::view(bar_o, xt::all(), e, xt::all()) =
                        xt::nanmean(o_k_masked_sampled, -1);
        // Determine forecast probability of threshold(s) exceedance to occur.
        // \require q_prd:
        //     Streamflow predictions.
        //     shape: (members, time)
        // \require q_thr:
        //     Streamflow exceedance threshold(s).