deterministic.hpp

#ifndef EVALHYD_DETERMINISTIC_HPP
#define EVALHYD_DETERMINISTIC_HPP

#include <vector>
#include <xtensor/xexpression.hpp>
#include <xtensor/xmath.hpp>

namespace evalhyd
{
    /// Compute the Nash-Sutcliffe Efficiency (NSE).
    ///
    /// If multi-dimensional arrays are provided, the arrays of simulations
    /// and observations must feature the same number of dimensions, they must
    /// be broadcastable, and their temporal dimensions must be along their
    /// last axis.
    ///
    /// \tparam A:
    ///     The type of data structures (e.g. `xt::xarray`, `xt::xtensor`).
    /// \param [in] sim:
    ///     Array of streamflow simulations.
    ///     shape: (..., time)
    /// \param [in] obs:
    ///     Array of streamflow observations.
    ///     shape: (1, time)
    /// \return
    ///     Array of computed efficiencies.
    ///     shape: (...)
    template <class A>
    inline A nse(
            const xt::xexpression<A>& sim,
            const xt::xexpression<A>& obs
    )
    {
        const A& q_sim = sim.derived_cast();
        const A& q_obs = obs.derived_cast();

        // check that data dimensions are compatible
        if (q_sim.dimension() != q_obs.dimension())
        {
            throw std::runtime_error(
                    "number of dimensions of 'sim' and 'obs' must be identical"
            );
        }

        // compute average observed flow
        A q_avg = xt::mean(q_obs, -1, xt::keep_dims);

        // compute squared errors operands
        A f_num = xt::sum(xt::square(q_obs - q_sim), -1, xt::keep_dims);
        A f_den = xt::sum(xt::square(q_obs - q_avg), -1, xt::keep_dims);

        // return computed NSE
        return 1 - (f_num / f_den);
    }
}

#endif //EVALHYD_DETERMINISTIC_HPP