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 <algorithm>
       8             : #include <numeric>
       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             : 
      17             : namespace scipp::core::element {
      18             : 
      19             : namespace {
      20     8059817 : constexpr auto iadd = [](const auto &x1, const scipp::index i1, const auto &x2,
      21             :                          const scipp::index i2) {
      22             :   using V = std::decay_t<decltype(x1)>;
      23             :   if constexpr (is_ValueAndVariance_v<V>) {
      24        9888 :     x1.value[i1] += x2.value[i2];
      25        9888 :     x1.variance[i1] += x2.variance[i2];
      26             :   } else {
      27     8049929 :     x1[i1] += x2[i2];
      28             :   }
      29     8059817 : };
      30             : } // namespace
      31             : 
      32             : namespace histogram_detail {
      33             : template <class Coord, class Weight>
      34             : using args = std::tuple<scipp::span<Weight>, scipp::span<const Coord>,
      35             :                         scipp::span<const Weight>, scipp::span<const Coord>>;
      36             : }
      37             : 
      38             : static constexpr auto histogram = overloaded{
      39             :     element::arg_list<histogram_detail::args<double, double>,
      40             :                       histogram_detail::args<double, float>,
      41             :                       histogram_detail::args<double, int64_t>,
      42             :                       histogram_detail::args<double, int32_t>,
      43             :                       histogram_detail::args<float, double>,
      44             :                       histogram_detail::args<float, float>,
      45             :                       histogram_detail::args<float, int64_t>,
      46             :                       histogram_detail::args<float, int32_t>,
      47             :                       histogram_detail::args<int32_t, double>,
      48             :                       histogram_detail::args<int32_t, float>,
      49             :                       histogram_detail::args<int32_t, int64_t>,
      50             :                       histogram_detail::args<int32_t, int32_t>,
      51             :                       histogram_detail::args<int64_t, double>,
      52             :                       histogram_detail::args<int64_t, float>,
      53             :                       histogram_detail::args<int64_t, int64_t>,
      54             :                       histogram_detail::args<int64_t, int32_t>,
      55             :                       histogram_detail::args<time_point, double>,
      56             :                       histogram_detail::args<time_point, float>,
      57             :                       histogram_detail::args<time_point, int64_t>,
      58             :                       histogram_detail::args<time_point, int32_t>>,
      59       54391 :     [](const auto &data, const auto &events, const auto &weights,
      60             :        const auto &edges) {
      61       54391 :       zero(data);
      62             :       // Special implementation for linear bins. Gives a 1x to 20x speedup
      63             :       // for few and many events per histogram, respectively.
      64       54391 :       if (scipp::numeric::islinspace(edges)) {
      65       53071 :         const auto params = core::linear_edge_params(edges);
      66     8120787 :         for (scipp::index i = 0; i < scipp::size(events); ++i) {
      67     8067716 :           const auto x = events[i];
      68     8067716 :           if (const auto bin = get_bin<scipp::index>(x, edges, params);
      69             :               bin >= 0)
      70     8051773 :             iadd(data, bin, weights, i);
      71             :         }
      72             :       } else {
      73        1320 :         core::expect::histogram::sorted_edges(edges);
      74        9373 :         for (scipp::index i = 0; i < scipp::size(events); ++i) {
      75        8054 :           const auto x = events[i];
      76        8054 :           auto it = std::upper_bound(edges.begin(), edges.end(), x);
      77        8054 :           if (it != edges.end() && it != edges.begin())
      78        8044 :             iadd(data, --it - edges.begin(), weights, i);
      79             :         }
      80             :       }
      81       54390 :     },
      82       47552 :     [](const units::Unit &events_unit, const units::Unit &weights_unit,
      83             :        const units::Unit &edge_unit) {
      84       47552 :       if (events_unit != edge_unit)
      85           0 :         throw except::UnitError(
      86           0 :             "Bin edges must have same unit as the input coordinate.");
      87       47552 :       return weights_unit;
      88             :     },
      89             :     transform_flags::expect_in_variance_if_out_variance,
      90             :     transform_flags::expect_no_variance_arg<1>,
      91             :     transform_flags::expect_no_variance_arg<3>};
      92             : 
      93             : } // namespace scipp::core::element