boost/histogram/histogram.hpp

   1 // Copyright 2015-2018 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_HISTOGRAM_HPP
   8 #define BOOST_HISTOGRAM_HISTOGRAM_HPP
   9
  10 #include <boost/histogram/detail/accumulator_traits.hpp>
  11 #include <boost/histogram/detail/argument_traits.hpp>
  12 #include <boost/histogram/detail/at.hpp>
  13 #include <boost/histogram/detail/axes.hpp>
  14 #include <boost/histogram/detail/common_type.hpp>
  15 #include <boost/histogram/detail/fill.hpp>
  16 #include <boost/histogram/detail/fill_n.hpp>
  17 #include <boost/histogram/detail/mutex_base.hpp>
  18 #include <boost/histogram/detail/non_member_container_access.hpp>
  19 #include <boost/histogram/detail/span.hpp>
  20 #include <boost/histogram/fwd.hpp>
  21 #include <boost/histogram/sample.hpp>
  22 #include <boost/histogram/storage_adaptor.hpp>
  23 #include <boost/histogram/unsafe_access.hpp>
  24 #include <boost/histogram/weight.hpp>
  25 #include <boost/mp11/integral.hpp>
  26 #include <boost/mp11/list.hpp>
  27 #include <boost/mp11/tuple.hpp>
  28 #include <boost/throw_exception.hpp>
  29 #include <mutex>
  30 #include <stdexcept>
  31 #include <tuple>
  32 #include <type_traits>
  33 #include <utility>
  34 #include <vector>
  35
  36 namespace boost {
  37 namespace histogram {
  38
  39 /** Central class of the histogram library.
  40
  41   Histogram uses the call operator to insert data, like the
  42   [Boost.Accumulators](https://www.boost.org/doc/libs/develop/doc/html/accumulators.html).
  43
  44   Use factory functions (see
  45   [make_histogram.hpp](histogram/reference.html#header.boost.histogram.make_histogram_hpp)
  46   and
  47   [make_profile.hpp](histogram/reference.html#header.boost.histogram.make_profile_hpp)) to
  48   conveniently create histograms rather than calling the ctors directly.
  49
  50   Use the [indexed](boost/histogram/indexed.html) range generator to iterate over filled
  51   histograms, which is convenient and faster than hand-written loops for multi-dimensional
  52   histograms.
  53
  54   @tparam Axes std::tuple of axis types OR std::vector of an axis type or axis::variant
  55   @tparam Storage class that implements the storage interface
  56  */
  57 template <class Axes, class Storage>
  58 class histogram : detail::mutex_base<Axes, Storage> {
  59   static_assert(std::is_same<std::decay_t<Storage>, Storage>::value,
  60                 "Storage may not be a reference or const or volatile");
  61   static_assert(mp11::mp_size<Axes>::value > 0, "at least one axis required");
  62
  63 public:
  64   using axes_type = Axes;
  65   using storage_type = Storage;
  66   using value_type = typename storage_type::value_type;
  67   // typedefs for boost::range_iterator
  68   using iterator = typename storage_type::iterator;
  69   using const_iterator = typename storage_type::const_iterator;
  70
  71 private:
  72   using mutex_base = typename detail::mutex_base<axes_type, storage_type>;
  73
  74 public:
  75   histogram() = default;
  76
  77   template <class A, class S>
  78   explicit histogram(histogram<A, S>&& rhs)
  79       : storage_(std::move(unsafe_access::storage(rhs)))
  80       , offset_(unsafe_access::offset(rhs)) {
  81     detail::axes_assign(axes_, std::move(unsafe_access::axes(rhs)));
  82     detail::throw_if_axes_is_too_large(axes_);
  83   }
  84
  85   template <class A, class S>
  86   explicit histogram(const histogram<A, S>& rhs)
  87       : storage_(unsafe_access::storage(rhs)), offset_(unsafe_access::offset(rhs)) {
  88     detail::axes_assign(axes_, unsafe_access::axes(rhs));
  89     detail::throw_if_axes_is_too_large(axes_);
  90   }
  91
  92   template <class A, class S>
  93   histogram& operator=(histogram<A, S>&& rhs) {
  94     detail::axes_assign(axes_, std::move(unsafe_access::axes(rhs)));
  95     detail::throw_if_axes_is_too_large(axes_);
  96     storage_ = std::move(unsafe_access::storage(rhs));
  97     offset_ = unsafe_access::offset(rhs);
  98     return *this;
  99   }
 100
 101   template <class A, class S>
 102   histogram& operator=(const histogram<A, S>& rhs) {
 103     detail::axes_assign(axes_, unsafe_access::axes(rhs));
 104     detail::throw_if_axes_is_too_large(axes_);
 105     storage_ = unsafe_access::storage(rhs);
 106     offset_ = unsafe_access::offset(rhs);
 107     return *this;
 108   }
 109
 110   template <class A, class = detail::requires_axes<A>>
 111   histogram(A&& a, Storage s)
 112       : axes_(std::forward<A>(a))
 113       , storage_(std::move(s))
 114       , offset_(detail::offset(axes_)) {
 115     detail::throw_if_axes_is_too_large(axes_);
 116     storage_.reset(detail::bincount(axes_));
 117   }
 118
 119   explicit histogram(Axes axes) : histogram(axes, storage_type()) {}
 120
 121   template <class... As, class = detail::requires_axes<std::tuple<std::decay_t<As>...>>>
 122   explicit histogram(As&&... as)
 123       : histogram(std::tuple<std::decay_t<As>...>(std::forward<As>(as)...),
 124                   storage_type()) {}
 125
 126   /// Number of axes (dimensions).
 127   constexpr unsigned rank() const noexcept { return detail::axes_rank(axes_); }
 128
 129   /// Total number of bins (including underflow/overflow).
 130   std::size_t size() const noexcept { return storage_.size(); }
 131
 132   /// Reset all bins to default initialized values.
 133   void reset() { storage_.reset(size()); }
 134
 135   /// Get N-th axis using a compile-time number.
 136   /// This version is more efficient than the one accepting a run-time number.
 137   template <unsigned N = 0>
 138   decltype(auto) axis(std::integral_constant<unsigned, N> = {}) const {
 139     detail::axis_index_is_valid(axes_, N);
 140     return detail::axis_get<N>(axes_);
 141   }
 142
 143   /// Get N-th axis with run-time number.
 144   /// Prefer the version that accepts a compile-time number, if you can use it.
 145   decltype(auto) axis(unsigned i) const {
 146     detail::axis_index_is_valid(axes_, i);
 147     return detail::axis_get(axes_, i);
 148   }
 149
 150   /// Apply unary functor/function to each axis.
 151   template <class Unary>
 152   auto for_each_axis(Unary&& unary) const {
 153     return detail::for_each_axis(axes_, std::forward<Unary>(unary));
 154   }
 155
 156   /** Fill histogram with values, an optional weight, and/or a sample.
 157
 158     Arguments are passed in order to the axis objects. Passing an argument type that is
 159     not convertible to the value type accepted by the axis or passing the wrong number
 160     of arguments causes a throw of `std::invalid_argument`.
 161
 162     __Optional weight__
 163
 164     An optional weight can be passed as the first or last argument
 165     with the [weight](boost/histogram/weight.html) helper function. Compilation fails if
 166     the storage elements do not support weights.
 167
 168     __Samples__
 169
 170     If the storage elements accept samples, pass them with the sample helper function
 171     in addition to the axis arguments, which can be the first or last argument. The
 172     [sample](boost/histogram/sample.html) helper function can pass one or more arguments
 173     to the storage element. If samples and weights are used together, they can be passed
 174     in any order at the beginning or end of the argument list.
 175
 176     __Axis with multiple arguments__
 177
 178     If the histogram contains an axis which accepts a `std::tuple` of arguments, the
 179     arguments for that axis need to passed as a `std::tuple`, for example,
 180     `std::make_tuple(1.2, 2.3)`. If the histogram contains only this axis and no other,
 181     the arguments can be passed directly.
 182   */
 183   template <class Arg0, class... Args>
 184   std::enable_if_t<(detail::is_tuple<Arg0>::value == false || sizeof...(Args) > 0),
 185                    iterator>
 186   operator()(const Arg0& arg0, const Args&... args) {
 187     return operator()(std::forward_as_tuple(arg0, args...));
 188   }
 189
 190   /// Fill histogram with values, an optional weight, and/or a sample from a `std::tuple`.
 191   template <class... Ts>
 192   iterator operator()(const std::tuple<Ts...>& args) {
 193     using arg_traits = detail::argument_traits<std::decay_t<Ts>...>;
 194     using acc_traits = detail::accumulator_traits<value_type>;
 195     constexpr bool weight_valid =
 196         arg_traits::wpos::value == -1 || acc_traits::wsupport::value;
 197     static_assert(weight_valid, "error: accumulator does not support weights");
 198     detail::sample_args_passed_vs_expected<typename arg_traits::sargs,
 199                                            typename acc_traits::args>();
 200     constexpr bool sample_valid =
 201         std::is_convertible<typename arg_traits::sargs, typename acc_traits::args>::value;
 202     std::lock_guard<typename mutex_base::type> guard{mutex_base::get()};
 203     return detail::fill(mp11::mp_bool<(weight_valid && sample_valid)>{}, arg_traits{},
 204                         offset_, storage_, axes_, args);
 205   }
 206
 207   /** Fill histogram with several values at once.
 208
 209     The argument must be an iterable with a size that matches the
 210     rank of the histogram. The element of an iterable may be 1) a value or 2) an iterable
 211     with contiguous storage over values or 3) a variant of 1) and 2). Sub-iterables must
 212     have the same length.
 213
 214     Values are passed to the corresponding histogram axis in order. If a single value is
 215     passed together with an iterable of values, the single value is treated like an
 216     iterable with matching length of copies of this value.
 217
 218     If the histogram has only one axis, an iterable of values may be passed directly.
 219
 220     @param args iterable as explained in the long description.
 221   */
 222   template <class Iterable, class = detail::requires_iterable<Iterable>>
 223   void fill(const Iterable& args) {
 224     using acc_traits = detail::accumulator_traits<value_type>;
 225     constexpr unsigned n_sample_args_expected =
 226         std::tuple_size<typename acc_traits::args>::value;
 227     static_assert(n_sample_args_expected == 0,
 228                   "sample argument is missing but required by accumulator");
 229     std::lock_guard<typename mutex_base::type> guard{mutex_base::get()};
 230     detail::fill_n(mp11::mp_bool<(n_sample_args_expected == 0)>{}, offset_, storage_,
 231                    axes_, detail::make_span(args));
 232   }
 233
 234   /** Fill histogram with several values and weights at once.
 235
 236     @param args iterable of values.
 237     @param weights single weight or an iterable of weights.
 238   */
 239   template <class Iterable, class T, class = detail::requires_iterable<Iterable>>
 240   void fill(const Iterable& args, const weight_type<T>& weights) {
 241     using acc_traits = detail::accumulator_traits<value_type>;
 242     constexpr bool weight_valid = acc_traits::wsupport::value;
 243     static_assert(weight_valid, "error: accumulator does not support weights");
 244     detail::sample_args_passed_vs_expected<std::tuple<>, typename acc_traits::args>();
 245     constexpr bool sample_valid =
 246         std::is_convertible<std::tuple<>, typename acc_traits::args>::value;
 247     std::lock_guard<typename mutex_base::type> guard{mutex_base::get()};
 248     detail::fill_n(mp11::mp_bool<(weight_valid && sample_valid)>{}, offset_, storage_,
 249                    axes_, detail::make_span(args),
 250                    weight(detail::to_ptr_size(weights.value)));
 251   }
 252
 253   /** Fill histogram with several values and weights at once.
 254
 255     @param weights single weight or an iterable of weights.
 256     @param args iterable of values.
 257   */
 258   template <class Iterable, class T, class = detail::requires_iterable<Iterable>>
 259   void fill(const weight_type<T>& weights, const Iterable& args) {
 260     fill(args, weights);
 261   }
 262
 263   /** Fill histogram with several values and samples at once.
 264
 265     @param args iterable of values.
 266     @param samples single sample or an iterable of samples.
 267   */
 268   template <class Iterable, class... Ts, class = detail::requires_iterable<Iterable>>
 269   void fill(const Iterable& args, const sample_type<std::tuple<Ts...>>& samples) {
 270     using acc_traits = detail::accumulator_traits<value_type>;
 271     using sample_args_passed =
 272         std::tuple<decltype(*detail::to_ptr_size(std::declval<Ts>()).first)...>;
 273     detail::sample_args_passed_vs_expected<sample_args_passed,
 274                                            typename acc_traits::args>();
 275     std::lock_guard<typename mutex_base::type> guard{mutex_base::get()};
 276     mp11::tuple_apply(
 277         [&](const auto&... sargs) {
 278           constexpr bool sample_valid =
 279               std::is_convertible<sample_args_passed, typename acc_traits::args>::value;
 280           detail::fill_n(mp11::mp_bool<(sample_valid)>{}, offset_, storage_, axes_,
 281                          detail::make_span(args), detail::to_ptr_size(sargs)...);
 282         },
 283         samples.value);
 284   }
 285
 286   /** Fill histogram with several values and samples at once.
 287
 288     @param samples single sample or an iterable of samples.
 289     @param args iterable of values.
 290   */
 291   template <class Iterable, class T, class = detail::requires_iterable<Iterable>>
 292   void fill(const sample_type<T>& samples, const Iterable& args) {
 293     fill(args, samples);
 294   }
 295
 296   template <class Iterable, class T, class... Ts,
 297             class = detail::requires_iterable<Iterable>>
 298   void fill(const Iterable& args, const weight_type<T>& weights,
 299             const sample_type<std::tuple<Ts...>>& samples) {
 300     using acc_traits = detail::accumulator_traits<value_type>;
 301     using sample_args_passed =
 302         std::tuple<decltype(*detail::to_ptr_size(std::declval<Ts>()).first)...>;
 303     detail::sample_args_passed_vs_expected<sample_args_passed,
 304                                            typename acc_traits::args>();
 305     std::lock_guard<typename mutex_base::type> guard{mutex_base::get()};
 306     mp11::tuple_apply(
 307         [&](const auto&... sargs) {
 308           constexpr bool weight_valid = acc_traits::wsupport::value;
 309           static_assert(weight_valid, "error: accumulator does not support weights");
 310           constexpr bool sample_valid =
 311               std::is_convertible<sample_args_passed, typename acc_traits::args>::value;
 312           detail::fill_n(mp11::mp_bool<(weight_valid && sample_valid)>{}, offset_,
 313                          storage_, axes_, detail::make_span(args),
 314                          weight(detail::to_ptr_size(weights.value)),
 315                          detail::to_ptr_size(sargs)...);
 316         },
 317         samples.value);
 318   }
 319
 320   template <class Iterable, class T, class U, class = detail::requires_iterable<Iterable>>
 321   void fill(const sample_type<T>& samples, const weight_type<U>& weights,
 322             const Iterable& args) {
 323     fill(args, weights, samples);
 324   }
 325
 326   template <class Iterable, class T, class U, class = detail::requires_iterable<Iterable>>
 327   void fill(const weight_type<T>& weights, const sample_type<U>& samples,
 328             const Iterable& args) {
 329     fill(args, weights, samples);
 330   }
 331
 332   template <class Iterable, class T, class U, class = detail::requires_iterable<Iterable>>
 333   void fill(const Iterable& args, const sample_type<T>& samples,
 334             const weight_type<U>& weights) {
 335     fill(args, weights, samples);
 336   }
 337
 338   /** Access cell value at integral indices.
 339
 340     You can pass indices as individual arguments, as a std::tuple of integers, or as an
 341     interable range of integers. Passing the wrong number of arguments causes a throw of
 342     std::invalid_argument. Passing an index which is out of bounds causes a throw of
 343     std::out_of_range.
 344
 345     @param i index of first axis.
 346     @param is indices of second, third, ... axes.
 347     @returns reference to cell value.
 348   */
 349   template <class... Indices>
 350   decltype(auto) at(axis::index_type i, Indices... is) {
 351     return at(std::forward_as_tuple(i, is...));
 352   }
 353
 354   /// Access cell value at integral indices (read-only).
 355   template <class... Indices>
 356   decltype(auto) at(axis::index_type i, Indices... is) const {
 357     return at(std::forward_as_tuple(i, is...));
 358   }
 359
 360   /// Access cell value at integral indices stored in `std::tuple`.
 361   template <class... Indices>
 362   decltype(auto) at(const std::tuple<Indices...>& is) {
 363     if (rank() != sizeof...(Indices))
 364       BOOST_THROW_EXCEPTION(
 365           std::invalid_argument("number of arguments != histogram rank"));
 366     const auto idx = detail::at(axes_, is);
 367     if (!is_valid(idx))
 368       BOOST_THROW_EXCEPTION(std::out_of_range("at least one index out of bounds"));
 369     BOOST_ASSERT(idx < storage_.size());
 370     return storage_[idx];
 371   }
 372
 373   /// Access cell value at integral indices stored in `std::tuple` (read-only).
 374   template <typename... Indices>
 375   decltype(auto) at(const std::tuple<Indices...>& is) const {
 376     if (rank() != sizeof...(Indices))
 377       BOOST_THROW_EXCEPTION(
 378           std::invalid_argument("number of arguments != histogram rank"));
 379     const auto idx = detail::at(axes_, is);
 380     if (!is_valid(idx))
 381       BOOST_THROW_EXCEPTION(std::out_of_range("at least one index out of bounds"));
 382     BOOST_ASSERT(idx < storage_.size());
 383     return storage_[idx];
 384   }
 385
 386   /// Access cell value at integral indices stored in iterable.
 387   template <class Iterable, class = detail::requires_iterable<Iterable>>
 388   decltype(auto) at(const Iterable& is) {
 389     if (rank() != detail::axes_rank(is))
 390       BOOST_THROW_EXCEPTION(
 391           std::invalid_argument("number of arguments != histogram rank"));
 392     const auto idx = detail::at(axes_, is);
 393     if (!is_valid(idx))
 394       BOOST_THROW_EXCEPTION(std::out_of_range("at least one index out of bounds"));
 395     BOOST_ASSERT(idx < storage_.size());
 396     return storage_[idx];
 397   }
 398
 399   /// Access cell value at integral indices stored in iterable (read-only).
 400   template <class Iterable, class = detail::requires_iterable<Iterable>>
 401   decltype(auto) at(const Iterable& is) const {
 402     if (rank() != detail::axes_rank(is))
 403       BOOST_THROW_EXCEPTION(
 404           std::invalid_argument("number of arguments != histogram rank"));
 405     const auto idx = detail::at(axes_, is);
 406     if (!is_valid(idx))
 407       BOOST_THROW_EXCEPTION(std::out_of_range("at least one index out of bounds"));
 408     BOOST_ASSERT(idx < storage_.size());
 409     return storage_[idx];
 410   }
 411
 412   /// Access value at index (number for rank = 1, else `std::tuple` or iterable).
 413   template <class Indices>
 414   decltype(auto) operator[](const Indices& is) {
 415     return at(is);
 416   }
 417
 418   /// Access value at index (read-only).
 419   template <class Indices>
 420   decltype(auto) operator[](const Indices& is) const {
 421     return at(is);
 422   }
 423
 424   /// Equality operator, tests equality for all axes and the storage.
 425   template <class A, class S>
 426   bool operator==(const histogram<A, S>& rhs) const noexcept {
 427     // testing offset is redundant, but offers fast return if it fails
 428     return offset_ == unsafe_access::offset(rhs) &&
 429            detail::axes_equal(axes_, unsafe_access::axes(rhs)) &&
 430            storage_ == unsafe_access::storage(rhs);
 431   }
 432
 433   /// Negation of the equality operator.
 434   template <class A, class S>
 435   bool operator!=(const histogram<A, S>& rhs) const noexcept {
 436     return !operator==(rhs);
 437   }
 438
 439   /// Add values of another histogram.
 440   template <class A, class S>
 441   std::enable_if_t<
 442       detail::has_operator_radd<value_type, typename histogram<A, S>::value_type>::value,
 443       histogram&>
 444   operator+=(const histogram<A, S>& rhs) {
 445     if (!detail::axes_equal(axes_, unsafe_access::axes(rhs)))
 446       BOOST_THROW_EXCEPTION(std::invalid_argument("axes of histograms differ"));
 447     auto rit = unsafe_access::storage(rhs).begin();
 448     for (auto&& x : storage_) x += *rit++;
 449     return *this;
 450   }
 451
 452   /// Subtract values of another histogram.
 453   template <class A, class S>
 454   std::enable_if_t<
 455       detail::has_operator_rsub<value_type, typename histogram<A, S>::value_type>::value,
 456       histogram&>
 457   operator-=(const histogram<A, S>& rhs) {
 458     if (!detail::axes_equal(axes_, unsafe_access::axes(rhs)))
 459       BOOST_THROW_EXCEPTION(std::invalid_argument("axes of histograms differ"));
 460     auto rit = unsafe_access::storage(rhs).begin();
 461     for (auto&& x : storage_) x -= *rit++;
 462     return *this;
 463   }
 464
 465   /// Multiply by values of another histogram.
 466   template <class A, class S>
 467   std::enable_if_t<
 468       detail::has_operator_rmul<value_type, typename histogram<A, S>::value_type>::value,
 469       histogram&>
 470   operator*=(const histogram<A, S>& rhs) {
 471     if (!detail::axes_equal(axes_, unsafe_access::axes(rhs)))
 472       BOOST_THROW_EXCEPTION(std::invalid_argument("axes of histograms differ"));
 473     auto rit = unsafe_access::storage(rhs).begin();
 474     for (auto&& x : storage_) x *= *rit++;
 475     return *this;
 476   }
 477
 478   /// Divide by values of another histogram.
 479   template <class A, class S>
 480   std::enable_if_t<
 481       detail::has_operator_rdiv<value_type, typename histogram<A, S>::value_type>::value,
 482       histogram&>
 483   operator/=(const histogram<A, S>& rhs) {
 484     if (!detail::axes_equal(axes_, unsafe_access::axes(rhs)))
 485       BOOST_THROW_EXCEPTION(std::invalid_argument("axes of histograms differ"));
 486     auto rit = unsafe_access::storage(rhs).begin();
 487     for (auto&& x : storage_) x /= *rit++;
 488     return *this;
 489   }
 490
 491   /// Multiply all values with a scalar.
 492   template <class V = value_type>
 493   std::enable_if_t<(detail::has_operator_rmul<V, double>::value &&
 494                     detail::has_operator_rmul<storage_type, double>::value == true),
 495                    histogram&>
 496   operator*=(const double x) {
 497     // use special implementation of scaling if available
 498     storage_ *= x;
 499     return *this;
 500   }
 501
 502   /// Multiply all values with a scalar.
 503   template <class V = value_type>
 504   std::enable_if_t<(detail::has_operator_rmul<V, double>::value &&
 505                     detail::has_operator_rmul<storage_type, double>::value == false),
 506                    histogram&>
 507   operator*=(const double x) {
 508     // generic implementation of scaling
 509     for (auto&& si : storage_) si *= x;
 510     return *this;
 511   }
 512
 513   /// Divide all values by a scalar.
 514   template <class V = value_type>
 515   std::enable_if_t<detail::has_operator_rmul<V, double>::value, histogram&> operator/=(
 516       const double x) {
 517     return operator*=(1.0 / x);
 518   }
 519
 520   /// Return value iterator to the beginning of the histogram.
 521   iterator begin() noexcept { return storage_.begin(); }
 522
 523   /// Return value iterator to the end in the histogram.
 524   iterator end() noexcept { return storage_.end(); }
 525
 526   /// Return value iterator to the beginning of the histogram (read-only).
 527   const_iterator begin() const noexcept { return storage_.begin(); }
 528
 529   /// Return value iterator to the end in the histogram (read-only).
 530   const_iterator end() const noexcept { return storage_.end(); }
 531
 532   /// Return value iterator to the beginning of the histogram (read-only).
 533   const_iterator cbegin() const noexcept { return begin(); }
 534
 535   /// Return value iterator to the end in the histogram (read-only).
 536   const_iterator cend() const noexcept { return end(); }
 537
 538   template <class Archive>
 539   void serialize(Archive& ar, unsigned /* version */) {
 540     detail::axes_serialize(ar, axes_);
 541     ar& make_nvp("storage", storage_);
 542     if (Archive::is_loading::value) {
 543       offset_ = detail::offset(axes_);
 544       detail::throw_if_axes_is_too_large(axes_);
 545     }
 546   }
 547
 548 private:
 549   axes_type axes_;
 550   storage_type storage_;
 551   std::size_t offset_ = 0;
 552
 553   friend struct unsafe_access;
 554 };
 555
 556 /**
 557   Pairwise add cells of two histograms and return histogram with the sum.
 558
 559   The returned histogram type is the most efficient and safest one constructible from the
 560   inputs, if they are not the same type. If one histogram has a tuple axis, the result has
 561   a tuple axis. The chosen storage is the one with the larger dynamic range.
 562 */
 563 template <class A1, class S1, class A2, class S2>
 564 auto operator+(const histogram<A1, S1>& a, const histogram<A2, S2>& b) {
 565   auto r = histogram<detail::common_axes<A1, A2>, detail::common_storage<S1, S2>>(a);
 566   return r += b;
 567 }
 568
 569 /** Pairwise multiply cells of two histograms and return histogram with the product.
 570
 571   For notes on the returned histogram type, see operator+.
 572 */
 573 template <class A1, class S1, class A2, class S2>
 574 auto operator*(const histogram<A1, S1>& a, const histogram<A2, S2>& b) {
 575   auto r = histogram<detail::common_axes<A1, A2>, detail::common_storage<S1, S2>>(a);
 576   return r *= b;
 577 }
 578
 579 /** Pairwise subtract cells of two histograms and return histogram with the difference.
 580
 581   For notes on the returned histogram type, see operator+.
 582 */
 583 template <class A1, class S1, class A2, class S2>
 584 auto operator-(const histogram<A1, S1>& a, const histogram<A2, S2>& b) {
 585   auto r = histogram<detail::common_axes<A1, A2>, detail::common_storage<S1, S2>>(a);
 586   return r -= b;
 587 }
 588
 589 /** Pairwise divide cells of two histograms and return histogram with the quotient.
 590
 591   For notes on the returned histogram type, see operator+.
 592 */
 593 template <class A1, class S1, class A2, class S2>
 594 auto operator/(const histogram<A1, S1>& a, const histogram<A2, S2>& b) {
 595   auto r = histogram<detail::common_axes<A1, A2>, detail::common_storage<S1, S2>>(a);
 596   return r /= b;
 597 }
 598
 599 /** Multiply all cells of the histogram by a number and return a new histogram.
 600
 601   If the original histogram has integer cells, the result has double cells.
 602 */
 603 template <class A, class S>
 604 auto operator*(const histogram<A, S>& h, double x) {
 605   auto r = histogram<A, detail::common_storage<S, dense_storage<double>>>(h);
 606   return r *= x;
 607 }
 608
 609 /** Multiply all cells of the histogram by a number and return a new histogram.
 610
 611   If the original histogram has integer cells, the result has double cells.
 612 */
 613 template <class A, class S>
 614 auto operator*(double x, const histogram<A, S>& h) {
 615   return h * x;
 616 }
 617
 618 /** Divide all cells of the histogram by a number and return a new histogram.
 619
 620   If the original histogram has integer cells, the result has double cells.
 621 */
 622 template <class A, class S>
 623 auto operator/(const histogram<A, S>& h, double x) {
 624   return h * (1.0 / x);
 625 }
 626
 627 #if __cpp_deduction_guides >= 201606
 628
 629 template <class... Axes, class = detail::requires_axes<std::tuple<std::decay_t<Axes>...>>>
 630 histogram(Axes...)->histogram<std::tuple<std::decay_t<Axes>...>>;
 631
 632 template <class... Axes, class S, class = detail::requires_storage_or_adaptible<S>>
 633 histogram(std::tuple<Axes...>, S)
 634     ->histogram<std::tuple<Axes...>, std::conditional_t<detail::is_adaptible<S>::value,
 635                                                         storage_adaptor<S>, S>>;
 636
 637 template <class Iterable, class = detail::requires_iterable<Iterable>,
 638           class = detail::requires_any_axis<typename Iterable::value_type>>
 639 histogram(Iterable)->histogram<std::vector<typename Iterable::value_type>>;
 640
 641 template <class Iterable, class S, class = detail::requires_iterable<Iterable>,
 642           class = detail::requires_any_axis<typename Iterable::value_type>,
 643           class = detail::requires_storage_or_adaptible<S>>
 644 histogram(Iterable, S)
 645     ->histogram<
 646         std::vector<typename Iterable::value_type>,
 647         std::conditional_t<detail::is_adaptible<S>::value, storage_adaptor<S>, S>>;
 648
 649 #endif
 650
 651 } // namespace histogram
 652 } // namespace boost
 653
 654 #endif