diff --git a/include/boost/histogram/basic_histogram.hpp b/include/boost/histogram/basic_histogram.hpp new file mode 100644 index 00000000..c2d2d305 --- /dev/null +++ b/include/boost/histogram/basic_histogram.hpp @@ -0,0 +1,125 @@ +#ifndef _BOOST_HISTOGRAM_BASE_HPP_ +#define _BOOST_HISTOGRAM_BASE_HPP_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#define BOOST_HISTOGRAM_AXIS_LIMIT 16 + +namespace boost { +namespace histogram { + +// holds collection of axis instances and computes the internal index +class basic_histogram { +public: + typedef container::static_vector axes_type; + typedef uintptr_t size_type; + + basic_histogram(const basic_histogram&); + basic_histogram& operator=(const basic_histogram&); + ~basic_histogram() {} + + unsigned dim() const { return axes_.size(); } + + template + T& axis(unsigned i) + { if (is_same::value) return axes_[i]; + return boost::get(axes_[i]); } + + template + const T& axis(unsigned i) const + { if (is_same::value) return axes_[i]; + return boost::get(axes_[i]); } + + unsigned bins(unsigned i) const { return size_[i]; } + unsigned shape(unsigned i) const { return size_[i] + 2 * uoflow_[i]; } + +protected: + basic_histogram() {} + explicit basic_histogram(const axes_type& axes); + +#define BOOST_HISTOGRAM_BASE_APPEND(z, n, unused) axes_.push_back(a ## n); +#define BOOST_HISTOGRAM_BASE_CTOR(z, n, unused) \ + basic_histogram( BOOST_PP_ENUM_PARAMS_Z(z, n, const axis_type& a) ) \ + { \ + axes_.reserve(n); \ + BOOST_PP_REPEAT(n, BOOST_HISTOGRAM_BASE_APPEND, unused) \ + update_buffers(); \ + } + +// generates constructors taking 1 to AXIS_LIMIT arguments +BOOST_PP_REPEAT_FROM_TO(1, BOOST_HISTOGRAM_AXIS_LIMIT, BOOST_HISTOGRAM_BASE_CTOR, nil) + + bool operator==(const basic_histogram&) const; + bool operator!=(const basic_histogram& o) const + { return !operator==(o); } + + template + inline + size_type pos(const Array& v) const { + int idx[BOOST_HISTOGRAM_AXIS_LIMIT]; + for (unsigned i = 0, n = axes_.size(); i < n; ++i) + idx[i] = apply_visitor(detail::index_visitor(v[i]), axes_[i]); + return linearize(idx); + } + + inline + size_type linearize(const int* idx) const { + size_type stride = 1, k = 0, i = axes_.size(); + while (i--) { + const int size = size_[i]; + const int range = size + 2 * uoflow_[i]; + int j = idx[i]; + // the following three lines work for any overflow setting + j += (j < 0) * (size + 2); // wrap around if j < 0 + if (j >= range) + return size_type(-1); // indicate out of range + k += j * stride; + stride *= range; + } + return k; + } + + // compute the number of fields needed for storage + size_type field_count() const; + +private: + axes_type axes_; + + // internal buffers + int32_t size_[BOOST_HISTOGRAM_AXIS_LIMIT]; + std::bitset uoflow_; + + void update_buffers(); ///< fills size_ and uoflow_ + + friend class serialization::access; + template + void serialize(Archive& ar, unsigned version) + { + using namespace serialization; + unsigned size = axes_.size(); + ar & size; + if (Archive::is_loading::value) { + axes_.resize(size); + ar & serialization::make_array(&axes_[0], size); + update_buffers(); + } else { + ar & serialization::make_array(&axes_[0], size); + } + } +}; + +} +} + +#endif diff --git a/src/basic_histogram.cpp b/src/basic_histogram.cpp new file mode 100644 index 00000000..abd00553 --- /dev/null +++ b/src/basic_histogram.cpp @@ -0,0 +1,63 @@ +#include +#include +#include + +namespace boost { +namespace histogram { + +basic_histogram::basic_histogram(const basic_histogram& o) : + axes_(o.axes_) +{ + update_buffers(); +} + +basic_histogram::basic_histogram(const axes_type& axes) : + axes_(axes) +{ + if (axes_.size() > BOOST_HISTOGRAM_AXIS_LIMIT) + throw std::invalid_argument("too many axes"); + update_buffers(); +} + +basic_histogram& +basic_histogram::operator=(const basic_histogram& o) +{ + axes_ = o.axes_; + update_buffers(); + return *this; +} + +bool +basic_histogram::operator==(const basic_histogram& o) + const +{ + if (axes_.size() != o.axes_.size()) + return false; + for (unsigned i = 0; i < axes_.size(); ++i) { + if (!apply_visitor(detail::cmp_visitor(), axes_[i], o.axes_[i])) + return false; + } + return true; +} + +basic_histogram::size_type +basic_histogram::field_count() + const +{ + size_type fc = 1; + for (unsigned i = 0, n = axes_.size(); i < n; ++i) + fc *= apply_visitor(detail::fields_visitor(), axes_[i]); + return fc; +} + +void +basic_histogram::update_buffers() +{ + for (unsigned i = 0, n = axes_.size(); i < n; ++i) { + size_[i] = apply_visitor(detail::bins_visitor(), axes_[i]); + uoflow_[i] = apply_visitor(detail::uoflow_visitor(), axes_[i]); + } +} + +} +} diff --git a/src/python/basic_histogram.cpp b/src/python/basic_histogram.cpp new file mode 100644 index 00000000..56da56c5 --- /dev/null +++ b/src/python/basic_histogram.cpp @@ -0,0 +1,325 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace histogram { + +namespace { + +python::object +variable_axis_init(python::tuple args, python::dict kwargs) { + using namespace python; + using python::tuple; + + object self = args[0]; + object pyinit = self.attr("__init__"); + + if (len(args) < 2) { + PyErr_SetString(PyExc_TypeError, "require at least two arguments"); + throw_error_already_set(); + } + + std::vector v; + for (int i = 1, n = len(args); i < n; ++i) { + v.push_back(extract(args[i])); + } + + std::string label; + bool uoflow = true; + while (len(kwargs) > 0) { + tuple kv = kwargs.popitem(); + std::string k = extract(kv[0]); + object v = kv[1]; + if (k == "label") + label = extract(v); + else if (k == "uoflow") + uoflow = extract(v); + else { + std::stringstream s; + s << "keyword " << k << " not recognized"; + PyErr_SetString(PyExc_KeyError, s.str().c_str()); + throw_error_already_set(); + } + } + return pyinit(v, label, uoflow); +} + +python::object +category_axis_init(python::tuple args, python::dict kwargs) { + using namespace python; + + object self = args[0]; + object pyinit = self.attr("__init__"); + + if (len(args) == 1) { + PyErr_SetString(PyExc_TypeError, "require at least one argument"); + throw_error_already_set(); + } + + if (len(kwargs) > 0) { + PyErr_SetString(PyExc_TypeError, "unknown keyword argument"); + throw_error_already_set(); + } + + if (len(args) == 2) { + extract es(args[1]); + if (es.check()) + pyinit(es); + else { + PyErr_SetString(PyExc_TypeError, "require one or several string arguments"); + throw_error_already_set(); + } + } + + std::vector c; + for (int i = 1, n = len(args); i < n; ++i) + c.push_back(extract(args[i])); + + return pyinit(c); +} + +// python::object +// regular_axis_data(const regular_axis& self) { +// #if HAVE_NUMPY +// py_intp dims[1] = { self.size() + 1 }; +// PyArrayObject* a = (PyArrayObject*)PyArray_SimpleNew(1, dims, NPY_DOUBLE); +// for (unsigned i = 0; i < self.size() + 1; ++i) { +// double* pv = (double*)PyArray_GETPTR1(a, i); +// *pv = self.left(i); +// } +// return python::object(handle<>((PyObject*)a)); +// #else +// list v; +// for (int i = 0; i < self.size() + 1; ++i) +// v.append(self.left(i)); +// return v; +// #endif +// } + +struct axis_visitor : public static_visitor +{ + template + python::object operator()(const T& t) const { return python::object(T(t)); } +}; + +template +unsigned +axis_len(const T& t) { + return t.bins() + 1; +} + +template <> +unsigned +axis_len(const category_axis& t) { + return t.bins(); +} + +template <> +unsigned +axis_len(const integer_axis& t) { + return t.bins(); +} + +template +typename T::value_type +axis_getitem(const T& t, int i) { + if (i == axis_len(t)) { + PyErr_SetString(PyExc_StopIteration, "no more"); + python::throw_error_already_set(); + } + return t[i]; +} + +template +struct has_index_method +{ + struct yes { char x[1]; }; + struct no { char x[2]; }; + template struct SFINAE {}; + template static yes test( SFINAE* ); + template static no test( ... ); + enum { value = sizeof(test(0)) == sizeof(yes) }; +}; + +std::string escape(const std::string& s) { + std::ostringstream os; + os << '\''; + for (unsigned i = 0; i < s.size(); ++i) { + const char c = s[i]; + if (c == '\'' && (i == 0 || s[i - 1] != '\\')) + os << "\\\'"; + else + os << c; + } + os << '\''; + return os.str(); +} + +std::string axis_repr(const regular_axis& a) { + std::stringstream s; + s << "regular_axis(" << a.bins() << ", " << a[0] << ", " << a[a.bins()]; + if (!a.label().empty()) + s << ", label=" << escape(a.label()); + if (!a.uoflow()) + s << ", uoflow=False"; + s << ")"; + return s.str(); +} + +std::string axis_repr(const polar_axis& a) { + std::stringstream s; + s << "polar_axis(" << a.bins(); + if (a[0] != 0.0) + s << ", " << a[0]; + if (!a.label().empty()) + s << ", label=" << escape(a.label()); + s << ")"; + return s.str(); +} + +std::string axis_repr(const variable_axis& a) { + std::stringstream s; + s << "variable_axis(" << a[0]; + for (int i = 1; i <= a.bins(); ++i) + s << ", " << a.left(i); + if (!a.label().empty()) + s << ", label=" << escape(a.label()); + if (!a.uoflow()) + s << ", uoflow=False"; + s << ")"; + return s.str(); +} + +std::string axis_repr(const category_axis& a) { + std::stringstream s; + s << "category_axis("; + for (int i = 0; i < a.bins(); ++i) + s << escape(a[i]) << (i == (a.bins() - 1)? ")" : ", "); + return s.str(); +} + +std::string axis_repr(const integer_axis& a) { + std::stringstream s; + s << "integer_axis(" << a[0] << ", " << a[a.bins() - 1]; + if (!a.label().empty()) + s << ", label=" << escape(a.label()); + if (!a.uoflow()) + s << ", uoflow=False"; + s << ")"; + return s.str(); +} + +template +struct axis_suite : public python::def_visitor > { + + template + static + typename enable_if_c::value, void>::type + add_axis_index(Class& cl) { + cl.def("index", &U::index); + } + + template + static + typename disable_if_c::value, void>::type + add_axis_index(Class& cl) {} + + template + static + typename enable_if, void>::type + add_axis_label(Class& cl) { + cl.add_property("label", + python::make_function((const std::string&(U::*)() const) &U::label, + python::return_value_policy()), + (void(U::*)(const std::string&)) &U::label); + } + + template + static + typename disable_if, void>::type + add_axis_label(Class& cl) {} + + template + static void + visit(Class& cl) + { + cl.add_property("bins", &T::bins); + add_axis_index(cl); + add_axis_label(cl); + cl.def("__len__", axis_len); + cl.def("__getitem__", axis_getitem); + cl.def("__repr__", (std::string (*)(const T&)) &axis_repr); + cl.def(python::self == python::self); + } +}; + +python::object +basic_histogram_axis(const basic_histogram& self, unsigned i) +{ + return apply_visitor(axis_visitor(), self.axis(i)); +} + +} // namespace + +void register_basic_histogram() { + using namespace python; + using python::arg; + + // used to pass arguments from raw python init to specialized C++ constructors + class_ >("vector_double", no_init); + class_ >("vector_string", no_init); + class_("axes", no_init); + + class_("regular_axis", no_init) + .def(init( + (arg("bin"), arg("min"), arg("max"), + arg("label") = std::string(), + arg("uoflow") = true))) + .def(axis_suite()) + ; + + class_("polar_axis", no_init) + .def(init( + (arg("bin"), arg("start") = 0.0, + arg("label") = std::string()))) + .def(axis_suite()) + ; + + class_("variable_axis", no_init) + .def("__init__", raw_function(variable_axis_init)) + .def(init, std::string, bool>()) + .def(axis_suite()) + ; + + class_("category_axis", no_init) + .def("__init__", raw_function(category_axis_init)) + .def(init()) + .def(init >()) + .def(axis_suite()) + ; + + class_("integer_axis", no_init) + .def(init( + (arg("min"), arg("max"), + arg("label") = std::string(), + arg("uoflow") = true))) + .def(axis_suite()) + ; + + class_("basic_histogram", no_init) + .add_property("dim", &basic_histogram::dim) + .def("bins", &basic_histogram::bins) + .def("shape", &basic_histogram::shape) + .def("axis", basic_histogram_axis) + ; +} + +} +}