Skip to content
GitLab
Select Git revision
  • master default
Find file BlameHistoryPermalink
evald.cpp 9.10 KiB
1 
#include <unordered_map>
2 
#include <vector>
3 
#include <array>
4 
#include <stdexcept>
5 
#include <xtensor/xexpression.hpp>
6 
#include <xtensor/xarray.hpp>
7 
#include <xtensor/xscalar.hpp>
8
9 
#include "evalhyd/evald.hpp"
10
11 
#include "utils.hpp"
12 
#include "masks.hpp"
13 
#include "maths.hpp"
14 
#include "uncertainty.hpp"
15 
#include "determinist/evaluator.hpp"
16
17 
namespace eh = evalhyd;
18
19 
namespace evalhyd
20 
{
21 
    std::vector<xt::xarray<double>> evald(
22 
            const xt::xtensor<double, 2>& q_obs,
23 
            const xt::xtensor<double, 2>& q_prd,
24 
            const std::vector<std::string>& metrics,
25 
            const std::string& transform,
26 
            const double exponent,
27 
            double epsilon,
28 
            const xt::xtensor<bool, 2>& t_msk,
29 
            const xt::xtensor<std::array<char, 32>, 1>& m_cdt,
30 
            const std::unordered_map<std::string, int>& bootstrap,
31 
            const std::vector<std::string>& dts
32 
    )
33 
    {
34 
        // check that the metrics to be computed are valid
35 
        utils::check_metrics(
36 
                metrics,
37 
                {"RMSE", "NSE", "KGE", "KGEPRIME"}
38 
        );
39
40 
        // check that optional parameters are valid
41 
        eh::utils::check_bootstrap(bootstrap);
42
43 
        // check that data dimensions are compatible
44 
        // > time
45 
        if (q_obs.shape(1) != q_prd.shape(1))
46 
            throw std::runtime_error(
47 
                    "observations and predictions feature different "
48 
                    "temporal lengths"
49 
            );
50 
        if (t_msk.size() > 0)
51 
            if (q_obs.shape(1) != t_msk.shape(1))
52 
                throw std::runtime_error(
53 
                        "observations and masks feature different "
54 
                        "temporal lengths"
55 
                );
56 
        if (!dts.empty())
57 
            if (q_obs.shape(1) != dts.size())
58 
                throw std::runtime_error(
59 
                        "observations and datetimes feature different "
60 
                        "temporal lengths"
61 
                );
62 
        // > series
63 
        if (q_obs.shape(0) != 1)
64 
            throw std::runtime_error(
65 
                    "observations contain more than one time series"
66 
            );
67
68 
        // retrieve dimensions
69 
        std::size_t n_tim = q_obs.shape(1);
70 
        std::size_t n_msk = t_msk.size() > 0 ? t_msk.shape(0) :
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
(m_cdt.size() > 0 ? m_cdt.shape(0) : 1); // initialise a mask if none provided // (corresponding to no temporal subset) auto gen_msk = [&]() { // if t_msk provided, it takes priority if (t_msk.size() > 0) return t_msk; // else if m_cdt provided, use them to generate t_msk else if (m_cdt.size() > 0) { xt::xtensor<bool, 2> c_msk = xt::zeros<bool>({n_msk, n_tim}); for (int m = 0; m < n_msk; m++) xt::view(c_msk, m) = eh::masks::generate_mask_from_conditions( m_cdt[0], xt::view(q_obs, 0), q_prd ); return c_msk; } // if neither t_msk nor m_cdt provided, generate dummy mask else return xt::xtensor<bool, 2>{xt::ones<bool>({std::size_t{1}, n_tim})}; }; auto msk = gen_msk(); // apply streamflow transformation if requested auto q_transform = [&](const xt::xtensor<double, 2>& q) { if ( transform == "none" || (transform == "pow" && exponent == 1)) { return q; } else if ( transform == "sqrt" ) { return xt::eval(xt::sqrt(q)); } else if ( transform == "inv" ) { if ( epsilon == -9 ) // determine an epsilon value to avoid zero divide epsilon = xt::mean(q_obs)() * 0.01; return xt::eval(1. / (q + epsilon)); } else if ( transform == "log" ) { if ( epsilon == -9 ) // determine an epsilon value to avoid log zero epsilon = xt::mean(q_obs)() * 0.01; return xt::eval(xt::log(q + epsilon)); } else if ( transform == "pow" ) { if ( exponent < 0 ) { if ( epsilon == -9 ) // determine an epsilon value to avoid zero divide epsilon = xt::mean(q_obs)() * 0.01; return xt::eval(xt::pow(q + epsilon, exponent)); } else { return xt::eval(xt::pow(q, exponent)); } }