regular.hpp 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431
  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. /**
  136. Axis for equidistant intervals on the real line.
  137. The most common binning strategy. Very fast. Binning is a O(1) operation.
  138. @tparam Value input value type, must be floating point.
  139. @tparam Transform builtin or user-defined transform type.
  140. @tparam MetaData type to store meta data.
  141. @tparam Options see boost::histogram::axis::option (all values allowed).
  142. */
  143. template <class Value, class Transform, class MetaData, class Options>
  144. class regular : public iterator_mixin<regular<Value, Transform, MetaData, Options>>,
  145. protected detail::replace_default<Transform, transform::id>,
  146. public metadata_base_t<MetaData> {
  147. // these must be private, so that they are not automatically inherited
  148. using value_type = Value;
  149. using transform_type = detail::replace_default<Transform, transform::id>;
  150. using metadata_base = metadata_base_t<MetaData>;
  151. using metadata_type = typename metadata_base::metadata_type;
  152. using options_type =
  153. detail::replace_default<Options, decltype(option::underflow | option::overflow)>;
  154. static_assert(std::is_nothrow_move_constructible<transform_type>::value,
  155. "transform must be no-throw move constructible");
  156. static_assert(std::is_nothrow_move_assignable<transform_type>::value,
  157. "transform must be no-throw move assignable");
  158. using unit_type = detail::get_unit_type<value_type>;
  159. using internal_value_type = detail::get_scale_type<value_type>;
  160. static_assert(std::is_floating_point<internal_value_type>::value,
  161. "regular axis requires floating point type");
  162. static_assert(
  163. (!options_type::test(option::circular) && !options_type::test(option::growth)) ||
  164. (options_type::test(option::circular) ^ options_type::test(option::growth)),
  165. "circular and growth options are mutually exclusive");
  166. public:
  167. constexpr regular() = default;
  168. /** Construct n bins over real transformed range [start, stop).
  169. *
  170. * @param trans transform instance to use.
  171. * @param n number of bins.
  172. * @param start low edge of first bin.
  173. * @param stop high edge of last bin.
  174. * @param meta description of the axis (optional).
  175. */
  176. regular(transform_type trans, unsigned n, value_type start, value_type stop,
  177. metadata_type meta = {})
  178. : transform_type(std::move(trans))
  179. , metadata_base(std::move(meta))
  180. , size_(static_cast<index_type>(n))
  181. , min_(this->forward(detail::get_scale(start)))
  182. , delta_(this->forward(detail::get_scale(stop)) - min_) {
  183. if (size() == 0) BOOST_THROW_EXCEPTION(std::invalid_argument("bins > 0 required"));
  184. if (!std::isfinite(min_) || !std::isfinite(delta_))
  185. BOOST_THROW_EXCEPTION(
  186. std::invalid_argument("forward transform of start or stop invalid"));
  187. if (delta_ == 0)
  188. BOOST_THROW_EXCEPTION(std::invalid_argument("range of axis is zero"));
  189. }
  190. /** Construct n bins over real range [start, stop).
  191. *
  192. * @param n number of bins.
  193. * @param start low edge of first bin.
  194. * @param stop high edge of last bin.
  195. * @param meta description of the axis (optional).
  196. */
  197. regular(unsigned n, value_type start, value_type stop, metadata_type meta = {})
  198. : regular({}, n, start, stop, std::move(meta)) {}
  199. /** Construct bins with the given step size over real transformed range
  200. * [start, stop).
  201. *
  202. * @param trans transform instance to use.
  203. * @param step width of a single bin.
  204. * @param start low edge of first bin.
  205. * @param stop upper limit of high edge of last bin (see below).
  206. * @param meta description of the axis (optional).
  207. *
  208. * The axis computes the number of bins as n = abs(stop - start) / step,
  209. * rounded down. This means that stop is an upper limit to the actual value
  210. * (start + n * step).
  211. */
  212. template <class T>
  213. regular(transform_type trans, step_type<T> step, value_type start, value_type stop,
  214. metadata_type meta = {})
  215. : regular(trans, static_cast<index_type>(std::abs(stop - start) / step.value),
  216. start,
  217. start + static_cast<index_type>(std::abs(stop - start) / step.value) *
  218. step.value,
  219. std::move(meta)) {}
  220. /** Construct bins with the given step size over real range [start, stop).
  221. *
  222. * @param step width of a single bin.
  223. * @param start low edge of first bin.
  224. * @param stop upper limit of high edge of last bin (see below).
  225. * @param meta description of the axis (optional).
  226. *
  227. * The axis computes the number of bins as n = abs(stop - start) / step,
  228. * rounded down. This means that stop is an upper limit to the actual value
  229. * (start + n * step).
  230. */
  231. template <class T>
  232. regular(step_type<T> step, value_type start, value_type stop, metadata_type meta = {})
  233. : regular({}, step, start, stop, std::move(meta)) {}
  234. /// Constructor used by algorithm::reduce to shrink and rebin (not for users).
  235. regular(const regular& src, index_type begin, index_type end, unsigned merge)
  236. : regular(src.transform(), (end - begin) / merge, src.value(begin), src.value(end),
  237. src.metadata()) {
  238. assert((end - begin) % merge == 0);
  239. if (options_type::test(option::circular) && !(begin == 0 && end == src.size()))
  240. BOOST_THROW_EXCEPTION(std::invalid_argument("cannot shrink circular axis"));
  241. }
  242. /// Return instance of the transform type.
  243. const transform_type& transform() const noexcept { return *this; }
  244. /// Return index for value argument.
  245. index_type index(value_type x) const noexcept {
  246. // Runs in hot loop, please measure impact of changes
  247. auto z = (this->forward(x / unit_type{}) - min_) / delta_;
  248. if (options_type::test(option::circular)) {
  249. if (std::isfinite(z)) {
  250. z -= std::floor(z);
  251. return static_cast<index_type>(z * size());
  252. }
  253. } else {
  254. if (z < 1) {
  255. if (z >= 0)
  256. return static_cast<index_type>(z * size());
  257. else
  258. return -1;
  259. }
  260. }
  261. return size(); // also returned if x is NaN
  262. }
  263. /// Returns index and shift (if axis has grown) for the passed argument.
  264. std::pair<index_type, index_type> update(value_type x) noexcept {
  265. assert(options_type::test(option::growth));
  266. const auto z = (this->forward(x / unit_type{}) - min_) / delta_;
  267. if (z < 1) { // don't use i here!
  268. if (z >= 0) {
  269. const auto i = static_cast<axis::index_type>(z * size());
  270. return {i, 0};
  271. }
  272. if (z != -std::numeric_limits<internal_value_type>::infinity()) {
  273. const auto stop = min_ + delta_;
  274. const auto i = static_cast<axis::index_type>(std::floor(z * size()));
  275. min_ += i * (delta_ / size());
  276. delta_ = stop - min_;
  277. size_ -= i;
  278. return {0, -i};
  279. }
  280. // z is -infinity
  281. return {-1, 0};
  282. }
  283. // z either beyond range, infinite, or NaN
  284. if (z < std::numeric_limits<internal_value_type>::infinity()) {
  285. const auto i = static_cast<axis::index_type>(z * size());
  286. const auto n = i - size() + 1;
  287. delta_ /= size();
  288. delta_ *= size() + n;
  289. size_ += n;
  290. return {i, -n};
  291. }
  292. // z either infinite or NaN
  293. return {size(), 0};
  294. }
  295. /// Return value for fractional index argument.
  296. value_type value(real_index_type i) const noexcept {
  297. auto z = i / size();
  298. if (!options_type::test(option::circular) && z < 0.0)
  299. z = -std::numeric_limits<internal_value_type>::infinity() * delta_;
  300. else if (options_type::test(option::circular) || z <= 1.0)
  301. z = (1.0 - z) * min_ + z * (min_ + delta_);
  302. else {
  303. z = std::numeric_limits<internal_value_type>::infinity() * delta_;
  304. }
  305. return static_cast<value_type>(this->inverse(z) * unit_type());
  306. }
  307. /// Return bin for index argument.
  308. decltype(auto) bin(index_type idx) const noexcept {
  309. return interval_view<regular>(*this, idx);
  310. }
  311. /// Returns the number of bins, without over- or underflow.
  312. index_type size() const noexcept { return size_; }
  313. /// Returns the options.
  314. static constexpr unsigned options() noexcept { return options_type::value; }
  315. template <class V, class T, class M, class O>
  316. bool operator==(const regular<V, T, M, O>& o) const noexcept {
  317. return detail::relaxed_equal{}(transform(), o.transform()) && size() == o.size() &&
  318. min_ == o.min_ && delta_ == o.delta_ &&
  319. detail::relaxed_equal{}(this->metadata(), o.metadata());
  320. }
  321. template <class V, class T, class M, class O>
  322. bool operator!=(const regular<V, T, M, O>& o) const noexcept {
  323. return !operator==(o);
  324. }
  325. template <class Archive>
  326. void serialize(Archive& ar, unsigned /* version */) {
  327. ar& make_nvp("transform", static_cast<transform_type&>(*this));
  328. ar& make_nvp("size", size_);
  329. ar& make_nvp("meta", this->metadata());
  330. ar& make_nvp("min", min_);
  331. ar& make_nvp("delta", delta_);
  332. }
  333. private:
  334. index_type size_{0};
  335. internal_value_type min_{0}, delta_{1};
  336. template <class V, class T, class M, class O>
  337. friend class regular;
  338. };
  339. #if __cpp_deduction_guides >= 201606
  340. template <class T>
  341. regular(unsigned, T, T)
  342. ->regular<detail::convert_integer<T, double>, transform::id, null_type>;
  343. template <class T, class M>
  344. regular(unsigned, T, T, M)
  345. ->regular<detail::convert_integer<T, double>, transform::id,
  346. detail::replace_cstring<std::decay_t<M>>>;
  347. template <class Tr, class T, class = detail::requires_transform<Tr, T>>
  348. regular(Tr, unsigned, T, T)->regular<detail::convert_integer<T, double>, Tr, null_type>;
  349. template <class Tr, class T, class M>
  350. regular(Tr, unsigned, T, T, M)
  351. ->regular<detail::convert_integer<T, double>, Tr,
  352. detail::replace_cstring<std::decay_t<M>>>;
  353. #endif
  354. /// Regular axis with circular option already set.
  355. template <class Value = double, class MetaData = use_default, class Options = use_default>
  356. #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
  357. using circular = regular<Value, transform::id, MetaData,
  358. decltype(detail::replace_default<Options, option::overflow_t>{} |
  359. option::circular)>;
  360. #else
  361. class circular;
  362. #endif
  363. } // namespace axis
  364. } // namespace histogram
  365. } // namespace boost
  366. #endif