From e12856b03d8336eeabf15cd150e84c4bed69d236 Mon Sep 17 00:00:00 2001 From: Hans Dembinski Date: Fri, 3 Nov 2017 00:13:43 +0100 Subject: [PATCH] replaced CApi numpy with boost numpy --- src/python/axis.cpp | 33 +++------ src/python/histogram.cpp | 139 +++++++++++++++----------------------- src/python/module.cpp | 20 ++---- test/python_suite_test.py | 22 ++++-- 4 files changed, 86 insertions(+), 128 deletions(-) diff --git a/src/python/axis.cpp b/src/python/axis.cpp index 36f03231..0fedf3c5 100644 --- a/src/python/axis.cpp +++ b/src/python/axis.cpp @@ -12,10 +12,7 @@ #include #include #ifdef HAVE_NUMPY -#define NO_IMPORT_ARRAY -#define PY_ARRAY_UNIQUE_SYMBOL boost_histogram_ARRAY_API -#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION -#include +#include #endif #include #include @@ -23,6 +20,8 @@ #include #include +namespace np = boost::python::numpy; + namespace boost { namespace histogram { @@ -132,22 +131,17 @@ template python::str axis_get_label(const T& t) { #ifdef HAVE_NUMPY template python::object axis_array_interface(const Axis& axis) { + using T = typename decay::type; python::dict d; auto shape = python::make_tuple(axis.size()+1); d["shape"] = shape; - // d["typestr"] = dtype_typestr(); - d["typestr"] = "|f8"; + d["typestr"] = python::dtype_typestr(); // make new array, and pass it to Python - auto dim = 1; - npy_intp shapes2[1] = { axis.size()+1 }; - auto *a = (PyArrayObject*)PyArray_SimpleNew(dim, shapes2, NPY_DOUBLE); - auto *buf = (double *)PyArray_DATA(a); - PyArray_CLEARFLAGS(a, NPY_ARRAY_WRITEABLE); - // auto a = python::numpy::empty(shape, python::numpy::dtype::get_builtin()); - // auto buf = reinterpret_cast(axis.get_data()); + auto a = np::empty(shape, np::dtype::get_builtin()); + auto buf = reinterpret_cast(a.get_data()); for (auto i = 0; i < axis.size()+1; ++i) buf[i] = axis[i].lower(); - d["data"] = python::object(python::handle<>((PyObject*)a)); + d["data"] = a; d["version"] = 3; return d; } @@ -158,16 +152,11 @@ template <> python::object axis_array_interface>(const axis::ca d["shape"] = shape; d["typestr"] = python::dtype_typestr(); // make new array, and pass it to Python - auto dim = 1; - npy_intp shapes2[1] = { axis.size() }; - auto *a = (PyArrayObject*)PyArray_SimpleNew(dim, shapes2, NPY_INT); - auto *buf = (int *)PyArray_DATA(a); - PyArray_CLEARFLAGS(a, NPY_ARRAY_WRITEABLE); - // auto a = python::numpy::empty(shape, python::numpy::dtype::get_builtin()); - // auto buf = reinterpret_cast(axis.get_data()); + auto a = np::empty(shape, np::dtype::get_builtin()); + auto buf = reinterpret_cast(a.get_data()); for (auto i = 0; i < axis.size(); ++i) buf[i] = axis[i]; - d["data"] = python::object(python::handle<>((PyObject*)a)); + d["data"] = a; d["version"] = 3; return d; } diff --git a/src/python/histogram.cpp b/src/python/histogram.cpp index d27ed3c1..d5dfbdcb 100644 --- a/src/python/histogram.cpp +++ b/src/python/histogram.cpp @@ -17,16 +17,16 @@ #include #include #ifdef HAVE_NUMPY -#define NO_IMPORT_ARRAY -#define PY_ARRAY_UNIQUE_SYMBOL boost_histogram_ARRAY_API -#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION -#include +#include #endif +#include #ifndef BOOST_HISTOGRAM_AXIS_LIMIT #define BOOST_HISTOGRAM_AXIS_LIMIT 32 #endif +namespace np = boost::python::numpy; + namespace boost { namespace histogram { @@ -78,38 +78,18 @@ public: object operator()(const array& b) const { // cannot pass non-existent memory to numpy; make new // zero-initialized uint8 array, and pass it - auto dim = len(shapes); - npy_intp shapes2[BOOST_HISTOGRAM_AXIS_LIMIT]; - for (auto i = 0; i < dim; ++i) { - shapes2[i] = extract(shapes[i]); - } - auto *a = (PyArrayObject*)PyArray_SimpleNew(dim, shapes2, NPY_UINT8); - for (auto i = 0; i < dim; ++i) { - PyArray_STRIDES(a)[i] = extract(strides[i]); - } - auto *buf = (uint8_t *)PyArray_DATA(a); - std::fill(buf, buf + b.size, uint8_t(0)); - PyArray_CLEARFLAGS(a, NPY_ARRAY_WRITEABLE); - return object(handle<>((PyObject*)a)); + return np::zeros(tuple(shapes), np::dtype::get_builtin()); } object operator()(const array& b) const { // cannot pass cpp_int to numpy; make new // double array, fill it and pass it - auto dim = len(shapes); - npy_intp shapes2[BOOST_HISTOGRAM_AXIS_LIMIT]; - for (auto i = 0; i < dim; ++i) { - shapes2[i] = extract(shapes[i]); - } - auto *a = (PyArrayObject*)PyArray_SimpleNew(dim, shapes2, NPY_DOUBLE); - for (auto i = 0; i < dim; ++i) { - PyArray_STRIDES(a)[i] = extract(strides[i]); - } - auto *buf = (double *)PyArray_DATA(a); - for (auto i = 0ul; i < b.size; ++i) { + auto a = np::empty(tuple(shapes), np::dtype::get_builtin()); + for (auto i = 0l, n = len(shapes); i < n; ++i) + const_cast(a.get_strides())[i] = python::extract(strides[i]); + auto *buf = (double *)a.get_data(); + for (auto i = 0ul; i < b.size; ++i) buf[i] = static_cast(b[i]); - } - PyArray_CLEARFLAGS(a, NPY_ARRAY_WRITEABLE); - return object(handle<>((PyObject*)a)); + return a; } }; @@ -204,51 +184,42 @@ python::object histogram_init(python::tuple args, python::dict kwargs) { } struct fetcher { - char type = 0; - union { -#ifdef HAVE_NUMPY - PyArrayObject* a; -#endif - double value; - }; - -#ifdef HAVE_NUMPY - ~fetcher() { - if (type == 2) Py_DECREF((PyObject*)a); - } -#endif - - long connect(python::object o) { - python::extract get_double(o); - if (get_double.check()) { - type = 1; - value = get_double(); - return 0; - } -#ifdef HAVE_NUMPY - a = (PyArrayObject*)PyArray_FROMANY(o.ptr(), NPY_DOUBLE, 1, 1, NPY_ARRAY_CARRAY); - if (a) { - type = 2; - return PyArray_SHAPE(a)[0]; - } -#endif - PyErr_SetString(PyExc_ValueError, "cannot convert argument"); - python::throw_error_already_set(); - return 0; - } - - bool empty() const { return type == 0; } - - double operator[](long i) const { -#ifdef HAVE_NUMPY - if (type == 2) { - return *static_cast(PyArray_GETPTR1(a, i)); - } -#endif - return value; - }; + virtual ~fetcher() {} + virtual double at(long) const = 0; + long n = 0; }; +#ifdef HAVE_NUMPY +struct fetcher_seq : public fetcher { + fetcher_seq(python::object o) + : array(np::from_object(o, np::dtype::get_builtin(), 1)) { + fetcher::n = array.shape(0); + } + ~fetcher_seq() {} + double at(long i) const { + return reinterpret_cast(array.get_data())[i]; + } + np::ndarray array; +}; +#endif + +struct fetcher_val : public fetcher { + fetcher_val(double val) + : value(val) {} + double at(long) const { return value; } + double value; +}; + +std::unique_ptr make_fetcher(python::object o) { + python::extract get_double(o); + if (get_double.check()) + return std::unique_ptr(new fetcher_val(get_double())); +#ifdef HAVE_NUMPY + return std::unique_ptr(new fetcher_seq(o)); +#endif + throw std::invalid_argument("python object is neither sequence nor number"); +} + python::object histogram_fill(python::tuple args, python::dict kwargs) { const auto nargs = python::len(args); dynamic_histogram &self = python::extract(args[0]); @@ -266,11 +237,11 @@ python::object histogram_fill(python::tuple args, python::dict kwargs) { python::throw_error_already_set(); } - fetcher fetch[BOOST_HISTOGRAM_AXIS_LIMIT]; + std::unique_ptr fetch[BOOST_HISTOGRAM_AXIS_LIMIT]; long n = 0; for (auto d = 0u; d < dim; ++d) { - python::object o = python::extract(args[1 + d]); - const auto on = fetch[d].connect(o); + fetch[d] = make_fetcher(args[1 + d]); + const auto on = fetch[d]->n; if (on > 0) { if (n && on != n) { PyErr_SetString(PyExc_ValueError, "lengths of sequences do not match"); @@ -280,15 +251,15 @@ python::object histogram_fill(python::tuple args, python::dict kwargs) { } } - fetcher fetch_weight; + std::unique_ptr fetch_weight; const auto nkwargs = python::len(kwargs); if (nkwargs > 0) { if (nkwargs > 1 || !kwargs.has_key("weight")) { PyErr_SetString(PyExc_RuntimeError, "only keyword weight allowed"); python::throw_error_already_set(); } - python::object o = kwargs.get("weight"); - const auto on = fetch_weight.connect(o); + fetch_weight = make_fetcher(kwargs.get("weight")); + const auto on = fetch_weight->n; if (on > 0) { if (n && on != n) { PyErr_SetString(PyExc_ValueError, "length of weight sequence does not match"); @@ -302,11 +273,11 @@ python::object histogram_fill(python::tuple args, python::dict kwargs) { if (!n) ++n; for (auto i = 0l; i < n; ++i) { for (auto d = 0u; d < dim; ++d) - v[d] = fetch[d][i]; - if (fetch_weight.empty()) { - self.fill(v, v + dim); + v[d] = fetch[d]->at(i); + if (fetch_weight) { + self.fill(v, v + dim, weight(fetch_weight->at(i))); } else { - self.fill(v, v + dim, weight(fetch_weight[i])); + self.fill(v, v + dim); } } diff --git a/src/python/module.cpp b/src/python/module.cpp index 5c30b7e9..0fc9adb4 100644 --- a/src/python/module.cpp +++ b/src/python/module.cpp @@ -8,19 +8,7 @@ #include #include #ifdef HAVE_NUMPY -#define PY_ARRAY_UNIQUE_SYMBOL boost_histogram_ARRAY_API -#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION -extern "C" { -#include -#if PY_MAJOR_VERSION >= 3 -static void *init_numpy() { - import_array(); - return NULL; -} -#else -static void init_numpy() { import_array(); } -#endif -} +#include #endif namespace boost { @@ -31,10 +19,10 @@ void register_histogram(); } BOOST_PYTHON_MODULE(histogram) { -#ifdef HAVE_NUMPY - init_numpy(); -#endif using namespace boost::python; +#ifdef HAVE_NUMPY + numpy::initialize(); +#endif scope current; object axis_module = object( borrowed(PyImport_AddModule("histogram.axis")) diff --git a/test/python_suite_test.py b/test/python_suite_test.py index 3a97b074..d8137351 100644 --- a/test/python_suite_test.py +++ b/test/python_suite_test.py @@ -701,15 +701,25 @@ class test_histogram(unittest.TestCase): @unittest.skipUnless(have_numpy, "requires build with numpy-support") def test_numpy_conversion_5(self): - a = histogram(integer(0, 2, uoflow=False)) - a.fill(0) - for i in xrange(80): + a = histogram(integer(0, 3, uoflow=False), + integer(0, 2, uoflow=False)) + a.fill(0, 0) + for i in range(80): a += a # a now holds a multiprecision type + a.fill(1, 0) + for i in range(2): a.fill(2, 0) + for i in range(3): a.fill(0, 1) + for i in range(4): a.fill(1, 1) + for i in range(5): a.fill(2, 1) a1 = numpy.asarray(a) - self.assertEqual(a1.shape, (2,)) - self.assertEqual(a1[0], float(2 ** 80)) - self.assertEqual(a1[1], 0) + self.assertEqual(a1.shape, (3, 2)) + self.assertEqual(a1[0, 0], float(2 ** 80)) + self.assertEqual(a1[1, 0], 1) + self.assertEqual(a1[2, 0], 2) + self.assertEqual(a1[0, 1], 3) + self.assertEqual(a1[1, 1], 4) + self.assertEqual(a1[2, 1], 5) @unittest.skipUnless(have_numpy, "requires build with numpy-support") def test_numpy_conversion_6(self):