boost/histogram/axis/regular.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_AXIS_REGULAR_HPP
   8 #define BOOST_HISTOGRAM_AXIS_REGULAR_HPP
   9
  10 #include <boost/assert.hpp>
  11 #include <boost/core/nvp.hpp>
  12 #include <boost/histogram/axis/interval_view.hpp>
  13 #include <boost/histogram/axis/iterator.hpp>
  14 #include <boost/histogram/axis/metadata_base.hpp>
  15 #include <boost/histogram/axis/option.hpp>
  16 #include <boost/histogram/detail/convert_integer.hpp>
  17 #include <boost/histogram/detail/relaxed_equal.hpp>
  18 #include <boost/histogram/detail/replace_type.hpp>
  19 #include <boost/histogram/fwd.hpp>
  20 #include <boost/mp11/utility.hpp>
  21 #include <boost/throw_exception.hpp>
  22 #include <cmath>
  23 #include <limits>
  24 #include <stdexcept>
  25 #include <string>
  26 #include <type_traits>
  27 #include <utility>
  28
  29 namespace boost {
  30 namespace histogram {
  31 namespace detail {
  32
  33 template <class T>
  34 using get_scale_type_helper = typename T::value_type;
  35
  36 template <class T>
  37 using get_scale_type = mp11::mp_eval_or<T, detail::get_scale_type_helper, T>;
  38
  39 struct one_unit {};
  40
  41 template <class T>
  42 T operator*(T&& t, const one_unit&) {
  43   return std::forward<T>(t);
  44 }
  45
  46 template <class T>
  47 T operator/(T&& t, const one_unit&) {
  48   return std::forward<T>(t);
  49 }
  50
  51 template <class T>
  52 using get_unit_type_helper = typename T::unit_type;
  53
  54 template <class T>
  55 using get_unit_type = mp11::mp_eval_or<one_unit, detail::get_unit_type_helper, T>;
  56
  57 template <class T, class R = get_scale_type<T>>
  58 R get_scale(const T& t) {
  59   return t / get_unit_type<T>();
  60 }
  61
  62 } // namespace detail
  63
  64 namespace axis {
  65
  66 namespace transform {
  67
  68 /// Identity transform for equidistant bins.
  69 struct id {
  70   /// Pass-through.
  71   template <class T>
  72   static T forward(T&& x) noexcept {
  73     return std::forward<T>(x);
  74   }
  75
  76   /// Pass-through.
  77   template <class T>
  78   static T inverse(T&& x) noexcept {
  79     return std::forward<T>(x);
  80   }
  81
  82   template <class Archive>
  83   void serialize(Archive&, unsigned /* version */) {}
  84 };
  85
  86 /// Log transform for equidistant bins in log-space.
  87 struct log {
  88   /// Returns log(x) of external value x.
  89   template <class T>
  90   static T forward(T x) {
  91     return std::log(x);
  92   }
  93
  94   /// Returns exp(x) for internal value x.
  95   template <class T>
  96   static T inverse(T x) {
  97     return std::exp(x);
  98   }
  99
 100   template <class Archive>
 101   void serialize(Archive&, unsigned /* version */) {}
 102 };
 103
 104 /// Sqrt transform for equidistant bins in sqrt-space.
 105 struct sqrt {
 106   /// Returns sqrt(x) of external value x.
 107   template <class T>
 108   static T forward(T x) {
 109     return std::sqrt(x);
 110   }
 111
 112   /// Returns x^2 of internal value x.
 113   template <class T>
 114   static T inverse(T x) {
 115     return x * x;
 116   }
 117
 118   template <class Archive>
 119   void serialize(Archive&, unsigned /* version */) {}
 120 };
 121
 122 /// Pow transform for equidistant bins in pow-space.
 123 struct pow {
 124   double power = 1; /**< power index */
 125
 126   /// Make transform with index p.
 127   explicit pow(double p) : power(p) {}
 128   pow() = default;
 129
 130   /// Returns pow(x, power) of external value x.
 131   template <class T>
 132   auto forward(T x) const {
 133     return std::pow(x, power);
 134   }
 135
 136   /// Returns pow(x, 1/power) of external value x.
 137   template <class T>
 138   auto inverse(T x) const {
 139     return std::pow(x, 1.0 / power);
 140   }
 141
 142   bool operator==(const pow& o) const noexcept { return power == o.power; }
 143
 144   template <class Archive>
 145   void serialize(Archive& ar, unsigned /* version */) {
 146     ar& make_nvp("power", power);
 147   }
 148 };
 149
 150 } // namespace transform
 151
 152 #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
 153 // Type envelope to mark value as step size
 154 template <class T>
 155 struct step_type {
 156   T value;
 157 };
 158 #endif
 159
 160 /**
 161   Helper function to mark argument as step size.
 162  */
 163 template <class T>
 164 step_type<T> step(T t) {
 165   return step_type<T>{t};
 166 }
 167
 168 /**
 169   Axis for equidistant intervals on the real line.
 170
 171   The most common binning strategy. Very fast. Binning is a O(1) operation.
 172
 173   @tparam Value input value type, must be floating point.
 174   @tparam Transform builtin or user-defined transform type.
 175   @tparam MetaData type to store meta data.
 176   @tparam Options see boost::histogram::axis::option (all values allowed).
 177  */
 178 template <class Value, class Transform, class MetaData, class Options>
 179 class regular : public iterator_mixin<regular<Value, Transform, MetaData, Options>>,
 180                 protected detail::replace_default<Transform, transform::id>,
 181                 public metadata_base<MetaData> {
 182   using value_type = Value;
 183   using transform_type = detail::replace_default<Transform, transform::id>;
 184   using metadata_type = typename metadata_base<MetaData>::metadata_type;
 185   using options_type =
 186       detail::replace_default<Options, decltype(option::underflow | option::overflow)>;
 187
 188   static_assert(std::is_nothrow_move_constructible<transform_type>::value,
 189                 "transform must be no-throw move constructible");
 190   static_assert(std::is_nothrow_move_assignable<transform_type>::value,
 191                 "transform must be no-throw move assignable");
 192
 193   using unit_type = detail::get_unit_type<value_type>;
 194   using internal_value_type = detail::get_scale_type<value_type>;
 195
 196   static_assert(std::is_floating_point<internal_value_type>::value,
 197                 "regular axis requires floating point type");
 198
 199   static_assert(
 200       (!options_type::test(option::circular) && !options_type::test(option::growth)) ||
 201           (options_type::test(option::circular) ^ options_type::test(option::growth)),
 202       "circular and growth options are mutually exclusive");
 203
 204 public:
 205   constexpr regular() = default;
 206
 207   /** Construct n bins over real transformed range [start, stop).
 208    *
 209    * @param trans    transform instance to use.
 210    * @param n        number of bins.
 211    * @param start    low edge of first bin.
 212    * @param stop     high edge of last bin.
 213    * @param meta     description of the axis (optional).
 214    */
 215   regular(transform_type trans, unsigned n, value_type start, value_type stop,
 216           metadata_type meta = {})
 217       : transform_type(std::move(trans))
 218       , metadata_base<MetaData>(std::move(meta))
 219       , size_(static_cast<index_type>(n))
 220       , min_(this->forward(detail::get_scale(start)))
 221       , delta_(this->forward(detail::get_scale(stop)) - min_) {
 222     if (size() == 0) BOOST_THROW_EXCEPTION(std::invalid_argument("bins > 0 required"));
 223     if (!std::isfinite(min_) || !std::isfinite(delta_))
 224       BOOST_THROW_EXCEPTION(
 225           std::invalid_argument("forward transform of start or stop invalid"));
 226     if (delta_ == 0)
 227       BOOST_THROW_EXCEPTION(std::invalid_argument("range of axis is zero"));
 228   }
 229
 230   /** Construct n bins over real range [start, stop).
 231    *
 232    * @param n        number of bins.
 233    * @param start    low edge of first bin.
 234    * @param stop     high edge of last bin.
 235    * @param meta     description of the axis (optional).
 236    */
 237   regular(unsigned n, value_type start, value_type stop, metadata_type meta = {})
 238       : regular({}, n, start, stop, std::move(meta)) {}
 239
 240   /** Construct bins with the given step size over real transformed range
 241    * [start, stop).
 242    *
 243    * @param trans   transform instance to use.
 244    * @param step    width of a single bin.
 245    * @param start   low edge of first bin.
 246    * @param stop    upper limit of high edge of last bin (see below).
 247    * @param meta    description of the axis (optional).
 248    *
 249    * The axis computes the number of bins as n = abs(stop - start) / step,
 250    * rounded down. This means that stop is an upper limit to the actual value
 251    * (start + n * step).
 252    */
 253   template <class T>
 254   regular(transform_type trans, step_type<T> step, value_type start, value_type stop,
 255           metadata_type meta = {})
 256       : regular(trans, static_cast<index_type>(std::abs(stop - start) / step.value),
 257                 start,
 258                 start + static_cast<index_type>(std::abs(stop - start) / step.value) *
 259                             step.value,
 260                 std::move(meta)) {}
 261
 262   /** Construct bins with the given step size over real range [start, stop).
 263    *
 264    * @param step    width of a single bin.
 265    * @param start   low edge of first bin.
 266    * @param stop    upper limit of high edge of last bin (see below).
 267    * @param meta    description of the axis (optional).
 268    *
 269    * The axis computes the number of bins as n = abs(stop - start) / step,
 270    * rounded down. This means that stop is an upper limit to the actual value
 271    * (start + n * step).
 272    */
 273   template <class T>
 274   regular(step_type<T> step, value_type start, value_type stop, metadata_type meta = {})
 275       : regular({}, step, start, stop, std::move(meta)) {}
 276
 277   /// Constructor used by algorithm::reduce to shrink and rebin (not for users).
 278   regular(const regular& src, index_type begin, index_type end, unsigned merge)
 279       : regular(src.transform(), (end - begin) / merge, src.value(begin), src.value(end),
 280                 src.metadata()) {
 281     BOOST_ASSERT((end - begin) % merge == 0);
 282     if (options_type::test(option::circular) && !(begin == 0 && end == src.size()))
 283       BOOST_THROW_EXCEPTION(std::invalid_argument("cannot shrink circular axis"));
 284   }
 285
 286   /// Return instance of the transform type.
 287   const transform_type& transform() const noexcept { return *this; }
 288
 289   /// Return index for value argument.
 290   index_type index(value_type x) const noexcept {
 291     // Runs in hot loop, please measure impact of changes
 292     auto z = (this->forward(x / unit_type{}) - min_) / delta_;
 293     if (options_type::test(option::circular)) {
 294       if (std::isfinite(z)) {
 295         z -= std::floor(z);
 296         return static_cast<index_type>(z * size());
 297       }
 298     } else {
 299       if (z < 1) {
 300         if (z >= 0)
 301           return static_cast<index_type>(z * size());
 302         else
 303           return -1;
 304       }
 305     }
 306     return size(); // also returned if x is NaN
 307   }
 308
 309   /// Returns index and shift (if axis has grown) for the passed argument.
 310   auto update(value_type x) noexcept {
 311     BOOST_ASSERT(options_type::test(option::growth));
 312     const auto z = (this->forward(x / unit_type{}) - min_) / delta_;
 313     if (z < 1) { // don't use i here!
 314       if (z >= 0) {
 315         const auto i = static_cast<axis::index_type>(z * size());
 316         return std::make_pair(i, 0);
 317       }
 318       if (z != -std::numeric_limits<internal_value_type>::infinity()) {
 319         const auto stop = min_ + delta_;
 320         const auto i = static_cast<axis::index_type>(std::floor(z * size()));
 321         min_ += i * (delta_ / size());
 322         delta_ = stop - min_;
 323         size_ -= i;
 324         return std::make_pair(0, -i);
 325       }
 326       // z is -infinity
 327       return std::make_pair(-1, 0);
 328     }
 329     // z either beyond range, infinite, or NaN
 330     if (z < std::numeric_limits<internal_value_type>::infinity()) {
 331       const auto i = static_cast<axis::index_type>(z * size());
 332       const auto n = i - size() + 1;
 333       delta_ /= size();
 334       delta_ *= size() + n;
 335       size_ += n;
 336       return std::make_pair(i, -n);
 337     }
 338     // z either infinite or NaN
 339     return std::make_pair(size(), 0);
 340   }
 341
 342   /// Return value for fractional index argument.
 343   value_type value(real_index_type i) const noexcept {
 344     auto z = i / size();
 345     if (!options_type::test(option::circular) && z < 0.0)
 346       z = -std::numeric_limits<internal_value_type>::infinity() * delta_;
 347     else if (options_type::test(option::circular) || z <= 1.0)
 348       z = (1.0 - z) * min_ + z * (min_ + delta_);
 349     else {
 350       z = std::numeric_limits<internal_value_type>::infinity() * delta_;
 351     }
 352     return static_cast<value_type>(this->inverse(z) * unit_type());
 353   }
 354
 355   /// Return bin for index argument.
 356   decltype(auto) bin(index_type idx) const noexcept {
 357     return interval_view<regular>(*this, idx);
 358   }
 359
 360   /// Returns the number of bins, without over- or underflow.
 361   index_type size() const noexcept { return size_; }
 362
 363   /// Returns the options.
 364   static constexpr unsigned options() noexcept { return options_type::value; }
 365
 366   template <class V, class T, class M, class O>
 367   bool operator==(const regular<V, T, M, O>& o) const noexcept {
 368     return detail::relaxed_equal(transform(), o.transform()) && size() == o.size() &&
 369            min_ == o.min_ && delta_ == o.delta_ && metadata_base<MetaData>::operator==(o);
 370   }
 371   template <class V, class T, class M, class O>
 372   bool operator!=(const regular<V, T, M, O>& o) const noexcept {
 373     return !operator==(o);
 374   }
 375
 376   template <class Archive>
 377   void serialize(Archive& ar, unsigned /* version */) {
 378     ar& make_nvp("transform", static_cast<transform_type&>(*this));
 379     ar& make_nvp("size", size_);
 380     ar& make_nvp("meta", this->metadata());
 381     ar& make_nvp("min", min_);
 382     ar& make_nvp("delta", delta_);
 383   }
 384
 385 private:
 386   index_type size_{0};
 387   internal_value_type min_{0}, delta_{1};
 388
 389   template <class V, class T, class M, class O>
 390   friend class regular;
 391 };
 392
 393 #if __cpp_deduction_guides >= 201606
 394
 395 template <class T>
 396 regular(unsigned, T, T)
 397     ->regular<detail::convert_integer<T, double>, transform::id, null_type>;
 398
 399 template <class T, class M>
 400 regular(unsigned, T, T, M)
 401     ->regular<detail::convert_integer<T, double>, transform::id,
 402               detail::replace_cstring<std::decay_t<M>>>;
 403
 404 template <class Tr, class T, class = detail::requires_transform<Tr, T>>
 405 regular(Tr, unsigned, T, T)->regular<detail::convert_integer<T, double>, Tr, null_type>;
 406
 407 template <class Tr, class T, class M>
 408 regular(Tr, unsigned, T, T, M)
 409     ->regular<detail::convert_integer<T, double>, Tr,
 410               detail::replace_cstring<std::decay_t<M>>>;
 411
 412 #endif
 413
 414 /// Regular axis with circular option already set.
 415 template <class Value = double, class MetaData = use_default, class Options = use_default>
 416 #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
 417 using circular = regular<Value, transform::id, MetaData,
 418                          decltype(detail::replace_default<Options, option::overflow_t>{} |
 419                                   option::circular)>;
 420 #else
 421 class circular;
 422 #endif
 423
 424 } // namespace axis
 425 } // namespace histogram
 426 } // namespace boost
 427
 428 #endif