From 7be6dd27f677e09cfd96bf9e1fe4e4c39ed55c37 Mon Sep 17 00:00:00 2001 From: Hans Dembinski Date: Fri, 11 Jan 2019 00:11:34 +0100 Subject: [PATCH] wip, support for growing axes --- CMakeLists.txt | 1 + include/boost/histogram/algorithm/project.hpp | 34 +- include/boost/histogram/algorithm/reduce.hpp | 37 +- .../histogram/axis/ostream_operators.hpp | 7 +- .../boost/histogram/axis/polymorphic_bin.hpp | 2 +- include/boost/histogram/axis/traits.hpp | 16 +- include/boost/histogram/axis/variable.hpp | 8 +- include/boost/histogram/axis/variant.hpp | 22 +- include/boost/histogram/detail/axes.hpp | 325 +----------------- include/boost/histogram/detail/linearize.hpp | 289 ++++++++++++++++ include/boost/histogram/detail/meta.hpp | 115 ++++--- include/boost/histogram/histogram.hpp | 15 +- include/boost/histogram/indexed.hpp | 19 +- include/boost/histogram/make_histogram.hpp | 8 +- test/axis_variant_test.cpp | 162 ++++----- test/detail_test.cpp | 26 -- test/histogram_growing_test.cpp | 74 ++++ test/meta_test.cpp | 46 +-- test/utility_histogram.hpp | 2 +- test/utility_meta.hpp | 15 +- 20 files changed, 635 insertions(+), 588 deletions(-) create mode 100644 include/boost/histogram/detail/linearize.hpp create mode 100644 test/histogram_growing_test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 8d0db76c..89bc4bf3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -93,6 +93,7 @@ compiled_test(test/axis_category_test.cpp) compiled_test(test/axis_variant_test.cpp) compiled_test(test/detail_test.cpp) compiled_test(test/histogram_dynamic_test.cpp) +compiled_test(test/histogram_growing_test.cpp) compiled_test(test/histogram_mixed_test.cpp) compiled_test(test/histogram_operators_test.cpp) compiled_test(test/histogram_test.cpp) diff --git a/include/boost/histogram/algorithm/project.hpp b/include/boost/histogram/algorithm/project.hpp index 21e89261..68d78fca 100644 --- a/include/boost/histogram/algorithm/project.hpp +++ b/include/boost/histogram/algorithm/project.hpp @@ -8,7 +8,6 @@ #define BOOST_HISTOGRAM_ALGORITHM_PROJECT_HPP #include -#include #include #include #include @@ -30,18 +29,27 @@ namespace algorithm { histogram is summed over the removed axes. */ template -auto project(const histogram& h, std::integral_constant n, Ns... ns) { - using LN = mp11::mp_list; +auto project(const histogram& h, std::integral_constant I, Ns... Is) { + using LN = mp11::mp_list; static_assert(mp11::mp_is_set::value, "indices must be unique"); - auto axes = detail::make_sub_axes(unsafe_access::axes(h), n, ns...); - auto result = histogram( - std::move(axes), detail::make_default(unsafe_access::storage(h))); + const auto& old_axes = unsafe_access::axes(h); + auto axes = detail::static_if>( + [&](const auto& old_axes) { + return std::make_tuple(std::get(old_axes), std::get(old_axes)...); + }, + [&](const auto& old_axes) { + return detail::naked({old_axes[I], old_axes[Is]...}); + }, + old_axes); - detail::axes_buffer idx(result.rank()); + const auto& old_storage = unsafe_access::storage(h); + using A2 = decltype(axes); + auto result = histogram(std::move(axes), detail::make_default(old_storage)); + auto idx = detail::make_stack_buffer(unsafe_access::axes(result)); for (auto x : indexed(h, true)) { auto i = idx.begin(); - mp11::mp_for_each([&i, &x](auto I) { *i++ = x[I]; }); + mp11::mp_for_each([&i, &x](auto J) { *i++ = x[J]; }); result.at(idx) += *x; } return result; @@ -61,22 +69,22 @@ auto project(const histogram& h, const Iterable& c) { const auto& old_axes = unsafe_access::axes(h); auto axes = detail::make_default(old_axes); axes.reserve(c.size()); - detail::axes_buffer seen(old_axes.size(), false); + auto seen = detail::make_stack_buffer(old_axes, false); for (auto d : c) { if (seen[d]) BOOST_THROW_EXCEPTION(std::invalid_argument("indices must be unique")); seen[d] = true; axes.emplace_back(old_axes[d]); } - auto result = - histogram(std::move(axes), detail::make_default(unsafe_access::storage(h))); - - detail::axes_buffer idx(result.rank()); + const auto& old_storage = unsafe_access::storage(h); + auto result = histogram(std::move(axes), detail::make_default(old_storage)); + auto idx = detail::make_stack_buffer(unsafe_access::axes(result)); for (auto x : indexed(h, true)) { auto i = idx.begin(); for (auto d : c) *i++ = x[d]; result.at(idx) += *x; } + return result; } diff --git a/include/boost/histogram/algorithm/reduce.hpp b/include/boost/histogram/algorithm/reduce.hpp index 59122ab2..4a1a36d2 100644 --- a/include/boost/histogram/algorithm/reduce.hpp +++ b/include/boost/histogram/algorithm/reduce.hpp @@ -67,33 +67,44 @@ inline reduce_option_type rebin(unsigned merge) { return rebin(0, merge); } template > decltype(auto) reduce(const Histogram& h, const C& options) { - using axes_type = typename Histogram::axes_type; + const auto& old_axes = unsafe_access::axes(h); struct option_item : reduce_option_type { int begin, end; - bool is_set() const noexcept { return reduce_option_type::merge > 0; } }; - auto options_internal = detail::axes_buffer(h.rank()); + auto opts = detail::make_stack_buffer(old_axes); for (const auto& o : options) { - auto& oi = options_internal[o.iaxis]; - if (oi.is_set()) // did we already set the option for this axis? + auto& oi = opts[o.iaxis]; + if (oi.merge > 0) // did we already set the option for this axis? BOOST_THROW_EXCEPTION(std::invalid_argument("indices must be unique")); oi.lower = o.lower; oi.upper = o.upper; oi.merge = o.merge; } - auto axes = detail::make_empty_axes(unsafe_access::axes(h)); + auto axes = detail::static_if>>( + [](const auto& c) { return detail::naked(); }, + [](const auto& c) { + using A = detail::naked; + auto axes = A(c.get_allocator()); + axes.reserve(c.size()); + detail::for_each_axis(c, [&axes](const auto& a) { + using U = detail::naked; + axes.emplace_back(U()); + }); + return axes; + }, + old_axes); unsigned iaxis = 0; h.for_each_axis([&](const auto& a) { - using T = detail::unqual; + using T = detail::naked; - auto& o = options_internal[iaxis]; + auto& o = opts[iaxis]; o.begin = 0; o.end = a.size(); - if (o.is_set()) { + if (o.merge > 0) { // option is set? if (o.lower < o.upper) { while (o.begin != o.end && a.value(o.begin) < o.lower) ++o.begin; while (o.end != o.begin && a.value(o.end - 1) >= o.upper) --o.end; @@ -112,13 +123,13 @@ decltype(auto) reduce(const Histogram& h, const C& options) { ++iaxis; }); - auto result = - Histogram(std::move(axes), detail::make_default(unsafe_access::storage(h))); + auto storage = detail::make_default(unsafe_access::storage(h)); + auto result = Histogram(std::move(axes), std::move(storage)); - detail::axes_buffer idx(h.rank()); + auto idx = detail::make_stack_buffer(unsafe_access::axes(result)); for (auto x : indexed(h, true)) { auto i = idx.begin(); - auto o = options_internal.begin(); + auto o = opts.begin(); for (auto j : x) { *i = (j - o->begin); if (*i <= -1) diff --git a/include/boost/histogram/axis/ostream_operators.hpp b/include/boost/histogram/axis/ostream_operators.hpp index 01affae5..8a1503ce 100644 --- a/include/boost/histogram/axis/ostream_operators.hpp +++ b/include/boost/histogram/axis/ostream_operators.hpp @@ -37,8 +37,7 @@ void stream_metadata(OStream& os, const T& t) { if (!oss.str().empty()) { os << ", metadata=" << std::quoted(oss.str()); } }, [&os](const auto&) { - using U = detail::unqual; - os << ", metadata=" << boost::core::demangled_name(BOOST_CORE_TYPEID(U)); + os << ", metadata=" << boost::core::demangled_name(BOOST_CORE_TYPEID(T)); }, t); } @@ -111,7 +110,7 @@ template std::basic_ostream& operator<<(std::basic_ostream& os, const polymorphic_bin& i) { if (i.is_discrete()) - os << i; + os << static_cast(i); else os << "[" << i.lower() << ", " << i.upper() << ")"; return os; @@ -169,7 +168,7 @@ std::basic_ostream& operator<<(std::basic_ostream& os, const variant& v) { visit( [&os](const auto& x) { - using A = detail::unqual; + using A = detail::naked; detail::static_if>( [&os](const auto& x) { os << x; }, [](const auto&) { diff --git a/include/boost/histogram/axis/polymorphic_bin.hpp b/include/boost/histogram/axis/polymorphic_bin.hpp index 3d854853..203d2602 100644 --- a/include/boost/histogram/axis/polymorphic_bin.hpp +++ b/include/boost/histogram/axis/polymorphic_bin.hpp @@ -38,7 +38,7 @@ public: polymorphic_bin(value_type lower, value_type upper) : lower_or_value_(lower), upper_(upper) {} - operator value_type() const noexcept { return lower_or_value_; } + operator const value_type&() const noexcept { return lower_or_value_; } value_type lower() const noexcept { return lower_or_value_; } value_type upper() const noexcept { return upper_; } diff --git a/include/boost/histogram/axis/traits.hpp b/include/boost/histogram/axis/traits.hpp index 965ba5e4..e94a9d05 100644 --- a/include/boost/histogram/axis/traits.hpp +++ b/include/boost/histogram/axis/traits.hpp @@ -30,15 +30,14 @@ axis::option_type options_impl(const T& t, std::true_type) { template decltype(auto) value_method_switch(FIntArg&& iarg, FDoubleArg&& darg, const T& t) { - using U = unqual; - return static_if>( + return static_if>( [](FIntArg&& iarg, FDoubleArg&& darg, const auto& t) { - using A = unqual; + using A = naked; return static_if, int>>( std::forward(iarg), std::forward(darg), t); }, [](FIntArg&&, FDoubleArg&&, const auto& t) { - using A = unqual; + using A = naked; BOOST_THROW_EXCEPTION(std::runtime_error(detail::cat( boost::core::demangled_name(BOOST_CORE_TYPEID(A)), " has no value method"))); return 0; @@ -48,16 +47,15 @@ decltype(auto) value_method_switch(FIntArg&& iarg, FDoubleArg&& darg, const T& t template R2 value_method_switch_with_return_type(FIntArg&& iarg, FDoubleArg&& darg, const T& t) { - using U = unqual; - return static_if>( + return static_if>( [](FIntArg&& iarg, FDoubleArg&& darg, const auto& t) -> R2 { - using A = unqual; + using A = naked; return static_if, int>>( std::forward(iarg), std::forward(darg), t); }, [](FIntArg&&, FDoubleArg&&, const auto&) -> R2 { BOOST_THROW_EXCEPTION(std::runtime_error( - detail::cat(boost::core::demangled_name(BOOST_CORE_TYPEID(U)), + detail::cat(boost::core::demangled_name(BOOST_CORE_TYPEID(T)), " has no value method or return type is not convertible to ", boost::core::demangled_name(BOOST_CORE_TYPEID(R1))))); }, @@ -69,7 +67,7 @@ namespace axis { namespace traits { template decltype(auto) metadata(T&& t) noexcept { - return detail::static_if>( + return detail::static_if>>( [](auto&& x) -> decltype(auto) { return x.metadata(); }, [](T &&) -> detail::copy_qualifiers { static axis::null_type m; diff --git a/include/boost/histogram/axis/variable.hpp b/include/boost/histogram/axis/variable.hpp index 8218057e..282bbf7c 100644 --- a/include/boost/histogram/axis/variable.hpp +++ b/include/boost/histogram/axis/variable.hpp @@ -116,7 +116,7 @@ public: BOOST_ASSERT((end - begin) % merge == 0); if (Options & option_type::circular && !(begin == 0 && end == src.size())) BOOST_THROW_EXCEPTION(std::invalid_argument("cannot shrink circular axis")); - using It = const detail::unqual*; + using It = const detail::naked*; struct skip_iterator { It it; unsigned skip; @@ -230,17 +230,17 @@ variable(std::initializer_list, M)->variable; template ()))>, double>> + detail::naked()))>, double>> variable(Iterable)->variable; template ()))>, double>> + detail::naked()))>, double>> variable(Iterable, const char*)->variable; template ()))>, double>> + detail::naked()))>, double>> variable(Iterable, M)->variable; #endif diff --git a/include/boost/histogram/axis/variant.hpp b/include/boost/histogram/axis/variant.hpp index 1b77b4fc..3431744c 100644 --- a/include/boost/histogram/axis/variant.hpp +++ b/include/boost/histogram/axis/variant.hpp @@ -45,16 +45,16 @@ namespace axis { template class variant : private boost::variant, public iterator_mixin> { using base_type = boost::variant; - using first_bounded_type = detail::unqual>; + using first_bounded_type = detail::naked>; - using types = mp11::mp_transform; + using types = mp11::mp_transform; template using requires_bounded_type = - mp11::mp_if>, void>; + mp11::mp_if>, void>; public: using metadata_type = - detail::unqual()))>; + detail::naked()))>; variant() = default; variant(const variant&) = default; @@ -80,7 +80,7 @@ public: variant& operator=(const variant& u) { visit( [this](const auto& u) { - using U = detail::unqual; + using U = detail::naked; detail::static_if>( [this](const auto& u) { this->operator=(u); }, [](const auto&) { @@ -146,7 +146,7 @@ public: int operator()(U&& x) const { return visit( [&x](const auto& a) { - using A = detail::unqual; + using A = detail::naked; using arg_t = detail::arg_type; return detail::static_if>( [&x](const auto& a) -> int { return a(x); }, @@ -272,21 +272,21 @@ const T* get(const variant* v) { } // pass-through version for generic programming, if U is axis instead of variant -template >> +template >> auto get(U&& u) -> detail::copy_qualifiers { return std::forward(u); } // pass-through version for generic programming, if U is axis instead of variant -template >> +template >> T* get(U* u) { - return std::is_same>::value ? reinterpret_cast(u) : nullptr; + return std::is_same>::value ? reinterpret_cast(u) : nullptr; } // pass-through version for generic programming, if U is axis instead of variant -template >> +template >> const T* get(const U* u) { - return std::is_same>::value ? reinterpret_cast(u) + return std::is_same>::value ? reinterpret_cast(u) : nullptr; } } // namespace axis diff --git a/include/boost/histogram/detail/axes.hpp b/include/boost/histogram/detail/axes.hpp index 3d388a2e..7aabd939 100644 --- a/include/boost/histogram/detail/axes.hpp +++ b/include/boost/histogram/detail/axes.hpp @@ -7,9 +7,7 @@ #ifndef BOOST_HISTOGRAM_DETAIL_AXES_HPP #define BOOST_HISTOGRAM_DETAIL_AXES_HPP -#include #include -#include #include #include #include @@ -19,86 +17,61 @@ #include #include #include -#include -#ifdef BOOST_HISTOGRAM_WITH_ACCUMULATORS_SUPPORT -#include -#endif - -/* Most of the histogram code is generic and works for any number of axes. Buffers with a - * fixed maximum capacity are used in some places, which have a size equal to the rank of - * a histogram. The buffers are statically allocated to improve performance, which means - * that they need a preset maximum capacity. 32 seems like a safe upper limit for the rank - * (you can nevertheless increase it here if necessary): the simplest non-trivial axis has - * 2 bins; even if counters are used which need only a byte of storage per bin, this still - * corresponds to 4 GB of storage. - */ -#ifndef BOOST_HISTOGRAM_DETAIL_AXES_LIMIT -#define BOOST_HISTOGRAM_DETAIL_AXES_LIMIT 32 -#endif namespace boost { namespace histogram { namespace detail { -template -struct is_accumulator_set : std::false_type {}; - -#ifdef BOOST_HISTOGRAM_WITH_ACCUMULATORS_SUPPORT -template -struct is_accumulator_set<::boost::accumulators::accumulator_set> - : std::true_type {}; -#endif - -template +template decltype(auto) axis_get(std::tuple& axes) { return std::get(axes); } -template +template decltype(auto) axis_get(const std::tuple& axes) { return std::get(axes); } -template +template decltype(auto) axis_get(T& axes) { return axes[N]; } -template +template decltype(auto) axis_get(const T& axes) { return axes[N]; } -template +template decltype(auto) axis_get(std::tuple& axes, unsigned i) { return mp11::mp_with_index( i, [&](auto I) { return axis::variant(std::get(axes)); }); } -template +template decltype(auto) axis_get(const std::tuple& axes, unsigned i) { return mp11::mp_with_index( i, [&](auto I) { return axis::variant(std::get(axes)); }); } -template +template decltype(auto) axis_get(T& axes, unsigned i) { return axes.at(i); } -template +template decltype(auto) axis_get(const T& axes, unsigned i) { return axes.at(i); } -template +template bool axes_equal(const std::tuple& t, const std::tuple& u) { return static_if, mp11::mp_list>>( [](const auto& a, const auto& b) { return relaxed_equal(a, b); }, [](const auto&, const auto&) { return false; }, t, u); } -template +template bool axes_equal(const std::tuple& t, const U& u) { if (sizeof...(Ts) != u.size()) return false; bool equal = true; @@ -110,18 +83,18 @@ bool axes_equal(const std::tuple& t, const U& u) { return equal; } -template +template bool axes_equal(const T& t, const std::tuple& u) { return axes_equal(u, t); } -template +template bool axes_equal(const T& t, const U& u) { if (t.size() != u.size()) return false; return std::equal(t.begin(), t.end(), u.begin()); } -template +template void axes_assign(std::tuple& t, const std::tuple& u) { static_if, mp11::mp_list>>( [](auto& a, const auto& b) { a = b; }, @@ -132,7 +105,7 @@ void axes_assign(std::tuple& t, const std::tuple& u) { t, u); } -template +template void axes_assign(std::tuple& t, const U& u) { mp11::mp_for_each>([&](auto I) { using T = mp11::mp_at_c, I>; @@ -140,7 +113,7 @@ void axes_assign(std::tuple& t, const U& u) { }); } -template +template void axes_assign(T& t, const std::tuple& u) { t.resize(sizeof...(Us)); mp11::mp_for_each>( @@ -153,18 +126,8 @@ void axes_assign(T& t, const U& u) { } template -constexpr std::size_t axes_size(const T& axes) noexcept { - return static_if>>( - [](const auto& a) { - using U = unqual; - return std::tuple_size::value; - }, - [&](const auto& a) { return a.size(); }, axes); -} - -template -void rank_check(const T& axes, const unsigned N) { - BOOST_ASSERT_MSG(N < axes_size(axes), "index out of range"); +void axis_index_is_valid(const T& axes, const unsigned N) { + BOOST_ASSERT_MSG(N < get_size(axes), "index out of range"); } template @@ -174,12 +137,12 @@ void for_each_axis_impl(std::true_type, const T& axes, F&& f) { template void for_each_axis_impl(std::false_type, const T& axes, F&& f) { - for (const auto& x : axes) f(x); + for (const auto& x : axes) std::forward(f)(x); } template void for_each_axis(const T& axes, F&& f) { - using U = mp11::mp_first>; + using U = mp11::mp_first; for_each_axis_impl(is_axis_variant(), axes, std::forward(f)); } @@ -188,28 +151,6 @@ void for_each_axis(const std::tuple& axes, F&& f) { mp11::tuple_for_each(axes, std::forward(f)); } -template -using axes_buffer = boost::container::static_vector< - T, mp11::mp_eval_if_c::value), - mp11::mp_size_t, - std::tuple_size, Axes>::value>; - -template -auto make_empty_axes(const T& t) { - auto r = T(t.get_allocator()); - r.reserve(t.size()); - for_each_axis(t, [&r](const auto& a) { - using U = unqual; - r.emplace_back(U()); - }); - return r; -} - -template -auto make_empty_axes(const std::tuple&) { - return std::tuple(); -} - template std::size_t bincount(const T& axes) { std::size_t n = 1; @@ -217,234 +158,6 @@ std::size_t bincount(const T& axes) { return n; } -template -auto make_sub_axes(const std::tuple& t, Ns... ns) { - return std::make_tuple(std::get(t)...); -} - -template -auto make_sub_axes(const T& t, Ns... ns) { - return T({t[ns]...}, t.get_allocator()); -} - -/// Index with an invalid state -struct optional_index { - std::size_t idx = 0; - std::size_t stride = 1; - operator bool() const { return stride > 0; } - std::size_t operator*() const { return idx; } -}; - -inline void linearize(optional_index& out, const int axis_shape, int j) noexcept { - // j is internal index, which is potentially shifted by +1 wrt external index - out.idx += j * out.stride; - // set stride to 0, if j is invalid - out.stride *= (0 <= j && j < axis_shape) * axis_shape; -} - -template -void linearize_value(std::true_type, optional_index& out, const A& axis, const U& u) { - const auto extend = axis::traits::extend(axis); - const auto opt = axis::traits::options(axis); - const auto j = axis(u) + (opt & axis::option_type::underflow); - return linearize(out, extend, j); -} - -template -void linearize_value(std::false_type, optional_index&, const A&, const U&) { - // protect against instantiation with wrong template argument - using arg_t = arg_type; - BOOST_THROW_EXCEPTION(std::invalid_argument( - detail::cat(boost::core::demangled_name(BOOST_CORE_TYPEID(A)), - ": cannot convert argument of type ", - boost::core::demangled_name(BOOST_CORE_TYPEID(U)), " to ", - boost::core::demangled_name(BOOST_CORE_TYPEID(arg_t))))); -} - -template -void linearize_value(optional_index& out, const A& axis, const U& u) { - // protect against instantiation with wrong template argument - using arg_t = arg_type; - return linearize_value(std::is_convertible(), out, axis, u); -} - -template -void linearize_value(optional_index& out, const axis::variant& axis, const U& u) { - axis::visit([&](const auto& a) { linearize_value(out, a, u); }, axis); -} - -template -void linearize_index(optional_index& out, const T& axis, const int j) { - const auto extend = axis::traits::extend(axis); - const auto opt = axis::traits::options(axis); - linearize(out, extend, j + (opt & axis::option_type::underflow)); -} - -// special case: if histogram::operator()(tuple(1, 2)) is called on 1d histogram with axis -// that accepts 2d tuple, this should not fail -// - solution is to forward tuples of size > 1 directly to axis for 1d histograms -// - has nice side-effect of making histogram::operator(1, 2) work as well -// - cannot detect call signature of axis at compile-time in all configurations -// (axis::variant provides generic call interface and hides concrete interface), -// so we throw at runtime if incompatible argument is passed (e.g. 3d tuple) -template -optional_index args_to_index(const std::tuple& axes, const U& args) { - optional_index idx; - if (N > 1) { - linearize_value(idx, std::get<0>(axes), sub_tuple(args)); - } else { - linearize_value(idx, std::get<0>(axes), std::get(args)); - } - return idx; -} - -template -optional_index args_to_index(const std::tuple& axes, const U& args) { - static_assert(sizeof...(Ts) + 2 == N, "number of arguments != histogram rank"); - optional_index idx; - mp11::mp_for_each>([&](auto I) { - linearize_value(idx, std::get(axes), std::get<(Offset + I)>(args)); - }); - return idx; -} - -// overload for dynamic axes -template -optional_index args_to_index(const T& axes, const U& args) { - const unsigned m = axes.size(); - optional_index idx; - if (m == 1 && N > 1) - linearize_value(idx, axes[0], sub_tuple(args)); - else { - if (m != N) - BOOST_THROW_EXCEPTION( - std::invalid_argument("number of arguments != histogram rank")); - mp11::mp_for_each>( - [&](auto I) { linearize_value(idx, axes[I], std::get<(Offset + I)>(args)); }); - } - return idx; -} - -template -constexpr std::pair weight_sample_indices() { - if (is_weight::value) return std::make_pair(0, -1); - if (is_sample::value) return std::make_pair(-1, 0); - return std::make_pair(-1, -1); -} - -template -constexpr std::pair weight_sample_indices() { - using L = mp11::mp_list; - const int n = sizeof...(Us) + 1; - if (is_weight>::value) { - if (is_sample>::value) return std::make_pair(0, 1); - if (is_sample>::value) return std::make_pair(0, n); - return std::make_pair(0, -1); - } - if (is_sample>::value) { - if (is_weight>::value) return std::make_pair(1, 0); - if (is_weight>::value) return std::make_pair(n, 0); - return std::make_pair(-1, 0); - } - if (is_weight>::value) { - // 0, n already covered - if (is_sample>::value) return std::make_pair(n, n - 1); - return std::make_pair(n, -1); - } - if (is_sample>::value) { - // n, 0 already covered - if (is_weight>::value) return std::make_pair(n - 1, n); - return std::make_pair(-1, n); - } - return std::make_pair(-1, -1); -} - -template -void fill_storage_impl(mp11::mp_int<-1>, mp11::mp_int<-1>, T&& t, U&&) { - static_if>([](auto&& t) { ++t; }, [](auto&& t) { t(); }, - std::forward(t)); -} - -template -void fill_storage_impl(IW, mp11::mp_int<-1>, T&& t, U&& args) { - static_if>( - [](auto&& t, const auto& w) { t += w; }, - [](auto&& t, const auto& w) { -#ifdef BOOST_HISTOGRAM_WITH_ACCUMULATORS_SUPPORT - static_if>>( - [w](auto&& t) { t(::boost::accumulators::weight = w); }, - [w](auto&& t) { t(w); }, t); -#else - t(w); -#endif - }, - std::forward(t), std::get(args).value); -} - -template -void fill_storage_impl(mp11::mp_int<-1>, IS, T&& t, U&& args) { - mp11::tuple_apply([&t](auto&&... args) { t(args...); }, - std::get(args).value); -} - -template -void fill_storage_impl(IW, IS, T&& t, U&& args) { -#ifdef BOOST_HISTOGRAM_WITH_ACCUMULATORS_SUPPORT - static_if>>( - [](auto&& t, const auto& w, const auto& s) { - mp11::tuple_apply( - [&](auto&&... args) { t(args..., ::boost::accumulators::weight = w); }, s); - }, - [](auto&& t, const auto& w, const auto& s) { - mp11::tuple_apply([&](auto&&... args) { t(w, args...); }, s); - }, - std::forward(t), std::get(args).value, - std::get(args).value); -#else - mp11::tuple_apply( - [&](auto&&... args2) { t(std::get(args).value, args2...); }, - std::get(args).value); -#endif -} - -template -void fill_impl(S& storage, const T& axes, const std::tuple& args) { - constexpr std::pair iws = weight_sample_indices(); - constexpr unsigned n = sizeof...(Us) - (iws.first > -1) - (iws.second > -1); - constexpr unsigned offset = (iws.first == 0 || iws.second == 0) - ? (iws.first == 1 || iws.second == 1 ? 2 : 1) - : 0; - optional_index idx = args_to_index(axes, args); - if (idx) { - fill_storage_impl(mp11::mp_int(), mp11::mp_int(), - storage[*idx], args); - } -} - -template -optional_index at_impl(const A& axes, const std::tuple& args) { - if (axes_size(axes) != sizeof...(Us)) - BOOST_THROW_EXCEPTION(std::invalid_argument("number of arguments != histogram rank")); - optional_index idx; - mp11::mp_for_each>([&](auto I) { - linearize_index(idx, axis_get(axes), static_cast(std::get(args))); - }); - return idx; -} - -template -optional_index at_impl(const A& axes, const U& args) { - if (axes_size(axes) != args.size()) - BOOST_THROW_EXCEPTION(std::invalid_argument("number of arguments != histogram rank")); - optional_index idx; - using std::begin; - auto it = begin(args); - for_each_axis(axes, - [&](const auto& a) { linearize_index(idx, a, static_cast(*it++)); }); - return idx; -} - } // namespace detail } // namespace histogram } // namespace boost diff --git a/include/boost/histogram/detail/linearize.hpp b/include/boost/histogram/detail/linearize.hpp new file mode 100644 index 00000000..e936d7a0 --- /dev/null +++ b/include/boost/histogram/detail/linearize.hpp @@ -0,0 +1,289 @@ +// Copyright 2015-2018 Hans Dembinski +// +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_HISTOGRAM_DETAIL_LINEARIZE_HPP +#define BOOST_HISTOGRAM_DETAIL_LINEARIZE_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#ifdef BOOST_HISTOGRAM_WITH_ACCUMULATORS_SUPPORT +#include +#endif + +namespace boost { +namespace histogram { +namespace detail { + +template +struct is_accumulator_set : std::false_type {}; + +#ifdef BOOST_HISTOGRAM_WITH_ACCUMULATORS_SUPPORT +template +struct is_accumulator_set<::boost::accumulators::accumulator_set> + : std::true_type {}; +#endif + +/// Index with an invalid state +struct optional_index { + std::size_t idx = 0; + std::size_t stride = 1; + operator bool() const { return stride > 0; } + std::size_t operator*() const { return idx; } +}; + +inline void linearize(optional_index& out, const int axis_shape, int j) noexcept { + // j is internal index, shifted by +1 wrt external index if axis has underflow bin + out.idx += j * out.stride; + // set stride to 0, if j is invalid + out.stride *= (0 <= j && j < axis_shape) * axis_shape; +} + +template +void linearize_value(optional_index& out, int& shift, A& axis, const V& value) { + using T = arg_type; + // turn instantiation with wrong template argument into runtime error + static_if>( + [&](auto& axis, const V& value) { + const auto opt = axis::traits::options(axis); + const auto j = static_if>>( + [&shift](auto& axis, const V& value) { + int j; + std::tie(j, shift) = axis.update(value); + return j; + }, + [](const auto& axis, const V& value) { return axis(value); }, + axis, value) + + (opt & axis::option_type::underflow); + linearize(out, axis::traits::extend(axis), j); + }, + [](A&, const V&) { + using T = arg_type; + BOOST_THROW_EXCEPTION(std::invalid_argument( + detail::cat(boost::core::demangled_name(BOOST_CORE_TYPEID(A)), + ": cannot convert argument of type ", + boost::core::demangled_name(BOOST_CORE_TYPEID(V)), " to ", + boost::core::demangled_name(BOOST_CORE_TYPEID(T))))); + }, + axis, value); +} + +template +void linearize_value(optional_index& out, int& s, axis::variant& axis, + const V& v) { + axis::visit([&](auto& a) { linearize_value(out, s, a, v); }, axis); +} + +template +void linearize_index(optional_index& out, const T& axis, const int j) { + const auto extend = axis::traits::extend(axis); + const auto opt = axis::traits::options(axis); + linearize(out, extend, j + (opt & axis::option_type::underflow)); +} + +template +void maybe_replace_storage(S& storage, const A& axes, const T& shifts) { + bool update_needed = false; + for (int s : shifts) update_needed |= s != 0; + if (!update_needed) return; + struct item { + int idx, size; + std::size_t stride; + }; + auto data = make_stack_buffer(axes); + auto sit = shifts.begin(); + auto dit = data.begin(); + std::size_t s = 1; + for_each_axis(axes, [&](const auto& a) { + const auto n = axis::traits::extend(a); + *dit++ = {0, n, s}; + s *= n - std::abs(*sit++); + }); + auto new_storage = make_default(storage); + new_storage.reset(detail::bincount(axes)); + for (const auto& x : storage) { + auto ns = new_storage.begin(); + sit = shifts.begin(); + for (const auto& d : data) { ns += (d.idx + std::max(*sit++, 9)) * d.stride; } + auto dit = data.begin(); + const auto last = data.end() - 1; + ++dit->idx; + while (dit != last && dit->idx == dit->size) { + dit->idx = 0; + ++(++dit)->idx; + } + *ns = x; + } + storage = std::move(new_storage); +} + +template +struct size_or_zero : mp11::mp_size_t<0> {}; + +template +struct size_or_zero> : mp11::mp_size_t {}; + +// special case: if histogram::operator()(tuple(1, 2)) is called on 1d histogram with axis +// that accepts 2d tuple, this should not fail +// - solution is to forward tuples of size > 1 directly to axis for 1d histograms +// - has nice side-effect of making histogram::operator(1, 2) work as well +// - cannot detect call signature of axis at compile-time in all configurations +// (axis::variant provides generic call interface and hides concrete interface), +// so we throw at runtime if incompatible argument is passed (e.g. 3d tuple) +template +optional_index args_to_index(S& storage, T& axes, const U& args) { + optional_index idx; + auto shifts = make_stack_buffer(axes); + const auto rank = get_size(axes); + if (rank == 1 && N > 1) + linearize_value(idx, shifts[0], axis_get<0>(axes), tuple_slice(args)); + else { + if (rank != N) + BOOST_THROW_EXCEPTION( + std::invalid_argument("number of arguments != histogram rank")); + constexpr unsigned M = size_or_zero>::value; + mp11::mp_for_each>([&](auto J) { + linearize_value(idx, shifts[J], axis_get(axes), std::get<(J + I)>(args)); + }); + } + maybe_replace_storage(storage, axes, shifts); + return idx; +} + +template +constexpr auto weight_sample_indices() { + if (is_weight::value) return std::make_pair(0, -1); + if (is_sample::value) return std::make_pair(-1, 0); + return std::make_pair(-1, -1); +} + +template +constexpr auto weight_sample_indices() { + using L = mp11::mp_list; + const int n = sizeof...(Us) + 1; + if (is_weight>::value) { + if (is_sample>::value) return std::make_pair(0, 1); + if (is_sample>::value) return std::make_pair(0, n); + return std::make_pair(0, -1); + } + if (is_sample>::value) { + if (is_weight>::value) return std::make_pair(1, 0); + if (is_weight>::value) return std::make_pair(n, 0); + return std::make_pair(-1, 0); + } + if (is_weight>::value) { + // 0, n already covered + if (is_sample>::value) return std::make_pair(n, n - 1); + return std::make_pair(n, -1); + } + if (is_sample>::value) { + // n, 0 already covered + if (is_weight>::value) return std::make_pair(n - 1, n); + return std::make_pair(-1, n); + } + return std::make_pair(-1, -1); +} + +template +void fill_storage_impl(mp11::mp_int<-1>, mp11::mp_int<-1>, T&& t, U&&) { + static_if>>([](auto&& t) { ++t; }, [](auto&& t) { t(); }, + std::forward(t)); +} + +template +void fill_storage_impl(IW, mp11::mp_int<-1>, T&& t, U&& args) { + static_if>>( + [](auto&& t, const auto& w) { t += w; }, + [](auto&& t, const auto& w) { +#ifdef BOOST_HISTOGRAM_WITH_ACCUMULATORS_SUPPORT + static_if>>( + [w](auto&& t) { t(::boost::accumulators::weight = w); }, + [w](auto&& t) { t(w); }, t); +#else + t(w); +#endif + }, + std::forward(t), std::get(args).value); +} + +template +void fill_storage_impl(mp11::mp_int<-1>, IS, T&& t, U&& args) { + mp11::tuple_apply([&t](auto&&... args) { t(args...); }, + std::get(args).value); +} + +template +void fill_storage_impl(IW, IS, T&& t, U&& args) { +#ifdef BOOST_HISTOGRAM_WITH_ACCUMULATORS_SUPPORT + static_if>>( + [](auto&& t, const auto& w, const auto& s) { + mp11::tuple_apply( + [&](auto&&... args) { t(args..., ::boost::accumulators::weight = w); }, s); + }, + [](auto&& t, const auto& w, const auto& s) { + mp11::tuple_apply([&](auto&&... args) { t(w, args...); }, s); + }, + std::forward(t), std::get(args).value, + std::get(args).value); +#else + mp11::tuple_apply( + [&](auto&&... args2) { t(std::get(args).value, args2...); }, + std::get(args).value); +#endif +} + +template +void fill_impl(S& storage, A& axes, const std::tuple& args) { + constexpr std::pair iws = weight_sample_indices(); + constexpr unsigned n = sizeof...(Us) - (iws.first > -1) - (iws.second > -1); + constexpr unsigned i = (iws.first == 0 || iws.second == 0) + ? (iws.first == 1 || iws.second == 1 ? 2 : 1) + : 0; + optional_index idx = args_to_index(storage, axes, args); + if (idx) { + fill_storage_impl(mp11::mp_int(), mp11::mp_int(), + storage[*idx], args); + } +} + +template +optional_index at_impl(const A& axes, const std::tuple& args) { + if (get_size(axes) != sizeof...(Us)) + BOOST_THROW_EXCEPTION(std::invalid_argument("number of arguments != histogram rank")); + optional_index idx; + mp11::mp_for_each>([&](auto I) { + linearize_index(idx, axis_get(axes), static_cast(std::get(args))); + }); + return idx; +} + +template +optional_index at_impl(const A& axes, const U& args) { + if (get_size(axes) != args.size()) + BOOST_THROW_EXCEPTION(std::invalid_argument("number of arguments != histogram rank")); + optional_index idx; + using std::begin; + auto it = begin(args); + for_each_axis(axes, + [&](const auto& a) { linearize_index(idx, a, static_cast(*it++)); }); + return idx; +} + +} // namespace detail +} // namespace histogram +} // namespace boost + +#endif diff --git a/include/boost/histogram/detail/meta.hpp b/include/boost/histogram/detail/meta.hpp index 06f00e6a..5bbc8e3e 100644 --- a/include/boost/histogram/detail/meta.hpp +++ b/include/boost/histogram/detail/meta.hpp @@ -7,6 +7,18 @@ #ifndef BOOST_HISTOGRAM_DETAIL_META_HPP #define BOOST_HISTOGRAM_DETAIL_META_HPP +/* Most of the histogram code is generic and works for any number of axes. Buffers with a + * fixed maximum capacity are used in some places, which have a size equal to the rank of + * a histogram. The buffers are statically allocated to improve performance, which means + * that they need a preset maximum capacity. 32 seems like a safe upper limit for the rank + * (you can nevertheless increase it here if necessary): the simplest non-trivial axis has + * 2 bins; even if counters are used which need only a byte of storage per bin, this still + * corresponds to 4 GB of storage. + */ +#ifndef BOOST_HISTOGRAM_DETAIL_AXES_LIMIT +#define BOOST_HISTOGRAM_DETAIL_AXES_LIMIT 32 +#endif + #include #if BOOST_WORKAROUND(BOOST_GCC, >= 60000) #pragma GCC diagnostic push @@ -17,6 +29,7 @@ #if BOOST_WORKAROUND(BOOST_GCC, >= 60000) #pragma GCC diagnostic pop #endif +#include #include #include #include @@ -30,16 +43,10 @@ namespace histogram { namespace detail { template -using unqual = std::remove_cv_t>; +using naked = std::remove_cv_t>; template -using convert_integer = mp11::mp_if>, U, T>; - -template -using mp_size = mp11::mp_size>; - -template -using mp_at_c = mp11::mp_at_c, N>; +using convert_integer = mp11::mp_if>, U, T>; template