boost/histogram/axis/variable.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_VARIABLE_HPP
   8 #define BOOST_HISTOGRAM_AXIS_VARIABLE_HPP
   9
  10 #include <algorithm>
  11 #include <boost/assert.hpp>
  12 #include <boost/core/nvp.hpp>
  13 #include <boost/histogram/axis/interval_view.hpp>
  14 #include <boost/histogram/axis/iterator.hpp>
  15 #include <boost/histogram/axis/metadata_base.hpp>
  16 #include <boost/histogram/axis/option.hpp>
  17 #include <boost/histogram/detail/convert_integer.hpp>
  18 #include <boost/histogram/detail/detect.hpp>
  19 #include <boost/histogram/detail/limits.hpp>
  20 #include <boost/histogram/detail/replace_type.hpp>
  21 #include <boost/histogram/fwd.hpp>
  22 #include <boost/throw_exception.hpp>
  23 #include <cmath>
  24 #include <limits>
  25 #include <memory>
  26 #include <stdexcept>
  27 #include <string>
  28 #include <type_traits>
  29 #include <utility>
  30 #include <vector>
  31
  32 namespace boost {
  33 namespace histogram {
  34 namespace axis {
  35
  36 /**
  37   Axis for non-equidistant bins on the real line.
  38
  39   Binning is a O(log(N)) operation. If speed matters and the problem domain
  40   allows it, prefer a regular axis, possibly with a transform.
  41
  42   @tparam Value input value type, must be floating point.
  43   @tparam MetaData type to store meta data.
  44   @tparam Options see boost::histogram::axis::option (all values allowed).
  45   @tparam Allocator allocator to use for dynamic memory management.
  46  */
  47 template <class Value, class MetaData, class Options, class Allocator>
  48 class variable : public iterator_mixin<variable<Value, MetaData, Options, Allocator>>,
  49                  public metadata_base<MetaData> {
  50   using value_type = Value;
  51   using metadata_type = typename metadata_base<MetaData>::metadata_type;
  52   using options_type =
  53       detail::replace_default<Options, decltype(option::underflow | option::overflow)>;
  54   using allocator_type = Allocator;
  55   using vector_type = std::vector<Value, allocator_type>;
  56
  57   static_assert(
  58       std::is_floating_point<value_type>::value,
  59       "current version of variable axis requires floating point type; "
  60       "if you need a variable axis with an integral type, please submit an issue");
  61
  62   static_assert(
  63       (!options_type::test(option::circular) && !options_type::test(option::growth)) ||
  64           (options_type::test(option::circular) ^ options_type::test(option::growth)),
  65       "circular and growth options are mutually exclusive");
  66
  67 public:
  68   constexpr variable() = default;
  69   explicit variable(allocator_type alloc) : vec_(alloc) {}
  70
  71   /** Construct from iterator range of bin edges.
  72    *
  73    * \param begin begin of edge sequence.
  74    * \param end   end of edge sequence.
  75    * \param meta  description of the axis.
  76    * \param alloc allocator instance to use.
  77    */
  78   template <class It, class = detail::requires_iterator<It>>
  79   variable(It begin, It end, metadata_type meta = {}, allocator_type alloc = {})
  80       : metadata_base<MetaData>(std::move(meta)), vec_(std::move(alloc)) {
  81     if (std::distance(begin, end) < 2)
  82       BOOST_THROW_EXCEPTION(std::invalid_argument("bins > 0 required"));
  83
  84     vec_.reserve(std::distance(begin, end));
  85     vec_.emplace_back(*begin++);
  86     bool strictly_ascending = true;
  87     while (begin != end) {
  88       if (*begin <= vec_.back()) strictly_ascending = false;
  89       vec_.emplace_back(*begin++);
  90     }
  91     if (!strictly_ascending)
  92       BOOST_THROW_EXCEPTION(
  93           std::invalid_argument("input sequence must be strictly ascending"));
  94   }
  95
  96   /** Construct variable axis from iterable range of bin edges.
  97    *
  98    * \param iterable iterable range of bin edges.
  99    * \param meta     description of the axis.
 100    * \param alloc    allocator instance to use.
 101    */
 102   template <class U, class = detail::requires_iterable<U>>
 103   variable(const U& iterable, metadata_type meta = {}, allocator_type alloc = {})
 104       : variable(std::begin(iterable), std::end(iterable), std::move(meta),
 105                  std::move(alloc)) {}
 106
 107   /** Construct variable axis from initializer list of bin edges.
 108    *
 109    * @param list  `std::initializer_list` of bin edges.
 110    * @param meta  description of the axis.
 111    * @param alloc allocator instance to use.
 112    */
 113   template <class U>
 114   variable(std::initializer_list<U> list, metadata_type meta = {},
 115            allocator_type alloc = {})
 116       : variable(list.begin(), list.end(), std::move(meta), std::move(alloc)) {}
 117
 118   /// Constructor used by algorithm::reduce to shrink and rebin (not for users).
 119   variable(const variable& src, index_type begin, index_type end, unsigned merge)
 120       : metadata_base<MetaData>(src), vec_(src.get_allocator()) {
 121     BOOST_ASSERT((end - begin) % merge == 0);
 122     if (options_type::test(option::circular) && !(begin == 0 && end == src.size()))
 123       BOOST_THROW_EXCEPTION(std::invalid_argument("cannot shrink circular axis"));
 124     vec_.reserve((end - begin) / merge);
 125     const auto beg = src.vec_.begin();
 126     for (index_type i = begin; i <= end; i += merge) vec_.emplace_back(*(beg + i));
 127   }
 128
 129   /// Return index for value argument.
 130   index_type index(value_type x) const noexcept {
 131     if (options_type::test(option::circular)) {
 132       const auto a = vec_[0];
 133       const auto b = vec_[size()];
 134       x -= std::floor((x - a) / (b - a)) * (b - a);
 135     }
 136     return static_cast<index_type>(std::upper_bound(vec_.begin(), vec_.end(), x) -
 137                                    vec_.begin() - 1);
 138   }
 139
 140   auto update(value_type x) noexcept {
 141     const auto i = index(x);
 142     if (std::isfinite(x)) {
 143       if (0 <= i) {
 144         if (i < size()) return std::make_pair(i, 0);
 145         const auto d = value(size()) - value(size() - 0.5);
 146         x = std::nextafter(x, (std::numeric_limits<value_type>::max)());
 147         x = (std::max)(x, vec_.back() + d);
 148         vec_.push_back(x);
 149         return std::make_pair(i, -1);
 150       }
 151       const auto d = value(0.5) - value(0);
 152       x = (std::min)(x, value(0) - d);
 153       vec_.insert(vec_.begin(), x);
 154       return std::make_pair(0, -i);
 155     }
 156     return std::make_pair(x < 0 ? -1 : size(), 0);
 157   }
 158
 159   /// Return value for fractional index argument.
 160   value_type value(real_index_type i) const noexcept {
 161     if (options_type::test(option::circular)) {
 162       auto shift = std::floor(i / size());
 163       i -= shift * size();
 164       double z;
 165       const auto k = static_cast<index_type>(std::modf(i, &z));
 166       const auto a = vec_[0];
 167       const auto b = vec_[size()];
 168       return (1.0 - z) * vec_[k] + z * vec_[k + 1] + shift * (b - a);
 169     }
 170     if (i < 0) return detail::lowest<value_type>();
 171     if (i == size()) return vec_.back();
 172     if (i > size()) return detail::highest<value_type>();
 173     const auto k = static_cast<index_type>(i); // precond: i >= 0
 174     const real_index_type z = i - k;
 175     return (1.0 - z) * vec_[k] + z * vec_[k + 1];
 176   }
 177
 178   /// Return bin for index argument.
 179   auto bin(index_type idx) const noexcept { return interval_view<variable>(*this, idx); }
 180
 181   /// Returns the number of bins, without over- or underflow.
 182   index_type size() const noexcept { return static_cast<index_type>(vec_.size()) - 1; }
 183
 184   /// Returns the options.
 185   static constexpr unsigned options() noexcept { return options_type::value; }
 186
 187   template <class V, class M, class O, class A>
 188   bool operator==(const variable<V, M, O, A>& o) const noexcept {
 189     const auto& a = vec_;
 190     const auto& b = o.vec_;
 191     return std::equal(a.begin(), a.end(), b.begin(), b.end()) &&
 192            metadata_base<MetaData>::operator==(o);
 193   }
 194
 195   template <class V, class M, class O, class A>
 196   bool operator!=(const variable<V, M, O, A>& o) const noexcept {
 197     return !operator==(o);
 198   }
 199
 200   /// Return allocator instance.
 201   auto get_allocator() const { return vec_.get_allocator(); }
 202
 203   template <class Archive>
 204   void serialize(Archive& ar, unsigned /* version */) {
 205     ar& make_nvp("seq", vec_);
 206     ar& make_nvp("meta", this->metadata());
 207   }
 208
 209 private:
 210   vector_type vec_;
 211
 212   template <class V, class M, class O, class A>
 213   friend class variable;
 214 };
 215
 216 #if __cpp_deduction_guides >= 201606
 217
 218 template <class T>
 219 variable(std::initializer_list<T>)
 220     ->variable<detail::convert_integer<T, double>, null_type>;
 221
 222 template <class T, class M>
 223 variable(std::initializer_list<T>, M)
 224     ->variable<detail::convert_integer<T, double>,
 225                detail::replace_type<std::decay_t<M>, const char*, std::string>>;
 226
 227 template <class Iterable, class = detail::requires_iterable<Iterable>>
 228 variable(Iterable)
 229     ->variable<
 230         detail::convert_integer<
 231             std::decay_t<decltype(*std::begin(std::declval<Iterable&>()))>, double>,
 232         null_type>;
 233
 234 template <class Iterable, class M>
 235 variable(Iterable, M)
 236     ->variable<
 237         detail::convert_integer<
 238             std::decay_t<decltype(*std::begin(std::declval<Iterable&>()))>, double>,
 239         detail::replace_type<std::decay_t<M>, const char*, std::string>>;
 240
 241 #endif
 242
 243 } // namespace axis
 244 } // namespace histogram
 245 } // namespace boost
 246
 247 #endif