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             : #include <algorithm>
       6             : #include <numeric>
       7             : 
       8             : #include "scipp/common/except.h"
       9             : #include "scipp/core/element/arithmetic.h"
      10             : #include "scipp/core/element/cumulative.h"
      11             : #include "scipp/core/subbin_sizes.h"
      12             : 
      13             : namespace scipp::core {
      14             : 
      15      101155 : SubbinSizes::SubbinSizes(const scipp::index offset, container_type &&sizes)
      16      101155 :     : m_offset(offset), m_sizes(std::move(sizes)) {
      17      101155 :   if (offset < 0)
      18           0 :     throw std::logic_error("Bad offset in class SubbinSizes.");
      19      101155 : }
      20             : 
      21       79539 : void SubbinSizes::operator=(const scipp::index value) {
      22    15827025 :   for (auto &size : m_sizes)
      23    15747486 :     size = value;
      24       79539 : }
      25             : 
      26       50510 : SubbinSizes &SubbinSizes::operator+=(const SubbinSizes &other) {
      27       50510 :   if (other.offset() < offset()) // avoid realloc if possible
      28           2 :     return *this = *this + other;
      29       50508 :   scipp::index current = other.offset() - offset();
      30       50508 :   const auto length = current + scipp::size(other.sizes());
      31       50508 :   if (length > scipp::size(sizes()))
      32       14859 :     m_sizes.resize(length);
      33     4901038 :   for (const auto &x : other.sizes())
      34     4850530 :     m_sizes[current++] += x;
      35       50508 :   return *this;
      36             : }
      37             : 
      38           5 : SubbinSizes &SubbinSizes::operator-=(const SubbinSizes &other) {
      39           5 :   return *this = *this - other;
      40             : }
      41             : 
      42       29045 : SubbinSizes SubbinSizes::cumsum_exclusive() const {
      43       29045 :   auto out = sizes();
      44       29045 :   scipp::index sum = 0;
      45    10947262 :   for (auto &x : out)
      46    10918217 :     element::exclusive_scan(sum, x);
      47       58090 :   return {offset(), std::move(out)};
      48       29045 : }
      49             : 
      50       29045 : scipp::index SubbinSizes::sum() const {
      51       29045 :   return std::accumulate(sizes().begin(), sizes().end(), scipp::index{0});
      52             : }
      53             : 
      54      125274 : SubbinSizes &SubbinSizes::add_intersection(const SubbinSizes &other) {
      55      125274 :   scipp::index delta = other.offset() - offset();
      56      125274 :   auto i = std::max(scipp::index(0), delta);
      57      125274 :   auto j = std::max(scipp::index(0), -delta);
      58    20013956 :   for (; j < scipp::size(other.sizes()) && i < scipp::size(sizes()); ++j, ++i) {
      59    19888682 :     m_sizes[i] += other.sizes()[j];
      60             :   }
      61      125274 :   return *this;
      62             : }
      63             : 
      64          43 : bool operator==(const SubbinSizes &a, const SubbinSizes &b) {
      65          43 :   return (a.offset() == b.offset()) && (a.sizes() == b.sizes());
      66             : }
      67             : 
      68             : template <class Op>
      69          17 : SubbinSizes binary(const SubbinSizes &a, const SubbinSizes &b, Op op) {
      70          17 :   const auto begin = std::min(a.offset(), b.offset());
      71          17 :   const auto end =
      72          17 :       std::max(a.offset() + a.sizes().size(), b.offset() + b.sizes().size());
      73          17 :   typename SubbinSizes::container_type sizes(end - begin);
      74          17 :   scipp::index current = a.offset() - begin;
      75          81 :   for (const auto &size : a.sizes())
      76          64 :     sizes[current++] += size;
      77          17 :   current = b.offset() - begin;
      78          44 :   for (const auto &size : b.sizes())
      79          27 :     op(sizes[current++], size);
      80          34 :   return {begin, std::move(sizes)};
      81          17 : }
      82             : 
      83           7 : SubbinSizes operator+(const SubbinSizes &a, const SubbinSizes &b) {
      84           7 :   return binary(a, b, element::add_equals);
      85             : }
      86             : 
      87          10 : SubbinSizes operator-(const SubbinSizes &a, const SubbinSizes &b) {
      88          10 :   return binary(a, b, element::subtract_equals);
      89             : }
      90             : 
      91             : /** Perform a step of an exclusive scan.
      92             :  *
      93             :  * The instance is the accumulant, the argument is the next value to be added.
      94             :  * Trims the accumulant to the offset and size of the argument.
      95             :  *
      96             :  * Note that effective cache use matters here, so we do not implemented this
      97             :  * via suboperations but a single loop.
      98             :  */
      99       50507 : void SubbinSizes::exclusive_scan(SubbinSizes &x) {
     100       50507 :   const scipp::index osize = scipp::size(sizes());
     101       50507 :   const scipp::index size = scipp::size(x.sizes());
     102       50507 :   const auto delta = x.offset() - offset();
     103       50507 :   m_sizes.reserve(size);
     104       50507 :   m_offset = x.m_offset;
     105     4901038 :   for (scipp::index i = 0; i < size; ++i) {
     106     4850531 :     const auto prev = delta + i >= osize ? 0 : m_sizes[delta + i];
     107     4850531 :     if (i >= osize)
     108        3859 :       m_sizes.push_back(prev + x.m_sizes[i]);
     109             :     else
     110     4846672 :       m_sizes[i] = prev + x.m_sizes[i];
     111     4850531 :     x.m_sizes[i] = prev;
     112             :   }
     113       50507 :   m_sizes.resize(size);
     114       50507 : }
     115             : 
     116             : } // namespace scipp::core