diff --git a/src/python/histogram.cpp b/src/python/histogram.cpp index 9fcc9eb9..13468031 100644 --- a/src/python/histogram.cpp +++ b/src/python/histogram.cpp @@ -34,12 +34,6 @@ using dynamic_histogram = histogram>; namespace python { -#ifdef HAVE_NUMPY -auto array_cast = [](handle<>& h) { - return downcast(h.get()); -}; -#endif - #ifdef HAVE_NUMPY class access { public: @@ -87,38 +81,38 @@ public: object operator()(const array& b) const { // cannot pass non-existent memory to numpy; make new // zero-initialized uint8 array, and pass it - int dim = len(shapes); + auto dim = len(shapes); npy_intp shapes2[BOOST_HISTOGRAM_AXIS_LIMIT]; - for (int i = 0; i < dim; ++i) { + for (auto i = 0; i < dim; ++i) { shapes2[i] = extract(shapes[i]); } - handle<> a(PyArray_SimpleNew(dim, shapes2, NPY_UINT8)); - for (int i = 0; i < dim; ++i) { - PyArray_STRIDES(array_cast(a))[i] = extract(strides[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 = static_cast(PyArray_DATA(array_cast(a))); + auto *buf = (uint8_t *)PyArray_DATA(a); std::fill(buf, buf + b.size, uint8_t(0)); - PyArray_CLEARFLAGS(array_cast(a), NPY_ARRAY_WRITEABLE); - return object(a); + PyArray_CLEARFLAGS(a, NPY_ARRAY_WRITEABLE); + return object(handle<>((PyObject*)a)); } object operator()(const array& b) const { // cannot pass cpp_int to numpy; make new // double array, fill it and pass it - int dim = len(shapes); + auto dim = len(shapes); npy_intp shapes2[BOOST_HISTOGRAM_AXIS_LIMIT]; - for (int i = 0; i < dim; ++i) { + for (auto i = 0; i < dim; ++i) { shapes2[i] = extract(shapes[i]); } - handle<> a(PyArray_SimpleNew(dim, shapes2, NPY_DOUBLE)); - for (int i = 0; i < dim; ++i) { - PyArray_STRIDES(array_cast(a))[i] = extract(strides[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 = static_cast(PyArray_DATA(array_cast(a))); - for (std::size_t i = 0; i < b.size; ++i) { + auto *buf = (double *)PyArray_DATA(a); + for (auto i = 0ul; i < b.size; ++i) { buf[i] = static_cast(b[i]); } - PyArray_CLEARFLAGS(array_cast(a), NPY_ARRAY_WRITEABLE); - return object(a); + PyArray_CLEARFLAGS(a, NPY_ARRAY_WRITEABLE); + return object(handle<>((PyObject*)a)); } }; @@ -225,89 +219,55 @@ python::object histogram_init(python::tuple args, python::dict kwargs) { return pyinit(h); } -python::object histogram_fill(python::tuple args, python::dict kwargs) { - const unsigned nargs = python::len(args); - dynamic_histogram &self = python::extract(args[0]); +struct fetcher { + char type = 0; + union { + PyArrayObject* a; + double value; + }; - python::object ow; - if (kwargs) { - if (len(kwargs) > 1 || !kwargs.has_key("weight")) { - PyErr_SetString(PyExc_RuntimeError, "only keyword weight allowed"); - python::throw_error_already_set(); - } - ow = kwargs.get("weight"); + ~fetcher() { + if (type == 2) Py_DECREF((PyObject*)a); } + long connect(python::object o) { + python::extract get_double(o); + if (get_double.check()) { + type = 1; + value = get_double(); + return 0; + } #ifdef HAVE_NUMPY - if (nargs == 2) { - python::object o = args[1]; - if (PySequence_Check(o.ptr())) { - // exception is thrown automatically if - python::handle<> a(PyArray_FROM_OTF(o.ptr(), NPY_DOUBLE, NPY_ARRAY_IN_ARRAY)); - - npy_intp *dims = PyArray_DIMS(python::array_cast(a)); - switch (PyArray_NDIM(python::array_cast(a))) { - case 1: - if (self.dim() > 1) { - PyErr_SetString(PyExc_ValueError, "array has to be two-dimensional"); - python::throw_error_already_set(); - } - break; - case 2: - if (self.dim() != dims[1]) { - PyErr_SetString(PyExc_ValueError, - "size of second dimension does not match"); - python::throw_error_already_set(); - } - break; - default: - PyErr_SetString(PyExc_ValueError, "array has wrong dimension"); - python::throw_error_already_set(); - } - - if (!ow.is_none()) { - if (PySequence_Check(ow.ptr())) { - - // exception is thrown automatically if handle below receives null - python::handle<> aw( - PyArray_FROM_OTF(ow.ptr(), NPY_DOUBLE, NPY_ARRAY_IN_ARRAY)); - - if (PyArray_NDIM(python::array_cast(aw)) != 1) { - PyErr_SetString(PyExc_ValueError, - "array has to be one-dimensional"); - python::throw_error_already_set(); - } - - if (PyArray_DIMS(python::array_cast(aw))[0] != dims[0]) { - PyErr_SetString(PyExc_ValueError, "sizes do not match"); - python::throw_error_already_set(); - } - - for (unsigned i = 0; i < dims[0]; ++i) { - double *v = reinterpret_cast(PyArray_GETPTR1(python::array_cast(a), i)); - double *w = reinterpret_cast(PyArray_GETPTR1(python::array_cast(aw), i)); - self.fill(v, v + self.dim(), weight(*w)); - } - - } else { - PyErr_SetString(PyExc_ValueError, "weight is not a sequence"); - python::throw_error_already_set(); - } - } else { - for (unsigned i = 0; i < dims[0]; ++i) { - double *v = reinterpret_cast(PyArray_GETPTR1(python::array_cast(a), i)); - self.fill(v, v + self.dim()); - } - } - - return python::object(); + 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; } -#endif /* HAVE_NUMPY */ + + 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; + }; +}; + +python::object histogram_fill(python::tuple args, python::dict kwargs) { + const auto nargs = python::len(args); + dynamic_histogram &self = python::extract(args[0]); const unsigned dim = nargs - 1; if (dim != self.dim()) { - PyErr_SetString(PyExc_RuntimeError, "wrong number of arguments"); + PyErr_SetString(PyExc_ValueError, "number of arguments and dimension do not match"); python::throw_error_already_set(); } @@ -318,15 +278,48 @@ python::object histogram_fill(python::tuple args, python::dict kwargs) { python::throw_error_already_set(); } - double v[BOOST_HISTOGRAM_AXIS_LIMIT]; - for (unsigned i = 0; i < dim; ++i) - v[i] = python::extract(args[1 + i]); + fetcher 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); + if (on > 0) { + if (n && on != n) { + PyErr_SetString(PyExc_ValueError, "lengths of sequences do not match"); + python::throw_error_already_set(); + } + n = on; + } + } - if (ow.is_none()) { - self.fill(v, v + self.dim()); - } else { - const double w = python::extract(ow); - self.fill(v, v + self.dim(), weight(w)); + double v[BOOST_HISTOGRAM_AXIS_LIMIT]; + fetcher 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); + if (on > 0) { + if (n && on != n) { + PyErr_SetString(PyExc_ValueError, "lengths of argument sequences do not match"); + python::throw_error_already_set(); + } + n = on; + } + } + + 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); + } else { + self.fill(v, v + dim, weight(fetch_weight(i))); + } } return python::object(); diff --git a/test/python_suite_test.py b/test/python_suite_test.py index 6b3eb6a8..ef3d11cf 100644 --- a/test/python_suite_test.py +++ b/test/python_suite_test.py @@ -349,9 +349,9 @@ class test_histogram(unittest.TestCase): h0 = histogram(integer(-1, 1, uoflow=False)) h1 = histogram(integer(-1, 1, uoflow=True)) for h in (h0, h1): - with self.assertRaises(RuntimeError): + with self.assertRaises(ValueError): h.fill() - with self.assertRaises(RuntimeError): + with self.assertRaises(ValueError): h.fill(1, 2) h.fill(-10) h.fill(-1) @@ -704,8 +704,9 @@ class test_histogram(unittest.TestCase): @unittest.skipUnless(have_numpy, "requires build with numpy-support") def test_fill_with_numpy_array_0(self): + ar = lambda *args: numpy.array(args) a = histogram(integer(0, 2, uoflow=False)) - a.fill(numpy.array([-1, 0, 1, 2, 1, 4])) + a.fill(ar(-1, 0, 1, 2, 1, 4)) a.fill((-1, 0)) a.fill([1, 2]) self.assertEqual(a.value(0), 2) @@ -715,15 +716,14 @@ class test_histogram(unittest.TestCase): with self.assertRaises(ValueError): a.fill(numpy.empty((2, 2))) with self.assertRaises(ValueError): - a.fill(numpy.empty((1, 2, 2))) - with self.assertRaises(RuntimeError): - a.fill(numpy.empty((2, 1)), 1) + a.fill(numpy.empty(2), 1) with self.assertRaises(ValueError): a.fill("abc") a = histogram(integer(0, 1, uoflow=False), regular(2, 0, 2, uoflow=False)) - a.fill(numpy.array([[-1., -1.], [0., 1.], [1., 0.1]])) + a.fill(ar(-1, 0, 1), + ar(-1., 1., 0.1)) self.assertEqual(a.value(0, 0), 0) self.assertEqual(a.value(0, 1), 1) self.assertEqual(a.value(1, 0), 1) @@ -742,11 +742,12 @@ class test_histogram(unittest.TestCase): @unittest.skipUnless(have_numpy, "requires build with numpy-support") def test_fill_with_numpy_array_1(self): + ar = lambda *args: numpy.array(args) a = histogram(integer(0, 2, uoflow=True)) v = numpy.array([-1, 0, 1, 2, 3, 4]) w = numpy.array([ 2, 3, 4, 5, 6, 7]) a.fill(v, weight=w) - a.fill([0, 1], weight=[2.0, 3.0]) + a.fill([0, 1], weight=[2, 3]) self.assertEqual(a.value(-1), 2) self.assertEqual(a.value(0), 5) self.assertEqual(a.value(1), 7) @@ -755,6 +756,11 @@ class test_histogram(unittest.TestCase): self.assertEqual(a.variance(0), 13) self.assertEqual(a.variance(1), 25) self.assertEqual(a.variance(2), 25) + a.fill([1, 2], weight=1) + a.fill(0, weight=[1, 2]) + self.assertEqual(a.value(0), 8) + self.assertEqual(a.value(1), 8) + self.assertEqual(a.value(2), 6) with self.assertRaises(RuntimeError): a.fill([1, 2], foo=[1, 1]) @@ -762,8 +768,6 @@ class test_histogram(unittest.TestCase): a.fill([1, 2], weight=[1]) with self.assertRaises(ValueError): a.fill([1, 2], weight="ab") - with self.assertRaises(ValueError): - a.fill([1, 2], weight=1) with self.assertRaises(RuntimeError): a.fill([1, 2], weight=[1, 1], foo=1) with self.assertRaises(ValueError): @@ -771,7 +775,7 @@ class test_histogram(unittest.TestCase): a = histogram(integer(0, 1, uoflow=False), regular(2, 0, 2, uoflow=False)) - a.fill(numpy.array([[-1., -1.], [0., 1.], [1., 0.1]])) + a.fill(ar(-1, 0, 1), ar(-1, 1, 0.1)) self.assertEqual(a.value(0, 0), 0) self.assertEqual(a.value(0, 1), 1) self.assertEqual(a.value(1, 0), 1) diff --git a/test/speed_numpy.py b/test/speed_numpy.py index 1edb1db5..99e33dee 100755 --- a/test/speed_numpy.py +++ b/test/speed_numpy.py @@ -8,7 +8,8 @@ import numpy as np from timeit import default_timer as timer -from histogram import histogram, regular_axis +from histogram import histogram +from histogram.axis import regular def compare_1d(n, distrib): if distrib == 0: @@ -24,7 +25,7 @@ def compare_1d(n, distrib): t = timer() - t best_numpy = min(t, best_numpy) - h = histogram(regular_axis(100, 0, 1)) + h = histogram(regular(100, 0, 1)) t = timer() h.fill(r) t = timer() - t @@ -50,9 +51,9 @@ def compare_2d(n, distrib): t = timer() - t best_numpy = min(t, best_numpy) - h = histogram(regular_axis(100, 0, 1), regular_axis(100, 0, 1)) + h = histogram(regular(100, 0, 1), regular(100, 0, 1)) t = timer() - h.fill(r) + h.fill(r[:,0], r[:,1]) t = timer() - t best_boost = min(t, best_boost) assert(np.all(w == np.array(h)[:-2,:-2])) @@ -78,11 +79,11 @@ def compare_3d(n, distrib): t = timer() - t best_numpy = min(t, best_numpy) - h = histogram(regular_axis(100, 0, 1), - regular_axis(100, 0, 1), - regular_axis(100, 0, 1)) + h = histogram(regular(100, 0, 1), + regular(100, 0, 1), + regular(100, 0, 1)) t = timer() - h.fill(r) + h.fill(r[:,0], r[:,1], r[:,2]) t = timer() - t best_boost = min(t, best_boost) assert(np.all(w == np.array(h)[:-2,:-2,:-2])) @@ -112,14 +113,14 @@ def compare_6d(n, distrib): t = timer() - t best_numpy = min(t, best_numpy) - h = histogram(regular_axis(10, 0, 1), - regular_axis(10, 0, 1), - regular_axis(10, 0, 1), - regular_axis(10, 0, 1), - regular_axis(10, 0, 1), - regular_axis(10, 0, 1)) + h = histogram(regular(10, 0, 1), + regular(10, 0, 1), + regular(10, 0, 1), + regular(10, 0, 1), + regular(10, 0, 1), + regular(10, 0, 1)) t = timer() - h.fill(r) + h.fill(r[:,0], r[:,1], r[:,2], r[:,3], r[:,4], r[:,5]) t = timer() - t best_boost = min(t, best_boost) assert(np.all(w == np.array(h)[:-2,:-2,:-2,:-2,:-2,:-2]))