regular.hpp 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360
  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/histogram/axis/interval_view.hpp>
  10. #include <boost/histogram/axis/iterator.hpp>
  11. #include <boost/histogram/axis/option.hpp>
  12. #include <boost/histogram/detail/compressed_pair.hpp>
  13. #include <boost/histogram/detail/meta.hpp>
  14. #include <boost/histogram/fwd.hpp>
  15. #include <boost/mp11/utility.hpp>
  16. #include <boost/throw_exception.hpp>
  17. #include <cmath>
  18. #include <limits>
  19. #include <stdexcept>
  20. #include <type_traits>
  21. #include <utility>
  22. namespace boost {
  23. namespace histogram {
  24. namespace axis {
  25. namespace transform {
  26. /// Identity transform for equidistant bins.
  27. struct id {
  28. /// Pass-through.
  29. template <typename T>
  30. static T forward(T&& x) noexcept {
  31. return std::forward<T>(x);
  32. }
  33. /// Pass-through.
  34. template <typename T>
  35. static T inverse(T&& x) noexcept {
  36. return std::forward<T>(x);
  37. }
  38. };
  39. /// Log transform for equidistant bins in log-space.
  40. struct log {
  41. /// Returns log(x) of external value x.
  42. template <typename T>
  43. static T forward(T x) {
  44. return std::log(x);
  45. }
  46. /// Returns exp(x) for internal value x.
  47. template <typename T>
  48. static T inverse(T x) {
  49. return std::exp(x);
  50. }
  51. };
  52. /// Sqrt transform for equidistant bins in sqrt-space.
  53. struct sqrt {
  54. /// Returns sqrt(x) of external value x.
  55. template <typename T>
  56. static T forward(T x) {
  57. return std::sqrt(x);
  58. }
  59. /// Returns x^2 of internal value x.
  60. template <typename T>
  61. static T inverse(T x) {
  62. return x * x;
  63. }
  64. };
  65. /// Pow transform for equidistant bins in pow-space.
  66. struct pow {
  67. double power = 1; /**< power index */
  68. /// Make transform with index p.
  69. explicit pow(double p) : power(p) {}
  70. pow() = default;
  71. /// Returns pow(x, power) of external value x.
  72. template <typename T>
  73. auto forward(T x) const {
  74. return std::pow(x, power);
  75. }
  76. /// Returns pow(x, 1/power) of external value x.
  77. template <typename T>
  78. auto inverse(T x) const {
  79. return std::pow(x, 1.0 / power);
  80. }
  81. bool operator==(const pow& o) const noexcept { return power == o.power; }
  82. };
  83. } // namespace transform
  84. #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
  85. // Type envelope to mark value as step size
  86. template <typename T>
  87. struct step_type {
  88. T value;
  89. };
  90. #endif
  91. /**
  92. Helper function to mark argument as step size.
  93. */
  94. template <typename T>
  95. auto step(T&& t) {
  96. return step_type<T&&>{std::forward<T>(t)};
  97. }
  98. /**
  99. Axis for equidistant intervals on the real line.
  100. The most common binning strategy. Very fast. Binning is a O(1) operation.
  101. @tparam Value input value type, must be floating point.
  102. @tparam Transform builtin or user-defined transform type.
  103. @tparam MetaData type to store meta data.
  104. @tparam Options see boost::histogram::axis::option (all values allowed).
  105. */
  106. template <class Value, class Transform, class MetaData, class Options>
  107. class regular : public iterator_mixin<regular<Value, Transform, MetaData, Options>>,
  108. protected detail::replace_default<Transform, transform::id> {
  109. using value_type = Value;
  110. using transform_type = detail::replace_default<Transform, transform::id>;
  111. using metadata_type = detail::replace_default<MetaData, std::string>;
  112. using options_type =
  113. detail::replace_default<Options, decltype(option::underflow | option::overflow)>;
  114. using unit_type = detail::get_unit_type<value_type>;
  115. using internal_value_type = detail::get_scale_type<value_type>;
  116. static_assert(std::is_floating_point<internal_value_type>::value,
  117. "variable axis requires floating point type");
  118. public:
  119. constexpr regular() = default;
  120. /** Construct n bins over real transformed range [start, stop).
  121. *
  122. * @param trans transform instance to use.
  123. * @param n number of bins.
  124. * @param start low edge of first bin.
  125. * @param stop high edge of last bin.
  126. * @param meta description of the axis (optional).
  127. */
  128. regular(transform_type trans, unsigned n, value_type start, value_type stop,
  129. metadata_type meta = {})
  130. : transform_type(std::move(trans))
  131. , size_meta_(static_cast<index_type>(n), std::move(meta))
  132. , min_(this->forward(detail::get_scale(start)))
  133. , delta_(this->forward(detail::get_scale(stop)) - min_) {
  134. if (size() == 0) BOOST_THROW_EXCEPTION(std::invalid_argument("bins > 0 required"));
  135. if (!std::isfinite(min_) || !std::isfinite(delta_))
  136. BOOST_THROW_EXCEPTION(
  137. std::invalid_argument("forward transform of start or stop invalid"));
  138. if (delta_ == 0)
  139. BOOST_THROW_EXCEPTION(std::invalid_argument("range of axis is zero"));
  140. }
  141. /** Construct n bins over real range [start, stop).
  142. *
  143. * @param n number of bins.
  144. * @param start low edge of first bin.
  145. * @param stop high edge of last bin.
  146. * @param meta description of the axis (optional).
  147. */
  148. regular(unsigned n, value_type start, value_type stop, metadata_type meta = {})
  149. : regular({}, n, start, stop, std::move(meta)) {}
  150. /** Construct bins with the given step size over real transformed range
  151. * [start, stop).
  152. *
  153. * @param trans transform instance to use.
  154. * @param step width of a single bin.
  155. * @param start low edge of first bin.
  156. * @param stop upper limit of high edge of last bin (see below).
  157. * @param meta description of the axis (optional).
  158. *
  159. * The axis computes the number of bins as n = abs(stop - start) / step,
  160. * rounded down. This means that stop is an upper limit to the actual value
  161. * (start + n * step).
  162. */
  163. template <class T>
  164. regular(transform_type trans, const step_type<T>& step, value_type start,
  165. value_type stop, metadata_type meta = {})
  166. : regular(trans, static_cast<index_type>(std::abs(stop - start) / step.value),
  167. start,
  168. start + static_cast<index_type>(std::abs(stop - start) / step.value) *
  169. step.value,
  170. std::move(meta)) {}
  171. /** Construct bins with the given step size over real range [start, stop).
  172. *
  173. * @param step width of a single bin.
  174. * @param start low edge of first bin.
  175. * @param stop upper limit of high edge of last bin (see below).
  176. * @param meta description of the axis (optional).
  177. *
  178. * The axis computes the number of bins as n = abs(stop - start) / step,
  179. * rounded down. This means that stop is an upper limit to the actual value
  180. * (start + n * step).
  181. */
  182. template <class T>
  183. regular(const step_type<T>& step, value_type start, value_type stop,
  184. metadata_type meta = {})
  185. : regular({}, step, start, stop, std::move(meta)) {}
  186. /// Constructor used by algorithm::reduce to shrink and rebin (not for users).
  187. regular(const regular& src, index_type begin, index_type end, unsigned merge)
  188. : regular(src.transform(), (end - begin) / merge, src.value(begin), src.value(end),
  189. src.metadata()) {
  190. BOOST_ASSERT((end - begin) % merge == 0);
  191. if (options_type::test(option::circular) && !(begin == 0 && end == src.size()))
  192. BOOST_THROW_EXCEPTION(std::invalid_argument("cannot shrink circular axis"));
  193. }
  194. /// Return instance of the transform type.
  195. const transform_type& transform() const noexcept { return *this; }
  196. /// Return index for value argument.
  197. index_type index(value_type x) const noexcept {
  198. // Runs in hot loop, please measure impact of changes
  199. auto z = (this->forward(x / unit_type{}) - min_) / delta_;
  200. if (options_type::test(option::circular)) {
  201. if (std::isfinite(z)) {
  202. z -= std::floor(z);
  203. return static_cast<index_type>(z * size());
  204. }
  205. } else {
  206. if (z < 1) {
  207. if (z >= 0)
  208. return static_cast<index_type>(z * size());
  209. else
  210. return -1;
  211. }
  212. }
  213. return size(); // also returned if x is NaN
  214. }
  215. /// Returns index and shift (if axis has grown) for the passed argument.
  216. auto update(value_type x) noexcept {
  217. BOOST_ASSERT(options_type::test(option::growth));
  218. const auto z = (this->forward(x / unit_type{}) - min_) / delta_;
  219. if (z < 1) { // don't use i here!
  220. if (z >= 0) {
  221. const auto i = static_cast<axis::index_type>(z * size());
  222. return std::make_pair(i, 0);
  223. }
  224. if (z != -std::numeric_limits<internal_value_type>::infinity()) {
  225. const auto stop = min_ + delta_;
  226. const auto i = static_cast<axis::index_type>(std::floor(z * size()));
  227. min_ += i * (delta_ / size());
  228. delta_ = stop - min_;
  229. size_meta_.first() -= i;
  230. return std::make_pair(0, -i);
  231. }
  232. // z is -infinity
  233. return std::make_pair(-1, 0);
  234. }
  235. // z either beyond range, infinite, or NaN
  236. if (z < std::numeric_limits<internal_value_type>::infinity()) {
  237. const auto i = static_cast<axis::index_type>(z * size());
  238. const auto n = i - size() + 1;
  239. delta_ /= size();
  240. delta_ *= size() + n;
  241. size_meta_.first() += n;
  242. return std::make_pair(i, -n);
  243. }
  244. // z either infinite or NaN
  245. return std::make_pair(size(), 0);
  246. }
  247. /// Return value for fractional index argument.
  248. value_type value(real_index_type i) const noexcept {
  249. auto z = i / size();
  250. if (!options_type::test(option::circular) && z < 0.0)
  251. z = -std::numeric_limits<internal_value_type>::infinity() * delta_;
  252. else if (options_type::test(option::circular) || z <= 1.0)
  253. z = (1.0 - z) * min_ + z * (min_ + delta_);
  254. else {
  255. z = std::numeric_limits<internal_value_type>::infinity() * delta_;
  256. }
  257. return static_cast<value_type>(this->inverse(z) * unit_type());
  258. }
  259. /// Return bin for index argument.
  260. decltype(auto) bin(index_type idx) const noexcept {
  261. return interval_view<regular>(*this, idx);
  262. }
  263. /// Returns the number of bins, without over- or underflow.
  264. index_type size() const noexcept { return size_meta_.first(); }
  265. /// Returns the options.
  266. static constexpr unsigned options() noexcept { return options_type::value; }
  267. /// Returns reference to metadata.
  268. metadata_type& metadata() noexcept { return size_meta_.second(); }
  269. /// Returns reference to const metadata.
  270. const metadata_type& metadata() const noexcept { return size_meta_.second(); }
  271. template <class V, class T, class M, class O>
  272. bool operator==(const regular<V, T, M, O>& o) const noexcept {
  273. return detail::relaxed_equal(transform(), o.transform()) && size() == o.size() &&
  274. detail::relaxed_equal(metadata(), o.metadata()) && min_ == o.min_ &&
  275. delta_ == o.delta_;
  276. }
  277. template <class V, class T, class M, class O>
  278. bool operator!=(const regular<V, T, M, O>& o) const noexcept {
  279. return !operator==(o);
  280. }
  281. template <class Archive>
  282. void serialize(Archive&, unsigned);
  283. private:
  284. detail::compressed_pair<index_type, metadata_type> size_meta_{0};
  285. internal_value_type min_{0}, delta_{1};
  286. template <class V, class T, class M, class O>
  287. friend class regular;
  288. };
  289. #if __cpp_deduction_guides >= 201606
  290. template <class T>
  291. regular(unsigned, T, T)->regular<detail::convert_integer<T, double>>;
  292. template <class T>
  293. regular(unsigned, T, T, const char*)->regular<detail::convert_integer<T, double>>;
  294. template <class T, class M>
  295. regular(unsigned, T, T, M)->regular<detail::convert_integer<T, double>, transform::id, M>;
  296. template <class Tr, class T>
  297. regular(Tr, unsigned, T, T)->regular<detail::convert_integer<T, double>, Tr>;
  298. template <class Tr, class T>
  299. regular(Tr, unsigned, T, T, const char*)->regular<detail::convert_integer<T, double>, Tr>;
  300. template <class Tr, class T, class M>
  301. regular(Tr, unsigned, T, T, M)->regular<detail::convert_integer<T, double>, Tr, M>;
  302. #endif
  303. template <class Value = double, class MetaData = use_default, class Options = use_default>
  304. using circular = regular<Value, transform::id, MetaData,
  305. decltype(detail::replace_default<Options, option::overflow_t>{} |
  306. option::circular)>;
  307. } // namespace axis
  308. } // namespace histogram
  309. } // namespace boost
  310. #endif