diff --git a/bottleneck/include/bottleneck.h b/bottleneck/include/bottleneck.h deleted file mode 100644 index 2c40557b4..000000000 --- a/bottleneck/include/bottleneck.h +++ /dev/null @@ -1,201 +0,0 @@ -// Copyright 2010-2019 Keith Goodman -// Copyright 2019 Bottleneck Developers -#ifndef BOTTLENECK_H_ -#define BOTTLENECK_H_ - -#include -#define NPY_NO_DEPRECATED_API NPY_1_11_API_VERSION -#include -#include - -/* THREADS=1 releases the GIL but increases function call - * overhead. THREADS=0 does not release the GIL but keeps - * function call overhead low. Curly brackets are for C89 - * support. */ -#define THREADS 1 -#if THREADS - #define BN_BEGIN_ALLOW_THREADS Py_BEGIN_ALLOW_THREADS { - #define BN_END_ALLOW_THREADS \ - ; \ - } \ - Py_END_ALLOW_THREADS -#else - #define BN_BEGIN_ALLOW_THREADS { - #define BN_END_ALLOW_THREADS } -#endif - -/* for ease of dtype templating */ -#define NPY_float64 NPY_FLOAT64 -#define NPY_float32 NPY_FLOAT32 -#define NPY_int64 NPY_INT64 -#define NPY_int32 NPY_INT32 -#define NPY_intp NPY_INTP -#define NPY_long NPY_LONG -#define NPY_MAX_int64 NPY_MAX_INT64 -#define NPY_MAX_int32 NPY_MAX_INT32 -#define NPY_MIN_int64 NPY_MIN_INT64 -#define NPY_MIN_int32 NPY_MIN_INT32 - -#define VARKEY (((unsigned)METH_VARARGS) | ((unsigned)METH_KEYWORDS)) -#define error_converting(x) (((x) == -1) && PyErr_Occurred()) - -#define VALUE_ERR(text) PyErr_SetString(PyExc_ValueError, text) -#define TYPE_ERR(text) PyErr_SetString(PyExc_TypeError, text) -#define MEMORY_ERR(text) PyErr_SetString(PyExc_MemoryError, text) -#define RUNTIME_ERR(text) PyErr_SetString(PyExc_RuntimeError, text) - -/* `inline`, `opt_3`, and isnan copied from NumPy. */ -#if HAVE_ATTRIBUTE_OPTIMIZE_OPT_3 - #define BN_OPT_3 __attribute__((optimize("O3"))) -#else - #define BN_OPT_3 -#endif - -#if HAVE___BUILTIN_ISNAN - #define bn_isnan(x) __builtin_isnan(x) -#elif HAVE_ISNAN - #define bn_isnan(x) isnan(x) -#elif HAVE__ISNAN - #define bn_isnan(x) _isnan(x) -#else - #define bn_isnan(x) ((x) != (x)) -#endif - -/* - * NAN and INFINITY like macros (same behavior as glibc for NAN, same as C99 - * for INFINITY). Copied from NumPy. - */ -static inline float __bn_inff(void) { - const union { - npy_uint32 __i; - float __f; - } __bint = {0x7f800000UL}; - return __bint.__f; -} - -static inline float __bn_nanf(void) { - const union { - npy_uint32 __i; - float __f; - } __bint = {0x7fc00000UL}; - return __bint.__f; -} - -#define BN_INFINITYF __bn_inff() -#define BN_NANF __bn_nanf() -#define BN_INFINITY ((npy_double)BN_INFINITYF) -#define BN_NAN ((npy_double)BN_NANF) - -/* WIRTH ----------------------------------------------------------------- */ - -/* - WIRTH macro based on: - Fast median search: an ANSI C implementation - Nicolas Devillard - ndevilla AT free DOT fr - July 1998 - which, in turn, took the algorithm from - Wirth, Niklaus - Algorithms + data structures = programs, p. 366 - Englewood Cliffs: Prentice-Hall, 1976 - - Adapted for Bottleneck: - (C) 2016 Keith Goodman -*/ - -#define WIRTH(dtype) \ - npy_##dtype x = B(dtype, k); \ - npy_intp m = l; \ - j = r; \ - do { \ - while (B(dtype, m) < x) \ - m++; \ - while (x < B(dtype, j)) \ - j--; \ - if (m <= j) { \ - const npy_##dtype atmp = B(dtype, m); \ - B(dtype, m) = B(dtype, j); \ - B(dtype, j) = atmp; \ - m++; \ - j--; \ - } \ - } while (m <= j); \ - if (j < k) l = m; \ - if (k < m) r = j; - -/* partition ------------------------------------------------------------- */ - -#define PARTITION(dtype) \ - while (l < r) { \ - const npy_##dtype al = B(dtype, l); \ - const npy_##dtype ak = B(dtype, k); \ - const npy_##dtype ar = B(dtype, r); \ - if (al > ak) { \ - if (ak < ar) { \ - if (al < ar) { \ - B(dtype, k) = al; \ - B(dtype, l) = ak; \ - } else { \ - B(dtype, k) = ar; \ - B(dtype, r) = ak; \ - } \ - } \ - } else { \ - if (ak > ar) { \ - if (al > ar) { \ - B(dtype, k) = al; \ - B(dtype, l) = ak; \ - } else { \ - B(dtype, k) = ar; \ - B(dtype, r) = ak; \ - } \ - } \ - } \ - WIRTH(dtype) \ - } - -/* slow ------------------------------------------------------------------ */ - -static PyObject *slow_module = NULL; - -static PyObject * -slow(char *name, PyObject *args, PyObject *kwds) { - PyObject *func = NULL; - PyObject *out = NULL; - - if (slow_module == NULL) { - /* bottleneck.slow has not been imported during the current - * python session. Only import it once per session to save time */ - slow_module = PyImport_ImportModule("bottleneck.slow"); - if (slow_module == NULL) { - PyErr_SetString(PyExc_RuntimeError, - "Cannot import bottleneck.slow"); - return NULL; - } - } - - func = PyObject_GetAttrString(slow_module, name); - if (func == NULL) { - PyErr_Format(PyExc_RuntimeError, - "Cannot import %s from bottleneck.slow", - name); - return NULL; - } - if (PyCallable_Check(func)) { - out = PyObject_Call(func, args, kwds); - if (out == NULL) { - Py_XDECREF(func); - return NULL; - } - } else { - Py_XDECREF(func); - PyErr_Format(PyExc_RuntimeError, - "bottleneck.slow.%s is not callable", - name); - return NULL; - } - Py_XDECREF(func); - - return out; -} - -#endif // BOTTLENECK_H_ diff --git a/bottleneck/include/iterators.h b/bottleneck/include/iterators.h deleted file mode 100644 index 1f5120a18..000000000 --- a/bottleneck/include/iterators.h +++ /dev/null @@ -1,405 +0,0 @@ -// Copyright 2010-2019 Keith Goodman -// Copyright 2019 Bottleneck Developers -#ifndef ITERATORS_H_ -#define ITERATORS_H_ - -#define NPY_NO_DEPRECATED_API NPY_1_11_API_VERSION -#include - -/* - Bottleneck iterators are based on ideas from NumPy's PyArray_IterAllButAxis - and PyArray_ITER_NEXT. -*/ - -#define C_CONTIGUOUS PyArray_IS_C_CONTIGUOUS -#define F_CONTIGUOUS PyArray_IS_F_CONTIGUOUS - -/* one input array ------------------------------------------------------- */ - -/* these iterators are used mainly by reduce functions such as nansum */ - -struct _iter { - int ndim_m2; /* ndim - 2 */ - int axis; /* axis to not iterate over */ - Py_ssize_t length; /* a.shape[axis] */ - Py_ssize_t astride; /* a.strides[axis] */ - npy_intp stride; /* element-level stride to take in the array */ - npy_intp i; /* integer used by some macros */ - npy_intp its; /* number of iterations completed */ - npy_intp nits; /* number of iterations iterator plans to make */ - npy_intp indices[NPY_MAXDIMS]; /* current location of iterator */ - npy_intp astrides[NPY_MAXDIMS]; /* a.strides, a.strides[axis] removed */ - npy_intp shape[NPY_MAXDIMS]; /* a.shape, a.shape[axis] removed */ - char * pa; /* pointer to data corresponding to indices */ - PyArrayObject *a_ravel; /* NULL or pointer to ravelled input array */ -}; -typedef struct _iter iter; - -static inline void init_iter_one(iter *it, PyArrayObject *a, int axis) { - const int ndim = PyArray_NDIM(a); - const npy_intp *shape = PyArray_SHAPE(a); - const npy_intp *strides = PyArray_STRIDES(a); - const npy_intp item_size = PyArray_ITEMSIZE(a); - - it->axis = axis; - it->its = 0; - it->nits = 1; - it->pa = PyArray_BYTES(a); - - it->ndim_m2 = -1; - it->length = 1; - it->astride = 0; - - if (ndim != 0) { - it->ndim_m2 = ndim - 2; - int j = 0; - for (int i = 0; i < ndim; i++) { - if (i == axis) { - it->astride = strides[i]; - it->length = shape[i]; - } else { - it->indices[j] = 0; - it->astrides[j] = strides[i]; - it->shape[j] = shape[i]; - it->nits *= shape[i]; - j++; - } - } - } - it->stride = it->astride / item_size; -} - -/* - * If both ravel != 0 and it.a_ravel != NULL then you are responsible for - * calling Py_DECREF(it.a_ravel) after you are done with the iterator. - * See nanargmin for an example. - */ -static inline void init_iter_all(iter *it, PyArrayObject *a, int ravel, int anyorder) { - int i = 0; - int j = 0; - const int ndim = PyArray_NDIM(a); - const npy_intp *shape = PyArray_SHAPE(a); - const npy_intp *strides = PyArray_STRIDES(a); - const npy_intp item_size = PyArray_ITEMSIZE(a); - - it->axis = 0; - it->its = 0; - it->nits = 1; - it->a_ravel = NULL; - - /* The fix for relaxed strides checking in numpy and the fix for - * issue #183 has left this if..else tree in need of a refactor from the - * the ground up */ - if (ndim == 1) { - it->ndim_m2 = -1; - it->length = shape[0]; - it->astride = strides[0]; - } else if (ndim == 0) { - it->ndim_m2 = -1; - it->length = 1; - it->astride = 0; - } else if (C_CONTIGUOUS(a) && !F_CONTIGUOUS(a)) { - /* The &&! in the next two else ifs is to deal with relaxed - * stride checking introduced in numpy 1.12.0; see gh #161 */ - it->ndim_m2 = -1; - it->axis = ndim - 1; - it->length = PyArray_SIZE(a); - it->astride = 0; - for (i = ndim - 1; i > -1; i--) { - /* protect against length zero strides such as in - * np.ones((2, 2))[..., np.newaxis] */ - if (strides[i] == 0) { - continue; - } - it->astride = strides[i]; - break; - } - } else if (F_CONTIGUOUS(a) && !C_CONTIGUOUS(a)) { - if (anyorder || !ravel) { - it->ndim_m2 = -1; - it->length = PyArray_SIZE(a); - it->astride = 0; - for (i = 0; i < ndim; i++) { - /* protect against length zero strides such as in - * np.ones((2, 2), order='F')[np.newaxis, ...] */ - if (strides[i] == 0) { - continue; - } - it->astride = strides[i]; - break; - } - } else { - it->ndim_m2 = -1; - if (anyorder) { - a = (PyArrayObject *)PyArray_Ravel(a, NPY_ANYORDER); - } else { - a = (PyArrayObject *)PyArray_Ravel(a, NPY_CORDER); - } - it->a_ravel = a; - it->length = PyArray_DIM(a, 0); - it->astride = PyArray_STRIDE(a, 0); - } - } else if (ravel) { - it->ndim_m2 = -1; - if (anyorder) { - a = (PyArrayObject *)PyArray_Ravel(a, NPY_ANYORDER); - } else { - a = (PyArrayObject *)PyArray_Ravel(a, NPY_CORDER); - } - it->a_ravel = a; - it->length = PyArray_DIM(a, 0); - it->astride = PyArray_STRIDE(a, 0); - } else { - it->ndim_m2 = ndim - 2; - it->astride = strides[0]; - for (i = 1; i < ndim; i++) { - if (strides[i] < it->astride) { - it->astride = strides[i]; - it->axis = i; - } - } - it->length = shape[it->axis]; - for (i = 0; i < ndim; i++) { - if (i != it->axis) { - it->indices[j] = 0; - it->astrides[j] = strides[i]; - it->shape[j] = shape[i]; - it->nits *= shape[i]; - j++; - } - } - } - - it->stride = it->astride / item_size; - it->pa = PyArray_BYTES(a); -} - -#define NEXT \ - for (it.i = it.ndim_m2; it.i > -1; it.i--) { \ - if (it.indices[it.i] < it.shape[it.i] - 1) { \ - it.pa += it.astrides[it.i]; \ - it.indices[it.i]++; \ - break; \ - } \ - it.pa -= it.indices[it.i] * it.astrides[it.i]; \ - it.indices[it.i] = 0; \ - } \ - it.its++; - -/* two input arrays ------------------------------------------------------ */ - -/* this iterator is used mainly by moving window functions such as move_sum */ - -struct _iter2 { - int ndim_m2; - int axis; - Py_ssize_t length; - Py_ssize_t astride; - Py_ssize_t ystride; - npy_intp i; - npy_intp its; - npy_intp nits; - npy_intp indices[NPY_MAXDIMS]; - npy_intp astrides[NPY_MAXDIMS]; - npy_intp ystrides[NPY_MAXDIMS]; - npy_intp shape[NPY_MAXDIMS]; - char * pa; - char * py; -}; -typedef struct _iter2 iter2; - -static inline void init_iter2(iter2 *it, PyArrayObject *a, PyObject *y, int axis) { - int i = 0; - int j = 0; - const int ndim = PyArray_NDIM(a); - const npy_intp *shape = PyArray_SHAPE(a); - const npy_intp *astrides = PyArray_STRIDES(a); - const npy_intp *ystrides = PyArray_STRIDES((PyArrayObject *)y); - - /* to avoid compiler warning of uninitialized variables */ - it->length = 0; - it->astride = 0; - it->ystride = 0; - - it->ndim_m2 = ndim - 2; - it->axis = axis; - it->its = 0; - it->nits = 1; - it->pa = PyArray_BYTES(a); - it->py = PyArray_BYTES((PyArrayObject *)y); - - for (i = 0; i < ndim; i++) { - if (i == axis) { - it->astride = astrides[i]; - it->ystride = ystrides[i]; - it->length = shape[i]; - } else { - it->indices[j] = 0; - it->astrides[j] = astrides[i]; - it->ystrides[j] = ystrides[i]; - it->shape[j] = shape[i]; - it->nits *= shape[i]; - j++; - } - } -} - -#define NEXT2 \ - for (it.i = it.ndim_m2; it.i > -1; it.i--) { \ - if (it.indices[it.i] < it.shape[it.i] - 1) { \ - it.pa += it.astrides[it.i]; \ - it.py += it.ystrides[it.i]; \ - it.indices[it.i]++; \ - break; \ - } \ - it.pa -= it.indices[it.i] * it.astrides[it.i]; \ - it.py -= it.indices[it.i] * it.ystrides[it.i]; \ - it.indices[it.i] = 0; \ - } \ - it.its++; - -/* three input arrays ---------------------------------------------------- */ - -/* this iterator is used mainly by rankdata and nanrankdata */ - -struct _iter3 { - int ndim_m2; - int axis; - Py_ssize_t length; - Py_ssize_t astride; - Py_ssize_t ystride; - Py_ssize_t zstride; - npy_intp i; - npy_intp its; - npy_intp nits; - npy_intp indices[NPY_MAXDIMS]; - npy_intp astrides[NPY_MAXDIMS]; - npy_intp ystrides[NPY_MAXDIMS]; - npy_intp zstrides[NPY_MAXDIMS]; - npy_intp shape[NPY_MAXDIMS]; - char * pa; - char * py; - char * pz; -}; -typedef struct _iter3 iter3; - -static inline void init_iter3(iter3 *it, PyArrayObject *a, PyObject *y, PyObject *z, int axis) { - int i = 0; - int j = 0; - const int ndim = PyArray_NDIM(a); - const npy_intp *shape = PyArray_SHAPE(a); - const npy_intp *astrides = PyArray_STRIDES(a); - const npy_intp *ystrides = PyArray_STRIDES((PyArrayObject *)y); - const npy_intp *zstrides = PyArray_STRIDES((PyArrayObject *)z); - - /* to avoid compiler warning of uninitialized variables */ - it->length = 0; - it->astride = 0; - it->ystride = 0; - it->zstride = 0; - - it->ndim_m2 = ndim - 2; - it->axis = axis; - it->its = 0; - it->nits = 1; - it->pa = PyArray_BYTES(a); - it->py = PyArray_BYTES((PyArrayObject *)y); - it->pz = PyArray_BYTES((PyArrayObject *)z); - - for (i = 0; i < ndim; i++) { - if (i == axis) { - it->astride = astrides[i]; - it->ystride = ystrides[i]; - it->zstride = zstrides[i]; - it->length = shape[i]; - } else { - it->indices[j] = 0; - it->astrides[j] = astrides[i]; - it->ystrides[j] = ystrides[i]; - it->zstrides[j] = zstrides[i]; - it->shape[j] = shape[i]; - it->nits *= shape[i]; - j++; - } - } -} - -#define NEXT3 \ - for (it.i = it.ndim_m2; it.i > -1; it.i--) { \ - if (it.indices[it.i] < it.shape[it.i] - 1) { \ - it.pa += it.astrides[it.i]; \ - it.py += it.ystrides[it.i]; \ - it.pz += it.zstrides[it.i]; \ - it.indices[it.i]++; \ - break; \ - } \ - it.pa -= it.indices[it.i] * it.astrides[it.i]; \ - it.py -= it.indices[it.i] * it.ystrides[it.i]; \ - it.pz -= it.indices[it.i] * it.zstrides[it.i]; \ - it.indices[it.i] = 0; \ - } \ - it.its++; - -/* macros used with iterators -------------------------------------------- */ - -/* most of these macros assume iterator is named `it` */ - -#define NDIM (it.ndim_m2 + 2) -#define SHAPE it.shape -#define SIZE (it.nits * it.length) -#define LENGTH it.length -#define INDEX it.i - -#define WHILE while (it.its < it.nits) -#define WHILE0 \ - it.i = 0; \ - while (it.i < min_count - 1) -#define WHILE1 while (it.i < window) -#define WHILE2 while (it.i < it.length) - -#define FOR for (it.i = 0; it.i < it.length; it.i++) -#define FOR_REVERSE for (it.i = it.length - 1; it.i > -1; it.i--) - -#define RESET it.its = 0; - -#define PA(dtype) (npy_##dtype *)(it.pa) - -#define A0(dtype) *(npy_##dtype *)(it.pa) -#define AI(dtype) *(npy_##dtype *)(it.pa + it.i * it.astride) -#define AX(dtype, x) *(npy_##dtype *)(it.pa + (x)*it.astride) -#define AOLD(dtype) *(npy_##dtype *)(it.pa + (it.i - window) * it.astride) - -#define SI(pa) pa[(it.i * it.stride)] -#define SX(pa, x) pa[(x)*it.stride] - -#define YPP *py++ -#define YI(dtype) *(npy_##dtype *)(it.py + it.i++ * it.ystride) -#define YX(dtype, x) *(npy_##dtype *)(it.py + (x)*it.ystride) - -#define ZX(dtype, x) *(npy_##dtype *)(it.pz + (x)*it.zstride) - -#define FILL_Y(value) \ - npy_intp _i; \ - npy_intp size = PyArray_SIZE((PyArrayObject *)y); \ - for (_i = 0; _i < size; _i++) { \ - YPP = value; \ - } - -#define REDUCE_CONTIGUOUS \ - (it.stride == 1 && ((it.ndim_m2 < 0) || (C_CONTIGUOUS(a) || F_CONTIGUOUS(a)))) -#define ONE_CONTIGUOUS \ - (it.stride == 1 && \ - ((it.ndim_m2 >= 0) && \ - ((C_CONTIGUOUS(a) || F_CONTIGUOUS(a)) && (it.axis == axis)))) -#define ONE_TRANSPOSE(dtype) \ - ((it.ndim_m2 == 0) && it.astrides[it.ndim_m2] == sizeof(dtype) && \ - ((C_CONTIGUOUS(a) && axis == 0) || (F_CONTIGUOUS(a) && axis == 1))) - -#define REDUCE_SPECIALIZE(code) \ - if (REDUCE_CONTIGUOUS) { \ - code \ - } else { \ - code \ - } - -#endif // ITERATORS_H_