intervals.hpp

// Copyright (c) 2023, INRAE.
// Distributed under the terms of the GPL-3 Licence.
// The full licence is in the file LICENCE, distributed with this software.

#ifndef EVALHYD_PROBABILIST_INTERVALS_HPP
#define EVALHYD_PROBABILIST_INTERVALS_HPP

#include <limits>

#include <xtensor/xtensor.hpp>
#include <xtensor/xview.hpp>
#include <xtensor/xindex_view.hpp>
#include <xtensor/xsort.hpp>


namespace evalhyd
{
    namespace probabilist
    {
        namespace elements
        {
            /// Compute the bounds of the predictive intervals by computing
            /// the quantiles of the predictive distribution corresponding
            /// to the confidence intervals.
            ///
            /// \param q_prd
            ///     Streamflow predictions.
            ///     shape: (sites, lead times, members, time)
            /// \param c_lvl
            ///     Confidence levels for the predictive intervals.
            ///     shape: (intervals,)
            /// \param n_sit
            ///     Number of sites.
            /// \param n_ldt
            ///     Number of lead times.
            /// \param n_itv
            ///     Number of predictive intervals.
            /// \param n_tim
            ///     Number of time steps.
            /// \return
            ///     Bounds of the predictive intervals.
            ///     shape: (sites, lead times, intervals, bounds, time)
            template <class XD4>
            inline xt::xtensor<double, 5> calc_itv_bnds(
                    const XD4& q_prd,
                    const xt::xtensor<double, 1>& c_lvl,
                    const xt::xtensor<double, 1>& q_lvl,
                    std::size_t n_sit,
                    std::size_t n_ldt,
                    std::size_t n_itv,
                    std::size_t n_tim
            )
            {
                xt::xtensor<double, 5> itv_bnds =
                        xt::zeros<double>({n_sit, n_ldt, n_itv, std::size_t {2}, n_tim});

                // determine quantiles forming the predictive intervals
                // from the confidence levels
                xt::xtensor<double, 2> quantiles =
                        xt::zeros<double>({n_itv, std::size_t {2}});
                xt::col(quantiles, 0) = 0.5 - c_lvl / 200.;
                xt::col(quantiles, 1) = 0.5 + c_lvl / 200.;

                // get or compute predictive interval bounds from quantiles
                if (q_lvl.size() > 0)
                {   
                    for (std::size_t i = 0; i < n_itv; i++)
                    {
                    auto a = xt::broadcast(xt::view(quantiles, i), std::vector<std::size_t>({q_lvl.size(), 2}));
                    auto b = xt::broadcast(q_lvl, std::vector<std::size_t>({2, q_lvl.size()}));