diff --git a/include/boost/histogram/dynamic_histogram.hpp b/include/boost/histogram/dynamic_histogram.hpp index 918d421a..4af62ea8 100644 --- a/include/boost/histogram/dynamic_histogram.hpp +++ b/include/boost/histogram/dynamic_histogram.hpp @@ -358,6 +358,8 @@ private: template friend class dynamic_histogram; + friend struct python_access; + template friend void serialize(Archiv&, dynamic_histogram&, unsigned); }; diff --git a/include/boost/histogram/dynamic_storage.hpp b/include/boost/histogram/dynamic_storage.hpp index 9a665d34..e3fa3ba3 100644 --- a/include/boost/histogram/dynamic_storage.hpp +++ b/include/boost/histogram/dynamic_storage.hpp @@ -15,13 +15,11 @@ #include #include #include -#include -#include -#include // malloc/free -#include // for bad_alloc exception +#include #include #include #include +#include namespace boost { namespace histogram { @@ -30,6 +28,8 @@ namespace detail { static_assert(std::is_pod::value, "weight_t must be POD"); + using mp_int = multiprecision::cpp_int; + using type_to_int = mpl::map< mpl::pair>, @@ -41,7 +41,9 @@ namespace detail { mpl::pair>, mpl::pair>, mpl::pair>, - mpl::pair> + mpl::pair>, + + mpl::pair> >; using int_to_type = mpl::map< @@ -49,16 +51,15 @@ namespace detail { mpl::pair, uint8_t>, mpl::pair, uint16_t>, mpl::pair, uint32_t>, - mpl::pair, uint64_t> + mpl::pair, uint64_t>, + mpl::pair, mp_int> >; template using storage_type = typename mpl::at< detail::int_to_type, - typename mpl::at< - detail::type_to_int, T - >::type + typename mpl::at::type >::type; template @@ -79,23 +80,44 @@ namespace detail { buffer(const buffer& o) : type_(o.type_), size_(o.size_), ptr_(nullptr) { - realloc(type_.depth_); - std::copy(static_cast(o.ptr_), - static_cast(o.ptr_) + size_ * type_.depth_, - static_cast(ptr_)); + switch (type_.id_) { + case -1: ptr_ = create(); copy_from(o.ptr_); break; + case 0: break; + case 1: ptr_ = create(); copy_from(o.ptr_); break; + case 2: ptr_ = create(); copy_from(o.ptr_); break; + case 3: ptr_ = create(); copy_from(o.ptr_); break; + case 4: ptr_ = create(); copy_from(o.ptr_); break; + case 5: ptr_ = create(); copy_from(o.ptr_); break; + } } buffer& operator=(const buffer& o) { if (this != &o) { if (size_ != o.size_ || type_.id_ != o.type_.id_) { + destroy_any(); size_ = o.size_; type_ = o.type_; - realloc(type_.depth_); + switch (type_.id_) { + case -1: ptr_ = create(); copy_from(o.ptr_); break; + case 0: break; + case 1: ptr_ = create(); copy_from(o.ptr_); break; + case 2: ptr_ = create(); copy_from(o.ptr_); break; + case 3: ptr_ = create(); copy_from(o.ptr_); break; + case 4: ptr_ = create(); copy_from(o.ptr_); break; + case 5: ptr_ = create(); copy_from(o.ptr_); break; + } + } else { + switch (type_.id_) { + case -1: copy_from(o.ptr_); break; + case 0: break; + case 1: copy_from(o.ptr_); break; + case 2: copy_from(o.ptr_); break; + case 3: copy_from(o.ptr_); break; + case 4: copy_from(o.ptr_); break; + case 5: copy_from(o.ptr_); break; + } } - std::copy(static_cast(o.ptr_), - static_cast(o.ptr_) + size_ * type_.depth_, - static_cast(ptr_)); } return *this; } @@ -127,65 +149,68 @@ namespace detail { ~buffer() { destroy_any(); } template - void create() { - type_.set(); - ptr_ = std::malloc(size_ * sizeof(T)); - new (ptr_) T[size_]; + T* create() { + std::allocator a; + T* ptr = a.allocate(size_); + new (ptr) T[size_]; + return ptr; } template void destroy() { - std::free(ptr_); + std::allocator a; + a.deallocate(static_cast(ptr_), size_); ptr_ = nullptr; } void destroy_any() { switch (type_.id_) { case -1: destroy(); break; - case 0: /* do nothing */ break; + case 0: break; case 1: destroy(); break; case 2: destroy(); break; case 3: destroy(); break; case 4: destroy(); break; + case 5: destroy(); break; } } + template + void copy_from(const void* p) { + std::copy(static_cast(p), static_cast(p) + size_, + static_cast(ptr_)); + } + template > void grow() { - static_assert(sizeof(U) >= sizeof(T), "U must be as large or larger than T"); - realloc(sizeof(U)); - std::copy_backward(&at(0), &at(size_), &at(size_)); + U* u = create(); + std::copy(&at(0), &at(size_), u); + destroy_any(); type_.set(); + ptr_ = u; } void wconvert() { switch (type_.id_) { - case -1: /* do nothing */ break; + case -1: break; case 0: initialize(); break; case 1: grow (); break; case 2: grow(); break; case 3: grow(); break; case 4: grow(); break; + case 5: grow(); break; } } template void initialize() { type_.set(); - ptr_ = nullptr; - realloc(sizeof(T)); + ptr_ = create(); std::fill(&at(0), &at(size_), T(0)); } - void realloc(unsigned d) - { - ptr_ = std::realloc(ptr_, size_ * d); - if (!ptr_ && (size_ * d > 0)) - throw std::bad_alloc(); - } - template T& at(std::size_t i) { return static_cast(ptr_)[i]; } @@ -217,12 +242,8 @@ namespace detail { } template <> - void increase_impl(buffer& b, std::size_t i) { - auto& bi = b.at(i); - if (bi < std::numeric_limits::max()) - ++bi; - else - throw std::overflow_error("histogram overflow"); + void increase_impl(buffer& b, std::size_t i) { + ++b.at(i); } template @@ -230,7 +251,7 @@ namespace detail { static void apply(buffer& b, std::size_t i, const TO& o) { auto& bi = b.at(i); if (static_cast(std::numeric_limits::max() - bi) >= o) - bi += o; + bi += static_cast(o); else { b.grow(); add_one_impl, TO>::apply(b, i, o); @@ -239,27 +260,26 @@ namespace detail { }; template - struct add_one_impl { + struct add_one_impl { static void apply(buffer& b, std::size_t i, const TO& o) { - auto& bi = b.at(i); - if (static_cast(std::numeric_limits::max() - bi) >= o) - bi += o; - else - throw std::overflow_error("histogram overflow"); + b.at(i) += o; } }; template void add_impl(buffer& b, const buffer& o) { - for (decltype(b.size_) i = 0; i < b.size_; ++i) + for (decltype(b.size_) i = 0; i < b.size_; ++i) { + const auto oi = o.at(i); switch (b.type_.id_) { - case -1: b.at(i) += o.at(i); break; + case -1: b.at(i) += oi; break; case 0: b.initialize(); // fall through - case 1: add_one_impl::apply(b, i, o.at(i)); break; - case 2: add_one_impl::apply(b, i, o.at(i)); break; - case 3: add_one_impl::apply(b, i, o.at(i)); break; - case 4: add_one_impl::apply(b, i, o.at(i)); break; + case 1: add_one_impl::apply(b, i, oi); break; + case 2: add_one_impl::apply(b, i, oi); break; + case 3: add_one_impl::apply(b, i, oi); break; + case 4: add_one_impl::apply(b, i, oi); break; + case 5: add_one_impl::apply(b, i, oi); break; } + } } template <> @@ -295,8 +315,9 @@ public: buffer_(o.size()) { using U = detail::storage_type; - buffer_.create(); - std::copy(o.data(), o.data() + buffer_.size_, &buffer_.at(0)); + buffer_.ptr_ = buffer_.create(); + buffer_.type_.set(); + std::copy(o.data(), o.data() + o.size(), &buffer_.at(0)); } template ; - buffer_.create(); - std::copy(o.data(), o.data() + buffer_.size_, &buffer_.at(0)); + buffer_.ptr_ = buffer_.create(); + buffer_.type_.set(); + std::copy(o.data(), o.data() + o.size(), &buffer_.at(0)); return *this; } @@ -324,6 +346,8 @@ public: private: detail::buffer buffer_; + friend struct python_access; + template friend void serialize(Archive&, dynamic_storage&, unsigned); }; @@ -338,6 +362,7 @@ void dynamic_storage::increase(std::size_t i) case 2: detail::increase_impl(buffer_, i); break; case 3: detail::increase_impl(buffer_, i); break; case 4: detail::increase_impl(buffer_, i); break; + case 5: detail::increase_impl(buffer_, i); break; } } @@ -353,11 +378,12 @@ dynamic_storage& dynamic_storage::operator+=(const dynamic_storage& o) { switch (o.buffer_.type_.id_) { case -1: detail::add_impl(buffer_, o.buffer_); break; - case 0: /* do nothing */ break; + case 0: break; case 1: detail::add_impl(buffer_, o.buffer_); break; case 2: detail::add_impl(buffer_, o.buffer_); break; case 3: detail::add_impl(buffer_, o.buffer_); break; case 4: detail::add_impl(buffer_, o.buffer_); break; + case 5: detail::add_impl(buffer_, o.buffer_); break; } return *this; } @@ -367,11 +393,12 @@ dynamic_storage::value_t dynamic_storage::value(std::size_t i) const { switch (buffer_.type_.id_) { case -1: return buffer_.at(i).w; - case 0: /* do nothing */ break; + case 0: break; case 1: return buffer_.at (i); case 2: return buffer_.at(i); case 3: return buffer_.at(i); case 4: return buffer_.at(i); + case 5: return static_cast(buffer_.at(i)); } return 0.0; } @@ -381,11 +408,12 @@ dynamic_storage::variance_t dynamic_storage::variance(std::size_t i) const { switch (buffer_.type_.id_) { case -1: return buffer_.at(i).w2; - case 0: /* do nothing */ break; + case 0: break; case 1: return buffer_.at (i); case 2: return buffer_.at(i); case 3: return buffer_.at(i); case 4: return buffer_.at(i); + case 5: return static_cast(buffer_.at(i)); } return 0.0; } diff --git a/include/boost/histogram/serialization.hpp b/include/boost/histogram/serialization.hpp index f3fc69ab..5383d6ab 100644 --- a/include/boost/histogram/serialization.hpp +++ b/include/boost/histogram/serialization.hpp @@ -72,21 +72,23 @@ inline void serialize(Archive& ar, dynamic_storage & store, unsigned version) ar & b.type_.depth_; if (Archive::is_loading::value) { switch (b.type_.id_) { - case -1: b.create(); break; + case -1: b.ptr_ = b.create(); break; case 0: b.ptr_ = nullptr; break; - case 1: b.create(); break; - case 2: b.create(); break; - case 3: b.create(); break; - case 4: b.create(); break; + case 1: b.ptr_ = b.create(); break; + case 2: b.ptr_ = b.create(); break; + case 3: b.ptr_ = b.create(); break; + case 4: b.ptr_ = b.create(); break; + case 5: b.ptr_ = b.create(); break; } } switch (b.type_.id_) { case -1: ar & serialization::make_array(&b.at(0), b.size_); break; - case 0: /* nothing to do */ break; + case 0: break; case 1: ar & serialization::make_array(&b.at(0), b.size_); break; case 2: ar & serialization::make_array(&b.at(0), b.size_); break; case 3: ar & serialization::make_array(&b.at(0), b.size_); break; case 4: ar & serialization::make_array(&b.at(0), b.size_); break; + case 5: ar & serialization::make_array(&b.at(0), b.size_); break; } } diff --git a/src/python/histogram.cpp b/src/python/histogram.cpp index abe76dd7..4cc8210b 100644 --- a/src/python/histogram.cpp +++ b/src/python/histogram.cpp @@ -229,23 +229,28 @@ histogram_variance(python::tuple args, python::dict kwargs) { return object(self.variance(idx + 0, idx + self.dim())); } -class histogram_access { -public: +struct python_access { static python::dict - histogram_array_interface(dynamic_histogram<>& self) { + array_interface(dynamic_histogram<>& self) { + const auto& b = self.storage_.buffer_; + if (b.type_.id_ == 5) { + PyErr_SetString(PyExc_RuntimeError, "cannot convert multiprecision storage to numpy array"); + boost::python::throw_error_already_set(); + } + python::dict d; python::list shapes; python::list strides; std::size_t stride = 1; - if (self.depth() == sizeof(detail::weight_t)) { + if (b.type_.id_ == -1) { stride *= sizeof(double); d["typestr"] = python::str("|f") + python::str(stride); strides.append(stride); stride *= 2; shapes.append(2); } else { - stride *= self.depth(); + stride *= b.type_.depth_; d["typestr"] = python::str("|u") + python::str(stride); } for (unsigned i = 0; i < self.dim(); ++i) { @@ -254,7 +259,7 @@ public: stride *= shape(self.axis(i)); } d["shape"] = python::tuple(shapes); - d["data"] = python::make_tuple(reinterpret_cast(self.data()), false); + d["data"] = python::make_tuple(reinterpret_cast(b.ptr_), false); d["strides"] = python::tuple(strides); return d; } @@ -279,9 +284,9 @@ void register_histogram() // shadowed C++ ctors .def(init::axes_t&>()) .add_property("__array_interface__", - &histogram_access::histogram_array_interface) + &python_access::array_interface) .add_property("dim", &dynamic_histogram<>::dim, - "dimensions of the histogram (number of axes)") + "dimensions of the histogram (number of axes)") .def("axis", histogram_axis, ":param int i: index of the axis\n" ":returns: axis object for axis i", diff --git a/test/dynamic_storage_test.cpp b/test/dynamic_storage_test.cpp index 2cac9984..90d4d63d 100644 --- a/test/dynamic_storage_test.cpp +++ b/test/dynamic_storage_test.cpp @@ -52,23 +52,16 @@ void dynamic_storage_increase_and_grow_impl() auto n2 = n; BOOST_CHECK_EQUAL(n.value(0), double(tmax - 1)); BOOST_CHECK_EQUAL(n2.value(0), double(tmax - 1)); + n.increase(0); + n.increase(0); dynamic_storage x(1); x.increase(0); - n.increase(0); n2 += x; - if (sizeof(T) == sizeof(uint64_t)) { - BOOST_CHECK_THROW(n.increase(0), std::overflow_error); - BOOST_CHECK_THROW(n2 += x, std::overflow_error); - } else { - n.increase(0); - n2 += x; - double v = tmax; - ++v; - BOOST_CHECK_EQUAL(n.value(0), v); - BOOST_CHECK_EQUAL(n2.value(0), v); - BOOST_CHECK(!(n.value(0) == double(tmax))); - BOOST_CHECK(!(n2.value(0) == double(tmax))); - } + n2 += x; + double v = tmax; + ++v; + BOOST_CHECK_EQUAL(n.value(0), v); + BOOST_CHECK_EQUAL(n2.value(0), v); } BOOST_AUTO_TEST_CASE(dynamic_storage_increase_and_grow) @@ -89,7 +82,7 @@ BOOST_AUTO_TEST_CASE(dynamic_storage_add_and_grow) BOOST_CHECK_EQUAL(y.value(0), 0.0); a += y; BOOST_CHECK_EQUAL(a.value(0), x); - for (unsigned i = 0; i < 63; ++i) { + for (unsigned i = 0; i < 80; ++i) { a += a; x += x; dynamic_storage b(1);