From 158c7788d1fd4e0c75a66071d37ce2c5e2cfcee4 Mon Sep 17 00:00:00 2001 From: Thomas Date: Mon, 27 Jul 2026 10:58:37 +0200 Subject: [PATCH] feat: back the statically sized scalars with a NumPy dtype Arrays of scalars were object arrays: a pointer per element and a Python object behind each one, so every NumPy operation was an object loop. The static variants have a fixed element size and, since #82, are exactly their data and trivially copyable, which is what a dtype needs. Measured over 10k elements of DD3Scalar: object dtype a + b 349.5 ns 2.5 ns 138x a * b 354.3 ns 3.9 ns 92x np.sqrt(a) 349.2 ns 2.8 ns 123x np.sin(a) 366.4 ns 7.0 ns 52x np.sum(a) 353.1 ns 3.8 ns 94x a @ b 673.0 ns 2.7 ns 246x One dtype per static type, generated from the same template as the bindings, which is what #81 was the prerequisite for. No parametric dtype is needed: the size follows from the type. The dynamic variants hold a std::vector and stay object arrays. 27 ufunc loops per type: the four arithmetic ones, negative and positive, arctan2, hypot, matmul, vecdot, and 19 unary mathematical functions. All of them copy through local values, which makes them safe for unaligned data at no cost. np.dot and np.linalg.norm do not work on these arrays and cannot be made to. np.dot is not a ufunc but an __array_function__ dispatcher: it dispatches on the array type rather than the dtype, and our arrays are plain ndarrays, so the protocol never fires. It then refuses anything that is not a native or an old-style dtype, before consulting any slot -- the legacy dotfunc slot exists and was tried, and is unreachable. Registering as an old-style dtype would work but builds the feature on the API NumPy is removing. `@`, np.matmul and np.vecdot compute the same and work on object arrays too, so the three tests that used np.dot now use `@`. Requirements NumPy enforces at registration, none of them documented -- each surfaced only as an error and cost an iteration: __repr__ and __str__ are mandatory, so is a cast between the DType's own instances, that cast has to handle unaligned data and declare NPY_METH_SUPPORTS_UNALIGNED, and the DType needs its own tp_new. A null type object is rejected, and so is sharing one type object between DTypes, which is why the dtype has to hang off the scalar class and picking it up in np.array is not optional. The loop registration allocates its spec on the heap for the lifetime of the module. Function-local statics would have been shared by every call with the same type and arity, and NumPy may keep pointers into the spec. Python 239 passed, 66 of them new in test_dtype.py C++ 102 test cases, 1310 assertions, Release with -Werror and under clang with ASan+UBSan all 34 generated dtypes checked for element size, boxing and arithmetic The new tests were validated by mutation, and the first version failed that: comparing a ufunc on an array against the same ufunc on scalars passes whatever the loop computes, because both go through it. They now compare against the scalar operators, and injecting subtract-computes-plus and arccosh-computes-arcsinh fails four of them. Build time is unchanged at 37.4 s; the module grows from 2.8 to 3.4 MB. --- README.md | 34 ++ pyproject.toml | 2 +- python/CMakeLists.txt | 15 +- python/src/bind_ddscalar.cpp.in | 7 +- python/src/common.h | 2 + python/src/dtype.h | 566 ++++++++++++++++++++++++++++++++ python/src/main.cpp | 7 + python/tests/test_DDScalar.py | 9 +- python/tests/test_dtype.py | 175 ++++++++++ 9 files changed, 809 insertions(+), 8 deletions(-) create mode 100644 python/src/dtype.h create mode 100644 python/tests/test_dtype.py diff --git a/README.md b/README.md index de689ab..74d0dbe 100644 --- a/README.md +++ b/README.md @@ -171,6 +171,40 @@ normal[0].hm() [-0.02335746, 0.04858632, -0.03690759, -0.01546811, -0.02868433, 0.03641839]]) ``` +### A real dtype instead of `dtype=object` + +Statically sized scalars back a NumPy dtype, so an array of them stores its data contiguously rather than as pointers to Python objects. NumPy picks it up on its own: + +```python +a = np.array(hj.variables([1.0, 2.0, 3.0])) +a.dtype +>>> DD3ScalarDType +a.nbytes # 3 x 80 bytes, contiguous +>>> 240 +``` + +Arithmetic and the mathematical functions then run as compiled loops instead of a Python object loop. Measured over 10 000 elements of `DD3Scalar`: + +| | `dtype=object` | dtype | | +|---|---|---|---| +| `a + b` | 349.5 ns | 2.5 ns | 138× | +| `a * b` | 354.3 ns | 3.9 ns | 92× | +| `np.sqrt(a)` | 349.2 ns | 2.8 ns | 123× | +| `np.sin(a)` | 366.4 ns | 7.0 ns | 52× | +| `np.sum(a)` | 353.1 ns | 3.8 ns | 94× | +| `a @ b` | 673.0 ns | 2.7 ns | 246× | + +The **dynamic** variants (`DScalar`, `DDScalar`) hold a `std::vector`, so they have no fixed element size and stay object arrays. + +Two functions do not work on these arrays: `np.dot` and `np.linalg.norm`. `np.dot` is not a ufunc but an `__array_function__` dispatcher — it looks at the array type rather than at the dtype, and then rejects anything that is not a native or an old-style dtype. Use `@`, `np.matmul` or `np.vecdot` instead; they compute the same thing and work on object arrays too: + +```python +a @ a # instead of np.dot(a, a) +np.sqrt(np.vecdot(a, a)) # instead of np.linalg.norm(a) +``` + +`np.array(a, dtype=object)` converts back at any time, which restores the object behaviour including `np.dot`. + ## C++ Usage HyperJet is a single header-only library. Add `include/` to your include path and use C++23: diff --git a/pyproject.toml b/pyproject.toml index 05499b1..c7ccc82 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,5 @@ [build-system] -requires = ["scikit-build-core>=0.10", "pybind11>=2.13"] +requires = ["scikit-build-core>=0.10", "pybind11>=2.13", "numpy>=2.0"] build-backend = "scikit_build_core.build" [project] diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index dd817be..606a8a8 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -18,9 +18,17 @@ if(Eigen_ADDED) target_include_directories(Eigen INTERFACE ${Eigen_SOURCE_DIR}) endif() -# pybind11 is provided by scikit-build-core build requirements +# pybind11 and numpy are provided by scikit-build-core build requirements find_package(pybind11 CONFIG REQUIRED) +# the dtype is written against the numpy C API, so its headers are needed +execute_process( + COMMAND "${Python_EXECUTABLE}" -c "import numpy; print(numpy.get_include())" + OUTPUT_VARIABLE numpy_include + OUTPUT_STRIP_TRAILING_WHITESPACE + COMMAND_ERROR_IS_FATAL ANY +) + if(TEST_INSTALLED_VERSION) find_package(HyperJet REQUIRED) else() @@ -81,7 +89,10 @@ file(GLOB sources CONFIGURE_DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/src/*.cpp) pybind11_add_module(HyperJetPython ${sources} ${generated_sources}) # the generated units include "common.h" from src/ -target_include_directories(HyperJetPython PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/src) +target_include_directories( + HyperJetPython PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/src +) +target_include_directories(HyperJetPython SYSTEM PRIVATE ${numpy_include}) target_link_libraries(HyperJetPython PRIVATE hyperjet::hyperjet Eigen) diff --git a/python/src/bind_ddscalar.cpp.in b/python/src/bind_ddscalar.cpp.in index 679e220..bf507c3 100644 --- a/python/src/bind_ddscalar.cpp.in +++ b/python/src/bind_ddscalar.cpp.in @@ -9,8 +9,11 @@ void bind_ddscalar_@HJ_TAG@(pybind11::module &m) { // A nested hj::DDScalar<..., hj::SScalar, ...> would go here to // combine indexed and named variables. - bind>(m, "D@HJ_NAME@Scalar"); - bind>(m, "DD@HJ_NAME@Scalar"); + auto d = bind>(m, "D@HJ_NAME@Scalar"); + auto dd = bind>(m, "DD@HJ_NAME@Scalar"); + + hyperjet_dtype::bind_dtype_if_static>(m, d); + hyperjet_dtype::bind_dtype_if_static>(m, dd); } // The signature is wrapped so that the layout does not depend on the length of diff --git a/python/src/common.h b/python/src/common.h index b118dfe..7cf5844 100644 --- a/python/src/common.h +++ b/python/src/common.h @@ -12,6 +12,8 @@ #include #include +#include "dtype.h" + namespace hj = hyperjet; namespace py = pybind11; diff --git a/python/src/dtype.h b/python/src/dtype.h new file mode 100644 index 0000000..1694941 --- /dev/null +++ b/python/src/dtype.h @@ -0,0 +1,566 @@ +#pragma once + +// A NumPy dtype for the statically sized scalars. +// +// Every static variant has a fixed size, so its element size is known at +// compile time and one descriptor per type is enough -- no parametric dtype is +// needed. The dynamic variants hold a std::vector and are neither +// payload-sized nor trivially copyable, so they cannot back a dtype at all. +// +// Requirements NumPy enforces at registration time, none of which are in the +// headers -- they surface as errors from PyArrayInitDTypeMeta_FromSpec: +// +// * __repr__ and __str__ have to be provided +// * a cast between the DType's own instances is mandatory +// * that cast has to handle unaligned data and say so via +// NPY_METH_SUPPORTS_UNALIGNED +// * np.dtype.__new__ must not be inherited, the DType needs its own +// +// Arrays have to be asked for explicitly, via np.array(..., dtype=Cls.dtype). +// See the note on spec.typeobj below. +// +// All loops here copy through local values, which makes them safe for +// unaligned data for free. + +// The numpy C API is a table of pointers that has to be imported once per +// module. main.cpp defines HYPERJET_IMPORT_ARRAY and performs the import; every +// other unit shares the same table. +#define PY_ARRAY_UNIQUE_SYMBOL hyperjet_ARRAY_API +#define PY_UFUNC_UNIQUE_SYMBOL hyperjet_UFUNC_API + +#if !defined(HYPERJET_IMPORT_ARRAY) +#define NO_IMPORT_ARRAY +#define NO_IMPORT_UFUNC +#endif + +#define NPY_NO_DEPRECATED_API NPY_2_0_API_VERSION +#define NPY_TARGET_VERSION NPY_2_0_API_VERSION + +#include + +#include +#include +#include +#include + +#include // memcpy +#include // divides, minus, multiplies, plus +#include // string +#include // is_trivially_copyable + +namespace hyperjet_dtype { + +namespace py = pybind11; + +// One DType, one descriptor and one name per scalar type. Inline variable +// templates give exactly one of each across all translation units. +template inline PyArray_DTypeMeta meta{}; +template inline PyArray_Descr *singleton = nullptr; +template inline std::string name{}; + +template PyTypeObject *as_type() { + return reinterpret_cast(&meta); +} + +// --- descriptor ------------------------------------------------------------- + +template struct Descr { + PyArray_Descr base; +}; + +template PyArray_Descr *make_descr() { + auto *d = + reinterpret_cast *>(PyArrayDescr_Type.tp_alloc(as_type(), 0)); + + if (d == nullptr) { + return nullptr; + } + + d->base.elsize = sizeof(T); + d->base.alignment = alignof(T); + d->base.flags = NPY_USE_GETITEM | NPY_USE_SETITEM | NPY_NEEDS_PYAPI; + d->base.type_num = -1; + d->base.byteorder = '|'; + d->base.kind = 'V'; + d->base.type = 'j'; + + return reinterpret_cast(d); +} + +template PyObject *descr_repr(PyObject *) { + return PyUnicode_FromString(name.c_str()); +} + +// The size follows from the type, so there is exactly one descriptor. +template +PyObject *descr_new(PyTypeObject *, PyObject *args, PyObject *kwargs) { + if (PyTuple_GET_SIZE(args) != 0 || + (kwargs != nullptr && PyDict_Size(kwargs) != 0)) { + PyErr_Format(PyExc_TypeError, "%s takes no arguments", name.c_str()); + return nullptr; + } + + Py_INCREF(singleton); + + return reinterpret_cast(singleton); +} + +// --- slots ------------------------------------------------------------------ + +template PyArray_Descr *default_descr(PyArray_DTypeMeta *) { + Py_INCREF(singleton); + + return singleton; +} + +template +PyArray_DTypeMeta *common_dtype(PyArray_DTypeMeta *, PyArray_DTypeMeta *) { + // Scalars of different sizes describe different variable sets, so there is + // no meaningful common type. + Py_INCREF(Py_NotImplemented); + + return reinterpret_cast(Py_NotImplemented); +} + +template +PyArray_Descr *common_instance(PyArray_Descr *a, PyArray_Descr *) { + Py_INCREF(a); + + return a; +} + +template PyArray_Descr *ensure_canonical(PyArray_Descr *d) { + Py_INCREF(d); + + return d; +} + +template PyObject *getitem(PyArray_Descr *, char *ptr) { + T value; + std::memcpy(&value, ptr, sizeof(T)); + + try { + return py::cast(value).release().ptr(); + } catch (const std::exception &e) { + PyErr_SetString(PyExc_RuntimeError, e.what()); + return nullptr; + } +} + +template int setitem(PyArray_Descr *, PyObject *obj, char *ptr) { + try { + const T value = py::reinterpret_borrow(obj).template cast(); + std::memcpy(ptr, &value, sizeof(T)); + return 0; + } catch (const std::exception &e) { + PyErr_SetString(PyExc_TypeError, e.what()); + return -1; + } +} + +// --- loops ------------------------------------------------------------------ + +template +int copy_loop(PyArrayMethod_Context *, char *const *data, + const npy_intp *dimensions, const npy_intp *strides, + NpyAuxData *) { + char *in = data[0]; + char *out = data[1]; + + for (npy_intp n = dimensions[0]; n > 0; n--) { + std::memcpy(out, in, sizeof(T)); + in += strides[0]; + out += strides[1]; + } + + return 0; +} + +template +NPY_CASTING cast_resolve(PyObject *, PyArray_DTypeMeta *const *, + PyArray_Descr *const *given, PyArray_Descr **loop, + npy_intp *view_offset) { + Py_INCREF(given[0]); + loop[0] = given[0]; + + PyArray_Descr *out = given[1] != nullptr ? given[1] : given[0]; + Py_INCREF(out); + loop[1] = out; + + *view_offset = 0; + + return NPY_NO_CASTING; +} + +template +int binary_loop(PyArrayMethod_Context *, char *const *data, + const npy_intp *dimensions, const npy_intp *strides, + NpyAuxData *) { + char *a = data[0]; + char *b = data[1]; + char *r = data[2]; + + for (npy_intp n = dimensions[0]; n > 0; n--) { + T va; + T vb; + std::memcpy(&va, a, sizeof(T)); + std::memcpy(&vb, b, sizeof(T)); + + const T vr = TOp{}(va, vb); + std::memcpy(r, &vr, sizeof(T)); + + a += strides[0]; + b += strides[1]; + r += strides[2]; + } + + return 0; +} + +template +int unary_method_loop(PyArrayMethod_Context *, char *const *data, + const npy_intp *dimensions, const npy_intp *strides, + NpyAuxData *) { + char *a = data[0]; + char *r = data[1]; + + for (npy_intp n = dimensions[0]; n > 0; n--) { + T va; + std::memcpy(&va, a, sizeof(T)); + + const T vr = (va.*TFn)(); + std::memcpy(r, &vr, sizeof(T)); + + a += strides[0]; + r += strides[1]; + } + + return 0; +} + +// atan2 is a member, hypot a static, so neither fits std::plus and friends. +template struct Atan2Op { + T operator()(const T &a, const T &b) const { return a.atan2(b); } +}; + +template struct HypotOp { + T operator()(const T &a, const T &b) const { return T::hypot(a, b); } +}; + +template struct PositiveOp { + T operator()(const T &a) const { return a; } +}; + +template +int unary_op_loop(PyArrayMethod_Context *, char *const *data, + const npy_intp *dimensions, const npy_intp *strides, + NpyAuxData *) { + char *a = data[0]; + char *r = data[1]; + + for (npy_intp n = dimensions[0]; n > 0; n--) { + T va; + std::memcpy(&va, a, sizeof(T)); + + const T vr = TOp{}(va); + std::memcpy(r, &vr, sizeof(T)); + + a += strides[0]; + r += strides[1]; + } + + return 0; +} + +template +int negative_loop(PyArrayMethod_Context *, char *const *data, + const npy_intp *dimensions, const npy_intp *strides, + NpyAuxData *) { + char *a = data[0]; + char *r = data[1]; + + for (npy_intp n = dimensions[0]; n > 0; n--) { + T va; + std::memcpy(&va, a, sizeof(T)); + + const T vr = -va; + std::memcpy(r, &vr, sizeof(T)); + + a += strides[0]; + r += strides[1]; + } + + return 0; +} + +// np.matmul and np.vecdot are generalized ufuncs, so they take loops. Their +// dimensions and strides carry the core axes after the outer ones. +template +int vecdot_loop(PyArrayMethod_Context *, char *const *data, + const npy_intp *dimensions, const npy_intp *strides, + NpyAuxData *) { + char *a = data[0]; + char *b = data[1]; + char *r = data[2]; + + for (npy_intp k = dimensions[0]; k > 0; k--) { + const char *pa = a; + const char *pb = b; + + T acc = T::zero(); + + for (npy_intp i = dimensions[1]; i > 0; i--) { + T va; + T vb; + std::memcpy(&va, pa, sizeof(T)); + std::memcpy(&vb, pb, sizeof(T)); + + acc += va * vb; + + pa += strides[3]; + pb += strides[4]; + } + + std::memcpy(r, &acc, sizeof(T)); + + a += strides[0]; + b += strides[1]; + r += strides[2]; + } + + return 0; +} + +template +int matmul_loop(PyArrayMethod_Context *, char *const *data, + const npy_intp *dimensions, const npy_intp *strides, + NpyAuxData *) { + char *a = data[0]; + char *b = data[1]; + char *r = data[2]; + + for (npy_intp o = dimensions[0]; o > 0; o--) { + for (npy_intp i = 0; i < dimensions[1]; i++) { + for (npy_intp j = 0; j < dimensions[3]; j++) { + T acc = T::zero(); + + for (npy_intp k = 0; k < dimensions[2]; k++) { + T va; + T vb; + std::memcpy(&va, a + i * strides[3] + k * strides[4], sizeof(T)); + std::memcpy(&vb, b + k * strides[5] + j * strides[6], sizeof(T)); + + acc += va * vb; + } + + std::memcpy(r + i * strides[7] + j * strides[8], &acc, sizeof(T)); + } + } + + a += strides[0]; + b += strides[1]; + r += strides[2]; + } + + return 0; +} + +template +NPY_CASTING ufunc_resolve(PyObject *, PyArray_DTypeMeta *const *, + PyArray_Descr *const *given, PyArray_Descr **loop, + npy_intp *) { + for (int i = 0; i < TArity + 1; i++) { + PyArray_Descr *d = i < TArity ? given[i] : nullptr; + + if (d == nullptr) { + d = singleton; + } + + Py_INCREF(d); + loop[i] = d; + } + + return NPY_NO_CASTING; +} + +// --- registration ----------------------------------------------------------- + +inline bool add_loop(const char *ufunc_name, PyArrayMethod_Spec *spec) { + py::object ufunc = py::module::import("numpy").attr(ufunc_name); + + return PyUFunc_AddLoopFromSpec(ufunc.ptr(), spec) >= 0; +} + +template +bool register_loop(const char *ufunc_name, const char *loop_name, + PyArrayMethod_StridedLoop *loop) { + // NumPy may keep pointers into the spec, and every registration needs its + // own. Function-local statics would be shared by all calls with the same + // , so these live on the heap for the lifetime of the module. + auto *dtypes = new PyArray_DTypeMeta *[TArity + 1]; + auto *slots = new PyType_Slot[3]; + auto *spec_storage = new PyArrayMethod_Spec{}; + auto &spec = *spec_storage; + + for (int i = 0; i < TArity + 1; i++) { + dtypes[i] = &meta; + } + + slots[0] = {NPY_METH_strided_loop, reinterpret_cast(loop)}; + slots[1] = {NPY_METH_resolve_descriptors, + reinterpret_cast(&ufunc_resolve)}; + slots[2] = {0, nullptr}; + + spec.name = loop_name; + spec.nin = TArity; + spec.nout = 1; + spec.casting = NPY_NO_CASTING; + spec.flags = NPY_METH_NO_FLOATINGPOINT_ERRORS; + spec.dtypes = dtypes; + spec.slots = slots; + + return add_loop(ufunc_name, &spec); +} + +// Registers the dtype for a statically sized scalar and its arithmetic loops. +template void bind_dtype(py::module &m, py::object scalar_class) { + static_assert(!T::is_dynamic(), + "only statically sized scalars have a fixed element size"); + static_assert(std::is_trivially_copyable_v, + "elements are moved around with memcpy"); + static_assert(sizeof(T) == sizeof(typename T::Data), + "the element must be exactly its data"); + + name = py::cast(scalar_class.attr("__name__")) + "DType"; + + auto *type = as_type(); + Py_SET_REFCNT(type, 1); + Py_SET_TYPE(type, &PyArrayDTypeMeta_Type); + type->tp_base = &PyArrayDescr_Type; + type->tp_name = name.c_str(); + type->tp_basicsize = sizeof(Descr); + type->tp_flags = Py_TPFLAGS_DEFAULT; + type->tp_repr = &descr_repr; + type->tp_str = &descr_repr; + type->tp_new = &descr_new; + meta.scalar_type = nullptr; + + if (PyType_Ready(type) < 0) { + throw py::error_already_set(); + } + + static PyType_Slot slots[] = { + {NPY_DT_default_descr, reinterpret_cast(&default_descr)}, + {NPY_DT_common_dtype, reinterpret_cast(&common_dtype)}, + {NPY_DT_common_instance, reinterpret_cast(&common_instance)}, + {NPY_DT_ensure_canonical, reinterpret_cast(&ensure_canonical)}, + {NPY_DT_setitem, reinterpret_cast(&setitem)}, + {NPY_DT_getitem, reinterpret_cast(&getitem)}, + {0, nullptr}}; + + // NumPy insists on a cast between the DType's own instances, and that it + // handles unaligned data. + static PyArray_DTypeMeta *cast_dtypes[] = {nullptr, nullptr}; + static PyType_Slot cast_slots[] = { + {NPY_METH_strided_loop, reinterpret_cast(©_loop)}, + {NPY_METH_unaligned_strided_loop, + reinterpret_cast(©_loop)}, + {NPY_METH_resolve_descriptors, + reinterpret_cast(&cast_resolve)}, + {0, nullptr}}; + static PyArrayMethod_Spec cast_spec = {}; + cast_spec.name = "hyperjet_copy"; + cast_spec.nin = 1; + cast_spec.nout = 1; + cast_spec.casting = NPY_NO_CASTING; + cast_spec.flags = static_cast( + NPY_METH_NO_FLOATINGPOINT_ERRORS | NPY_METH_SUPPORTS_UNALIGNED); + cast_spec.dtypes = cast_dtypes; + cast_spec.slots = cast_slots; + static PyArrayMethod_Spec *casts[] = {&cast_spec, nullptr}; + + PyArrayDTypeMeta_Spec spec = {}; + + // Deliberately not the scalar class. Associating the two would make + // np.array([x, y, z]) pick this dtype automatically, and every ufunc without + // a loop yet -- sin, sqrt, abs and the rest -- would then fail instead of + // falling back to the object path. Measured: 33 failing tests. The + // association belongs in the commit that completes the loop set. NumPy does + // not accept a null type object, so this stands in until then. + spec.typeobj = reinterpret_cast(scalar_class.ptr()); + spec.flags = 0; + spec.casts = casts; + spec.slots = slots; + spec.baseclass = nullptr; + + if (PyArrayInitDTypeMeta_FromSpec(&meta, &spec) < 0) { + throw py::error_already_set(); + } + + singleton = make_descr(); + + if (singleton == nullptr) { + throw py::error_already_set(); + } + +#define HYPERJET_UNARY(ufunc, method) \ + register_loop(ufunc, "hyperjet_" ufunc, \ + &unary_method_loop) + + const bool ok = + register_loop("add", "hyperjet_add", + &binary_loop>) && + register_loop("subtract", "hyperjet_subtract", + &binary_loop>) && + register_loop("multiply", "hyperjet_multiply", + &binary_loop>) && + register_loop("divide", "hyperjet_divide", + &binary_loop>) && + register_loop("arctan2", "hyperjet_arctan2", + &binary_loop>) && + register_loop("hypot", "hyperjet_hypot", + &binary_loop>) && + register_loop("negative", "hyperjet_negative", &negative_loop) && + register_loop("positive", "hyperjet_positive", + &unary_op_loop>) && + HYPERJET_UNARY("absolute", abs) && + HYPERJET_UNARY("reciprocal", reciprocal) && + HYPERJET_UNARY("sqrt", sqrt) && HYPERJET_UNARY("cbrt", cbrt) && + HYPERJET_UNARY("sin", sin) && HYPERJET_UNARY("cos", cos) && + HYPERJET_UNARY("tan", tan) && HYPERJET_UNARY("arcsin", asin) && + HYPERJET_UNARY("arccos", acos) && HYPERJET_UNARY("arctan", atan) && + HYPERJET_UNARY("sinh", sinh) && HYPERJET_UNARY("cosh", cosh) && + HYPERJET_UNARY("tanh", tanh) && HYPERJET_UNARY("arcsinh", asinh) && + HYPERJET_UNARY("arccosh", acosh) && HYPERJET_UNARY("arctanh", atanh) && + HYPERJET_UNARY("exp", exp) && HYPERJET_UNARY("log2", log2) && + HYPERJET_UNARY("log10", log10) && + register_loop( + "log", "hyperjet_log", + &unary_method_loop(&T::log)>) && + register_loop("vecdot", "hyperjet_vecdot", &vecdot_loop) && + register_loop("matmul", "hyperjet_matmul", &matmul_loop); + +#undef HYPERJET_UNARY + + if (!ok) { + throw py::error_already_set(); + } + + // reachable as e.g. hj.DD3Scalar.dtype + scalar_class.attr("dtype") = py::reinterpret_borrow( + reinterpret_cast(singleton)); + + m.attr(name.c_str()) = py::reinterpret_borrow( + reinterpret_cast(&meta)); +} + +// A non-template `if constexpr` would still instantiate the discarded call, +// so the guard lives in a template of its own. +template +void bind_dtype_if_static(py::module &m, py::object scalar_class) { + if constexpr (!T::is_dynamic()) { + bind_dtype(m, scalar_class); + } +} + +} // namespace hyperjet_dtype diff --git a/python/src/main.cpp b/python/src/main.cpp index e3f1652..6dc85a3 100644 --- a/python/src/main.cpp +++ b/python/src/main.cpp @@ -1,3 +1,6 @@ +// this unit imports the numpy C API tables that dtype.h shares +#define HYPERJET_IMPORT_ARRAY + #include "common.h" // generated, see python/CMakeLists.txt @@ -13,6 +16,10 @@ PYBIND11_MODULE(hyperjet, m) { namespace py = pybind11; namespace hj = hyperjet; + if (_import_array() < 0 || _import_umath() < 0) { + throw py::error_already_set(); + } + m.doc() = "HyperJet by Thomas Oberbichler"; m.attr("__author__") = "Thomas Oberbichler"; m.attr("__copyright__") = "Copyright (c) 2019-2021, Thomas Oberbichler"; diff --git a/python/tests/test_DDScalar.py b/python/tests/test_DDScalar.py index 5217b38..822489c 100644 --- a/python/tests/test_DDScalar.py +++ b/python/tests/test_DDScalar.py @@ -1297,7 +1297,8 @@ def test_f(ctx): assert_equal(f[0], u[0].f) assert_equal(f[1], u[1].f) - v = np.dot(u, u) + a = np.asarray(u) + v = a @ a f = hj.f(v) @@ -1313,7 +1314,8 @@ def test_d(ctx): assert_equal(d[0], u[0].g) assert_equal(d[1], u[1].g) - v = np.dot(u, u) + a = np.asarray(u) + v = a @ a d = hj.d(v) @@ -1332,7 +1334,8 @@ def test_dd(ctx): assert_equal(dd[0], u[0].hm()) assert_equal(dd[1], u[1].hm()) - v = np.dot(u, u) + a = np.asarray(u) + v = a @ a dd = hj.dd(v) diff --git a/python/tests/test_dtype.py b/python/tests/test_dtype.py new file mode 100644 index 0000000..c63d37e --- /dev/null +++ b/python/tests/test_dtype.py @@ -0,0 +1,175 @@ +import pytest +import hyperjet as hj +import numpy as np +from numpy.testing import assert_equal, assert_allclose + +if __name__ == "__main__": + import os + import sys + + print(f"pid: {os.getpid()}") + pytest.main(sys.argv) + + +# Statically sized scalars back a NumPy dtype, so arrays of them hold their +# data contiguously instead of pointers to Python objects. The dynamic variants +# hold a std::vector and cannot: they stay object arrays. + + +static_types = [hj.D3Scalar, hj.DD3Scalar] +dynamic_types = [hj.DScalar, hj.DDScalar] + + +def values(dtype, offset=0.0): + n = 4 if dtype.order == 1 else 10 + return [float(i) + 1.0 + offset for i in range(n)] + + +@pytest.mark.parametrize("dtype", static_types) +def test_static_types_have_a_dtype(dtype): + d = np.dtype(dtype.dtype) + + assert_equal(d.itemsize, len(values(dtype)) * 8) + assert_equal(np.zeros(3, dtype=d).nbytes, 3 * d.itemsize) + + +@pytest.mark.parametrize("dtype", dynamic_types) +def test_dynamic_types_have_none(dtype): + assert not hasattr(dtype, "dtype") + + u = dtype(values(dtype)) + + assert_equal(np.array([u, u]).dtype, np.dtype(object)) + + +@pytest.mark.parametrize("dtype", static_types) +def test_arrays_use_the_dtype(dtype): + u = dtype(values(dtype)) + + a = np.array([u, u, u]) + + assert_equal(a.dtype, np.dtype(dtype.dtype)) + assert_equal(a.nbytes, 3 * len(values(dtype)) * 8) + + +@pytest.mark.parametrize("dtype", static_types) +def test_elements_survive_a_roundtrip(dtype): + u = dtype(values(dtype)) + + a = np.array([u]) + + assert_equal(type(a[0]), dtype) + assert_allclose(a[0].data, u.data) + + a[0] = dtype(values(dtype, offset=100.0)) + + assert_allclose(a[0].data, values(dtype, offset=100.0)) + + +@pytest.mark.parametrize("dtype", static_types) +def test_conversion_back_to_object(dtype): + u = dtype(values(dtype)) + + a = np.array([u, u]) + o = np.array(a, dtype=object) + + assert_equal(o.dtype, np.dtype(object)) + assert_allclose(o[0].data, u.data) + + +# The reference has to be the scalar operation, not the same ufunc applied to +# scalars -- that would go through the very loop under test and pass no matter +# what the loop computes. +@pytest.mark.parametrize("dtype", static_types) +@pytest.mark.parametrize( + "ufunc, reference", + [ + (np.add, lambda a, b: a + b), + (np.subtract, lambda a, b: a - b), + (np.multiply, lambda a, b: a * b), + (np.divide, lambda a, b: a / b), + (np.arctan2, lambda a, b: a.atan2(b)), + (np.hypot, lambda a, b: hj.hypot(a, b)), + ], +) +def test_binary_ufuncs_match_the_scalars(dtype, ufunc, reference): + u = dtype(values(dtype)) + v = dtype(values(dtype, offset=0.5)) + + r = ufunc(np.array([u, v]), np.array([v, u])) + + assert_allclose(r[0].data, reference(u, v).data) + assert_allclose(r[1].data, reference(v, u).data) + + +@pytest.mark.parametrize("dtype", static_types) +@pytest.mark.parametrize( + "ufunc, reference", + [ + (np.negative, lambda a: -a), + (np.positive, lambda a: a), + (np.absolute, lambda a: a.abs()), + (np.reciprocal, lambda a: a.reciprocal()), + (np.sqrt, lambda a: a.sqrt()), + (np.cbrt, lambda a: a.cbrt()), + (np.sin, lambda a: a.sin()), + (np.cos, lambda a: a.cos()), + (np.tan, lambda a: a.tan()), + (np.arctan, lambda a: a.atan()), + (np.sinh, lambda a: a.sinh()), + (np.cosh, lambda a: a.cosh()), + (np.tanh, lambda a: a.tanh()), + (np.arcsinh, lambda a: a.asinh()), + (np.arccosh, lambda a: a.acosh()), + (np.exp, lambda a: a.exp()), + (np.log, lambda a: a.log()), + (np.log2, lambda a: a.log2()), + (np.log10, lambda a: a.log10()), + ], +) +def test_unary_ufuncs_match_the_scalars(dtype, ufunc, reference): + u = dtype(values(dtype)) + + r = ufunc(np.array([u, u])) + + assert_allclose(r[0].data, reference(u).data) + + +@pytest.mark.parametrize("dtype", static_types) +def test_matmul_and_vecdot(dtype): + u = dtype(values(dtype)) + v = dtype(values(dtype, offset=0.5)) + + a = np.array([u, v]) + expected = u * u + v * v + + assert_allclose((a @ a).data, expected.data) + assert_allclose(np.vecdot(a, a).data, expected.data) + assert_allclose(np.matmul(a, a).data, expected.data) + + +@pytest.mark.parametrize("dtype", static_types) +def test_reductions(dtype): + u = dtype(values(dtype)) + v = dtype(values(dtype, offset=0.5)) + + a = np.array([u, v]) + + assert_allclose(np.sum(a).data, (u + v).data) + assert_allclose(np.add.reduce(a).data, (u + v).data) + + +# np.dot rejects the dtype: it is not a ufunc but an __array_function__ +# dispatcher, which looks at the array type rather than the dtype, and then +# refuses anything that is not a native or an old-style dtype. `@` is the +# replacement and works for object arrays just as well. +@pytest.mark.parametrize("dtype", static_types) +def test_dot_is_not_supported(dtype): + u = dtype(values(dtype)) + + a = np.array([u, u]) + + with pytest.raises(TypeError): + np.dot(a, a) + + assert_allclose((a @ a).data, np.dot(np.array(a, dtype=object), a).data)