regular.hpp 16 KB

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