evald.cpp

#include <unordered_map>
#include <vector>
#include <array>
#include <stdexcept>
#include <xtensor/xexpression.hpp>
#include <xtensor/xarray.hpp>
#include <xtensor/xscalar.hpp>

#include "evalhyd/evald.hpp"

#include "utils.hpp"
#include "masks.hpp"
#include "maths.hpp"
#include "uncertainty.hpp"
#include "determinist/evaluator.hpp"

namespace eh = evalhyd;

namespace evalhyd
{
    std::vector<xt::xarray<double>> evald(
            const xt::xtensor<double, 2>& q_obs,
            const xt::xtensor<double, 2>& q_prd,
            const std::vector<std::string>& metrics,
            const std::string& transform,
            const double exponent,
            double epsilon,
            const xt::xtensor<bool, 2>& t_msk,
            const xt::xtensor<std::array<char, 32>, 1>& m_cdt,
            const std::unordered_map<std::string, int>& bootstrap,
            const std::vector<std::string>& dts
    )
    {
        // check that the metrics to be computed are valid
        utils::check_metrics(
                metrics,
                {"RMSE", "NSE", "KGE", "KGEPRIME"}
        );

        // check that optional parameters are valid
        eh::utils::check_bootstrap(bootstrap);

        // check that data dimensions are compatible
        // > time
        if (q_obs.shape(1) != q_prd.shape(1))
            throw std::runtime_error(
                    "observations and predictions feature different "
                    "temporal lengths"
            );
        if (t_msk.size() > 0)
            if (q_obs.shape(1) != t_msk.shape(1))
                throw std::runtime_error(
                        "observations and masks feature different "
                        "temporal lengths"
                );
        if (!dts.empty())
            if (q_obs.shape(1) != dts.size())
                throw std::runtime_error(
                        "observations and datetimes feature different "
                        "temporal lengths"
                );
        // > series
        if (q_obs.shape(0) != 1)
            throw std::runtime_error(
                    "observations contain more than one time series"
            );

        // retrieve dimensions
        std::size_t n_tim = q_obs.shape(1);
        std::size_t n_msk = t_msk.size() > 0 ? t_msk.shape(0) :