From 645511027846c82fed21e56938f25a00c97c067e Mon Sep 17 00:00:00 2001
From: fbourgin <francois.bourgin@inrae.fr>
Date: Tue, 8 Apr 2025 13:13:10 +0200
Subject: [PATCH] add a new variable q_lvl to get the level of quantiles

---
 .../evalhyd/detail/probabilist/evaluator.hpp  |  8 ++++--
 .../evalhyd/detail/probabilist/intervals.hpp  | 28 +++++++++++++++++--
 .../evalhyd/detail/probabilist/quantiles.hpp  | 15 ++++++++--
 include/evalhyd/evalp.hpp                     |  4 ++-
 4 files changed, 46 insertions(+), 9 deletions(-)

diff --git a/include/evalhyd/detail/probabilist/evaluator.hpp b/include/evalhyd/detail/probabilist/evaluator.hpp
index 3f94fd6..619104e 100644
--- a/include/evalhyd/detail/probabilist/evaluator.hpp
+++ b/include/evalhyd/detail/probabilist/evaluator.hpp
@@ -36,6 +36,7 @@ namespace evalhyd
             // members for optional input data
             const XD2& _q_thr;
             const xt::xtensor<double, 1>& _c_lvl;
+            const xt::xtensor<double, 1>& _q_lvl;
             xtl::xoptional<const std::string, bool> _events;
             xt::xtensor<bool, 4> t_msk;
             const std::vector<xt::xkeep_slice<int>>& b_exp;
@@ -336,7 +337,7 @@ namespace evalhyd
                 if (!itv_bnds.has_value())
                 {
                     itv_bnds = elements::calc_itv_bnds(
-                            q_prd, get_c_lvl(),
+                            q_prd, get_c_lvl(), _q_lvl,
                             n_sit, n_ldt, n_itv, n_tim
                     );
                 }
@@ -393,7 +394,7 @@ namespace evalhyd
                 if (!qs.has_value())
                 {
                     qs = intermediate::calc_qs(
-                            q_obs, get_q_qnt(), n_mbr
+                            q_obs, get_q_qnt(), n_mbr, _q_lvl
                     );
                 }
                 return qs.value();
@@ -482,12 +483,13 @@ namespace evalhyd
                       const XD4& prd,
                       const XD2& thr,
                       const xt::xtensor<double, 1>& lvl,
+                      const xt::xtensor<double, 1>& qlvl,
                       xtl::xoptional<const std::string&, bool> events,
                       const XB4& msk,
                       const std::vector<xt::xkeep_slice<int>>& exp,
                       const long int seed) :
                     q_obs{obs}, q_prd{prd},
