boost/histogram/indexed.hpp

   1 // Copyright 2015-2016 Hans Dembinski
   2 //
   3 // Distributed under the Boost Software License, Version 1.0.
   4 // (See accompanying file LICENSE_1_0.txt
   5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_HISTOGRAM_INDEXED_HPP
   8 #define BOOST_HISTOGRAM_INDEXED_HPP
   9
  10 #include <array>
  11 #include <boost/config.hpp>
  12 #include <boost/histogram/axis/traits.hpp>
  13 #include <boost/histogram/detail/axes.hpp>
  14 #include <boost/histogram/detail/iterator_adaptor.hpp>
  15 #include <boost/histogram/detail/operators.hpp>
  16 #include <boost/histogram/fwd.hpp>
  17 #include <iterator>
  18 #include <type_traits>
  19 #include <utility>
  20
  21 namespace boost {
  22 namespace histogram {
  23
  24 /** Coverage mode of the indexed range generator.
  25
  26   Defines options for the iteration strategy.
  27 */
  28 enum class coverage {
  29   inner, /*!< iterate over inner bins, exclude underflow and overflow */
  30   all,   /*!< iterate over all bins, including underflow and overflow */
  31 };
  32
  33 /** Input iterator range over histogram bins with multi-dimensional index.
  34
  35   The iterator returned by begin() can only be incremented. begin() may only be called
  36   once, calling it a second time returns the end() iterator. If several copies of the
  37   input iterators exist, the other copies become invalid if one of them is incremented.
  38 */
  39 template <class Histogram>
  40 class BOOST_ATTRIBUTE_NODISCARD indexed_range {
  41 private:
  42   using histogram_type = Histogram;
  43   static constexpr std::size_t buffer_size =
  44       detail::buffer_size<typename std::remove_const_t<histogram_type>::axes_type>::value;
  45
  46 public:
  47   using value_iterator = std::conditional_t<std::is_const<histogram_type>::value,
  48                                             typename histogram_type::const_iterator,
  49                                             typename histogram_type::iterator>;
  50   using value_reference = typename std::iterator_traits<value_iterator>::reference;
  51   using value_type = typename std::iterator_traits<value_iterator>::value_type;
  52
  53   class iterator;
  54   using range_iterator = iterator; ///< deprecated
  55
  56   /** Lightweight view to access value and index of current cell.
  57
  58     The methods provide access to the current cell indices and bins. It acts like a
  59     pointer to the cell value, and in a limited way also like a reference. To interoperate
  60     with the algorithms of the standard library, the accessor is implicitly convertible to
  61     a cell value. Assignments and comparisons are passed through to the cell. An accessor
  62     is coupled to its parent indexed_range::iterator. Moving the parent iterator
  63     forward also updates the linked accessor. Accessors are not copyable. They cannot be
  64     stored in containers, but indexed_range::iterator can be stored.
  65   */
  66   class BOOST_ATTRIBUTE_NODISCARD accessor : detail::mirrored<accessor, void> {
  67   public:
  68     /// Array-like view into the current multi-dimensional index.
  69     class index_view {
  70       using index_pointer = const typename iterator::index_data*;
  71
  72     public:
  73       using const_reference = const axis::index_type&;
  74       using reference = const_reference; ///< deprecated
  75
  76       /// implementation detail
  77       class const_iterator
  78           : public detail::iterator_adaptor<const_iterator, index_pointer,
  79                                             const_reference> {
  80       public:
  81         const_reference operator*() const noexcept { return const_iterator::base()->idx; }
  82
  83       private:
  84         explicit const_iterator(index_pointer i) noexcept
  85             : const_iterator::iterator_adaptor_(i) {}
  86
  87         friend class index_view;
  88       };
  89
  90       const_iterator begin() const noexcept { return const_iterator{begin_}; }
  91       const_iterator end() const noexcept { return const_iterator{end_}; }
  92       std::size_t size() const noexcept {
  93         return static_cast<std::size_t>(end_ - begin_);
  94       }
  95       const_reference operator[](unsigned d) const noexcept { return begin_[d].idx; }
  96       const_reference at(unsigned d) const { return begin_[d].idx; }
  97
  98     private:
  99       /// implementation detail
 100       index_view(index_pointer b, index_pointer e) : begin_(b), end_(e) {}
 101
 102       index_pointer begin_, end_;
 103       friend class accessor;
 104     };
 105
 106     // assignment is pass-through
 107     accessor& operator=(const accessor& o) {
 108       get() = o.get();
 109       return *this;
 110     }
 111
 112     // assignment is pass-through
 113     template <class T>
 114     accessor& operator=(const T& x) {
 115       get() = x;
 116       return *this;
 117     }
 118
 119     /// Returns the cell reference.
 120     value_reference get() const noexcept { return *(iter_.iter_); }
 121     /// @copydoc get()
 122     value_reference operator*() const noexcept { return get(); }
 123     /// Access fields and methods of the cell object.
 124     value_iterator operator->() const noexcept { return iter_.iter_; }
 125
 126     /// Access current index.
 127     /// @param d axis dimension.
 128     axis::index_type index(unsigned d = 0) const noexcept {
 129       return iter_.indices_[d].idx;
 130     }
 131
 132     /// Access indices as an iterable range.
 133     index_view indices() const noexcept {
 134       BOOST_ASSERT(iter_.indices_.hist_);
 135       return {iter_.indices_.data(),
 136               iter_.indices_.data() + iter_.indices_.hist_->rank()};
 137     }
 138
 139     /// Access current bin.
 140     /// @tparam N axis dimension.
 141     template <unsigned N = 0>
 142     decltype(auto) bin(std::integral_constant<unsigned, N> = {}) const {
 143       BOOST_ASSERT(iter_.indices_.hist_);
 144       return iter_.indices_.hist_->axis(std::integral_constant<unsigned, N>())
 145           .bin(index(N));
 146     }
 147
 148     /// Access current bin.
 149     /// @param d axis dimension.
 150     decltype(auto) bin(unsigned d) const {
 151       BOOST_ASSERT(iter_.indices_.hist_);
 152       return iter_.indices_.hist_->axis(d).bin(index(d));
 153     }
 154
 155     /** Computes density in current cell.
 156
 157       The density is computed as the cell value divided by the product of bin widths. Axes
 158       without bin widths, like axis::category, are treated as having unit bin with.
 159     */
 160     double density() const {
 161       BOOST_ASSERT(iter_.indices_.hist_);
 162       double x = 1;
 163       unsigned d = 0;
 164       iter_.indices_.hist_->for_each_axis([&](const auto& a) {
 165         const auto w = axis::traits::width_as<double>(a, this->index(d++));
 166         x *= w ? w : 1;
 167       });
 168       return get() / x;
 169     }
 170
 171     // forward all comparison operators to the value
 172     bool operator<(const accessor& o) noexcept { return get() < o.get(); }
 173     bool operator>(const accessor& o) noexcept { return get() > o.get(); }
 174     bool operator==(const accessor& o) noexcept { return get() == o.get(); }
 175     bool operator!=(const accessor& o) noexcept { return get() != o.get(); }
 176     bool operator<=(const accessor& o) noexcept { return get() <= o.get(); }
 177     bool operator>=(const accessor& o) noexcept { return get() >= o.get(); }
 178
 179     template <class U>
 180     bool operator<(const U& o) const noexcept {
 181       return get() < o;
 182     }
 183
 184     template <class U>
 185     bool operator>(const U& o) const noexcept {
 186       return get() > o;
 187     }
 188
 189     template <class U>
 190     bool operator==(const U& o) const noexcept {
 191       return get() == o;
 192     }
 193
 194     template <class U>
 195     bool operator!=(const U& o) const noexcept {
 196       return get() != o;
 197     }
 198
 199     template <class U>
 200     bool operator<=(const U& o) const noexcept {
 201       return get() <= o;
 202     }
 203
 204     template <class U>
 205     bool operator>=(const U& o) const noexcept {
 206       return get() >= o;
 207     }
 208
 209     operator value_type() const noexcept { return get(); }
 210
 211   private:
 212     accessor(iterator& i) noexcept : iter_(i) {}
 213
 214     accessor(const accessor&) = default; // only callable by indexed_range::iterator
 215
 216     iterator& iter_;
 217
 218     friend class iterator;
 219   };
 220
 221   /// implementation detail
 222   class iterator {
 223   public:
 224     using value_type = typename indexed_range::value_type;
 225     using reference = accessor;
 226
 227   private:
 228     struct pointer_proxy {
 229       reference* operator->() noexcept { return std::addressof(ref_); }
 230       reference ref_;
 231     };
 232
 233   public:
 234     using pointer = pointer_proxy;
 235     using difference_type = std::ptrdiff_t;
 236     using iterator_category = std::forward_iterator_tag;
 237
 238     reference operator*() noexcept { return *this; }
 239     pointer operator->() noexcept { return pointer_proxy{operator*()}; }
 240
 241     iterator& operator++() {
 242       BOOST_ASSERT(indices_.hist_);
 243       std::size_t stride = 1;
 244       auto c = indices_.begin();
 245       ++c->idx;
 246       ++iter_;
 247       while (c->idx == c->end && (c != (indices_.end() - 1))) {
 248         c->idx = c->begin;
 249         iter_ -= (c->end - c->begin) * stride;
 250         stride *= c->extent;
 251         ++c;
 252         ++c->idx;
 253         iter_ += stride;
 254       }
 255       return *this;
 256     }
 257
 258     iterator operator++(int) {
 259       auto prev = *this;
 260       operator++();
 261       return prev;
 262     }
 263
 264     bool operator==(const iterator& x) const noexcept { return iter_ == x.iter_; }
 265     bool operator!=(const iterator& x) const noexcept { return !operator==(x); }
 266
 267   private:
 268     iterator(histogram_type* h) : iter_(h->begin()), indices_(h) {}
 269
 270     value_iterator iter_;
 271
 272     struct index_data {
 273       axis::index_type idx, begin, end, extent;
 274     };
 275
 276     struct indices_t : std::array<index_data, buffer_size> {
 277       using base_type = std::array<index_data, buffer_size>;
 278
 279       indices_t(histogram_type* h) noexcept : hist_{h} {}
 280
 281       constexpr auto end() noexcept { return base_type::begin() + hist_->rank(); }
 282       constexpr auto end() const noexcept { return base_type::begin() + hist_->rank(); }
 283       histogram_type* hist_;
 284     } indices_;
 285
 286     friend class indexed_range;
 287   };
 288
 289   indexed_range(histogram_type& hist, coverage cov) : begin_(&hist), end_(&hist) {
 290     std::size_t stride = 1;
 291     auto ca = begin_.indices_.begin();
 292     const auto clast = ca + begin_.indices_.hist_->rank() - 1;
 293     begin_.indices_.hist_->for_each_axis(
 294         [ca, clast, cov, &stride, this](const auto& a) mutable {
 295           using opt = axis::traits::static_options<std::decay_t<decltype(a)>>;
 296           constexpr int under = opt::test(axis::option::underflow);
 297           constexpr int over = opt::test(axis::option::overflow);
 298           const auto size = a.size();
 299
 300           ca->extent = size + under + over;
 301           // -1 if underflow and cover all, else 0
 302           ca->begin = cov == coverage::all ? -under : 0;
 303           // size + 1 if overflow and cover all, else size
 304           ca->end = cov == coverage::all ? size + over : size;
 305           ca->idx = ca->begin;
 306
 307           begin_.iter_ += (ca->begin + under) * stride;
 308           end_.iter_ += ((ca < clast ? ca->begin : ca->end) + under) * stride;
 309
 310           stride *= ca->extent;
 311           ++ca;
 312         });
 313   }
 314
 315   iterator begin() noexcept { return begin_; }
 316   iterator end() noexcept { return end_; }
 317
 318 private:
 319   iterator begin_, end_;
 320 };
 321
 322 /** Generates an indexed range of <a
 323   href="https://en.cppreference.com/w/cpp/named_req/ForwardIterator">forward iterators</a>
 324   over the histogram cells.
 325
 326   Use this in a range-based for loop:
 327
 328   ```
 329   for (auto&& x : indexed(hist)) { ... }
 330   ```
 331
 332   This generates an optimized loop which is nearly always faster than a hand-written loop
 333   over the histogram cells. The iterators dereference to an indexed_range::accessor, which
 334   has methods to query the current indices and bins and acts like a pointer to the cell
 335   value. The returned iterators are forward iterators. They can be stored in a container,
 336   but may not be used after the life-time of the histogram ends.
 337
 338   @returns indexed_range
 339
 340   @param hist Reference to the histogram.
 341   @param cov  Iterate over all or only inner bins (optional, default: inner).
 342  */
 343 template <class Histogram>
 344 auto indexed(Histogram&& hist, coverage cov = coverage::inner) {
 345   return indexed_range<std::remove_reference_t<Histogram>>{std::forward<Histogram>(hist),
 346                                                            cov};
 347 }
 348
 349 } // namespace histogram
 350 } // namespace boost
 351
 352 #endif