variable.hpp 9.0 KB

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