From 4f861238803aa9ab913a1011eabed3505a3f703d Mon Sep 17 00:00:00 2001 From: Julian Taylor Date: Wed, 21 Jan 2015 19:54:12 +0100 Subject: MAINT: merge _compiled_base module into multiarray Allows access to internal functions for the file. --- numpy/add_newdocs.py | 18 +- numpy/bento.info | 2 +- numpy/core/bscript | 1 + numpy/core/setup.py | 1 + numpy/core/src/multiarray/compiled_base.c | 1585 ++++++++++++++++++ numpy/core/src/multiarray/compiled_base.h | 24 + numpy/core/src/multiarray/multiarraymodule.c | 23 +- .../core/src/multiarray/multiarraymodule_onefile.c | 1 + numpy/core/src/umath/umathmodule.c | 45 + numpy/lib/bento.info | 6 - numpy/lib/bscript | 11 - numpy/lib/function_base.py | 6 +- numpy/lib/index_tricks.py | 2 +- numpy/lib/npyio.py | 2 +- numpy/lib/setup.py | 22 - numpy/lib/src/_compiled_base.c | 1711 -------------------- 16 files changed, 1693 insertions(+), 1767 deletions(-) create mode 100644 numpy/core/src/multiarray/compiled_base.c create mode 100644 numpy/core/src/multiarray/compiled_base.h delete mode 100644 numpy/lib/bento.info delete mode 100644 numpy/lib/bscript delete mode 100644 numpy/lib/setup.py delete mode 100644 numpy/lib/src/_compiled_base.c diff --git a/numpy/add_newdocs.py b/numpy/add_newdocs.py index 73efdb6a9..f5809a105 100644 --- a/numpy/add_newdocs.py +++ b/numpy/add_newdocs.py @@ -4832,11 +4832,11 @@ add_newdoc('numpy.core.umath', 'seterrobj', ############################################################################## # -# lib._compiled_base functions +# compiled_base functions # ############################################################################## -add_newdoc('numpy.lib._compiled_base', 'digitize', +add_newdoc('numpy.core.multiarray', 'digitize', """ digitize(x, bins, right=False) @@ -4916,7 +4916,7 @@ add_newdoc('numpy.lib._compiled_base', 'digitize', array([1, 3, 3, 4, 5]) """) -add_newdoc('numpy.lib._compiled_base', 'bincount', +add_newdoc('numpy.core.multiarray', 'bincount', """ bincount(x, weights=None, minlength=None) @@ -4989,7 +4989,7 @@ add_newdoc('numpy.lib._compiled_base', 'bincount', """) -add_newdoc('numpy.lib._compiled_base', 'ravel_multi_index', +add_newdoc('numpy.core.multiarray', 'ravel_multi_index', """ ravel_multi_index(multi_index, dims, mode='raise', order='C') @@ -5046,7 +5046,7 @@ add_newdoc('numpy.lib._compiled_base', 'ravel_multi_index', 1621 """) -add_newdoc('numpy.lib._compiled_base', 'unravel_index', +add_newdoc('numpy.core.multiarray', 'unravel_index', """ unravel_index(indices, dims, order='C') @@ -5089,7 +5089,7 @@ add_newdoc('numpy.lib._compiled_base', 'unravel_index', """) -add_newdoc('numpy.lib._compiled_base', 'add_docstring', +add_newdoc('numpy.core.multiarray', 'add_docstring', """ add_docstring(obj, docstring) @@ -5099,7 +5099,7 @@ add_newdoc('numpy.lib._compiled_base', 'add_docstring', raise a TypeError """) -add_newdoc('numpy.lib._compiled_base', 'add_newdoc_ufunc', +add_newdoc('numpy.core.umath', '_add_newdoc_ufunc', """ add_ufunc_docstring(ufunc, new_docstring) @@ -5125,7 +5125,7 @@ add_newdoc('numpy.lib._compiled_base', 'add_newdoc_ufunc', and then throwing away the ufunc. """) -add_newdoc('numpy.lib._compiled_base', 'packbits', +add_newdoc('numpy.core.multiarray', 'packbits', """ packbits(myarray, axis=None) @@ -5169,7 +5169,7 @@ add_newdoc('numpy.lib._compiled_base', 'packbits', """) -add_newdoc('numpy.lib._compiled_base', 'unpackbits', +add_newdoc('numpy.core.multiarray', 'unpackbits', """ unpackbits(myarray, axis=None) diff --git a/numpy/bento.info b/numpy/bento.info index 1bed9c711..9b91fe3ab 100644 --- a/numpy/bento.info +++ b/numpy/bento.info @@ -1,5 +1,5 @@ Recurse: - core, fft, lib, linalg, random + core, fft, linalg, random Library: Packages: diff --git a/numpy/core/bscript b/numpy/core/bscript index 0f0f6816a..48d18b885 100644 --- a/numpy/core/bscript +++ b/numpy/core/bscript @@ -440,6 +440,7 @@ def pre_build(context): pjoin('src', 'multiarray', 'buffer.c'), pjoin('src', 'multiarray', 'calculation.c'), pjoin('src', 'multiarray', 'common.c'), + pjoin('src', 'multiarray', 'compiled_base.c'), pjoin('src', 'multiarray', 'conversion_utils.c'), pjoin('src', 'multiarray', 'convert.c'), pjoin('src', 'multiarray', 'convert_datatype.c'), diff --git a/numpy/core/setup.py b/numpy/core/setup.py index f178c4032..75d64d81b 100644 --- a/numpy/core/setup.py +++ b/numpy/core/setup.py @@ -797,6 +797,7 @@ def configuration(parent_package='',top_path=None): join('src', 'multiarray', 'array_assign_array.c'), join('src', 'multiarray', 'buffer.c'), join('src', 'multiarray', 'calculation.c'), + join('src', 'multiarray', 'compiled_base.c'), join('src', 'multiarray', 'common.c'), join('src', 'multiarray', 'convert.c'), join('src', 'multiarray', 'convert_datatype.c'), diff --git a/numpy/core/src/multiarray/compiled_base.c b/numpy/core/src/multiarray/compiled_base.c new file mode 100644 index 000000000..d258d1ad8 --- /dev/null +++ b/numpy/core/src/multiarray/compiled_base.c @@ -0,0 +1,1585 @@ +#define NPY_NO_DEPRECATED_API NPY_API_VERSION +#include +#include +#include + +#define _MULTIARRAYMODULE +#include "numpy/arrayobject.h" +#include "numpy/npy_3kcompat.h" +#include "numpy/npy_math.h" +#include "npy_config.h" + + +/* + * Returns -1 if the array is monotonic decreasing, + * +1 if the array is monotonic increasing, + * and 0 if the array is not monotonic. + */ +static int +check_array_monotonic(const double *a, npy_int lena) +{ + npy_intp i; + double next; + double last = a[0]; + + /* Skip repeated values at the beginning of the array */ + for (i = 1; (i < lena) && (a[i] == last); i++); + + if (i == lena) { + /* all bin edges hold the same value */ + return 1; + } + + next = a[i]; + if (last < next) { + /* Possibly monotonic increasing */ + for (i += 1; i < lena; i++) { + last = next; + next = a[i]; + if (last > next) { + return 0; + } + } + return 1; + } + else { + /* last > next, possibly monotonic decreasing */ + for (i += 1; i < lena; i++) { + last = next; + next = a[i]; + if (last < next) { + return 0; + } + } + return -1; + } +} + +/* Find the minimum and maximum of an integer array */ +static void +minmax(const npy_intp *data, npy_intp data_len, npy_intp *mn, npy_intp *mx) +{ + npy_intp min = *data; + npy_intp max = *data; + + while (--data_len) { + const npy_intp val = *(++data); + if (val < min) { + min = val; + } + else if (val > max) { + max = val; + } + } + + *mn = min; + *mx = max; +} + +/* + * arr_bincount is registered as bincount. + * + * bincount accepts one, two or three arguments. The first is an array of + * non-negative integers The second, if present, is an array of weights, + * which must be promotable to double. Call these arguments list and + * weight. Both must be one-dimensional with len(weight) == len(list). If + * weight is not present then bincount(list)[i] is the number of occurrences + * of i in list. If weight is present then bincount(self,list, weight)[i] + * is the sum of all weight[j] where list [j] == i. Self is not used. + * The third argument, if present, is a minimum length desired for the + * output array. + */ +NPY_NO_EXPORT PyObject * +arr_bincount(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds) +{ + PyArray_Descr *type; + PyObject *list = NULL, *weight=Py_None, *mlength=Py_None; + PyArrayObject *lst=NULL, *ans=NULL, *wts=NULL; + npy_intp *numbers, *ians, len , mx, mn, ans_size, minlength; + int i; + double *weights , *dans; + static char *kwlist[] = {"list", "weights", "minlength", NULL}; + + if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|OO", + kwlist, &list, &weight, &mlength)) { + goto fail; + } + + lst = (PyArrayObject *)PyArray_ContiguousFromAny(list, NPY_INTP, 1, 1); + if (lst == NULL) { + goto fail; + } + len = PyArray_SIZE(lst); + type = PyArray_DescrFromType(NPY_INTP); + + if (mlength == Py_None) { + minlength = 0; + } + else { + minlength = PyArray_PyIntAsIntp(mlength); + if (minlength <= 0) { + if (!PyErr_Occurred()) { + PyErr_SetString(PyExc_ValueError, + "minlength must be positive"); + } + goto fail; + } + } + + /* handle empty list */ + if (len == 0) { + if (!(ans = (PyArrayObject *)PyArray_Zeros(1, &minlength, type, 0))){ + goto fail; + } + Py_DECREF(lst); + return (PyObject *)ans; + } + + numbers = (npy_intp *) PyArray_DATA(lst); + minmax(numbers, len, &mn, &mx); + if (mn < 0) { + PyErr_SetString(PyExc_ValueError, + "The first argument of bincount must be non-negative"); + goto fail; + } + ans_size = mx + 1; + if (mlength != Py_None) { + if (ans_size < minlength) { + ans_size = minlength; + } + } + if (weight == Py_None) { + ans = (PyArrayObject *)PyArray_Zeros(1, &ans_size, type, 0); + if (ans == NULL) { + goto fail; + } + ians = (npy_intp *)(PyArray_DATA(ans)); + NPY_BEGIN_ALLOW_THREADS; + for (i = 0; i < len; i++) + ians [numbers [i]] += 1; + NPY_END_ALLOW_THREADS; + Py_DECREF(lst); + } + else { + wts = (PyArrayObject *)PyArray_ContiguousFromAny( + weight, NPY_DOUBLE, 1, 1); + if (wts == NULL) { + goto fail; + } + weights = (double *)PyArray_DATA (wts); + if (PyArray_SIZE(wts) != len) { + PyErr_SetString(PyExc_ValueError, + "The weights and list don't have the same length."); + goto fail; + } + type = PyArray_DescrFromType(NPY_DOUBLE); + ans = (PyArrayObject *)PyArray_Zeros(1, &ans_size, type, 0); + if (ans == NULL) { + goto fail; + } + dans = (double *)PyArray_DATA(ans); + NPY_BEGIN_ALLOW_THREADS; + for (i = 0; i < len; i++) { + dans[numbers[i]] += weights[i]; + } + NPY_END_ALLOW_THREADS; + Py_DECREF(lst); + Py_DECREF(wts); + } + return (PyObject *)ans; + +fail: + Py_XDECREF(lst); + Py_XDECREF(wts); + Py_XDECREF(ans); + return NULL; +} + +/* + * digitize(x, bins, right=False) returns an array of integers the same length + * as x. The values i returned are such that bins[i - 1] <= x < bins[i] if + * bins is monotonically increasing, or bins[i - 1] > x >= bins[i] if bins + * is monotonically decreasing. Beyond the bounds of bins, returns either + * i = 0 or i = len(bins) as appropriate. If right == True the comparison + * is bins [i - 1] < x <= bins[i] or bins [i - 1] >= x > bins[i] + */ +NPY_NO_EXPORT PyObject * +arr_digitize(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds) +{ + PyObject *obj_x = NULL; + PyObject *obj_bins = NULL; + PyArrayObject *arr_x = NULL; + PyArrayObject *arr_bins = NULL; + PyObject *ret = NULL; + npy_intp len_bins; + int monotonic, right = 0; + NPY_BEGIN_THREADS_DEF + + static char *kwlist[] = {"x", "bins", "right", NULL}; + + if (!PyArg_ParseTupleAndKeywords(args, kwds, "OO|i", kwlist, + &obj_x, &obj_bins, &right)) { + goto fail; + } + + /* PyArray_SearchSorted will make `x` contiguous even if we don't */ + arr_x = (PyArrayObject *)PyArray_FROMANY(obj_x, NPY_DOUBLE, 0, 0, + NPY_ARRAY_CARRAY_RO); + if (arr_x == NULL) { + goto fail; + } + + /* TODO: `bins` could be strided, needs change to check_array_monotonic */ + arr_bins = (PyArrayObject *)PyArray_FROMANY(obj_bins, NPY_DOUBLE, 1, 1, + NPY_ARRAY_CARRAY_RO); + if (arr_bins == NULL) { + goto fail; + } + + len_bins = PyArray_SIZE(arr_bins); + if (len_bins == 0) { + PyErr_SetString(PyExc_ValueError, "bins must have non-zero length"); + goto fail; + } + + NPY_BEGIN_THREADS_THRESHOLDED(len_bins) + monotonic = check_array_monotonic((const double *)PyArray_DATA(arr_bins), + len_bins); + NPY_END_THREADS + + if (monotonic == 0) { + PyErr_SetString(PyExc_ValueError, + "bins must be monotonically increasing or decreasing"); + goto fail; + } + + /* PyArray_SearchSorted needs an increasing array */ + if (monotonic == - 1) { + PyArrayObject *arr_tmp = NULL; + npy_intp shape = PyArray_DIM(arr_bins, 0); + npy_intp stride = -PyArray_STRIDE(arr_bins, 0); + void *data = (void *)(PyArray_BYTES(arr_bins) - stride * (shape - 1)); + + arr_tmp = (PyArrayObject *)PyArray_New(&PyArray_Type, 1, &shape, + NPY_DOUBLE, &stride, data, 0, + PyArray_FLAGS(arr_bins), NULL); + if (!arr_tmp) { + goto fail; + } + + if (PyArray_SetBaseObject(arr_tmp, (PyObject *)arr_bins) < 0) { + + Py_DECREF(arr_tmp); + goto fail; + } + arr_bins = arr_tmp; + } + + ret = PyArray_SearchSorted(arr_bins, (PyObject *)arr_x, + right ? NPY_SEARCHLEFT : NPY_SEARCHRIGHT, NULL); + if (!ret) { + goto fail; + } + + /* If bins is decreasing, ret has bins from end, not start */ + if (monotonic == -1) { + npy_intp *ret_data = + (npy_intp *)PyArray_DATA((PyArrayObject *)ret); + npy_intp len_ret = PyArray_SIZE((PyArrayObject *)ret); + + NPY_BEGIN_THREADS_THRESHOLDED(len_ret) + while (len_ret--) { + *ret_data = len_bins - *ret_data; + ret_data++; + } + NPY_END_THREADS + } + + fail: + Py_XDECREF(arr_x); + Py_XDECREF(arr_bins); + return ret; +} + +/* + * Insert values from an input array into an output array, at positions + * indicated by a mask. If the arrays are of dtype object (indicated by + * the objarray flag), take care of reference counting. + * + * This function implements the copying logic of arr_insert() defined + * below. + */ +static void +arr_insert_loop(char *mptr, char *vptr, char *input_data, char *zero, + char *avals_data, int melsize, int delsize, int objarray, + int totmask, int numvals, int nd, npy_intp *instrides, + npy_intp *inshape) +{ + int mindx, rem_indx, indx, i, copied; + + /* + * Walk through mask array, when non-zero is encountered + * copy next value in the vals array to the input array. + * If we get through the value array, repeat it as necessary. + */ + copied = 0; + for (mindx = 0; mindx < totmask; mindx++) { + if (memcmp(mptr,zero,melsize) != 0) { + /* compute indx into input array */ + rem_indx = mindx; + indx = 0; + for (i = nd - 1; i > 0; --i) { + indx += (rem_indx % inshape[i]) * instrides[i]; + rem_indx /= inshape[i]; + } + indx += rem_indx * instrides[0]; + /* fprintf(stderr, "mindx = %d, indx=%d\n", mindx, indx); */ + /* Copy value element over to input array */ + memcpy(input_data+indx,vptr,delsize); + if (objarray) { + Py_INCREF(*((PyObject **)vptr)); + } + vptr += delsize; + copied += 1; + /* If we move past value data. Reset */ + if (copied >= numvals) { + vptr = avals_data; + } + } + mptr += melsize; + } +} + +/* + * Returns input array with values inserted sequentially into places + * indicated by the mask + */ +NPY_NO_EXPORT PyObject * +arr_insert(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict) +{ + PyObject *mask = NULL, *vals = NULL; + PyArrayObject *ainput = NULL, *amask = NULL, *avals = NULL, *tmp = NULL; + int numvals, totmask, sameshape; + char *input_data, *mptr, *vptr, *zero = NULL; + int melsize, delsize, nd, objarray, k; + npy_intp *instrides, *inshape; + + static char *kwlist[] = {"input", "mask", "vals", NULL}; + + if (!PyArg_ParseTupleAndKeywords(args, kwdict, "O&OO", kwlist, + PyArray_Converter, &ainput, + &mask, &vals)) { + goto fail; + } + + amask = (PyArrayObject *)PyArray_FROM_OF(mask, NPY_ARRAY_CARRAY); + if (amask == NULL) { + goto fail; + } + /* Cast an object array */ + if (PyArray_DESCR(amask)->type_num == NPY_OBJECT) { + tmp = (PyArrayObject *)PyArray_Cast(amask, NPY_INTP); + if (tmp == NULL) { + goto fail; + } + Py_DECREF(amask); + amask = tmp; + } + + sameshape = 1; + if (PyArray_NDIM(amask) == PyArray_NDIM(ainput)) { + for (k = 0; k < PyArray_NDIM(amask); k++) { + if (PyArray_DIMS(amask)[k] != PyArray_DIMS(ainput)[k]) { + sameshape = 0; + } + } + } + else { + /* Test to see if amask is 1d */ + if (PyArray_NDIM(amask) != 1) { + sameshape = 0; + } + else if ((PyArray_SIZE(ainput)) != PyArray_SIZE(amask)) { + sameshape = 0; + } + } + if (!sameshape) { + PyErr_SetString(PyExc_TypeError, + "mask array must be 1-d or same shape as input array"); + goto fail; + } + + avals = (PyArrayObject *)PyArray_FromObject(vals, + PyArray_DESCR(ainput)->type_num, 0, 1); + if (avals == NULL) { + goto fail; + } + numvals = PyArray_SIZE(avals); + nd = PyArray_NDIM(ainput); + input_data = PyArray_DATA(ainput); + mptr = PyArray_DATA(amask); + melsize = PyArray_DESCR(amask)->elsize; + vptr = PyArray_DATA(avals); + delsize = PyArray_DESCR(avals)->elsize; + zero = PyArray_Zero(amask); + if (zero == NULL) { + goto fail; + } + objarray = (PyArray_DESCR(ainput)->type_num == NPY_OBJECT); + + /* Handle zero-dimensional case separately */ + if (nd == 0) { + if (memcmp(mptr,zero,melsize) != 0) { + /* Copy value element over to input array */ + memcpy(input_data,vptr,delsize); + if (objarray) { + Py_INCREF(*((PyObject **)vptr)); + } + } + Py_DECREF(amask); + Py_DECREF(avals); + PyDataMem_FREE(zero); + Py_DECREF(ainput); + Py_INCREF(Py_None); + return Py_None; + } + + totmask = (int) PyArray_SIZE(amask); + instrides = PyArray_STRIDES(ainput); + inshape = PyArray_DIMS(ainput); + if (objarray) { + /* object array, need to refcount, can't release the GIL */ + arr_insert_loop(mptr, vptr, input_data, zero, PyArray_DATA(avals), + melsize, delsize, objarray, totmask, numvals, nd, + instrides, inshape); + } + else { + /* No increfs take place in arr_insert_loop, so release the GIL */ + NPY_BEGIN_ALLOW_THREADS; + arr_insert_loop(mptr, vptr, input_data, zero, PyArray_DATA(avals), + melsize, delsize, objarray, totmask, numvals, nd, + instrides, inshape); + NPY_END_ALLOW_THREADS; + } + + Py_DECREF(amask); + Py_DECREF(avals); + PyDataMem_FREE(zero); + Py_DECREF(ainput); + Py_INCREF(Py_None); + return Py_None; + +fail: + PyDataMem_FREE(zero); + Py_XDECREF(ainput); + Py_XDECREF(amask); + Py_XDECREF(avals); + return NULL; +} + +/** @brief Use bisection on a sorted array to find first entry > key. + * + * Use bisection to find an index i s.t. arr[i] <= key < arr[i + 1]. If there is + * no such i the error returns are: + * key < arr[0] -- -1 + * key == arr[len - 1] -- len - 1 + * key > arr[len - 1] -- len + * The array is assumed contiguous and sorted in ascending order. + * + * @param key key value. + * @param arr contiguous sorted array to be searched. + * @param len length of the array. + * @return index + */ +static npy_intp +binary_search(double key, double arr [], npy_intp len) +{ + npy_intp imin = 0; + npy_intp imax = len; + + if (key > arr[len - 1]) { + return len; + } + while (imin < imax) { + npy_intp imid = imin + ((imax - imin) >> 1); + if (key >= arr[imid]) { + imin = imid + 1; + } + else { + imax = imid; + } + } + return imin - 1; +} + +NPY_NO_EXPORT PyObject * +arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict) +{ + + PyObject *fp, *xp, *x; + PyObject *left = NULL, *right = NULL; + PyArrayObject *afp = NULL, *axp = NULL, *ax = NULL, *af = NULL; + npy_intp i, lenx, lenxp; + double lval, rval; + double *dy, *dx, *dz, *dres, *slopes; + + static char *kwlist[] = {"x", "xp", "fp", "left", "right", NULL}; + + if (!PyArg_ParseTupleAndKeywords(args, kwdict, "OOO|OO", kwlist, + &x, &xp, &fp, &left, &right)) { + return NULL; + } + + afp = (PyArrayObject *)PyArray_ContiguousFromAny(fp, NPY_DOUBLE, 1, 1); + if (afp == NULL) { + return NULL; + } + axp = (PyArrayObject *)PyArray_ContiguousFromAny(xp, NPY_DOUBLE, 1, 1); + if (axp == NULL) { + goto fail; + } + ax = (PyArrayObject *)PyArray_ContiguousFromAny(x, NPY_DOUBLE, 1, 0); + if (ax == NULL) { + goto fail; + } + lenxp = PyArray_DIMS(axp)[0]; + if (lenxp == 0) { + PyErr_SetString(PyExc_ValueError, + "array of sample points is empty"); + goto fail; + } + if (PyArray_DIMS(afp)[0] != lenxp) { + PyErr_SetString(PyExc_ValueError, + "fp and xp are not of the same length."); + goto fail; + } + + af = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(ax), + PyArray_DIMS(ax), NPY_DOUBLE); + if (af == NULL) { + goto fail; + } + lenx = PyArray_SIZE(ax); + + dy = (double *)PyArray_DATA(afp); + dx = (double *)PyArray_DATA(axp); + dz = (double *)PyArray_DATA(ax); + dres = (double *)PyArray_DATA(af); + + /* Get left and right fill values. */ + if ((left == NULL) || (left == Py_None)) { + lval = dy[0]; + } + else { + lval = PyFloat_AsDouble(left); + if ((lval == -1) && PyErr_Occurred()) { + goto fail; + } + } + if ((right == NULL) || (right == Py_None)) { + rval = dy[lenxp-1]; + } + else { + rval = PyFloat_AsDouble(right); + if ((rval == -1) && PyErr_Occurred()) { + goto fail; + } + } + + /* only pre-calculate slopes if there are relatively few of them. */ + if (lenxp <= lenx) { + slopes = (double *) PyArray_malloc((lenxp - 1)*sizeof(double)); + if (! slopes) { + goto fail; + } + NPY_BEGIN_ALLOW_THREADS; + for (i = 0; i < lenxp - 1; i++) { + slopes[i] = (dy[i + 1] - dy[i])/(dx[i + 1] - dx[i]); + } + for (i = 0; i < lenx; i++) { + const double x = dz[i]; + npy_intp j; + + if (npy_isnan(x)) { + dres[i] = x; + continue; + } + + j = binary_search(x, dx, lenxp); + if (j == -1) { + dres[i] = lval; + } + else if (j == lenxp - 1) { + dres[i] = dy[j]; + } + else if (j == lenxp) { + dres[i] = rval; + } + else { + dres[i] = slopes[j]*(x - dx[j]) + dy[j]; + } + } + NPY_END_ALLOW_THREADS; + PyArray_free(slopes); + } + else { + NPY_BEGIN_ALLOW_THREADS; + for (i = 0; i < lenx; i++) { + const double x = dz[i]; + npy_intp j; + + if (npy_isnan(x)) { + dres[i] = x; + continue; + } + + j = binary_search(x, dx, lenxp); + if (j == -1) { + dres[i] = lval; + } + else if (j == lenxp - 1) { + dres[i] = dy[j]; + } + else if (j == lenxp) { + dres[i] = rval; + } + else { + const double slope = (dy[j + 1] - dy[j])/(dx[j + 1] - dx[j]); + dres[i] = slope*(x - dx[j]) + dy[j]; + } + } + NPY_END_ALLOW_THREADS; + } + + Py_DECREF(afp); + Py_DECREF(axp); + Py_DECREF(ax); + return (PyObject *)af; + +fail: + Py_XDECREF(afp); + Py_XDECREF(axp); + Py_XDECREF(ax); + Py_XDECREF(af); + return NULL; +} + +/* + * Converts a Python sequence into 'count' PyArrayObjects + * + * seq - Input Python object, usually a tuple but any sequence works. + * op - Where the arrays are placed. + * count - How many arrays there should be (errors if it doesn't match). + * paramname - The name of the parameter that produced 'seq'. + */ +static int sequence_to_arrays(PyObject *seq, + PyArrayObject **op, int count, + char *paramname) +{ + int i; + + if (!PySequence_Check(seq) || PySequence_Size(seq) != count) { + PyErr_Format(PyExc_ValueError, + "parameter %s must be a sequence of length %d", + paramname, count); + return -1; + } + + for (i = 0; i < count; ++i) { + PyObject *item = PySequence_GetItem(seq, i); + if (item == NULL) { + while (--i >= 0) { + Py_DECREF(op[i]); + op[i] = NULL; + } + return -1; + } + + op[i] = (PyArrayObject *)PyArray_FromAny(item, NULL, 0, 0, 0, NULL); + if (op[i] == NULL) { + while (--i >= 0) { + Py_DECREF(op[i]); + op[i] = NULL; + } + Py_DECREF(item); + return -1; + } + + Py_DECREF(item); + } + + return 0; +} + +/* Inner loop for unravel_index */ +static int +ravel_multi_index_loop(int ravel_ndim, npy_intp *ravel_dims, + npy_intp *ravel_strides, + npy_intp count, + NPY_CLIPMODE *modes, + char **coords, npy_intp *coords_strides) +{ + int i; + char invalid; + npy_intp j, m; + + NPY_BEGIN_ALLOW_THREADS; + invalid = 0; + while (count--) { + npy_intp raveled = 0; + for (i = 0; i < ravel_ndim; ++i) { + m = ravel_dims[i]; + j = *(npy_intp *)coords[i]; + switch (modes[i]) { + case NPY_RAISE: + if (j < 0 || j >= m) { + invalid = 1; + goto end_while; + } + break; + case NPY_WRAP: + if (j < 0) { + j += m; + if (j < 0) { + j = j % m; + if (j != 0) { + j += m; + } + } + } + else if (j >= m) { + j -= m; + if (j >= m) { + j = j % m; + } + } + break; + case NPY_CLIP: + if (j < 0) { + j = 0; + } + else if (j >= m) { + j = m - 1; + } + break; + + } + raveled += j * ravel_strides[i]; + + coords[i] += coords_strides[i]; + } + *(npy_intp *)coords[ravel_ndim] = raveled; + coords[ravel_ndim] += coords_strides[ravel_ndim]; + } +end_while: + NPY_END_ALLOW_THREADS; + if (invalid) { + PyErr_SetString(PyExc_ValueError, + "invalid entry in coordinates array"); + return NPY_FAIL; + } + return NPY_SUCCEED; +} + +/* ravel_multi_index implementation - see add_newdocs.py */ +NPY_NO_EXPORT PyObject * +arr_ravel_multi_index(PyObject *self, PyObject *args, PyObject *kwds) +{ + int i, s; + PyObject *mode0=NULL, *coords0=NULL; + PyArrayObject *ret = NULL; + PyArray_Dims dimensions={0,0}; + npy_intp ravel_strides[NPY_MAXDIMS]; + NPY_ORDER order = NPY_CORDER; + NPY_CLIPMODE modes[NPY_MAXDIMS]; + + PyArrayObject *op[NPY_MAXARGS]; + PyArray_Descr *dtype[NPY_MAXARGS]; + npy_uint32 op_flags[NPY_MAXARGS]; + + NpyIter *iter = NULL; + + char *kwlist[] = {"multi_index", "dims", "mode", "order", NULL}; + + memset(op, 0, sizeof(op)); + dtype[0] = NULL; + + if (!PyArg_ParseTupleAndKeywords(args, kwds, + "OO&|OO&:ravel_multi_index", kwlist, + &coords0, + PyArray_IntpConverter, &dimensions, + &mode0, + PyArray_OrderConverter, &order)) { + goto fail; + } + + if (dimensions.len+1 > NPY_MAXARGS) { + PyErr_SetString(PyExc_ValueError, + "too many dimensions passed to ravel_multi_index"); + goto fail; + } + + if (!PyArray_ConvertClipmodeSequence(mode0, modes, dimensions.len)) { + goto fail; + } + + switch (order) { + case NPY_CORDER: + s = 1; + for (i = dimensions.len-1; i >= 0; --i) { + ravel_strides[i] = s; + s *= dimensions.ptr[i]; + } + break; + case NPY_FORTRANORDER: + s = 1; + for (i = 0; i < dimensions.len; ++i) { + ravel_strides[i] = s; + s *= dimensions.ptr[i]; + } + break; + default: + PyErr_SetString(PyExc_ValueError, + "only 'C' or 'F' order is permitted"); + goto fail; + } + + /* Get the multi_index into op */ + if (sequence_to_arrays(coords0, op, dimensions.len, "multi_index") < 0) { + goto fail; + } + + + for (i = 0; i < dimensions.len; ++i) { + op_flags[i] = NPY_ITER_READONLY| + NPY_ITER_ALIGNED; + } + op_flags[dimensions.len] = NPY_ITER_WRITEONLY| + NPY_ITER_ALIGNED| + NPY_ITER_ALLOCATE; + dtype[0] = PyArray_DescrFromType(NPY_INTP); + for (i = 1; i <= dimensions.len; ++i) { + dtype[i] = dtype[0]; + } + + iter = NpyIter_MultiNew(dimensions.len+1, op, NPY_ITER_BUFFERED| + NPY_ITER_EXTERNAL_LOOP| + NPY_ITER_ZEROSIZE_OK, + NPY_KEEPORDER, + NPY_SAME_KIND_CASTING, + op_flags, dtype); + if (iter == NULL) { + goto fail; + } + + if (NpyIter_GetIterSize(iter) != 0) { + NpyIter_IterNextFunc *iternext; + char **dataptr; + npy_intp *strides; + npy_intp *countptr; + + iternext = NpyIter_GetIterNext(iter, NULL); + if (iternext == NULL) { + goto fail; + } + dataptr = NpyIter_GetDataPtrArray(iter); + strides = NpyIter_GetInnerStrideArray(iter); + countptr = NpyIter_GetInnerLoopSizePtr(iter); + + do { + if (ravel_multi_index_loop(dimensions.len, dimensions.ptr, + ravel_strides, *countptr, modes, + dataptr, strides) != NPY_SUCCEED) { + goto fail; + } + } while(iternext(iter)); + } + + ret = NpyIter_GetOperandArray(iter)[dimensions.len]; + Py_INCREF(ret); + + Py_DECREF(dtype[0]); + for (i = 0; i < dimensions.len; ++i) { + Py_XDECREF(op[i]); + } + PyDimMem_FREE(dimensions.ptr); + NpyIter_Deallocate(iter); + return PyArray_Return(ret); + +fail: + Py_XDECREF(dtype[0]); + for (i = 0; i < dimensions.len; ++i) { + Py_XDECREF(op[i]); + } + PyDimMem_FREE(dimensions.ptr); + NpyIter_Deallocate(iter); + return NULL; +} + +/* C-order inner loop for unravel_index */ +static int +unravel_index_loop_corder(int unravel_ndim, npy_intp *unravel_dims, + npy_intp unravel_size, npy_intp count, + char *indices, npy_intp indices_stride, + npy_intp *coords) +{ + int i; + char invalid; + npy_intp val; + + NPY_BEGIN_ALLOW_THREADS; + invalid = 0; + while (count--) { + val = *(npy_intp *)indices; + if (val < 0 || val >= unravel_size) { + invalid = 1; + break; + } + for (i = unravel_ndim-1; i >= 0; --i) { + coords[i] = val % unravel_dims[i]; + val /= unravel_dims[i]; + } + coords += unravel_ndim; + indices += indices_stride; + } + NPY_END_ALLOW_THREADS; + if (invalid) { + PyErr_SetString(PyExc_ValueError, + "invalid entry in index array"); + return NPY_FAIL; + } + return NPY_SUCCEED; +} + +/* Fortran-order inner loop for unravel_index */ +static int +unravel_index_loop_forder(int unravel_ndim, npy_intp *unravel_dims, + npy_intp unravel_size, npy_intp count, + char *indices, npy_intp indices_stride, + npy_intp *coords) +{ + int i; + char invalid; + npy_intp val; + + NPY_BEGIN_ALLOW_THREADS; + invalid = 0; + while (count--) { + val = *(npy_intp *)indices; + if (val < 0 || val >= unravel_size) { + invalid = 1; + break; + } + for (i = 0; i < unravel_ndim; ++i) { + *coords++ = val % unravel_dims[i]; + val /= unravel_dims[i]; + } + indices += indices_stride; + } + NPY_END_ALLOW_THREADS; + if (invalid) { + PyErr_SetString(PyExc_ValueError, + "invalid entry in index array"); + return NPY_FAIL; + } + return NPY_SUCCEED; +} + +/* unravel_index implementation - see add_newdocs.py */ +NPY_NO_EXPORT PyObject * +arr_unravel_index(PyObject *self, PyObject *args, PyObject *kwds) +{ + PyObject *indices0 = NULL, *ret_tuple = NULL; + PyArrayObject *ret_arr = NULL; + PyArrayObject *indices = NULL; + PyArray_Descr *dtype = NULL; + PyArray_Dims dimensions={0,0}; + NPY_ORDER order = NPY_CORDER; + npy_intp unravel_size; + + NpyIter *iter = NULL; + int i, ret_ndim; + npy_intp ret_dims[NPY_MAXDIMS], ret_strides[NPY_MAXDIMS]; + + char *kwlist[] = {"indices", "dims", "order", NULL}; + + if (!PyArg_ParseTupleAndKeywords(args, kwds, "OO&|O&:unravel_index", + kwlist, + &indices0, + PyArray_IntpConverter, &dimensions, + PyArray_OrderConverter, &order)) { + goto fail; + } + + if (dimensions.len == 0) { + PyErr_SetString(PyExc_ValueError, + "dims must have at least one value"); + goto fail; + } + + unravel_size = PyArray_MultiplyList(dimensions.ptr, dimensions.len); + + if (!PyArray_Check(indices0)) { + indices = (PyArrayObject*)PyArray_FromAny(indices0, + NULL, 0, 0, 0, NULL); + if (indices == NULL) { + goto fail; + } + } + else { + indices = (PyArrayObject *)indices0; + Py_INCREF(indices); + } + + dtype = PyArray_DescrFromType(NPY_INTP); + if (dtype == NULL) { + goto fail; + } + + iter = NpyIter_New(indices, NPY_ITER_READONLY| + NPY_ITER_ALIGNED| + NPY_ITER_BUFFERED| + NPY_ITER_ZEROSIZE_OK| + NPY_ITER_DONT_NEGATE_STRIDES| + NPY_ITER_MULTI_INDEX, + NPY_KEEPORDER, NPY_SAME_KIND_CASTING, + dtype); + if (iter == NULL) { + goto fail; + } + + /* + * Create the return array with a layout compatible with the indices + * and with a dimension added to the end for the multi-index + */ + ret_ndim = PyArray_NDIM(indices) + 1; + if (NpyIter_GetShape(iter, ret_dims) != NPY_SUCCEED) { + goto fail; + } + ret_dims[ret_ndim-1] = dimensions.len; + if (NpyIter_CreateCompatibleStrides(iter, + dimensions.len*sizeof(npy_intp), ret_strides) != NPY_SUCCEED) { + goto fail; + } + ret_strides[ret_ndim-1] = sizeof(npy_intp); + + /* Remove the multi-index and inner loop */ + if (NpyIter_RemoveMultiIndex(iter) != NPY_SUCCEED) { + goto fail; + } + if (NpyIter_EnableExternalLoop(iter) != NPY_SUCCEED) { + goto fail; + } + + ret_arr = (PyArrayObject *)PyArray_NewFromDescr(&PyArray_Type, dtype, + ret_ndim, ret_dims, ret_strides, NULL, 0, NULL); + dtype = NULL; + if (ret_arr == NULL) { + goto fail; + } + + if (order == NPY_CORDER) { + if (NpyIter_GetIterSize(iter) != 0) { + NpyIter_IterNextFunc *iternext; + char **dataptr; + npy_intp *strides; + npy_intp *countptr, count; + npy_intp *coordsptr = (npy_intp *)PyArray_DATA(ret_arr); + + iternext = NpyIter_GetIterNext(iter, NULL); + if (iternext == NULL) { + goto fail; + } + dataptr = NpyIter_GetDataPtrArray(iter); + strides = NpyIter_GetInnerStrideArray(iter); + countptr = NpyIter_GetInnerLoopSizePtr(iter); + + do { + count = *countptr; + if (unravel_index_loop_corder(dimensions.len, dimensions.ptr, + unravel_size, count, *dataptr, *strides, + coordsptr) != NPY_SUCCEED) { + goto fail; + } + coordsptr += count*dimensions.len; + } while(iternext(iter)); + } + } + else if (order == NPY_FORTRANORDER) { + if (NpyIter_GetIterSize(iter) != 0) { + NpyIter_IterNextFunc *iternext; + char **dataptr; + npy_intp *strides; + npy_intp *countptr, count; + npy_intp *coordsptr = (npy_intp *)PyArray_DATA(ret_arr); + + iternext = NpyIter_GetIterNext(iter, NULL); + if (iternext == NULL) { + goto fail; + } + dataptr = NpyIter_GetDataPtrArray(iter); + strides = NpyIter_GetInnerStrideArray(iter); + countptr = NpyIter_GetInnerLoopSizePtr(iter); + + do { + count = *countptr; + if (unravel_index_loop_forder(dimensions.len, dimensions.ptr, + unravel_size, count, *dataptr, *strides, + coordsptr) != NPY_SUCCEED) { + goto fail; + } + coordsptr += count*dimensions.len; + } while(iternext(iter)); + } + } + else { + PyErr_SetString(PyExc_ValueError, + "only 'C' or 'F' order is permitted"); + goto fail; + } + + /* Now make a tuple of views, one per index */ + ret_tuple = PyTuple_New(dimensions.len); + if (ret_tuple == NULL) { + goto fail; + } + for (i = 0; i < dimensions.len; ++i) { + PyArrayObject *view; + + view = (PyArrayObject *)PyArray_New(&PyArray_Type, ret_ndim-1, + ret_dims, NPY_INTP, + ret_strides, + PyArray_BYTES(ret_arr) + i*sizeof(npy_intp), + 0, 0, NULL); + if (view == NULL) { + goto fail; + } + Py_INCREF(ret_arr); + if (PyArray_SetBaseObject(view, (PyObject *)ret_arr) < 0) { + Py_DECREF(view); + goto fail; + } + PyTuple_SET_ITEM(ret_tuple, i, PyArray_Return(view)); + } + + Py_DECREF(ret_arr); + Py_XDECREF(indices); + PyDimMem_FREE(dimensions.ptr); + NpyIter_Deallocate(iter); + + return ret_tuple; + +fail: + Py_XDECREF(ret_tuple); + Py_XDECREF(ret_arr); + Py_XDECREF(dtype); + Py_XDECREF(indices); + PyDimMem_FREE(dimensions.ptr); + NpyIter_Deallocate(iter); + return NULL; +} + + +/* Can only be called if doc is currently NULL */ +NPY_NO_EXPORT PyObject * +arr_add_docstring(PyObject *NPY_UNUSED(dummy), PyObject *args) +{ + PyObject *obj; + PyObject *str; + char *docstr; + static char *msg = "already has a docstring"; + PyObject *tp_dict; + PyObject *myobj; + static PyTypeObject *PyMemberDescr_TypePtr = NULL; + static PyTypeObject *PyGetSetDescr_TypePtr = NULL; + static PyTypeObject *PyMethodDescr_TypePtr = NULL; + + /* Don't add docstrings */ + if (Py_OptimizeFlag > 1) { + Py_INCREF(Py_None); + return Py_None; + } + + if (PyGetSetDescr_TypePtr == NULL) { + tp_dict = PyArrayDescr_Type.tp_dict; + /* Get "subdescr" */ + myobj = PyDict_GetItemString(tp_dict, "fields"); + if (myobj != NULL) { + PyGetSetDescr_TypePtr = Py_TYPE(myobj); + } + } + if (PyMemberDescr_TypePtr == NULL) { + myobj = PyDict_GetItemString(tp_dict, "alignment"); + if (myobj != NULL) { + PyMemberDescr_TypePtr = Py_TYPE(myobj); + } + } + if (PyMethodDescr_TypePtr == NULL) { + myobj = PyDict_GetItemString(tp_dict, "newbyteorder"); + if (myobj != NULL) { + PyMethodDescr_TypePtr = Py_TYPE(myobj); + } + } + +#if defined(NPY_PY3K) + if (!PyArg_ParseTuple(args, "OO!", &obj, &PyUnicode_Type, &str)) { + return NULL; + } + + docstr = PyBytes_AS_STRING(PyUnicode_AsUTF8String(str)); +#else + if (!PyArg_ParseTuple(args, "OO!", &obj, &PyString_Type, &str)) { + return NULL; + } + + docstr = PyString_AS_STRING(str); +#endif + +#define _TESTDOC1(typebase) (Py_TYPE(obj) == &Py##typebase##_Type) +#define _TESTDOC2(typebase) (Py_TYPE(obj) == Py##typebase##_TypePtr) +#define _ADDDOC(typebase, doc, name) do { \ + Py##typebase##Object *new = (Py##typebase##Object *)obj; \ + if (!(doc)) { \ + doc = docstr; \ + } \ + else { \ + PyErr_Format(PyExc_RuntimeError, "%s method %s", name, msg); \ + return NULL; \ + } \ + } while (0) + + if (_TESTDOC1(CFunction)) { + _ADDDOC(CFunction, new->m_ml->ml_doc, new->m_ml->ml_name); + } + else if (_TESTDOC1(Type)) { + _ADDDOC(Type, new->tp_doc, new->tp_name); + } + else if (_TESTDOC2(MemberDescr)) { + _ADDDOC(MemberDescr, new->d_member->doc, new->d_member->name); + } + else if (_TESTDOC2(GetSetDescr)) { + _ADDDOC(GetSetDescr, new->d_getset->doc, new->d_getset->name); + } + else if (_TESTDOC2(MethodDescr)) { + _ADDDOC(MethodDescr, new->d_method->ml_doc, new->d_method->ml_name); + } + else { + PyObject *doc_attr; + + doc_attr = PyObject_GetAttrString(obj, "__doc__"); + if (doc_attr != NULL && doc_attr != Py_None) { + PyErr_Format(PyExc_RuntimeError, "object %s", msg); + return NULL; + } + Py_XDECREF(doc_attr); + + if (PyObject_SetAttrString(obj, "__doc__", str) < 0) { + PyErr_SetString(PyExc_TypeError, + "Cannot set a docstring for that object"); + return NULL; + } + Py_INCREF(Py_None); + return Py_None; + } + +#undef _TESTDOC1 +#undef _TESTDOC2 +#undef _ADDDOC + + Py_INCREF(str); + Py_INCREF(Py_None); + return Py_None; +} + + +/* + * This function packs boolean values in the input array into the bits of a + * byte array. Truth values are determined as usual: 0 is false, everything + * else is true. + */ +static NPY_INLINE void +pack_inner(const char *inptr, + npy_intp element_size, /* in bytes */ + npy_intp n_in, + npy_intp in_stride, + char *outptr, + npy_intp n_out, + npy_intp out_stride) +{ + /* + * Loop through the elements of inptr. + * Determine whether or not it is nonzero. + * Yes: set corresponding bit (and adjust build value) + * No: move on + * Every 8th value, set the value of build and increment the outptr + */ + npy_intp index; + int remain = n_in % 8; /* uneven bits */ + + if (remain == 0) { /* assumes n_in > 0 */ + remain = 8; + } + for (index = 0; index < n_out; index++) { + char build = 0; + int i, maxi; + npy_intp j; + + maxi = (index == n_out - 1) ? remain : 8; + for (i = 0; i < maxi; i++) { + build <<= 1; + for (j = 0; j < element_size; j++) { + build |= (inptr[j] != 0); + } + inptr += in_stride; + } + if (index == n_out - 1) { + build <<= 8 - remain; + } + *outptr = build; + outptr += out_stride; + } +} + +static PyObject * +pack_bits(PyObject *input, int axis) +{ + PyArrayObject *inp; + PyArrayObject *new = NULL; + PyArrayObject *out = NULL; + npy_intp outdims[NPY_MAXDIMS]; + int i; + PyArrayIterObject *it, *ot; + NPY_BEGIN_THREADS_DEF; + + inp = (PyArrayObject *)PyArray_FROM_O(input); + + if (inp == NULL) { + return NULL; + } + if (!PyArray_ISBOOL(inp) && !PyArray_ISINTEGER(inp)) { + PyErr_SetString(PyExc_TypeError, + "Expected an input array of integer or boolean data type"); + goto fail; + } + + new = (PyArrayObject *)PyArray_CheckAxis(inp, &axis, 0); + Py_DECREF(inp); + if (new == NULL) { + return NULL; + } + /* Handle empty array separately */ + if (PyArray_SIZE(new) == 0) { + return PyArray_Copy(new); + } + + if (PyArray_NDIM(new) == 0) { + char *optr, *iptr; + + out = (PyArrayObject *)PyArray_New(Py_TYPE(new), 0, NULL, NPY_UBYTE, + NULL, NULL, 0, 0, NULL); + if (out == NULL) { + goto fail; + } + optr = PyArray_DATA(out); + iptr = PyArray_DATA(new); + *optr = 0; + for (i = 0; i < PyArray_ITEMSIZE(new); i++) { + if (*iptr != 0) { + *optr = 1; + break; + } + iptr++; + } + goto finish; + } + + + /* Setup output shape */ + for (i = 0; i < PyArray_NDIM(new); i++) { + outdims[i] = PyArray_DIM(new, i); + } + + /* + * Divide axis dimension by 8 + * 8 -> 1, 9 -> 2, 16 -> 2, 17 -> 3 etc.. + */ + outdims[axis] = ((outdims[axis] - 1) >> 3) + 1; + + /* Create output array */ + out = (PyArrayObject *)PyArray_New(Py_TYPE(new), + PyArray_NDIM(new), outdims, NPY_UBYTE, + NULL, NULL, 0, PyArray_ISFORTRAN(new), NULL); + if (out == NULL) { + goto fail; + } + /* Setup iterators to iterate over all but given axis */ + it = (PyArrayIterObject *)PyArray_IterAllButAxis((PyObject *)new, &axis); + ot = (PyArrayIterObject *)PyArray_IterAllButAxis((PyObject *)out, &axis); + if (it == NULL || ot == NULL) { + Py_XDECREF(it); + Py_XDECREF(ot); + goto fail; + } + + NPY_BEGIN_THREADS_THRESHOLDED(PyArray_DIM(out, axis)); + while (PyArray_ITER_NOTDONE(it)) { + pack_inner(PyArray_ITER_DATA(it), PyArray_ITEMSIZE(new), + PyArray_DIM(new, axis), PyArray_STRIDE(new, axis), + PyArray_ITER_DATA(ot), PyArray_DIM(out, axis), + PyArray_STRIDE(out, axis)); + PyArray_ITER_NEXT(it); + PyArray_ITER_NEXT(ot); + } + NPY_END_THREADS; + + Py_DECREF(it); + Py_DECREF(ot); + +finish: + Py_DECREF(new); + return (PyObject *)out; + +fail: + Py_XDECREF(new); + Py_XDECREF(out); + return NULL; +} + +static PyObject * +unpack_bits(PyObject *input, int axis) +{ + PyArrayObject *inp; + PyArrayObject *new = NULL; + PyArrayObject *out = NULL; + npy_intp outdims[NPY_MAXDIMS]; + int i; + PyArrayIterObject *it, *ot; + npy_intp n_in, in_stride, out_stride; + NPY_BEGIN_THREADS_DEF; + + inp = (PyArrayObject *)PyArray_FROM_O(input); + + if (inp == NULL) { + return NULL; + } + if (PyArray_TYPE(inp) != NPY_UBYTE) { + PyErr_SetString(PyExc_TypeError, + "Expected an input array of unsigned byte data type"); + goto fail; + } + + new = (PyArrayObject *)PyArray_CheckAxis(inp, &axis, 0); + Py_DECREF(inp); + if (new == NULL) { + return NULL; + } + /* Handle zero-dim array separately */ + if (PyArray_SIZE(new) == 0) { + return PyArray_Copy(new); + } + + if (PyArray_NDIM(new) == 0) { + /* Handle 0-d array by converting it to a 1-d array */ + PyArrayObject *temp; + PyArray_Dims newdim = {NULL, 1}; + npy_intp shape = 1; + + newdim.ptr = &shape; + temp = (PyArrayObject *)PyArray_Newshape(new, &newdim, NPY_CORDER); + if (temp == NULL) { + goto fail; + } + Py_DECREF(new); + new = temp; + } + + /* Setup output shape */ + for (i=0; i>= 1; + } + inptr += in_stride; + } + PyArray_ITER_NEXT(it); + PyArray_ITER_NEXT(ot); + } + NPY_END_THREADS; + + Py_DECREF(it); + Py_DECREF(ot); + + Py_DECREF(new); + return (PyObject *)out; + +fail: + Py_XDECREF(new); + Py_XDECREF(out); + return NULL; +} + + +NPY_NO_EXPORT PyObject * +io_pack(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds) +{ + PyObject *obj; + int axis = NPY_MAXDIMS; + static char *kwlist[] = {"in", "axis", NULL}; + + if (!PyArg_ParseTupleAndKeywords( args, kwds, "O|O&" , kwlist, + &obj, PyArray_AxisConverter, &axis)) { + return NULL; + } + return pack_bits(obj, axis); +} + +NPY_NO_EXPORT PyObject * +io_unpack(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds) +{ + PyObject *obj; + int axis = NPY_MAXDIMS; + static char *kwlist[] = {"in", "axis", NULL}; + + if (!PyArg_ParseTupleAndKeywords( args, kwds, "O|O&" , kwlist, + &obj, PyArray_AxisConverter, &axis)) { + return NULL; + } + return unpack_bits(obj, axis); +} diff --git a/numpy/core/src/multiarray/compiled_base.h b/numpy/core/src/multiarray/compiled_base.h new file mode 100644 index 000000000..19e3778ad --- /dev/null +++ b/numpy/core/src/multiarray/compiled_base.h @@ -0,0 +1,24 @@ +#ifndef _NPY_PRIVATE__COMPILED_BASE_H_ +#define _NPY_PRIVATE__COMPILED_BASE_H_ +#include + +NPY_NO_EXPORT PyObject * +arr_insert(PyObject *, PyObject *, PyObject *); +NPY_NO_EXPORT PyObject * +arr_bincount(PyObject *, PyObject *, PyObject *); +NPY_NO_EXPORT PyObject * +arr_digitize(PyObject *, PyObject *, PyObject *kwds); +NPY_NO_EXPORT PyObject * +arr_interp(PyObject *, PyObject *, PyObject *); +NPY_NO_EXPORT PyObject * +arr_ravel_multi_index(PyObject *, PyObject *, PyObject *); +NPY_NO_EXPORT PyObject * +arr_unravel_index(PyObject *, PyObject *, PyObject *); +NPY_NO_EXPORT PyObject * +arr_add_docstring(PyObject *, PyObject *); +NPY_NO_EXPORT PyObject * +io_pack(PyObject *, PyObject *, PyObject *); +NPY_NO_EXPORT PyObject * +io_unpack(PyObject *, PyObject *, PyObject *); + +#endif diff --git a/numpy/core/src/multiarray/multiarraymodule.c b/numpy/core/src/multiarray/multiarraymodule.c index e6b317ed1..0c48e1130 100644 --- a/numpy/core/src/multiarray/multiarraymodule.c +++ b/numpy/core/src/multiarray/multiarraymodule.c @@ -59,6 +59,7 @@ NPY_NO_EXPORT int NPY_NUMUSERTYPES = 0; #include "cblasfuncs.h" #include "vdot.h" #include "templ_common.h" /* for npy_mul_with_overflow_intp */ +#include "compiled_base.h" /* Only here for API compatibility */ NPY_NO_EXPORT PyTypeObject PyBigArray_Type; @@ -3881,8 +3882,6 @@ array_may_share_memory(PyObject *NPY_UNUSED(ignored), PyObject *args) } } - - static struct PyMethodDef array_module_methods[] = { {"_get_ndarray_c_version", (PyCFunction)array__get_ndarray_c_version, @@ -4029,6 +4028,26 @@ static struct PyMethodDef array_module_methods[] = { {"test_interrupt", (PyCFunction)test_interrupt, METH_VARARGS, NULL}, + {"_insert", (PyCFunction)arr_insert, + METH_VARARGS | METH_KEYWORDS, + "Insert vals sequentially into equivalent 1-d positions " + "indicated by mask."}, + {"bincount", (PyCFunction)arr_bincount, + METH_VARARGS | METH_KEYWORDS, NULL}, + {"digitize", (PyCFunction)arr_digitize, + METH_VARARGS | METH_KEYWORDS, NULL}, + {"interp", (PyCFunction)arr_interp, + METH_VARARGS | METH_KEYWORDS, NULL}, + {"ravel_multi_index", (PyCFunction)arr_ravel_multi_index, + METH_VARARGS | METH_KEYWORDS, NULL}, + {"unravel_index", (PyCFunction)arr_unravel_index, + METH_VARARGS | METH_KEYWORDS, NULL}, + {"add_docstring", (PyCFunction)arr_add_docstring, + METH_VARARGS, NULL}, + {"packbits", (PyCFunction)io_pack, + METH_VARARGS | METH_KEYWORDS, NULL}, + {"unpackbits", (PyCFunction)io_unpack, + METH_VARARGS | METH_KEYWORDS, NULL}, {NULL, NULL, 0, NULL} /* sentinel */ }; diff --git a/numpy/core/src/multiarray/multiarraymodule_onefile.c b/numpy/core/src/multiarray/multiarraymodule_onefile.c index 7adfdf9da..3940d009b 100644 --- a/numpy/core/src/multiarray/multiarraymodule_onefile.c +++ b/numpy/core/src/multiarray/multiarraymodule_onefile.c @@ -54,6 +54,7 @@ #include "arrayobject.c" #include "numpymemoryview.c" #include "multiarraymodule.c" +#include "compiled_base.c" #if defined(HAVE_CBLAS) #include "python_xerbla.c" diff --git a/numpy/core/src/umath/umathmodule.c b/numpy/core/src/umath/umathmodule.c index 2c84027f9..d792e8b24 100644 --- a/numpy/core/src/umath/umathmodule.c +++ b/numpy/core/src/umath/umathmodule.c @@ -197,6 +197,49 @@ ufunc_frompyfunc(PyObject *NPY_UNUSED(dummy), PyObject *args, PyObject *NPY_UNUS return (PyObject *)self; } +/* docstring in numpy.add_newdocs.py */ +static PyObject * +add_newdoc_ufunc(PyObject *NPY_UNUSED(dummy), PyObject *args) +{ + PyUFuncObject *ufunc; + PyObject *str; + char *docstr, *newdocstr; + +#if defined(NPY_PY3K) + if (!PyArg_ParseTuple(args, "O!O!", &PyUFunc_Type, &ufunc, + &PyUnicode_Type, &str)) { + return NULL; + } + docstr = PyBytes_AS_STRING(PyUnicode_AsUTF8String(str)); +#else + if (!PyArg_ParseTuple(args, "O!O!", &PyUFunc_Type, &ufunc, + &PyString_Type, &str)) { + return NULL; + } + docstr = PyString_AS_STRING(str); +#endif + + if (NULL != ufunc->doc) { + PyErr_SetString(PyExc_ValueError, + "Cannot change docstring of ufunc with non-NULL docstring"); + return NULL; + } + + /* + * This introduces a memory leak, as the memory allocated for the doc + * will not be freed even if the ufunc itself is deleted. In practice + * this should not be a problem since the user would have to + * repeatedly create, document, and throw away ufuncs. + */ + newdocstr = malloc(strlen(docstr) + 1); + strcpy(newdocstr, docstr); + ufunc->doc = newdocstr; + + Py_INCREF(Py_None); + return Py_None; +} + + /* ***************************************************************************** ** SETUP UFUNCS ** @@ -241,6 +284,8 @@ static struct PyMethodDef methods[] = { {"geterrobj", (PyCFunction) ufunc_geterr, METH_VARARGS, NULL}, + {"_add_newdoc_ufunc", (PyCFunction)add_newdoc_ufunc, + METH_VARARGS, NULL}, {NULL, NULL, 0, NULL} /* sentinel */ }; diff --git a/numpy/lib/bento.info b/numpy/lib/bento.info deleted file mode 100644 index 9f4fa6f0f..000000000 --- a/numpy/lib/bento.info +++ /dev/null @@ -1,6 +0,0 @@ -HookFile: bscript - -Library: - Extension: _compiled_base - Sources: - src/_compiled_base.c diff --git a/numpy/lib/bscript b/numpy/lib/bscript deleted file mode 100644 index 61debfe41..000000000 --- a/numpy/lib/bscript +++ /dev/null @@ -1,11 +0,0 @@ -from bento.commands import hooks - -@hooks.pre_build -def build(context): - context.tweak_extension("_compiled_base", - includes=["../core/include", "../core/include/numpy", "../core", - "../core/src/private"], - defines=['_FILE_OFFSET_BITS=64', - '_LARGEFILE_SOURCE=1', - '_LARGEFILE64_SOURCE=1'] - ) diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py index 135053e43..2baf83830 100644 --- a/numpy/lib/function_base.py +++ b/numpy/lib/function_base.py @@ -23,9 +23,9 @@ from numpy.core.fromnumeric import ( from numpy.core.numerictypes import typecodes, number from numpy.lib.twodim_base import diag from .utils import deprecate -from ._compiled_base import _insert, add_docstring -from ._compiled_base import digitize, bincount, interp as compiled_interp -from ._compiled_base import add_newdoc_ufunc +from numpy.core.multiarray import _insert, add_docstring +from numpy.core.multiarray import digitize, bincount, interp as compiled_interp +from numpy.core.umath import _add_newdoc_ufunc as add_newdoc_ufunc from numpy.compat import long # Force range to be a generator, for np.delete's usage. diff --git a/numpy/lib/index_tricks.py b/numpy/lib/index_tricks.py index f83024961..eb9aad6ad 100644 --- a/numpy/lib/index_tricks.py +++ b/numpy/lib/index_tricks.py @@ -12,7 +12,7 @@ from numpy.core.numerictypes import find_common_type from . import function_base import numpy.matrixlib as matrix from .function_base import diff -from numpy.lib._compiled_base import ravel_multi_index, unravel_index +from numpy.core.multiarray import ravel_multi_index, unravel_index from numpy.lib.stride_tricks import as_strided makemat = matrix.matrix diff --git a/numpy/lib/npyio.py b/numpy/lib/npyio.py index 5f274f27c..40a0b6c80 100644 --- a/numpy/lib/npyio.py +++ b/numpy/lib/npyio.py @@ -11,7 +11,7 @@ from operator import itemgetter import numpy as np from . import format from ._datasource import DataSource -from ._compiled_base import packbits, unpackbits +from numpy.core.multiarray import packbits, unpackbits from ._iotools import ( LineSplitter, NameValidator, StringConverter, ConverterError, ConverterLockError, ConversionWarning, _is_string_like, has_nested_fields, diff --git a/numpy/lib/setup.py b/numpy/lib/setup.py deleted file mode 100644 index 62d1dfbb8..000000000 --- a/numpy/lib/setup.py +++ /dev/null @@ -1,22 +0,0 @@ -from __future__ import division, print_function - -from os.path import join - -def configuration(parent_package='',top_path=None): - from numpy.distutils.misc_util import Configuration - - config = Configuration('lib', parent_package, top_path) - - config.add_include_dirs(join('..', 'core', 'include')) - - config.add_extension('_compiled_base', - sources=[join('src', '_compiled_base.c')] - ) - - config.add_data_dir('tests') - - return config - -if __name__ == '__main__': - from numpy.distutils.core import setup - setup(configuration=configuration) diff --git a/numpy/lib/src/_compiled_base.c b/numpy/lib/src/_compiled_base.c deleted file mode 100644 index 9d9344526..000000000 --- a/numpy/lib/src/_compiled_base.c +++ /dev/null @@ -1,1711 +0,0 @@ -#define NPY_NO_DEPRECATED_API NPY_API_VERSION -#include -#include -#include - -#include "numpy/arrayobject.h" -#include "numpy/npy_3kcompat.h" -#include "npy_config.h" -#include "numpy/ufuncobject.h" - - -/* - * Returns -1 if the array is monotonic decreasing, - * +1 if the array is monotonic increasing, - * and 0 if the array is not monotonic. - */ -static int -check_array_monotonic(const double *a, npy_int lena) -{ - npy_intp i; - double next; - double last = a[0]; - - /* Skip repeated values at the beginning of the array */ - for (i = 1; (i < lena) && (a[i] == last); i++); - - if (i == lena) { - /* all bin edges hold the same value */ - return 1; - } - - next = a[i]; - if (last < next) { - /* Possibly monotonic increasing */ - for (i += 1; i < lena; i++) { - last = next; - next = a[i]; - if (last > next) { - return 0; - } - } - return 1; - } - else { - /* last > next, possibly monotonic decreasing */ - for (i += 1; i < lena; i++) { - last = next; - next = a[i]; - if (last < next) { - return 0; - } - } - return -1; - } -} - -/* Find the minimum and maximum of an integer array */ -static void -minmax(const npy_intp *data, npy_intp data_len, npy_intp *mn, npy_intp *mx) -{ - npy_intp min = *data; - npy_intp max = *data; - - while (--data_len) { - const npy_intp val = *(++data); - if (val < min) { - min = val; - } - else if (val > max) { - max = val; - } - } - - *mn = min; - *mx = max; -} - -/* - * arr_bincount is registered as bincount. - * - * bincount accepts one, two or three arguments. The first is an array of - * non-negative integers The second, if present, is an array of weights, - * which must be promotable to double. Call these arguments list and - * weight. Both must be one-dimensional with len(weight) == len(list). If - * weight is not present then bincount(list)[i] is the number of occurrences - * of i in list. If weight is present then bincount(self,list, weight)[i] - * is the sum of all weight[j] where list [j] == i. Self is not used. - * The third argument, if present, is a minimum length desired for the - * output array. - */ -static PyObject * -arr_bincount(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds) -{ - PyArray_Descr *type; - PyObject *list = NULL, *weight=Py_None, *mlength=Py_None; - PyArrayObject *lst=NULL, *ans=NULL, *wts=NULL; - npy_intp *numbers, *ians, len , mx, mn, ans_size, minlength; - int i; - double *weights , *dans; - static char *kwlist[] = {"list", "weights", "minlength", NULL}; - - if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|OO", - kwlist, &list, &weight, &mlength)) { - goto fail; - } - - lst = (PyArrayObject *)PyArray_ContiguousFromAny(list, NPY_INTP, 1, 1); - if (lst == NULL) { - goto fail; - } - len = PyArray_SIZE(lst); - type = PyArray_DescrFromType(NPY_INTP); - - if (mlength == Py_None) { - minlength = 0; - } - else { - minlength = PyArray_PyIntAsIntp(mlength); - if (minlength <= 0) { - if (!PyErr_Occurred()) { - PyErr_SetString(PyExc_ValueError, - "minlength must be positive"); - } - goto fail; - } - } - - /* handle empty list */ - if (len == 0) { - if (!(ans = (PyArrayObject *)PyArray_Zeros(1, &minlength, type, 0))){ - goto fail; - } - Py_DECREF(lst); - return (PyObject *)ans; - } - - numbers = (npy_intp *) PyArray_DATA(lst); - minmax(numbers, len, &mn, &mx); - if (mn < 0) { - PyErr_SetString(PyExc_ValueError, - "The first argument of bincount must be non-negative"); - goto fail; - } - ans_size = mx + 1; - if (mlength != Py_None) { - if (ans_size < minlength) { - ans_size = minlength; - } - } - if (weight == Py_None) { - ans = (PyArrayObject *)PyArray_Zeros(1, &ans_size, type, 0); - if (ans == NULL) { - goto fail; - } - ians = (npy_intp *)(PyArray_DATA(ans)); - NPY_BEGIN_ALLOW_THREADS; - for (i = 0; i < len; i++) - ians [numbers [i]] += 1; - NPY_END_ALLOW_THREADS; - Py_DECREF(lst); - } - else { - wts = (PyArrayObject *)PyArray_ContiguousFromAny( - weight, NPY_DOUBLE, 1, 1); - if (wts == NULL) { - goto fail; - } - weights = (double *)PyArray_DATA (wts); - if (PyArray_SIZE(wts) != len) { - PyErr_SetString(PyExc_ValueError, - "The weights and list don't have the same length."); - goto fail; - } - type = PyArray_DescrFromType(NPY_DOUBLE); - ans = (PyArrayObject *)PyArray_Zeros(1, &ans_size, type, 0); - if (ans == NULL) { - goto fail; - } - dans = (double *)PyArray_DATA(ans); - NPY_BEGIN_ALLOW_THREADS; - for (i = 0; i < len; i++) { - dans[numbers[i]] += weights[i]; - } - NPY_END_ALLOW_THREADS; - Py_DECREF(lst); - Py_DECREF(wts); - } - return (PyObject *)ans; - -fail: - Py_XDECREF(lst); - Py_XDECREF(wts); - Py_XDECREF(ans); - return NULL; -} - -/* - * digitize(x, bins, right=False) returns an array of integers the same length - * as x. The values i returned are such that bins[i - 1] <= x < bins[i] if - * bins is monotonically increasing, or bins[i - 1] > x >= bins[i] if bins - * is monotonically decreasing. Beyond the bounds of bins, returns either - * i = 0 or i = len(bins) as appropriate. If right == True the comparison - * is bins [i - 1] < x <= bins[i] or bins [i - 1] >= x > bins[i] - */ -static PyObject * -arr_digitize(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds) -{ - PyObject *obj_x = NULL; - PyObject *obj_bins = NULL; - PyArrayObject *arr_x = NULL; - PyArrayObject *arr_bins = NULL; - PyObject *ret = NULL; - npy_intp len_bins; - int monotonic, right = 0; - NPY_BEGIN_THREADS_DEF - - static char *kwlist[] = {"x", "bins", "right", NULL}; - - if (!PyArg_ParseTupleAndKeywords(args, kwds, "OO|i", kwlist, - &obj_x, &obj_bins, &right)) { - goto fail; - } - - /* PyArray_SearchSorted will make `x` contiguous even if we don't */ - arr_x = (PyArrayObject *)PyArray_FROMANY(obj_x, NPY_DOUBLE, 0, 0, - NPY_ARRAY_CARRAY_RO); - if (arr_x == NULL) { - goto fail; - } - - /* TODO: `bins` could be strided, needs change to check_array_monotonic */ - arr_bins = (PyArrayObject *)PyArray_FROMANY(obj_bins, NPY_DOUBLE, 1, 1, - NPY_ARRAY_CARRAY_RO); - if (arr_bins == NULL) { - goto fail; - } - - len_bins = PyArray_SIZE(arr_bins); - if (len_bins == 0) { - PyErr_SetString(PyExc_ValueError, "bins must have non-zero length"); - goto fail; - } - - NPY_BEGIN_THREADS_THRESHOLDED(len_bins) - monotonic = check_array_monotonic((const double *)PyArray_DATA(arr_bins), - len_bins); - NPY_END_THREADS - - if (monotonic == 0) { - PyErr_SetString(PyExc_ValueError, - "bins must be monotonically increasing or decreasing"); - goto fail; - } - - /* PyArray_SearchSorted needs an increasing array */ - if (monotonic == - 1) { - PyArrayObject *arr_tmp = NULL; - npy_intp shape = PyArray_DIM(arr_bins, 0); - npy_intp stride = -PyArray_STRIDE(arr_bins, 0); - void *data = (void *)(PyArray_BYTES(arr_bins) - stride * (shape - 1)); - - arr_tmp = (PyArrayObject *)PyArray_New(&PyArray_Type, 1, &shape, - NPY_DOUBLE, &stride, data, 0, - PyArray_FLAGS(arr_bins), NULL); - if (!arr_tmp) { - goto fail; - } - - if (PyArray_SetBaseObject(arr_tmp, (PyObject *)arr_bins) < 0) { - - Py_DECREF(arr_tmp); - goto fail; - } - arr_bins = arr_tmp; - } - - ret = PyArray_SearchSorted(arr_bins, (PyObject *)arr_x, - right ? NPY_SEARCHLEFT : NPY_SEARCHRIGHT, NULL); - if (!ret) { - goto fail; - } - - /* If bins is decreasing, ret has bins from end, not start */ - if (monotonic == -1) { - npy_intp *ret_data = - (npy_intp *)PyArray_DATA((PyArrayObject *)ret); - npy_intp len_ret = PyArray_SIZE((PyArrayObject *)ret); - - NPY_BEGIN_THREADS_THRESHOLDED(len_ret) - while (len_ret--) { - *ret_data = len_bins - *ret_data; - ret_data++; - } - NPY_END_THREADS - } - - fail: - Py_XDECREF(arr_x); - Py_XDECREF(arr_bins); - return ret; -} - -static char arr_insert__doc__[] = "Insert vals sequentially into equivalent 1-d positions indicated by mask."; - -/* - * Insert values from an input array into an output array, at positions - * indicated by a mask. If the arrays are of dtype object (indicated by - * the objarray flag), take care of reference counting. - * - * This function implements the copying logic of arr_insert() defined - * below. - */ -static void -arr_insert_loop(char *mptr, char *vptr, char *input_data, char *zero, - char *avals_data, int melsize, int delsize, int objarray, - int totmask, int numvals, int nd, npy_intp *instrides, - npy_intp *inshape) -{ - int mindx, rem_indx, indx, i, copied; - - /* - * Walk through mask array, when non-zero is encountered - * copy next value in the vals array to the input array. - * If we get through the value array, repeat it as necessary. - */ - copied = 0; - for (mindx = 0; mindx < totmask; mindx++) { - if (memcmp(mptr,zero,melsize) != 0) { - /* compute indx into input array */ - rem_indx = mindx; - indx = 0; - for (i = nd - 1; i > 0; --i) { - indx += (rem_indx % inshape[i]) * instrides[i]; - rem_indx /= inshape[i]; - } - indx += rem_indx * instrides[0]; - /* fprintf(stderr, "mindx = %d, indx=%d\n", mindx, indx); */ - /* Copy value element over to input array */ - memcpy(input_data+indx,vptr,delsize); - if (objarray) { - Py_INCREF(*((PyObject **)vptr)); - } - vptr += delsize; - copied += 1; - /* If we move past value data. Reset */ - if (copied >= numvals) { - vptr = avals_data; - } - } - mptr += melsize; - } -} - -/* - * Returns input array with values inserted sequentially into places - * indicated by the mask - */ -static PyObject * -arr_insert(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict) -{ - PyObject *mask = NULL, *vals = NULL; - PyArrayObject *ainput = NULL, *amask = NULL, *avals = NULL, *tmp = NULL; - int numvals, totmask, sameshape; - char *input_data, *mptr, *vptr, *zero = NULL; - int melsize, delsize, nd, objarray, k; - npy_intp *instrides, *inshape; - - static char *kwlist[] = {"input", "mask", "vals", NULL}; - - if (!PyArg_ParseTupleAndKeywords(args, kwdict, "O&OO", kwlist, - PyArray_Converter, &ainput, - &mask, &vals)) { - goto fail; - } - - amask = (PyArrayObject *)PyArray_FROM_OF(mask, NPY_ARRAY_CARRAY); - if (amask == NULL) { - goto fail; - } - /* Cast an object array */ - if (PyArray_DESCR(amask)->type_num == NPY_OBJECT) { - tmp = (PyArrayObject *)PyArray_Cast(amask, NPY_INTP); - if (tmp == NULL) { - goto fail; - } - Py_DECREF(amask); - amask = tmp; - } - - sameshape = 1; - if (PyArray_NDIM(amask) == PyArray_NDIM(ainput)) { - for (k = 0; k < PyArray_NDIM(amask); k++) { - if (PyArray_DIMS(amask)[k] != PyArray_DIMS(ainput)[k]) { - sameshape = 0; - } - } - } - else { - /* Test to see if amask is 1d */ - if (PyArray_NDIM(amask) != 1) { - sameshape = 0; - } - else if ((PyArray_SIZE(ainput)) != PyArray_SIZE(amask)) { - sameshape = 0; - } - } - if (!sameshape) { - PyErr_SetString(PyExc_TypeError, - "mask array must be 1-d or same shape as input array"); - goto fail; - } - - avals = (PyArrayObject *)PyArray_FromObject(vals, - PyArray_DESCR(ainput)->type_num, 0, 1); - if (avals == NULL) { - goto fail; - } - numvals = PyArray_SIZE(avals); - nd = PyArray_NDIM(ainput); - input_data = PyArray_DATA(ainput); - mptr = PyArray_DATA(amask); - melsize = PyArray_DESCR(amask)->elsize; - vptr = PyArray_DATA(avals); - delsize = PyArray_DESCR(avals)->elsize; - zero = PyArray_Zero(amask); - if (zero == NULL) { - goto fail; - } - objarray = (PyArray_DESCR(ainput)->type_num == NPY_OBJECT); - - /* Handle zero-dimensional case separately */ - if (nd == 0) { - if (memcmp(mptr,zero,melsize) != 0) { - /* Copy value element over to input array */ - memcpy(input_data,vptr,delsize); - if (objarray) { - Py_INCREF(*((PyObject **)vptr)); - } - } - Py_DECREF(amask); - Py_DECREF(avals); - PyDataMem_FREE(zero); - Py_DECREF(ainput); - Py_INCREF(Py_None); - return Py_None; - } - - totmask = (int) PyArray_SIZE(amask); - instrides = PyArray_STRIDES(ainput); - inshape = PyArray_DIMS(ainput); - if (objarray) { - /* object array, need to refcount, can't release the GIL */ - arr_insert_loop(mptr, vptr, input_data, zero, PyArray_DATA(avals), - melsize, delsize, objarray, totmask, numvals, nd, - instrides, inshape); - } - else { - /* No increfs take place in arr_insert_loop, so release the GIL */ - NPY_BEGIN_ALLOW_THREADS; - arr_insert_loop(mptr, vptr, input_data, zero, PyArray_DATA(avals), - melsize, delsize, objarray, totmask, numvals, nd, - instrides, inshape); - NPY_END_ALLOW_THREADS; - } - - Py_DECREF(amask); - Py_DECREF(avals); - PyDataMem_FREE(zero); - Py_DECREF(ainput); - Py_INCREF(Py_None); - return Py_None; - -fail: - PyDataMem_FREE(zero); - Py_XDECREF(ainput); - Py_XDECREF(amask); - Py_XDECREF(avals); - return NULL; -} - -/** @brief Use bisection on a sorted array to find first entry > key. - * - * Use bisection to find an index i s.t. arr[i] <= key < arr[i + 1]. If there is - * no such i the error returns are: - * key < arr[0] -- -1 - * key == arr[len - 1] -- len - 1 - * key > arr[len - 1] -- len - * The array is assumed contiguous and sorted in ascending order. - * - * @param key key value. - * @param arr contiguous sorted array to be searched. - * @param len length of the array. - * @return index - */ -static npy_intp -binary_search(double key, double arr [], npy_intp len) -{ - npy_intp imin = 0; - npy_intp imax = len; - - if (key > arr[len - 1]) { - return len; - } - while (imin < imax) { - npy_intp imid = imin + ((imax - imin) >> 1); - if (key >= arr[imid]) { - imin = imid + 1; - } - else { - imax = imid; - } - } - return imin - 1; -} - -static PyObject * -arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict) -{ - - PyObject *fp, *xp, *x; - PyObject *left = NULL, *right = NULL; - PyArrayObject *afp = NULL, *axp = NULL, *ax = NULL, *af = NULL; - npy_intp i, lenx, lenxp; - double lval, rval; - double *dy, *dx, *dz, *dres, *slopes; - - static char *kwlist[] = {"x", "xp", "fp", "left", "right", NULL}; - - if (!PyArg_ParseTupleAndKeywords(args, kwdict, "OOO|OO", kwlist, - &x, &xp, &fp, &left, &right)) { - return NULL; - } - - afp = (PyArrayObject *)PyArray_ContiguousFromAny(fp, NPY_DOUBLE, 1, 1); - if (afp == NULL) { - return NULL; - } - axp = (PyArrayObject *)PyArray_ContiguousFromAny(xp, NPY_DOUBLE, 1, 1); - if (axp == NULL) { - goto fail; - } - ax = (PyArrayObject *)PyArray_ContiguousFromAny(x, NPY_DOUBLE, 1, 0); - if (ax == NULL) { - goto fail; - } - lenxp = PyArray_DIMS(axp)[0]; - if (lenxp == 0) { - PyErr_SetString(PyExc_ValueError, - "array of sample points is empty"); - goto fail; - } - if (PyArray_DIMS(afp)[0] != lenxp) { - PyErr_SetString(PyExc_ValueError, - "fp and xp are not of the same length."); - goto fail; - } - - af = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(ax), - PyArray_DIMS(ax), NPY_DOUBLE); - if (af == NULL) { - goto fail; - } - lenx = PyArray_SIZE(ax); - - dy = (double *)PyArray_DATA(afp); - dx = (double *)PyArray_DATA(axp); - dz = (double *)PyArray_DATA(ax); - dres = (double *)PyArray_DATA(af); - - /* Get left and right fill values. */ - if ((left == NULL) || (left == Py_None)) { - lval = dy[0]; - } - else { - lval = PyFloat_AsDouble(left); - if ((lval == -1) && PyErr_Occurred()) { - goto fail; - } - } - if ((right == NULL) || (right == Py_None)) { - rval = dy[lenxp-1]; - } - else { - rval = PyFloat_AsDouble(right); - if ((rval == -1) && PyErr_Occurred()) { - goto fail; - } - } - - /* only pre-calculate slopes if there are relatively few of them. */ - if (lenxp <= lenx) { - slopes = (double *) PyArray_malloc((lenxp - 1)*sizeof(double)); - if (! slopes) { - goto fail; - } - NPY_BEGIN_ALLOW_THREADS; - for (i = 0; i < lenxp - 1; i++) { - slopes[i] = (dy[i + 1] - dy[i])/(dx[i + 1] - dx[i]); - } - for (i = 0; i < lenx; i++) { - const double x = dz[i]; - npy_intp j; - - if (npy_isnan(x)) { - dres[i] = x; - continue; - } - - j = binary_search(x, dx, lenxp); - if (j == -1) { - dres[i] = lval; - } - else if (j == lenxp - 1) { - dres[i] = dy[j]; - } - else if (j == lenxp) { - dres[i] = rval; - } - else { - dres[i] = slopes[j]*(x - dx[j]) + dy[j]; - } - } - NPY_END_ALLOW_THREADS; - PyArray_free(slopes); - } - else { - NPY_BEGIN_ALLOW_THREADS; - for (i = 0; i < lenx; i++) { - const double x = dz[i]; - npy_intp j; - - if (npy_isnan(x)) { - dres[i] = x; - continue; - } - - j = binary_search(x, dx, lenxp); - if (j == -1) { - dres[i] = lval; - } - else if (j == lenxp - 1) { - dres[i] = dy[j]; - } - else if (j == lenxp) { - dres[i] = rval; - } - else { - const double slope = (dy[j + 1] - dy[j])/(dx[j + 1] - dx[j]); - dres[i] = slope*(x - dx[j]) + dy[j]; - } - } - NPY_END_ALLOW_THREADS; - } - - Py_DECREF(afp); - Py_DECREF(axp); - Py_DECREF(ax); - return (PyObject *)af; - -fail: - Py_XDECREF(afp); - Py_XDECREF(axp); - Py_XDECREF(ax); - Py_XDECREF(af); - return NULL; -} - -/* - * Converts a Python sequence into 'count' PyArrayObjects - * - * seq - Input Python object, usually a tuple but any sequence works. - * op - Where the arrays are placed. - * count - How many arrays there should be (errors if it doesn't match). - * paramname - The name of the parameter that produced 'seq'. - */ -static int sequence_to_arrays(PyObject *seq, - PyArrayObject **op, int count, - char *paramname) -{ - int i; - - if (!PySequence_Check(seq) || PySequence_Size(seq) != count) { - PyErr_Format(PyExc_ValueError, - "parameter %s must be a sequence of length %d", - paramname, count); - return -1; - } - - for (i = 0; i < count; ++i) { - PyObject *item = PySequence_GetItem(seq, i); - if (item == NULL) { - while (--i >= 0) { - Py_DECREF(op[i]); - op[i] = NULL; - } - return -1; - } - - op[i] = (PyArrayObject *)PyArray_FromAny(item, NULL, 0, 0, 0, NULL); - if (op[i] == NULL) { - while (--i >= 0) { - Py_DECREF(op[i]); - op[i] = NULL; - } - Py_DECREF(item); - return -1; - } - - Py_DECREF(item); - } - - return 0; -} - -/* Inner loop for unravel_index */ -static int -ravel_multi_index_loop(int ravel_ndim, npy_intp *ravel_dims, - npy_intp *ravel_strides, - npy_intp count, - NPY_CLIPMODE *modes, - char **coords, npy_intp *coords_strides) -{ - int i; - char invalid; - npy_intp j, m; - - NPY_BEGIN_ALLOW_THREADS; - invalid = 0; - while (count--) { - npy_intp raveled = 0; - for (i = 0; i < ravel_ndim; ++i) { - m = ravel_dims[i]; - j = *(npy_intp *)coords[i]; - switch (modes[i]) { - case NPY_RAISE: - if (j < 0 || j >= m) { - invalid = 1; - goto end_while; - } - break; - case NPY_WRAP: - if (j < 0) { - j += m; - if (j < 0) { - j = j % m; - if (j != 0) { - j += m; - } - } - } - else if (j >= m) { - j -= m; - if (j >= m) { - j = j % m; - } - } - break; - case NPY_CLIP: - if (j < 0) { - j = 0; - } - else if (j >= m) { - j = m - 1; - } - break; - - } - raveled += j * ravel_strides[i]; - - coords[i] += coords_strides[i]; - } - *(npy_intp *)coords[ravel_ndim] = raveled; - coords[ravel_ndim] += coords_strides[ravel_ndim]; - } -end_while: - NPY_END_ALLOW_THREADS; - if (invalid) { - PyErr_SetString(PyExc_ValueError, - "invalid entry in coordinates array"); - return NPY_FAIL; - } - return NPY_SUCCEED; -} - -/* ravel_multi_index implementation - see add_newdocs.py */ -static PyObject * -arr_ravel_multi_index(PyObject *self, PyObject *args, PyObject *kwds) -{ - int i, s; - PyObject *mode0=NULL, *coords0=NULL; - PyArrayObject *ret = NULL; - PyArray_Dims dimensions={0,0}; - npy_intp ravel_strides[NPY_MAXDIMS]; - NPY_ORDER order = NPY_CORDER; - NPY_CLIPMODE modes[NPY_MAXDIMS]; - - PyArrayObject *op[NPY_MAXARGS]; - PyArray_Descr *dtype[NPY_MAXARGS]; - npy_uint32 op_flags[NPY_MAXARGS]; - - NpyIter *iter = NULL; - - char *kwlist[] = {"multi_index", "dims", "mode", "order", NULL}; - - memset(op, 0, sizeof(op)); - dtype[0] = NULL; - - if (!PyArg_ParseTupleAndKeywords(args, kwds, - "OO&|OO&:ravel_multi_index", kwlist, - &coords0, - PyArray_IntpConverter, &dimensions, - &mode0, - PyArray_OrderConverter, &order)) { - goto fail; - } - - if (dimensions.len+1 > NPY_MAXARGS) { - PyErr_SetString(PyExc_ValueError, - "too many dimensions passed to ravel_multi_index"); - goto fail; - } - - if (!PyArray_ConvertClipmodeSequence(mode0, modes, dimensions.len)) { - goto fail; - } - - switch (order) { - case NPY_CORDER: - s = 1; - for (i = dimensions.len-1; i >= 0; --i) { - ravel_strides[i] = s; - s *= dimensions.ptr[i]; - } - break; - case NPY_FORTRANORDER: - s = 1; - for (i = 0; i < dimensions.len; ++i) { - ravel_strides[i] = s; - s *= dimensions.ptr[i]; - } - break; - default: - PyErr_SetString(PyExc_ValueError, - "only 'C' or 'F' order is permitted"); - goto fail; - } - - /* Get the multi_index into op */ - if (sequence_to_arrays(coords0, op, dimensions.len, "multi_index") < 0) { - goto fail; - } - - - for (i = 0; i < dimensions.len; ++i) { - op_flags[i] = NPY_ITER_READONLY| - NPY_ITER_ALIGNED; - } - op_flags[dimensions.len] = NPY_ITER_WRITEONLY| - NPY_ITER_ALIGNED| - NPY_ITER_ALLOCATE; - dtype[0] = PyArray_DescrFromType(NPY_INTP); - for (i = 1; i <= dimensions.len; ++i) { - dtype[i] = dtype[0]; - } - - iter = NpyIter_MultiNew(dimensions.len+1, op, NPY_ITER_BUFFERED| - NPY_ITER_EXTERNAL_LOOP| - NPY_ITER_ZEROSIZE_OK, - NPY_KEEPORDER, - NPY_SAME_KIND_CASTING, - op_flags, dtype); - if (iter == NULL) { - goto fail; - } - - if (NpyIter_GetIterSize(iter) != 0) { - NpyIter_IterNextFunc *iternext; - char **dataptr; - npy_intp *strides; - npy_intp *countptr; - - iternext = NpyIter_GetIterNext(iter, NULL); - if (iternext == NULL) { - goto fail; - } - dataptr = NpyIter_GetDataPtrArray(iter); - strides = NpyIter_GetInnerStrideArray(iter); - countptr = NpyIter_GetInnerLoopSizePtr(iter); - - do { - if (ravel_multi_index_loop(dimensions.len, dimensions.ptr, - ravel_strides, *countptr, modes, - dataptr, strides) != NPY_SUCCEED) { - goto fail; - } - } while(iternext(iter)); - } - - ret = NpyIter_GetOperandArray(iter)[dimensions.len]; - Py_INCREF(ret); - - Py_DECREF(dtype[0]); - for (i = 0; i < dimensions.len; ++i) { - Py_XDECREF(op[i]); - } - PyDimMem_FREE(dimensions.ptr); - NpyIter_Deallocate(iter); - return PyArray_Return(ret); - -fail: - Py_XDECREF(dtype[0]); - for (i = 0; i < dimensions.len; ++i) { - Py_XDECREF(op[i]); - } - PyDimMem_FREE(dimensions.ptr); - NpyIter_Deallocate(iter); - return NULL; -} - -/* C-order inner loop for unravel_index */ -static int -unravel_index_loop_corder(int unravel_ndim, npy_intp *unravel_dims, - npy_intp unravel_size, npy_intp count, - char *indices, npy_intp indices_stride, - npy_intp *coords) -{ - int i; - char invalid; - npy_intp val; - - NPY_BEGIN_ALLOW_THREADS; - invalid = 0; - while (count--) { - val = *(npy_intp *)indices; - if (val < 0 || val >= unravel_size) { - invalid = 1; - break; - } - for (i = unravel_ndim-1; i >= 0; --i) { - coords[i] = val % unravel_dims[i]; - val /= unravel_dims[i]; - } - coords += unravel_ndim; - indices += indices_stride; - } - NPY_END_ALLOW_THREADS; - if (invalid) { - PyErr_SetString(PyExc_ValueError, - "invalid entry in index array"); - return NPY_FAIL; - } - return NPY_SUCCEED; -} - -/* Fortran-order inner loop for unravel_index */ -static int -unravel_index_loop_forder(int unravel_ndim, npy_intp *unravel_dims, - npy_intp unravel_size, npy_intp count, - char *indices, npy_intp indices_stride, - npy_intp *coords) -{ - int i; - char invalid; - npy_intp val; - - NPY_BEGIN_ALLOW_THREADS; - invalid = 0; - while (count--) { - val = *(npy_intp *)indices; - if (val < 0 || val >= unravel_size) { - invalid = 1; - break; - } - for (i = 0; i < unravel_ndim; ++i) { - *coords++ = val % unravel_dims[i]; - val /= unravel_dims[i]; - } - indices += indices_stride; - } - NPY_END_ALLOW_THREADS; - if (invalid) { - PyErr_SetString(PyExc_ValueError, - "invalid entry in index array"); - return NPY_FAIL; - } - return NPY_SUCCEED; -} - -/* unravel_index implementation - see add_newdocs.py */ -static PyObject * -arr_unravel_index(PyObject *self, PyObject *args, PyObject *kwds) -{ - PyObject *indices0 = NULL, *ret_tuple = NULL; - PyArrayObject *ret_arr = NULL; - PyArrayObject *indices = NULL; - PyArray_Descr *dtype = NULL; - PyArray_Dims dimensions={0,0}; - NPY_ORDER order = NPY_CORDER; - npy_intp unravel_size; - - NpyIter *iter = NULL; - int i, ret_ndim; - npy_intp ret_dims[NPY_MAXDIMS], ret_strides[NPY_MAXDIMS]; - - char *kwlist[] = {"indices", "dims", "order", NULL}; - - if (!PyArg_ParseTupleAndKeywords(args, kwds, "OO&|O&:unravel_index", - kwlist, - &indices0, - PyArray_IntpConverter, &dimensions, - PyArray_OrderConverter, &order)) { - goto fail; - } - - if (dimensions.len == 0) { - PyErr_SetString(PyExc_ValueError, - "dims must have at least one value"); - goto fail; - } - - unravel_size = PyArray_MultiplyList(dimensions.ptr, dimensions.len); - - if (!PyArray_Check(indices0)) { - indices = (PyArrayObject*)PyArray_FromAny(indices0, - NULL, 0, 0, 0, NULL); - if (indices == NULL) { - goto fail; - } - } - else { - indices = (PyArrayObject *)indices0; - Py_INCREF(indices); - } - - dtype = PyArray_DescrFromType(NPY_INTP); - if (dtype == NULL) { - goto fail; - } - - iter = NpyIter_New(indices, NPY_ITER_READONLY| - NPY_ITER_ALIGNED| - NPY_ITER_BUFFERED| - NPY_ITER_ZEROSIZE_OK| - NPY_ITER_DONT_NEGATE_STRIDES| - NPY_ITER_MULTI_INDEX, - NPY_KEEPORDER, NPY_SAME_KIND_CASTING, - dtype); - if (iter == NULL) { - goto fail; - } - - /* - * Create the return array with a layout compatible with the indices - * and with a dimension added to the end for the multi-index - */ - ret_ndim = PyArray_NDIM(indices) + 1; - if (NpyIter_GetShape(iter, ret_dims) != NPY_SUCCEED) { - goto fail; - } - ret_dims[ret_ndim-1] = dimensions.len; - if (NpyIter_CreateCompatibleStrides(iter, - dimensions.len*sizeof(npy_intp), ret_strides) != NPY_SUCCEED) { - goto fail; - } - ret_strides[ret_ndim-1] = sizeof(npy_intp); - - /* Remove the multi-index and inner loop */ - if (NpyIter_RemoveMultiIndex(iter) != NPY_SUCCEED) { - goto fail; - } - if (NpyIter_EnableExternalLoop(iter) != NPY_SUCCEED) { - goto fail; - } - - ret_arr = (PyArrayObject *)PyArray_NewFromDescr(&PyArray_Type, dtype, - ret_ndim, ret_dims, ret_strides, NULL, 0, NULL); - dtype = NULL; - if (ret_arr == NULL) { - goto fail; - } - - if (order == NPY_CORDER) { - if (NpyIter_GetIterSize(iter) != 0) { - NpyIter_IterNextFunc *iternext; - char **dataptr; - npy_intp *strides; - npy_intp *countptr, count; - npy_intp *coordsptr = (npy_intp *)PyArray_DATA(ret_arr); - - iternext = NpyIter_GetIterNext(iter, NULL); - if (iternext == NULL) { - goto fail; - } - dataptr = NpyIter_GetDataPtrArray(iter); - strides = NpyIter_GetInnerStrideArray(iter); - countptr = NpyIter_GetInnerLoopSizePtr(iter); - - do { - count = *countptr; - if (unravel_index_loop_corder(dimensions.len, dimensions.ptr, - unravel_size, count, *dataptr, *strides, - coordsptr) != NPY_SUCCEED) { - goto fail; - } - coordsptr += count*dimensions.len; - } while(iternext(iter)); - } - } - else if (order == NPY_FORTRANORDER) { - if (NpyIter_GetIterSize(iter) != 0) { - NpyIter_IterNextFunc *iternext; - char **dataptr; - npy_intp *strides; - npy_intp *countptr, count; - npy_intp *coordsptr = (npy_intp *)PyArray_DATA(ret_arr); - - iternext = NpyIter_GetIterNext(iter, NULL); - if (iternext == NULL) { - goto fail; - } - dataptr = NpyIter_GetDataPtrArray(iter); - strides = NpyIter_GetInnerStrideArray(iter); - countptr = NpyIter_GetInnerLoopSizePtr(iter); - - do { - count = *countptr; - if (unravel_index_loop_forder(dimensions.len, dimensions.ptr, - unravel_size, count, *dataptr, *strides, - coordsptr) != NPY_SUCCEED) { - goto fail; - } - coordsptr += count*dimensions.len; - } while(iternext(iter)); - } - } - else { - PyErr_SetString(PyExc_ValueError, - "only 'C' or 'F' order is permitted"); - goto fail; - } - - /* Now make a tuple of views, one per index */ - ret_tuple = PyTuple_New(dimensions.len); - if (ret_tuple == NULL) { - goto fail; - } - for (i = 0; i < dimensions.len; ++i) { - PyArrayObject *view; - - view = (PyArrayObject *)PyArray_New(&PyArray_Type, ret_ndim-1, - ret_dims, NPY_INTP, - ret_strides, - PyArray_BYTES(ret_arr) + i*sizeof(npy_intp), - 0, 0, NULL); - if (view == NULL) { - goto fail; - } - Py_INCREF(ret_arr); - if (PyArray_SetBaseObject(view, (PyObject *)ret_arr) < 0) { - Py_DECREF(view); - goto fail; - } - PyTuple_SET_ITEM(ret_tuple, i, PyArray_Return(view)); - } - - Py_DECREF(ret_arr); - Py_XDECREF(indices); - PyDimMem_FREE(dimensions.ptr); - NpyIter_Deallocate(iter); - - return ret_tuple; - -fail: - Py_XDECREF(ret_tuple); - Py_XDECREF(ret_arr); - Py_XDECREF(dtype); - Py_XDECREF(indices); - PyDimMem_FREE(dimensions.ptr); - NpyIter_Deallocate(iter); - return NULL; -} - - -static PyTypeObject *PyMemberDescr_TypePtr = NULL; -static PyTypeObject *PyGetSetDescr_TypePtr = NULL; -static PyTypeObject *PyMethodDescr_TypePtr = NULL; - -/* Can only be called if doc is currently NULL */ -static PyObject * -arr_add_docstring(PyObject *NPY_UNUSED(dummy), PyObject *args) -{ - PyObject *obj; - PyObject *str; - char *docstr; - static char *msg = "already has a docstring"; - - /* Don't add docstrings */ - if (Py_OptimizeFlag > 1) { - Py_INCREF(Py_None); - return Py_None; - } -#if defined(NPY_PY3K) - if (!PyArg_ParseTuple(args, "OO!", &obj, &PyUnicode_Type, &str)) { - return NULL; - } - - docstr = PyBytes_AS_STRING(PyUnicode_AsUTF8String(str)); -#else - if (!PyArg_ParseTuple(args, "OO!", &obj, &PyString_Type, &str)) { - return NULL; - } - - docstr = PyString_AS_STRING(str); -#endif - -#define _TESTDOC1(typebase) (Py_TYPE(obj) == &Py##typebase##_Type) -#define _TESTDOC2(typebase) (Py_TYPE(obj) == Py##typebase##_TypePtr) -#define _ADDDOC(typebase, doc, name) do { \ - Py##typebase##Object *new = (Py##typebase##Object *)obj; \ - if (!(doc)) { \ - doc = docstr; \ - } \ - else { \ - PyErr_Format(PyExc_RuntimeError, "%s method %s", name, msg); \ - return NULL; \ - } \ - } while (0) - - if (_TESTDOC1(CFunction)) { - _ADDDOC(CFunction, new->m_ml->ml_doc, new->m_ml->ml_name); - } - else if (_TESTDOC1(Type)) { - _ADDDOC(Type, new->tp_doc, new->tp_name); - } - else if (_TESTDOC2(MemberDescr)) { - _ADDDOC(MemberDescr, new->d_member->doc, new->d_member->name); - } - else if (_TESTDOC2(GetSetDescr)) { - _ADDDOC(GetSetDescr, new->d_getset->doc, new->d_getset->name); - } - else if (_TESTDOC2(MethodDescr)) { - _ADDDOC(MethodDescr, new->d_method->ml_doc, new->d_method->ml_name); - } - else { - PyObject *doc_attr; - - doc_attr = PyObject_GetAttrString(obj, "__doc__"); - if (doc_attr != NULL && doc_attr != Py_None) { - PyErr_Format(PyExc_RuntimeError, "object %s", msg); - return NULL; - } - Py_XDECREF(doc_attr); - - if (PyObject_SetAttrString(obj, "__doc__", str) < 0) { - PyErr_SetString(PyExc_TypeError, - "Cannot set a docstring for that object"); - return NULL; - } - Py_INCREF(Py_None); - return Py_None; - } - -#undef _TESTDOC1 -#undef _TESTDOC2 -#undef _ADDDOC - - Py_INCREF(str); - Py_INCREF(Py_None); - return Py_None; -} - - -/* docstring in numpy.add_newdocs.py */ -static PyObject * -add_newdoc_ufunc(PyObject *NPY_UNUSED(dummy), PyObject *args) -{ - PyUFuncObject *ufunc; - PyObject *str; - char *docstr, *newdocstr; - -#if defined(NPY_PY3K) - if (!PyArg_ParseTuple(args, "O!O!", &PyUFunc_Type, &ufunc, - &PyUnicode_Type, &str)) { - return NULL; - } - docstr = PyBytes_AS_STRING(PyUnicode_AsUTF8String(str)); -#else - if (!PyArg_ParseTuple(args, "O!O!", &PyUFunc_Type, &ufunc, - &PyString_Type, &str)) { - return NULL; - } - docstr = PyString_AS_STRING(str); -#endif - - if (NULL != ufunc->doc) { - PyErr_SetString(PyExc_ValueError, - "Cannot change docstring of ufunc with non-NULL docstring"); - return NULL; - } - - /* - * This introduces a memory leak, as the memory allocated for the doc - * will not be freed even if the ufunc itself is deleted. In practice - * this should not be a problem since the user would have to - * repeatedly create, document, and throw away ufuncs. - */ - newdocstr = malloc(strlen(docstr) + 1); - strcpy(newdocstr, docstr); - ufunc->doc = newdocstr; - - Py_INCREF(Py_None); - return Py_None; -} - -/* - * This function packs boolean values in the input array into the bits of a - * byte array. Truth values are determined as usual: 0 is false, everything - * else is true. - */ -static NPY_INLINE void -pack_inner(const char *inptr, - npy_intp element_size, /* in bytes */ - npy_intp n_in, - npy_intp in_stride, - char *outptr, - npy_intp n_out, - npy_intp out_stride) -{ - /* - * Loop through the elements of inptr. - * Determine whether or not it is nonzero. - * Yes: set corresponding bit (and adjust build value) - * No: move on - * Every 8th value, set the value of build and increment the outptr - */ - npy_intp index; - int remain = n_in % 8; /* uneven bits */ - - if (remain == 0) { /* assumes n_in > 0 */ - remain = 8; - } - for (index = 0; index < n_out; index++) { - char build = 0; - int i, maxi; - npy_intp j; - - maxi = (index == n_out - 1) ? remain : 8; - for (i = 0; i < maxi; i++) { - build <<= 1; - for (j = 0; j < element_size; j++) { - build |= (inptr[j] != 0); - } - inptr += in_stride; - } - if (index == n_out - 1) { - build <<= 8 - remain; - } - *outptr = build; - outptr += out_stride; - } -} - -static PyObject * -pack_bits(PyObject *input, int axis) -{ - PyArrayObject *inp; - PyArrayObject *new = NULL; - PyArrayObject *out = NULL; - npy_intp outdims[NPY_MAXDIMS]; - int i; - PyArrayIterObject *it, *ot; - NPY_BEGIN_THREADS_DEF; - - inp = (PyArrayObject *)PyArray_FROM_O(input); - - if (inp == NULL) { - return NULL; - } - if (!PyArray_ISBOOL(inp) && !PyArray_ISINTEGER(inp)) { - PyErr_SetString(PyExc_TypeError, - "Expected an input array of integer or boolean data type"); - goto fail; - } - - new = (PyArrayObject *)PyArray_CheckAxis(inp, &axis, 0); - Py_DECREF(inp); - if (new == NULL) { - return NULL; - } - /* Handle empty array separately */ - if (PyArray_SIZE(new) == 0) { - return PyArray_Copy(new); - } - - if (PyArray_NDIM(new) == 0) { - char *optr, *iptr; - - out = (PyArrayObject *)PyArray_New(Py_TYPE(new), 0, NULL, NPY_UBYTE, - NULL, NULL, 0, 0, NULL); - if (out == NULL) { - goto fail; - } - optr = PyArray_DATA(out); - iptr = PyArray_DATA(new); - *optr = 0; - for (i = 0; i < PyArray_ITEMSIZE(new); i++) { - if (*iptr != 0) { - *optr = 1; - break; - } - iptr++; - } - goto finish; - } - - - /* Setup output shape */ - for (i = 0; i < PyArray_NDIM(new); i++) { - outdims[i] = PyArray_DIM(new, i); - } - - /* - * Divide axis dimension by 8 - * 8 -> 1, 9 -> 2, 16 -> 2, 17 -> 3 etc.. - */ - outdims[axis] = ((outdims[axis] - 1) >> 3) + 1; - - /* Create output array */ - out = (PyArrayObject *)PyArray_New(Py_TYPE(new), - PyArray_NDIM(new), outdims, NPY_UBYTE, - NULL, NULL, 0, PyArray_ISFORTRAN(new), NULL); - if (out == NULL) { - goto fail; - } - /* Setup iterators to iterate over all but given axis */ - it = (PyArrayIterObject *)PyArray_IterAllButAxis((PyObject *)new, &axis); - ot = (PyArrayIterObject *)PyArray_IterAllButAxis((PyObject *)out, &axis); - if (it == NULL || ot == NULL) { - Py_XDECREF(it); - Py_XDECREF(ot); - goto fail; - } - - NPY_BEGIN_THREADS_THRESHOLDED(PyArray_DIM(out, axis)); - while (PyArray_ITER_NOTDONE(it)) { - pack_inner(PyArray_ITER_DATA(it), PyArray_ITEMSIZE(new), - PyArray_DIM(new, axis), PyArray_STRIDE(new, axis), - PyArray_ITER_DATA(ot), PyArray_DIM(out, axis), - PyArray_STRIDE(out, axis)); - PyArray_ITER_NEXT(it); - PyArray_ITER_NEXT(ot); - } - NPY_END_THREADS; - - Py_DECREF(it); - Py_DECREF(ot); - -finish: - Py_DECREF(new); - return (PyObject *)out; - -fail: - Py_XDECREF(new); - Py_XDECREF(out); - return NULL; -} - -static PyObject * -unpack_bits(PyObject *input, int axis) -{ - PyArrayObject *inp; - PyArrayObject *new = NULL; - PyArrayObject *out = NULL; - npy_intp outdims[NPY_MAXDIMS]; - int i; - PyArrayIterObject *it, *ot; - npy_intp n_in, in_stride, out_stride; - NPY_BEGIN_THREADS_DEF; - - inp = (PyArrayObject *)PyArray_FROM_O(input); - - if (inp == NULL) { - return NULL; - } - if (PyArray_TYPE(inp) != NPY_UBYTE) { - PyErr_SetString(PyExc_TypeError, - "Expected an input array of unsigned byte data type"); - goto fail; - } - - new = (PyArrayObject *)PyArray_CheckAxis(inp, &axis, 0); - Py_DECREF(inp); - if (new == NULL) { - return NULL; - } - /* Handle zero-dim array separately */ - if (PyArray_SIZE(new) == 0) { - return PyArray_Copy(new); - } - - if (PyArray_NDIM(new) == 0) { - /* Handle 0-d array by converting it to a 1-d array */ - PyArrayObject *temp; - PyArray_Dims newdim = {NULL, 1}; - npy_intp shape = 1; - - newdim.ptr = &shape; - temp = (PyArrayObject *)PyArray_Newshape(new, &newdim, NPY_CORDER); - if (temp == NULL) { - goto fail; - } - Py_DECREF(new); - new = temp; - } - - /* Setup output shape */ - for (i=0; i>= 1; - } - inptr += in_stride; - } - PyArray_ITER_NEXT(it); - PyArray_ITER_NEXT(ot); - } - NPY_END_THREADS; - - Py_DECREF(it); - Py_DECREF(ot); - - Py_DECREF(new); - return (PyObject *)out; - -fail: - Py_XDECREF(new); - Py_XDECREF(out); - return NULL; -} - - -static PyObject * -io_pack(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds) -{ - PyObject *obj; - int axis = NPY_MAXDIMS; - static char *kwlist[] = {"in", "axis", NULL}; - - if (!PyArg_ParseTupleAndKeywords( args, kwds, "O|O&" , kwlist, - &obj, PyArray_AxisConverter, &axis)) { - return NULL; - } - return pack_bits(obj, axis); -} - -static PyObject * -io_unpack(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds) -{ - PyObject *obj; - int axis = NPY_MAXDIMS; - static char *kwlist[] = {"in", "axis", NULL}; - - if (!PyArg_ParseTupleAndKeywords( args, kwds, "O|O&" , kwlist, - &obj, PyArray_AxisConverter, &axis)) { - return NULL; - } - return unpack_bits(obj, axis); -} - -/* The docstrings for many of these methods are in add_newdocs.py. */ -static struct PyMethodDef methods[] = { - {"_insert", (PyCFunction)arr_insert, - METH_VARARGS | METH_KEYWORDS, arr_insert__doc__}, - {"bincount", (PyCFunction)arr_bincount, - METH_VARARGS | METH_KEYWORDS, NULL}, - {"digitize", (PyCFunction)arr_digitize, - METH_VARARGS | METH_KEYWORDS, NULL}, - {"interp", (PyCFunction)arr_interp, - METH_VARARGS | METH_KEYWORDS, NULL}, - {"ravel_multi_index", (PyCFunction)arr_ravel_multi_index, - METH_VARARGS | METH_KEYWORDS, NULL}, - {"unravel_index", (PyCFunction)arr_unravel_index, - METH_VARARGS | METH_KEYWORDS, NULL}, - {"add_docstring", (PyCFunction)arr_add_docstring, - METH_VARARGS, NULL}, - {"add_newdoc_ufunc", (PyCFunction)add_newdoc_ufunc, - METH_VARARGS, NULL}, - {"packbits", (PyCFunction)io_pack, - METH_VARARGS | METH_KEYWORDS, NULL}, - {"unpackbits", (PyCFunction)io_unpack, - METH_VARARGS | METH_KEYWORDS, NULL}, - {NULL, NULL, 0, NULL} /* sentinel */ -}; - -static void -define_types(void) -{ - PyObject *tp_dict; - PyObject *myobj; - - tp_dict = PyArrayDescr_Type.tp_dict; - /* Get "subdescr" */ - myobj = PyDict_GetItemString(tp_dict, "fields"); - if (myobj == NULL) { - return; - } - PyGetSetDescr_TypePtr = Py_TYPE(myobj); - myobj = PyDict_GetItemString(tp_dict, "alignment"); - if (myobj == NULL) { - return; - } - PyMemberDescr_TypePtr = Py_TYPE(myobj); - myobj = PyDict_GetItemString(tp_dict, "newbyteorder"); - if (myobj == NULL) { - return; - } - PyMethodDescr_TypePtr = Py_TYPE(myobj); - return; -} - -#if defined(NPY_PY3K) -static struct PyModuleDef moduledef = { - PyModuleDef_HEAD_INIT, - "_compiled_base", - NULL, - -1, - methods, - NULL, - NULL, - NULL, - NULL -}; -#endif - -#if defined(NPY_PY3K) -#define RETVAL m -PyMODINIT_FUNC PyInit__compiled_base(void) -#else -#define RETVAL -PyMODINIT_FUNC -init_compiled_base(void) -#endif -{ - PyObject *m, *d; - -#if defined(NPY_PY3K) - m = PyModule_Create(&moduledef); -#else - m = Py_InitModule("_compiled_base", methods); -#endif - if (!m) { - return RETVAL; - } - - /* Import the array objects */ - import_array(); - import_umath(); - - /* Add some symbolic constants to the module */ - d = PyModule_GetDict(m); - - /* - * PyExc_Exception should catch all the standard errors that are - * now raised instead of the string exception "numpy.lib.error". - * This is for backward compatibility with existing code. - */ - PyDict_SetItemString(d, "error", PyExc_Exception); - - - /* define PyGetSetDescr_Type and PyMemberDescr_Type */ - define_types(); - - return RETVAL; -} -- cgit v1.2.1 From 90199336e34891448ce4052f322228906a72c774 Mon Sep 17 00:00:00 2001 From: Julian Taylor Date: Wed, 21 Jan 2015 20:47:46 +0100 Subject: MAINT: update .gitignore for new template files --- .gitignore | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index 3952be1fe..8387f4a03 100644 --- a/.gitignore +++ b/.gitignore @@ -126,8 +126,8 @@ numpy/core/src/npysort/selection.c numpy/core/src/npysort/sort.c numpy/core/src/private/npy_binsearch.h numpy/core/src/private/npy_partition.h -numpy/core/src/private/scalarmathmodule.h -numpy/core/src/scalarmathmodule.c +numpy/core/src/private/templ_common.h +numpy/core/src/umath/scalarmath.c numpy/core/src/umath/funcs.inc numpy/core/src/umath/loops.[ch] numpy/core/src/umath/operand_flag_tests.c -- cgit v1.2.1