histogram/axis/regular.hpp

0001 // Copyright 2015-2018 Hans Dembinski
0002 //
0003 // Distributed under the Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_HISTOGRAM_AXIS_REGULAR_HPP
0008 #define BOOST_HISTOGRAM_AXIS_REGULAR_HPP
0009
0010 #include <boost/core/nvp.hpp>
0011 #include <boost/histogram/axis/interval_view.hpp>
0012 #include <boost/histogram/axis/iterator.hpp>
0013 #include <boost/histogram/axis/metadata_base.hpp>
0014 #include <boost/histogram/axis/option.hpp>
0015 #include <boost/histogram/detail/convert_integer.hpp>
0016 #include <boost/histogram/detail/relaxed_equal.hpp>
0017 #include <boost/histogram/detail/replace_type.hpp>
0018 #include <boost/histogram/fwd.hpp>
0019 #include <boost/mp11/utility.hpp>
0020 #include <boost/throw_exception.hpp>
0021 #include <cassert>
0022 #include <cmath>
0023 #include <limits>
0024 #include <stdexcept>
0025 #include <string>
0026 #include <type_traits>
0027 #include <utility>
0028
0029 namespace boost {
0030 namespace histogram {
0031 namespace detail {
0032
0033 template <class T>
0034 using get_scale_type_helper = typename T::value_type;
0035
0036 template <class T>
0037 using get_scale_type = mp11::mp_eval_or<T, detail::get_scale_type_helper, T>;
0038
0039 struct one_unit {};
0040
0041 template <class T>
0042 T operator*(T&& t, const one_unit&) {
0043   return std::forward<T>(t);
0044 }
0045
0046 template <class T>
0047 T operator/(T&& t, const one_unit&) {
0048   return std::forward<T>(t);
0049 }
0050
0051 template <class T>
0052 using get_unit_type_helper = typename T::unit_type;
0053
0054 template <class T>
0055 using get_unit_type = mp11::mp_eval_or<one_unit, detail::get_unit_type_helper, T>;
0056
0057 template <class T, class R = get_scale_type<T>>
0058 R get_scale(const T& t) {
0059   return t / get_unit_type<T>();
0060 }
0061
0062 } // namespace detail
0063
0064 namespace axis {
0065
0066 namespace transform {
0067
0068 /// Identity transform for equidistant bins.
0069 struct id {
0070   /// Pass-through.
0071   template <class T>
0072   static T forward(T&& x) noexcept {
0073     return std::forward<T>(x);
0074   }
0075
0076   /// Pass-through.
0077   template <class T>
0078   static T inverse(T&& x) noexcept {
0079     return std::forward<T>(x);
0080   }
0081
0082   template <class Archive>
0083   void serialize(Archive&, unsigned /* version */) {}
0084 };
0085
0086 /// Log transform for equidistant bins in log-space.
0087 struct log {
0088   /// Returns log(x) of external value x.
0089   template <class T>
0090   static T forward(T x) {
0091     return std::log(x);
0092   }
0093
0094   /// Returns exp(x) for internal value x.
0095   template <class T>
0096   static T inverse(T x) {
0097     return std::exp(x);
0098   }
0099
0100   template <class Archive>
0101   void serialize(Archive&, unsigned /* version */) {}
0102 };
0103
0104 /// Sqrt transform for equidistant bins in sqrt-space.
0105 struct sqrt {
0106   /// Returns sqrt(x) of external value x.
0107   template <class T>
0108   static T forward(T x) {
0109     return std::sqrt(x);
0110   }
0111
0112   /// Returns x^2 of internal value x.
0113   template <class T>
0114   static T inverse(T x) {
0115     return x * x;
0116   }
0117
0118   template <class Archive>
0119   void serialize(Archive&, unsigned /* version */) {}
0120 };
0121
0122 /// Pow transform for equidistant bins in pow-space.
0123 struct pow {
0124   double power = 1; /**< power index */
0125
0126   /// Make transform with index p.
0127   explicit pow(double p) : power(p) {}
0128   pow() = default;
0129
0130   /// Returns pow(x, power) of external value x.
0131   template <class T>
0132   auto forward(T x) const {
0133     return std::pow(x, power);
0134   }
0135
0136   /// Returns pow(x, 1/power) of external value x.
0137   template <class T>
0138   auto inverse(T x) const {
0139     return std::pow(x, 1.0 / power);
0140   }
0141
0142   bool operator==(const pow& o) const noexcept { return power == o.power; }
0143
0144   template <class Archive>
0145   void serialize(Archive& ar, unsigned /* version */) {
0146     ar& make_nvp("power", power);
0147   }
0148 };
0149
0150 } // namespace transform
0151
0152 #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
0153 // Type envelope to mark value as step size
0154 template <class T>
0155 struct step_type {
0156   T value;
0157 };
0158 #endif
0159
0160 /**
0161   Helper function to mark argument as step size.
0162  */
0163 template <class T>
0164 step_type<T> step(T t) {
0165   return step_type<T>{t};
0166 }
0167
0168 /** Axis for equidistant intervals on the real line.
0169
0170   The most common binning strategy. Very fast. Binning is a O(1) operation.
0171
0172   If the axis has an overflow bin (the default), a value on the upper edge of the last
0173   bin is put in the overflow bin. The axis range represents a semi-open interval.
0174
0175   If the overflow bin is deactivated, then a value on the upper edge of the last bin is
0176   still counted towards the last bin. The axis range represents a closed interval.
0177
0178   The options `growth` and `circular` are mutually exclusive.
0179
0180   @tparam Value input value type, must be floating point.
0181   @tparam Transform builtin or user-defined transform type.
0182   @tparam MetaData type to store meta data.
0183   @tparam Options see boost::histogram::axis::option.
0184  */
0185 template <class Value, class Transform, class MetaData, class Options>
0186 class regular : public iterator_mixin<regular<Value, Transform, MetaData, Options>>,
0187                 protected detail::replace_default<Transform, transform::id>,
0188                 public metadata_base_t<MetaData> {
0189   // these must be private, so that they are not automatically inherited
0190   using value_type = Value;
0191   using transform_type = detail::replace_default<Transform, transform::id>;
0192   using metadata_base = metadata_base_t<MetaData>;
0193   using metadata_type = typename metadata_base::metadata_type;
0194   using options_type =
0195       detail::replace_default<Options, decltype(option::underflow | option::overflow)>;
0196
0197   using unit_type = detail::get_unit_type<value_type>;
0198   using internal_value_type = detail::get_scale_type<value_type>;
0199
0200 public:
0201   constexpr regular() = default;
0202
0203   /** Construct n bins over real transformed range [start, stop).
0204
0205     @param trans    transform instance to use.
0206     @param n        number of bins.
0207     @param start    low edge of first bin.
0208     @param stop     high edge of last bin.
0209     @param meta     description of the axis (optional).
0210     @param options  see boost::histogram::axis::option (optional).
0211
0212     The constructor throws `std::invalid_argument` if n is zero, or if start and stop
0213     produce an invalid range after transformation.
0214
0215     The arguments meta and alloc are passed by value. If you move either of them into the
0216     axis and the constructor throws, their values are lost. Do not move if you cannot
0217     guarantee that the bin description is not valid.
0218    */
0219   regular(transform_type trans, unsigned n, value_type start, value_type stop,
0220           metadata_type meta = {}, options_type options = {})
0221       : transform_type(std::move(trans))
0222       , metadata_base(std::move(meta))
0223       , size_(static_cast<index_type>(n))
0224       , min_(this->forward(detail::get_scale(start)))
0225       , delta_(this->forward(detail::get_scale(stop)) - min_) {
0226     // static_asserts were moved here from class scope to satisfy deduction in gcc>=11
0227     static_assert(std::is_nothrow_move_constructible<transform_type>::value,
0228                   "transform must be no-throw move constructible");
0229     static_assert(std::is_nothrow_move_assignable<transform_type>::value,
0230                   "transform must be no-throw move assignable");
0231     static_assert(std::is_floating_point<internal_value_type>::value,
0232                   "regular axis requires floating point type");
0233     static_assert(!(options.test(option::circular) && options.test(option::growth)),
0234                   "circular and growth options are mutually exclusive");
0235     if (size() <= 0) BOOST_THROW_EXCEPTION(std::invalid_argument("bins > 0 required"));
0236     if (!std::isfinite(min_) || !std::isfinite(delta_))
0237       BOOST_THROW_EXCEPTION(
0238           std::invalid_argument("forward transform of start or stop invalid"));
0239     if (delta_ == 0)
0240       BOOST_THROW_EXCEPTION(std::invalid_argument("range of axis is zero"));
0241   }
0242
0243   /** Construct n bins over real range [start, stop).
0244
0245      @param n        number of bins.
0246      @param start    low edge of first bin.
0247      @param stop     high edge of last bin.
0248      @param meta     description of the axis (optional).
0249      @param options  see boost::histogram::axis::option (optional).
0250    */
0251   explicit regular(unsigned n, value_type start, value_type stop, metadata_type meta = {},
0252                    options_type options = {})
0253       : regular({}, n, start, stop, std::move(meta), options) {}
0254
0255   /** Construct bins with the given step size over real transformed range
0256      [start, stop).
0257
0258      @param trans    transform instance to use.
0259      @param step     width of a single bin.
0260      @param start    low edge of first bin.
0261      @param stop     upper limit of high edge of last bin (see below).
0262      @param meta     description of the axis (optional).
0263      @param options  see boost::histogram::axis::option (optional).
0264
0265      The axis computes the number of bins as n = abs(stop - start) / step,
0266      rounded down. This means that stop is an upper limit to the actual value
0267      (start + n * step).
0268    */
0269   template <class T>
0270   explicit regular(transform_type trans, step_type<T> step, value_type start,
0271                    value_type stop, metadata_type meta = {}, options_type options = {})
0272       : regular(trans, static_cast<index_type>(std::abs(stop - start) / step.value),
0273                 start,
0274                 start + static_cast<index_type>(std::abs(stop - start) / step.value) *
0275                             step.value,
0276                 std::move(meta), options) {}
0277
0278   /** Construct bins with the given step size over real range [start, stop).
0279
0280      @param step     width of a single bin.
0281      @param start    low edge of first bin.
0282      @param stop     upper limit of high edge of last bin (see below).
0283      @param meta     description of the axis (optional).
0284      @param options  see boost::histogram::axis::option (optional).
0285
0286      The axis computes the number of bins as n = abs(stop - start) / step,
0287      rounded down. This means that stop is an upper limit to the actual value
0288      (start + n * step).
0289    */
0290   template <class T>
0291   explicit regular(step_type<T> step, value_type start, value_type stop,
0292                    metadata_type meta = {}, options_type options = {})
0293       : regular({}, step, start, stop, std::move(meta), options) {}
0294
0295   /// Constructor used by algorithm::reduce to shrink and rebin (not for users).
0296   regular(const regular& src, index_type begin, index_type end, unsigned merge)
0297       : regular(src.transform(), (end - begin) / merge, src.value(begin), src.value(end),
0298                 src.metadata()) {
0299     assert((end - begin) % merge == 0);
0300     if (options_type::test(option::circular) && !(begin == 0 && end == src.size()))
0301       BOOST_THROW_EXCEPTION(std::invalid_argument("cannot shrink circular axis"));
0302   }
0303
0304   /// Return instance of the transform type.
0305   const transform_type& transform() const noexcept { return *this; }
0306
0307   /// Return index for value argument.
0308   index_type index(value_type x) const noexcept {
0309     // Runs in hot loop, please measure impact of changes
0310     auto z = (this->forward(x / unit_type{}) - min_) / delta_;
0311     if (options_type::test(option::circular)) {
0312       if (std::isfinite(z)) {
0313         z -= std::floor(z);
0314         return static_cast<index_type>(z * size());
0315       }
0316     } else {
0317       if (z < 1) {
0318         if (z >= 0)
0319           return static_cast<index_type>(z * size());
0320         else
0321           return -1;
0322       }
0323       // upper edge of last bin is inclusive if overflow bin is not present
0324       if (!options_type::test(option::overflow) && z == 1) return size() - 1;
0325     }
0326     return size(); // also returned if x is NaN
0327   }
0328
0329   /// Returns index and shift (if axis has grown) for the passed argument.
0330   std::pair<index_type, index_type> update(value_type x) noexcept {
0331     assert(options_type::test(option::growth));
0332     const auto z = (this->forward(x / unit_type{}) - min_) / delta_;
0333     if (z < 1) { // don't use i here!
0334       if (z >= 0) {
0335         const auto i = static_cast<axis::index_type>(z * size());
0336         return {i, 0};
0337       }
0338       if (z != -std::numeric_limits<internal_value_type>::infinity()) {
0339         const auto stop = min_ + delta_;
0340         const auto i = static_cast<axis::index_type>(std::floor(z * size()));
0341         min_ += i * (delta_ / size());
0342         delta_ = stop - min_;
0343         size_ -= i;
0344         return {0, -i};
0345       }
0346       // z is -infinity
0347       return {-1, 0};
0348     }
0349     // z either beyond range, infinite, or NaN
0350     if (z < std::numeric_limits<internal_value_type>::infinity()) {
0351       const auto i = static_cast<axis::index_type>(z * size());
0352       const auto n = i - size() + 1;
0353       delta_ /= size();
0354       delta_ *= size() + n;
0355       size_ += n;
0356       return {i, -n};
0357     }
0358     // z either infinite or NaN
0359     return {size(), 0};
0360   }
0361
0362   /// Return value for fractional index argument.
0363   value_type value(real_index_type i) const noexcept {
0364     auto z = i / size();
0365     if (!options_type::test(option::circular) && z < 0.0)
0366       z = -std::numeric_limits<internal_value_type>::infinity() * delta_;
0367     else if (options_type::test(option::circular) || z <= 1.0)
0368       z = (1.0 - z) * min_ + z * (min_ + delta_);
0369     else {
0370       z = std::numeric_limits<internal_value_type>::infinity() * delta_;
0371     }
0372     return static_cast<value_type>(this->inverse(z) * unit_type());
0373   }
0374
0375   /// Return bin for index argument.
0376   decltype(auto) bin(index_type idx) const noexcept {
0377     return interval_view<regular>(*this, idx);
0378   }
0379
0380   /// Returns the number of bins, without over- or underflow.
0381   index_type size() const noexcept { return size_; }
0382
0383   /// Returns the options.
0384   static constexpr unsigned options() noexcept { return options_type::value; }
0385
0386   template <class V, class T, class M, class O>
0387   bool operator==(const regular<V, T, M, O>& o) const noexcept {
0388     return detail::relaxed_equal{}(transform(), o.transform()) && size() == o.size() &&
0389            min_ == o.min_ && delta_ == o.delta_ &&
0390            detail::relaxed_equal{}(this->metadata(), o.metadata());
0391   }
0392   template <class V, class T, class M, class O>
0393   bool operator!=(const regular<V, T, M, O>& o) const noexcept {
0394     return !operator==(o);
0395   }
0396
0397   template <class Archive>
0398   void serialize(Archive& ar, unsigned /* version */) {
0399     ar& make_nvp("transform", static_cast<transform_type&>(*this));
0400     ar& make_nvp("size", size_);
0401     ar& make_nvp("meta", this->metadata());
0402     ar& make_nvp("min", min_);
0403     ar& make_nvp("delta", delta_);
0404   }
0405
0406 private:
0407   index_type size_{0};
0408   internal_value_type min_{0}, delta_{1};
0409
0410   template <class V, class T, class M, class O>
0411   friend class regular;
0412 };
0413
0414 #if __cpp_deduction_guides >= 201606
0415
0416 template <class T>
0417 regular(unsigned, T, T)
0418     -> regular<detail::convert_integer<T, double>, transform::id, null_type>;
0419
0420 template <class T, class M>
0421 regular(unsigned, T, T, M) -> regular<detail::convert_integer<T, double>, transform::id,
0422                                       detail::replace_cstring<std::decay_t<M>>>;
0423
0424 template <class T, class M, unsigned B>
0425 regular(unsigned, T, T, M, const option::bitset<B>&)
0426     -> regular<detail::convert_integer<T, double>, transform::id,
0427                detail::replace_cstring<std::decay_t<M>>, option::bitset<B>>;
0428
0429 template <class Tr, class T, class = detail::requires_transform<Tr, T>>
0430 regular(Tr, unsigned, T, T) -> regular<detail::convert_integer<T, double>, Tr, null_type>;
0431
0432 template <class Tr, class T, class M>
0433 regular(Tr, unsigned, T, T, M) -> regular<detail::convert_integer<T, double>, Tr,
0434                                           detail::replace_cstring<std::decay_t<M>>>;
0435
0436 template <class Tr, class T, class M, unsigned B>
0437 regular(Tr, unsigned, T, T, M, const option::bitset<B>&)
0438     -> regular<detail::convert_integer<T, double>, Tr,
0439                detail::replace_cstring<std::decay_t<M>>, option::bitset<B>>;
0440
0441 #endif
0442
0443 /// Regular axis with circular option already set.
0444 template <class Value = double, class MetaData = use_default, class Options = use_default>
0445 #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
0446 using circular = regular<Value, transform::id, MetaData,
0447                          decltype(detail::replace_default<Options, option::overflow_t>{} |
0448                                   option::circular)>;
0449 #else
0450 class circular;
0451 #endif
0452
0453 } // namespace axis
0454 } // namespace histogram
0455 } // namespace boost
0456
0457 #endif