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 158 : [](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 79 : zero(data_new); 29 79 : const auto oldSize = scipp::size(xold) - 1; 30 79 : const auto newSize = scipp::size(xnew) - 1; 31 79 : scipp::index iold = 0; 32 79 : scipp::index inew = 0; 33 852 : while ((iold < oldSize) && (inew < newSize)) { 34 773 : const T xo_low = xold[iold]; 35 773 : const T xo_high = xold[iold + 1]; 36 773 : const T xn_low = xnew[inew]; 37 773 : const T xn_high = xnew[inew + 1]; 38 773 : if (!less(xo_low, xn_high)) 39 5 : inew++; // old and new bins do not overlap 40 768 : 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 698 : std::abs(min(xn_high, xo_high, less) - max(xn_low, xo_low, less)); 48 698 : const auto owidth = std::abs(xo_high - xo_low); 49 698 : 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 690 : data_new[inew] += data_old[iold] * scale; 56 : } 57 698 : if (less(xo_high, xn_high)) { 58 536 : iold++; 59 : } else { 60 162 : inew++; 61 : } 62 : } 63 : } 64 : }, 65 57 : [](const units::Unit &target_edges, const units::Unit &data, 66 : const units::Unit &edges) { 67 57 : if (target_edges != edges) 68 0 : throw except::UnitError( 69 : "Input and output bin edges must have the same unit."); 70 57 : 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