Line data    Source code

       1             : // SPDX-License-Identifier: BSD-3-Clause
       2             : // Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
       3             : /// @file
       4             : /// @author Simon Heybrock, Igor Gudich
       5             : #include "scipp/core/element/rebin.h"
       6             : #include "scipp/core/parallel.h"
       7             : #include "scipp/units/except.h"
       8             : #include "scipp/variable/arithmetic.h"
       9             : #include "scipp/variable/astype.h"
      10             : #include "scipp/variable/rebin.h"
      11             : #include "scipp/variable/reduction.h"
      12             : #include "scipp/variable/shape.h"
      13             : #include "scipp/variable/transform_subspan.h"
      14             : #include "scipp/variable/util.h"
      15             : 
      16             : using namespace scipp::core::element;
      17             : 
      18             : namespace scipp::variable {
      19             : 
      20             : template <typename T, class Less>
      21          40 : void rebin_non_inner(const Dim dim, const Variable &oldT,
      22             :                      // cppcheck-suppress constParameter # bug in cppcheck
      23             :                      Variable &newT, const Variable &oldCoord,
      24             :                      const Variable &newCoord) {
      25          40 :   if (oldCoord.ndim() != 1 || newCoord.ndim() != 1)
      26           0 :     throw std::invalid_argument(
      27             :         "Internal error in rebin, this should be unreachable.");
      28          40 :   const auto oldSize = oldT.dims()[dim];
      29          40 :   const auto newSize = newT.dims()[dim];
      30             : 
      31          40 :   const auto *xold = oldCoord.values<T>().data();
      32          40 :   const auto *xnew = newCoord.values<T>().data();
      33             : 
      34         392 :   auto add_from_bin = [&](auto &&slice, const auto xn_low, const auto xn_high,
      35             :                           const scipp::index iold) {
      36          88 :     auto xo_low = xold[iold];
      37          88 :     auto xo_high = xold[iold + 1];
      38             :     // delta is the overlap of the bins on the x axis
      39         176 :     const auto delta = std::abs(std::min<double>(xn_high, xo_high, Less{}) -
      40          88 :                                 std::max<double>(xn_low, xo_low, Less{}));
      41          88 :     const auto owidth = std::abs(xo_high - xo_low);
      42          88 :     slice += oldT.slice({dim, iold}) *
      43         176 :              ((delta / owidth) * (slice.unit() / oldT.unit()));
      44             :   };
      45         373 :   auto accumulate_bin = [&](auto &&slice, const auto xn_low,
      46             :                             const auto xn_high) {
      47          56 :     scipp::index begin =
      48          56 :         std::upper_bound(xold, xold + oldSize + 1, xn_low, Less{}) - xold;
      49          56 :     scipp::index end =
      50          56 :         std::upper_bound(xold, xold + oldSize + 1, xn_high, Less{}) - xold;
      51          56 :     if (begin == oldSize + 1 || end == 0)
      52           0 :       return;
      53          56 :     begin = std::max(scipp::index(0), begin - 1);
      54          56 :     add_from_bin(slice, xn_low, xn_high, begin);
      55          56 :     if (begin + 1 < end - 1)
      56          27 :       sum_into(slice, oldT.slice({dim, begin + 1, end - 1}));
      57          56 :     if (begin != end - 1 && end < oldSize + 1)
      58          32 :       add_from_bin(slice, xn_low, xn_high, end - 1);
      59             :   };
      60          80 :   auto accumulate_bins = [&](const auto &range) {
      61          96 :     for (scipp::index inew = range.begin(); inew < range.end(); ++inew) {
      62          56 :       auto xn_low = xnew[inew];
      63          56 :       auto xn_high = xnew[inew + 1];
      64          56 :       accumulate_bin(newT.slice({dim, inew}), xn_low, xn_high);
      65             :     }
      66             :   };
      67          40 :   core::parallel::parallel_for(core::parallel::blocked_range(0, newSize),
      68             :                                accumulate_bins);
      69          40 : }
      70             : 
      71             : namespace {
      72             : template <class Out, class OutEdge, class In, class InEdge>
      73             : using args = std::tuple<scipp::span<Out>, scipp::span<const OutEdge>,
      74             :                         scipp::span<const In>, scipp::span<const InEdge>>;
      75             : 
      76             : struct Greater {
      77             :   template <class A, class B>
      78         402 :   constexpr bool operator()(const A a, const B b) const noexcept {
      79         402 :     return a > b;
      80             :   }
      81             : };
      82             : 
      83             : struct Less {
      84             :   template <class A, class B>
      85        3705 :   constexpr bool operator()(const A a, const B b) const noexcept {
      86        3705 :     return a < b;
      87             :   }
      88             : };
      89             : 
      90             : } // namespace
      91             : 
      92         100 : Variable rebin(const Variable &var, const Dim dim, const Variable &oldCoord,
      93             :                const Variable &newCoord) {
      94             :   // The code branch dealing with non-stride-1 data cannot handle non-1D edges.
      95             :   // This is likely a rare case in practice so a slow transpose of input and
      96             :   // output should be sufficient for now.
      97         100 :   if (var.stride(dim) != 1 && (oldCoord.ndim() != 1 || newCoord.ndim() != 1))
      98             :     // We *copy* the transpose to ensure that memory order of dims matches input
      99             :     return copy(
     100           4 :         transpose(rebin(as_contiguous(var, dim), dim, oldCoord, newCoord),
     101           4 :                   var.dims().labels()));
     102             :   // Rebin could also implemented for count-densities. However, it may be better
     103             :   // to avoid this since it increases complexity. Instead, densities could
     104             :   // always be computed on-the-fly for visualization, if required.
     105          98 :   if (!is_edges(var.dims(), oldCoord.dims(), dim))
     106           0 :     throw except::BinEdgeError(
     107           0 :         "The input does not have coordinates with bin-edges.");
     108             : 
     109          98 :   if (is_bins(var))
     110           1 :     throw except::TypeError("The input variable cannot be binned data. Use "
     111           2 :                             "`bin` or `histogram` instead of `rebin`.");
     112             : 
     113             :   using transform_args = std::tuple<
     114             :       args<double, double, int64_t, double>,
     115             :       args<double, double, int32_t, double>,
     116             :       args<double, core::time_point, double, core::time_point>,
     117             :       args<float, core::time_point, float, core::time_point>,
     118             :       args<double, core::time_point, int64_t, core::time_point>,
     119             :       args<double, core::time_point, int32_t, core::time_point>,
     120             :       args<double, core::time_point, bool, core::time_point>,
     121             :       args<double, double, double, double>, args<float, float, float, float>,
     122             :       args<float, double, float, double>, args<float, float, float, double>,
     123             :       args<double, double, bool, double>>;
     124             : 
     125         165 :   const bool ascending = allsorted(oldCoord, dim, SortOrder::Ascending) &&
     126          68 :                          allsorted(newCoord, dim, SortOrder::Ascending);
     127         126 :   if (!ascending && !(allsorted(oldCoord, dim, SortOrder::Descending) &&
     128          29 :                       allsorted(newCoord, dim, SortOrder::Descending)))
     129           0 :     throw except::BinEdgeError(
     130           0 :         "Rebin: The old or new bin edges are not sorted.");
     131         190 :   const auto out_type = (is_int(var.dtype()) || var.dtype() == dtype<bool>)
     132          97 :                             ? dtype<double>
     133          97 :                             : var.dtype();
     134             :   // Both code branches below require stride 1 for input and output edges.
     135          97 :   const auto oldEdges = as_contiguous(oldCoord, dim);
     136          97 :   const auto newEdges = as_contiguous(newCoord, dim);
     137          97 :   Variable rebinned;
     138          97 :   if (var.stride(dim) == 1) {
     139          57 :     if (ascending) {
     140          84 :       rebinned = transform_subspan<transform_args>(
     141          42 :           out_type, dim, newEdges.dims()[dim] - 1, newEdges, var, oldEdges,
     142          42 :           core::element::rebin<Less>, "rebin");
     143             :     } else {
     144          30 :       rebinned = transform_subspan<transform_args>(
     145          15 :           out_type, dim, newEdges.dims()[dim] - 1, newEdges, var, oldEdges,
     146          15 :           core::element::rebin<Greater>, "rebin");
     147             :     }
     148             :   } else {
     149          40 :     auto dims = var.dims();
     150          40 :     dims.resize(dim, newEdges.dims()[dim] - 1);
     151          40 :     rebinned = Variable(astype(Variable(var, Dimensions{}), out_type), dims);
     152          40 :     if (newEdges.dims().ndim() > 1)
     153           0 :       throw std::runtime_error(
     154           0 :           "Not inner rebin works only for 1d coordinates for now.");
     155          40 :     if (oldEdges.dtype() == dtype<double>) {
     156          40 :       if (ascending)
     157          26 :         rebin_non_inner<double, Less>(dim, var, rebinned, oldEdges, newEdges);
     158             :       else
     159          14 :         rebin_non_inner<double, Greater>(dim, var, rebinned, oldEdges,
     160             :                                          newEdges);
     161           0 :     } else if (oldEdges.dtype() == dtype<float>) {
     162           0 :       if (ascending)
     163           0 :         rebin_non_inner<float, Less>(dim, var, rebinned, oldEdges, newEdges);
     164             :       else
     165           0 :         rebin_non_inner<float, Greater>(dim, var, rebinned, oldEdges, newEdges);
     166             :     } else {
     167           0 :       throw except::TypeError("Rebinning is possible only for coords of types "
     168           0 :                               "`float64` or `float32`.");
     169             :     }
     170          40 :   }
     171             :   // If rebinned dimension has stride 1 but is not an inner dimension then we
     172             :   // need to transpose the output of transform_subspan to retain the input
     173             :   // dimension order.
     174          97 :   return transpose(rebinned, var.dims().labels());
     175          97 : }
     176             : 
     177             : } // namespace scipp::variable