regular.hpp 14 KB

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