histogram/axis/variable.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_VARIABLE_HPP
0008 #define BOOST_HISTOGRAM_AXIS_VARIABLE_HPP
0009
0010 #include <algorithm>
0011 #include <boost/core/nvp.hpp>
0012 #include <boost/histogram/axis/interval_view.hpp>
0013 #include <boost/histogram/axis/iterator.hpp>
0014 #include <boost/histogram/axis/metadata_base.hpp>
0015 #include <boost/histogram/axis/option.hpp>
0016 #include <boost/histogram/detail/convert_integer.hpp>
0017 #include <boost/histogram/detail/detect.hpp>
0018 #include <boost/histogram/detail/limits.hpp>
0019 #include <boost/histogram/detail/relaxed_equal.hpp>
0020 #include <boost/histogram/detail/replace_type.hpp>
0021 #include <boost/histogram/fwd.hpp>
0022 #include <boost/throw_exception.hpp>
0023 #include <cassert>
0024 #include <cmath>
0025 #include <limits>
0026 #include <memory>
0027 #include <stdexcept>
0028 #include <string>
0029 #include <type_traits>
0030 #include <utility>
0031 #include <vector>
0032
0033 namespace boost {
0034 namespace histogram {
0035 namespace axis {
0036
0037 /** Axis for non-equidistant bins on the real line.
0038
0039   Binning is a O(log(N)) operation. If speed matters and the problem domain
0040   allows it, prefer a regular axis, possibly with a transform.
0041
0042   If the axis has an overflow bin (the default), a value on the upper edge of the last
0043   bin is put in the overflow bin. The axis range represents a semi-open interval.
0044
0045   If the overflow bin is deactivated, then a value on the upper edge of the last bin is
0046   still counted towards the last bin. The axis range represents a closed interval. This
0047   is the desired behavior for random numbers drawn from a bounded interval, which is
0048   usually closed.
0049
0050   @tparam Value     input value type, must be floating point.
0051   @tparam MetaData  type to store meta data.
0052   @tparam Options   see boost::histogram::axis::option.
0053   @tparam Allocator allocator to use for dynamic memory management.
0054 */
0055 template <class Value, class MetaData, class Options, class Allocator>
0056 class variable : public iterator_mixin<variable<Value, MetaData, Options, Allocator>>,
0057                  public metadata_base_t<MetaData> {
0058   // these must be private, so that they are not automatically inherited
0059   using value_type = Value;
0060   using metadata_base = metadata_base_t<MetaData>;
0061   using metadata_type = typename metadata_base::metadata_type;
0062   using options_type =
0063       detail::replace_default<Options, decltype(option::underflow | option::overflow)>;
0064   using allocator_type = Allocator;
0065   using vector_type = std::vector<Value, allocator_type>;
0066
0067 public:
0068   constexpr variable() = default;
0069   explicit variable(allocator_type alloc) : vec_(alloc) {}
0070
0071   /** Construct from forward iterator range of bin edges.
0072
0073     @param begin   begin of edge sequence.
0074     @param end     end of edge sequence.
0075     @param meta    description of the axis (optional).
0076     @param options see boost::histogram::axis::option (optional).
0077     @param alloc   allocator instance to use (optional).
0078
0079     The constructor throws `std::invalid_argument` if iterator range is invalid, if less
0080     than two edges are provided or if bin edges are not in ascending order.
0081
0082     The arguments meta and alloc are passed by value. If you move either of them into the
0083     axis and the constructor throws, their values are lost. Do not move if you cannot
0084     guarantee that the bin description is not valid.
0085    */
0086   template <class It, class = detail::requires_iterator<It>>
0087   variable(It begin, It end, metadata_type meta = {}, options_type options = {},
0088            allocator_type alloc = {})
0089       : metadata_base(std::move(meta)), vec_(std::move(alloc)) {
0090     // static_asserts were moved here from class scope to satisfy deduction in gcc>=11
0091     static_assert(
0092         std::is_floating_point<value_type>::value,
0093         "current version of variable axis requires floating point type; "
0094         "if you need a variable axis with an integral type, please submit an issue");
0095     static_assert((!options.test(option::circular) && !options.test(option::growth)) ||
0096                       (options.test(option::circular) ^ options.test(option::growth)),
0097                   "circular and growth options are mutually exclusive");
0098
0099     const auto n = std::distance(begin, end);
0100     if (n < 0)
0101       BOOST_THROW_EXCEPTION(
0102           std::invalid_argument("end must be reachable by incrementing begin"));
0103
0104     if (n < 2) BOOST_THROW_EXCEPTION(std::invalid_argument("bins > 1 required"));
0105
0106     vec_.reserve(n);
0107     vec_.emplace_back(*begin++);
0108     bool strictly_ascending = true;
0109     for (; begin != end; ++begin) {
0110       strictly_ascending &= vec_.back() < *begin;
0111       vec_.emplace_back(*begin);
0112     }
0113
0114     if (!strictly_ascending)
0115       BOOST_THROW_EXCEPTION(
0116           std::invalid_argument("input sequence must be strictly ascending"));
0117   }
0118
0119   // kept for backward compatibility; requires_allocator is a workaround for deduction
0120   // guides in gcc>=11
0121   template <class It, class A, class = detail::requires_iterator<It>,
0122             class = detail::requires_allocator<A>>
0123   variable(It begin, It end, metadata_type meta, A alloc)
0124       : variable(begin, end, std::move(meta), {}, std::move(alloc)) {}
0125
0126   /** Construct variable axis from iterable range of bin edges.
0127
0128      @param iterable iterable range of bin edges.
0129      @param meta     description of the axis (optional).
0130      @param options  see boost::histogram::axis::option (optional).
0131      @param alloc    allocator instance to use (optional).
0132    */
0133   template <class U, class = detail::requires_iterable<U>>
0134   variable(const U& iterable, metadata_type meta = {}, options_type options = {},
0135            allocator_type alloc = {})
0136       : variable(std::begin(iterable), std::end(iterable), std::move(meta), options,
0137                  std::move(alloc)) {}
0138
0139   // kept for backward compatibility; requires_allocator is a workaround for deduction
0140   // guides in gcc>=11
0141   template <class U, class A, class = detail::requires_iterable<U>,
0142             class = detail::requires_allocator<A>>
0143   variable(const U& iterable, metadata_type meta, A alloc)
0144       : variable(std::begin(iterable), std::end(iterable), std::move(meta), {},
0145                  std::move(alloc)) {}
0146
0147   /** Construct variable axis from initializer list of bin edges.
0148
0149      @param list     `std::initializer_list` of bin edges.
0150      @param meta     description of the axis (optional).
0151      @param options  see boost::histogram::axis::option (optional).
0152      @param alloc    allocator instance to use (optional).
0153    */
0154   template <class U>
0155   variable(std::initializer_list<U> list, metadata_type meta = {},
0156            options_type options = {}, allocator_type alloc = {})
0157       : variable(list.begin(), list.end(), std::move(meta), options, std::move(alloc)) {}
0158
0159   // kept for backward compatibility; requires_allocator is a workaround for deduction
0160   // guides in gcc>=11
0161   template <class U, class A, class = detail::requires_allocator<A>>
0162   variable(std::initializer_list<U> list, metadata_type meta, A alloc)
0163       : variable(list.begin(), list.end(), std::move(meta), {}, std::move(alloc)) {}
0164
0165   /// Constructor used by algorithm::reduce to shrink and rebin (not for users).
0166   variable(const variable& src, index_type begin, index_type end, unsigned merge)
0167       : metadata_base(src), vec_(src.get_allocator()) {
0168     assert((end - begin) % merge == 0);
0169     if (options_type::test(option::circular) && !(begin == 0 && end == src.size()))
0170       BOOST_THROW_EXCEPTION(std::invalid_argument("cannot shrink circular axis"));
0171     vec_.reserve((end - begin) / merge);
0172     const auto beg = src.vec_.begin();
0173     for (index_type i = begin; i <= end; i += merge) vec_.emplace_back(*(beg + i));
0174   }
0175
0176   /// Return index for value argument.
0177   index_type index(value_type x) const noexcept {
0178     if (options_type::test(option::circular)) {
0179       const auto a = vec_[0];
0180       const auto b = vec_[size()];
0181       x -= std::floor((x - a) / (b - a)) * (b - a);
0182     }
0183     // upper edge of last bin is inclusive if overflow bin is not present
0184     if (!options_type::test(option::overflow) && x == vec_.back()) return size() - 1;
0185     return static_cast<index_type>(std::upper_bound(vec_.begin(), vec_.end(), x) -
0186                                    vec_.begin() - 1);
0187   }
0188
0189   std::pair<index_type, index_type> update(value_type x) noexcept {
0190     const auto i = index(x);
0191     if (std::isfinite(x)) {
0192       if (0 <= i) {
0193         if (i < size()) return std::make_pair(i, 0);
0194         const auto d = value(size()) - value(size() - 0.5);
0195         x = std::nextafter(x, (std::numeric_limits<value_type>::max)());
0196         x = (std::max)(x, vec_.back() + d);
0197         vec_.push_back(x);
0198         return {i, -1};
0199       }
0200       const auto d = value(0.5) - value(0);
0201       x = (std::min)(x, value(0) - d);
0202       vec_.insert(vec_.begin(), x);
0203       return {0, -i};
0204     }
0205     return {x < 0 ? -1 : size(), 0};
0206   }
0207
0208   /// Return value for fractional index argument.
0209   value_type value(real_index_type i) const noexcept {
0210     if (options_type::test(option::circular)) {
0211       auto shift = std::floor(i / size());
0212       i -= shift * size();
0213       double z;
0214       const auto k = static_cast<index_type>(std::modf(i, &z));
0215       const auto a = vec_[0];
0216       const auto b = vec_[size()];
0217       return (1.0 - z) * vec_[k] + z * vec_[k + 1] + shift * (b - a);
0218     }
0219     if (i < 0) return detail::lowest<value_type>();
0220     if (i == size()) return vec_.back();
0221     if (i > size()) return detail::highest<value_type>();
0222     const auto k = static_cast<index_type>(i); // precond: i >= 0
0223     const real_index_type z = i - k;
0224     // check z == 0 needed to avoid returning nan when vec_[k + 1] is infinity
0225     return (1.0 - z) * vec_[k] + (z == 0 ? 0 : z * vec_[k + 1]);
0226   }
0227
0228   /// Return bin for index argument.
0229   auto bin(index_type idx) const noexcept { return interval_view<variable>(*this, idx); }
0230
0231   /// Returns the number of bins, without over- or underflow.
0232   index_type size() const noexcept { return static_cast<index_type>(vec_.size()) - 1; }
0233
0234   /// Returns the options.
0235   static constexpr unsigned options() noexcept { return options_type::value; }
0236
0237   template <class V, class M, class O, class A>
0238   bool operator==(const variable<V, M, O, A>& o) const noexcept {
0239     const auto& a = vec_;
0240     const auto& b = o.vec_;
0241     return std::equal(a.begin(), a.end(), b.begin(), b.end()) &&
0242            detail::relaxed_equal{}(this->metadata(), o.metadata());
0243   }
0244
0245   template <class V, class M, class O, class A>
0246   bool operator!=(const variable<V, M, O, A>& o) const noexcept {
0247     return !operator==(o);
0248   }
0249
0250   /// Return allocator instance.
0251   auto get_allocator() const { return vec_.get_allocator(); }
0252
0253   template <class Archive>
0254   void serialize(Archive& ar, unsigned /* version */) {
0255     ar& make_nvp("seq", vec_);
0256     ar& make_nvp("meta", this->metadata());
0257   }
0258
0259 private:
0260   vector_type vec_;
0261
0262   template <class V, class M, class O, class A>
0263   friend class variable;
0264 };
0265
0266 #if __cpp_deduction_guides >= 201606
0267
0268 template <class T>
0269 variable(std::initializer_list<T>)
0270     -> variable<detail::convert_integer<T, double>, null_type>;
0271
0272 template <class T, class M>
0273 variable(std::initializer_list<T>, M)
0274     -> variable<detail::convert_integer<T, double>,
0275                 detail::replace_type<std::decay_t<M>, const char*, std::string>>;
0276
0277 template <class T, class M, unsigned B>
0278 variable(std::initializer_list<T>, M, const option::bitset<B>&)
0279     -> variable<detail::convert_integer<T, double>,
0280                 detail::replace_type<std::decay_t<M>, const char*, std::string>,
0281                 option::bitset<B>>;
0282
0283 template <class Iterable, class = detail::requires_iterable<Iterable>>
0284 variable(Iterable) -> variable<
0285     detail::convert_integer<
0286         std::decay_t<decltype(*std::begin(std::declval<Iterable&>()))>, double>,
0287     null_type>;
0288
0289 template <class Iterable, class M>
0290 variable(Iterable, M) -> variable<
0291     detail::convert_integer<
0292         std::decay_t<decltype(*std::begin(std::declval<Iterable&>()))>, double>,
0293     detail::replace_type<std::decay_t<M>, const char*, std::string>>;
0294
0295 template <class Iterable, class M, unsigned B>
0296 variable(Iterable, M, const option::bitset<B>&) -> variable<
0297     detail::convert_integer<
0298         std::decay_t<decltype(*std::begin(std::declval<Iterable&>()))>, double>,
0299     detail::replace_type<std::decay_t<M>, const char*, std::string>, option::bitset<B>>;
0300
0301 #endif
0302
0303 } // namespace axis
0304 } // namespace histogram
0305 } // namespace boost
0306
0307 #endif