fix bug for intervals climatology

include/evalhyd/detail/probabilist/evaluator.hpp

@@ -741,7 +741,7 @@ namespace evalhyd
                 {
                     WSS = metrics::calc_WSS(
                             q_obs, get_c_lvl(), get_clim_bnds(), get_WS(),
-                            n_sit, n_ldt, n_itv, n_msk, n_exp
+                            t_msk, b_exp, n_sit, n_ldt, n_itv, n_msk, n_exp
                     );
                 }
                 return WSS.value();

include/evalhyd/detail/probabilist/intervals.hpp

@@ -7,6 +7,7 @@
 
 #include <xtensor/xtensor.hpp>
 #include <xtensor/xview.hpp>
+#include <xtensor/xindex_view.hpp>
 #include <xtensor/xsort.hpp>
 
 #include "../maths.hpp"
@@ -220,21 +221,21 @@ namespace evalhyd
                             {
                                 for (std::size_t i = 0; i < n_itv; i++)
                                 {
+                                    auto obs = xt::view(q_obs_masked_sampled,
+                                                        s, l, xt::all());
+                                    auto obs_filtered =
+                                            xt::filter(obs, !xt::isnan(obs));
+
                                     // lower bound
                                     xt::view(clim_bnds, s, l, m, e, i, 0) =
                                             evalhyd::maths::quantile(
-                                                    xt::view(q_obs_masked_sampled,
-                                                             s, l, xt::all()),
-                                                    quantiles(i, 0)
+                                                    obs_filtered, quantiles(i, 0)
                                             );
                                     // upper bound
                                     xt::view(clim_bnds, s, l, m, e, i, 1) =
                                             evalhyd::maths::quantile(
-                                                    xt::view(q_obs_masked_sampled,
-                                                             s, l, xt::all()),
-                                                    quantiles(i, 1)
+                                                    obs_filtered, quantiles(i, 1)
                                             );
-
                                 }
                             }
                         }
@@ -571,6 +572,12 @@ namespace evalhyd
             /// \param WS
             ///     Winkler scores.
             ///     shape: (sites, lead times, subsets, samples, intervals)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (sites, lead times, subsets, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
             /// \param n_sit
             ///     Number of sites.
             /// \param n_ldt
@@ -590,6 +597,8 @@ namespace evalhyd
                     const xt::xtensor<double, 1>& c_lvl,
                     const xt::xtensor<double, 6>& clim_bnds,
                     const xt::xtensor<double, 5>& WS,
+                    const xt::xtensor<bool, 4>& t_msk,
+                    const std::vector<xt::xkeep_slice<int>>& b_exp,
                     std::size_t n_sit,
                     std::size_t n_ldt,
                     std::size_t n_itv,
@@ -601,21 +610,38 @@ namespace evalhyd
                 xt::xtensor<double, 5> WS_clim =
                         xt::zeros<double>({n_sit, n_ldt, n_msk, n_exp, n_itv});
 
-                for (std::size_t m = 0; m < n_msk; m++)
+                for (std::size_t l = 0; l < n_ldt; l++)
                 {
-                    for (std::size_t e = 0; e < n_exp; e++)
+                    for (std::size_t m = 0; m < n_msk; m++)
                     {
-                        auto ws_clim = intermediate::calc_ws(
-                                q_obs, c_lvl,
-                                xt::view(clim_bnds, xt::all(), xt::all(), m, e,
-                                         xt::all(), xt::all(), xt::newaxis())
-                        );
+                        for (std::size_t e = 0; e < n_exp; e++)
+                        {
+                            auto ws_clim = intermediate::calc_ws(
+                                    q_obs, c_lvl,
+                                    xt::view(clim_bnds, xt::all(), l, xt::newaxis(),
+                                             m, e, xt::all(), xt::all(),
+                                             xt::newaxis())
+                            );
+
+                            // apply the mask
+                            // (using NaN workaround until reducers work on masked_view)
+                            auto ws_clim_masked = xt::where(
+                                    xt::view(t_msk, xt::all(), l, m, xt::newaxis(), xt::all()),
+                                    xt::view(ws_clim, xt::all(), l, xt::all(), xt::all()),
+                                    NAN
+                            );
+
+                            // apply the bootstrap sampling
+                            auto ws_clim_masked_sampled =
+                                    xt::view(ws_clim_masked, xt::all(), xt::all(), b_exp[e]);
 
-                        xt::view(WS_clim, xt::all(), xt::all(), m, e, xt::all()) =
-                                xt::nanmean(ws_clim, -1);
+                            xt::view(WS_clim, xt::all(), l, m, e, xt::all()) =
+                                    xt::nanmean(ws_clim_masked_sampled, -1);
+                        }
                     }
                 }
 
+                // compute the Winkler skill score
                 return 1 - (WS / WS_clim);
             }
         }