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
       5             : #pragma once
       6             : 
       7             : #include <cmath>
       8             : #include <type_traits>
       9             : 
      10             : #include "scipp/common/numeric.h"
      11             : #include "scipp/common/overloaded.h"
      12             : #include "scipp/core/element/arg_list.h"
      13             : #include "scipp/core/element/util.h"
      14             : #include "scipp/core/histogram.h"
      15             : #include "scipp/core/transform_common.h"
      16             : #include "scipp/core/value_and_variance.h"
      17             : 
      18             : namespace scipp::core::element {
      19             : 
      20             : template <class Less>
      21             : static constexpr auto rebin = overloaded{
      22         150 :     [](const auto &data_new, const auto &xnew, const auto &data_old,
      23             :        const auto &xold) {
      24             :       using T = decltype(xold[0] + (xold[0] - xnew[0]));
      25             :       // Note: using const rather than constexpr here is required to
      26             :       // avoid an internal compiler error on Windows/MSVC
      27             :       const Less less;
      28          75 :       zero(data_new);
      29          75 :       const auto oldSize = scipp::size(xold) - 1;
      30          75 :       const auto newSize = scipp::size(xnew) - 1;
      31          75 :       scipp::index iold = 0;
      32          75 :       scipp::index inew = 0;
      33         416 :       while ((iold < oldSize) && (inew < newSize)) {
      34         341 :         const T xo_low = xold[iold];
      35         341 :         const T xo_high = xold[iold + 1];
      36         341 :         const T xn_low = xnew[inew];
      37         341 :         const T xn_high = xnew[inew + 1];
      38         341 :         if (!less(xo_low, xn_high))
      39           5 :           inew++; // old and new bins do not overlap
      40         336 :         else if (!less(xn_low, xo_high))
      41          70 :           iold++; // old and new bins do not overlap
      42             :         else {
      43             :           // delta is the overlap of the bins on the x axis
      44             :           using std::min;
      45             :           using std::max;
      46             :           const auto delta =
      47         266 :               std::abs(min(xn_high, xo_high, less) - max(xn_low, xo_low, less));
      48         266 :           const auto owidth = std::abs(xo_high - xo_low);
      49         266 :           const auto scale = delta / owidth;
      50             :           if constexpr (is_ValueAndVariance_v<
      51             :                             std::decay_t<decltype(data_old)>>) {
      52           8 :             data_new.value[inew] += data_old.value[iold] * scale;
      53           8 :             data_new.variance[inew] += data_old.variance[iold] * scale;
      54             :           } else {
      55         258 :             data_new[inew] += data_old[iold] * scale;
      56             :           }
      57         266 :           if (less(xo_high, xn_high)) {
      58         136 :             iold++;
      59             :           } else {
      60         130 :             inew++;
      61             :           }
      62             :         }
      63             :       }
      64             :     },
      65          53 :     [](const units::Unit &target_edges, const units::Unit &data,
      66             :        const units::Unit &edges) {
      67          53 :       if (target_edges != edges)
      68           0 :         throw except::UnitError(
      69             :             "Input and output bin edges must have the same unit.");
      70          53 :       return data;
      71             :     },
      72             :     transform_flags::expect_in_variance_if_out_variance,
      73             :     transform_flags::expect_no_variance_arg<1>,
      74             :     transform_flags::expect_no_variance_arg<3>};
      75             : 
      76             : } // namespace scipp::core::element