fix broadcasting problems with multi-sites/multi-leadtimes

include/evalhyd/detail/probabilist/ranks.hpp

@@ -107,7 +107,8 @@ namespace evalhyd
                 return ranks;
             }
 
-            /// Compute the number of observations for all possible rank values.
+            /// Compute the number of observations in each interval of the
+            /// rank diagram.
             ///
             /// \param r_k
             ///     Ranks of streamflow observations.
@@ -129,7 +130,7 @@ namespace evalhyd
             /// \param n_exp
             ///     Number of bootstrap samples.
             /// \return
-            ///     Tallies of streamflow observations for all possible ranks.
+            ///     Tallies of streamflow observations in each rank interval.
             ///     shape: (sites, lead times, subsets, samples, ranks)
             inline xt::xtensor<double, 5> calc_o_j(
                     const xt::xtensor<double, 3>& r_k,
@@ -235,7 +236,9 @@ namespace evalhyd
                                          m, b_exp[e]);
 
                         // calculate length of subset
-                        auto l = xt::sum(t_msk_sampled, -1);
+                        auto l = xt::eval(
+                                xt::sum(t_msk_sampled, -1, xt::keep_dims)
+                        );
 
                         // compute the rank diagram
                         xt::view(REL_DIAG, xt::all(), xt::all(), m, e, xt::all()) =
@@ -300,20 +303,29 @@ namespace evalhyd
                                          m, b_exp[e]);
 
                         // calculate length of subset
-                        auto l = xt::sum(t_msk_sampled, -1);
+                        auto l = xt::eval(
+                                xt::sum(t_msk_sampled, -1, xt::keep_dims)
+                        );
 
                         // compute the Delta score
                         // \Delta = \sum_{k=1}^{N+1} (r_k - \frac{M}{N+1})^2
+                        auto delta =  xt::nansum(
+                                xt::square(
+                                        xt::view(o_j, xt::all(), xt::all(), m, e, xt::all())
+                                        - (l / (n_mbr + 1))
+                                ),
+                                -1
+                        );
+
                         // \Delta_o = \frac{MN}{N+1}
+                        auto delta_o = (
+                                xt::view(l, xt::all(), xt::all(), 0)
+                                * n_mbr / (n_mbr + 1)
+                        );
+
                         // \delta = $\frac{\Delta}{\Delta_o}
                         xt::view(DS, xt::all(), xt::all(), m, e) =
-                                xt::nansum(
-                                        xt::square(
-                                                xt::view(o_j, xt::all(), xt::all(), m, e, xt::all())
-                                                - (l / (n_mbr + 1))
-                                        ),
-                                        -1
-                                ) / (l * n_mbr / (n_mbr + 1));
+                               delta / delta_o;
                     }
                 }