-                    _q_thr{thr}, _c_lvl{lvl}, _events{events},
+                    _q_thr{thr}, _c_lvl{lvl}, _q_lvl{qlvl}, _events{events},
                     t_msk(msk), b_exp(exp),
                     random_seed{seed}
             {
diff --git a/include/evalhyd/detail/probabilist/intervals.hpp b/include/evalhyd/detail/probabilist/intervals.hpp
index adbd487..21316f2 100644
--- a/include/evalhyd/detail/probabilist/intervals.hpp
+++ b/include/evalhyd/detail/probabilist/intervals.hpp
@@ -44,6 +44,7 @@ namespace evalhyd
             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,
@@ -60,15 +61,38 @@ namespace evalhyd
                 xt::col(quantiles, 0) = 0.5 - c_lvl / 200.;
                 xt::col(quantiles, 1) = 0.5 + c_lvl / 200.;
 
-                // compute predictive interval bounds from quantiles
-                for (std::size_t i = 0; i < n_itv; i++)
+                // get or compute predictive interval bounds from quantiles
+                if (q_lvl.size() > 0)
+                {   
+                    for (std::size_t i = 0; i < n_itv; i++)
+                    {
+                    auto res = xt::where(xt::equal(xt::view(quantiles, i), q_lvl));
+                    if (res.size() != 2) 
+                    {
+                        throw std::runtime_error(
+                                "interval-based metric requested, "
+                                "but *c_lvl* not matching *q_lvl"
+                        );
+                    } else 
+                    {
+                        xt::view(itv_bnds, xt::all(), xt::all(), i, 0, xt::all()) =
+                                xt::view(q_prd, res[0]);
+                        xt::view(itv_bnds, xt::all(), xt::all(), i, 1, xt::all()) = 
+                                xt::view(q_prd, res[1]);
+                    }
+                    }
+                }
+                else
                 {
+                    for (std::size_t i = 0; i < n_itv; i++)
+                    {
                     auto q =  xt::quantile(q_prd, xt::view(quantiles, i), 2);
 
                     xt::view(itv_bnds, xt::all(), xt::all(), i, 0, xt::all()) =
                             xt::view(q, 0);
                     xt::view(itv_bnds, xt::all(), xt::all(), i, 1, xt::all()) =
                             xt::view(q, 1);
+                    }
                 }
 
                 return itv_bnds;
diff --git a/include/evalhyd/detail/probabilist/quantiles.hpp b/include/evalhyd/detail/probabilist/quantiles.hpp
index 88bd528..bb13751 100644
--- a/include/evalhyd/detail/probabilist/quantiles.hpp
+++ b/include/evalhyd/detail/probabilist/quantiles.hpp
@@ -60,13 +60,22 @@ namespace evalhyd
             inline xt::xtensor<double, 4> calc_qs(
                     const XD2 &q_obs,
                     const xt::xtensor<double, 4>& q_qnt,
-                    std::size_t n_mbr
+                    std::size_t n_mbr,
+                    const xt::xtensor<double, 1>& q_lvl
             )
             {
-                // compute the quantile order $alpha$
-                xt::xtensor<double, 1> alpha =
+                xt::xtensor<double, 1> alpha;
+                // get or compute the quantile order $alpha$
+                if (q_lvl.size() > 0)
+                {
+                    alpha = xt::sort(q_lvl);
+                }
+                else
+                {
+                    alpha =
                         xt::arange<double>(1., double(n_mbr + 1))
                         / double(n_mbr + 1);
+                }
 
                 // calculate the difference
                 xt::xtensor<double, 4> diff =
diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp
index f7f4340..0039677 100644
--- a/include/evalhyd/evalp.hpp
+++ b/include/evalhyd/evalp.hpp
@@ -181,6 +181,7 @@ namespace evalhyd
             xtl::xoptional<const std::string, bool> events =
                     xtl::missing<const std::string>(),
             const std::vector<double>& c_lvl = {},
+            const std::vector<double>& q_lvl = {},
             const xt::xexpression<XB4>& t_msk = XB4({}),
             const xt::xexpression<XS2>& m_cdt = XS2({}),
             xtl::xoptional<const std::unordered_map<std::string, int>, bool> bootstrap =
@@ -233,6 +234,7 @@ namespace evalhyd
 
         // adapt vector to tensor
         const xt::xtensor<double, 1> c_lvl_ = xt::adapt(c_lvl);
+        const xt::xtensor<double, 1> q_lvl_ = xt::adapt(q_lvl);
 
         // check that the metrics/diagnostics to be computed are valid
         utils::check_metrics(
@@ -418,7 +420,7 @@ namespace evalhyd
 
         // instantiate determinist evaluator
         probabilist::Evaluator<XD2, XD4, XB4> evaluator(
-                q_obs_, q_prd_, q_thr_, c_lvl_, events,
+                q_obs_, q_prd_, q_thr_, c_lvl_, q_lvl_, events,
                 t_msk_.size() > 0 ? t_msk_: (m_cdt_.size() > 0 ? c_msk : t_msk_),
                 b_exp,
                 random_seed
-- 
GitLab