diff options
-rw-r--r-- | doc/release/2.0.0-notes.rst | 22 | ||||
-rw-r--r-- | numpy/add_newdocs.py | 2 | ||||
-rw-r--r-- | numpy/core/fromnumeric.py | 2 | ||||
-rw-r--r-- | numpy/core/src/multiarray/item_selection.c | 229 | ||||
-rw-r--r-- | numpy/core/src/multiarray/methods.c | 5 | ||||
-rw-r--r-- | numpy/core/tests/test_maskna.py | 30 | ||||
-rw-r--r-- | numpy/lib/twodim_base.py | 15 |
7 files changed, 176 insertions, 129 deletions
diff --git a/doc/release/2.0.0-notes.rst b/doc/release/2.0.0-notes.rst index 195a24d01..d5aa99639 100644 --- a/doc/release/2.0.0-notes.rst +++ b/doc/release/2.0.0-notes.rst @@ -11,11 +11,16 @@ Highlights New features ============ + +Mask-based NA missing values +---------------------------- + Support for NA missing values similar to those in R has been implemented. This was done by adding optional NA masks to the core array object. -Support for NA is not complete, here is a list of the things that do and -do not work with NA values: +While a significant amount of the NumPy functionality has been extended to +support NA masks, not everything is yet supported. Here is a list of things +that do and do not work with NA values: What works with NA: * Basic indexing and slicing, as well as full boolean mask indexing. @@ -23,9 +28,7 @@ What works with NA: * UFunc.reduce methods, with a new skipna parameter. * Array methods: + ndarray.clip, ndarray.min, ndarray.max, ndarray.sum, ndarray.prod, - - - ndarray.conjugate + ndarray.conjugate, ndarray.diagonal What doesn't work with NA: * Fancy indexing, such as with lists and partial boolean masks. @@ -46,6 +49,15 @@ Custom formatter for printing arrays Changes ======= +The default casting rule for UFunc out= parameters has been changed from +'unsafe' to 'same_kind'. Most usages which violate the 'same_kind' +rule are likely bugs, so this change may expose previously undetected +errors in projects that depend on NumPy. + +The functions np.diag, np.diagonal, and <ndarray>.diagonal now return a +view into the original array instead of making a copy. This makes these +functions more consistent with NumPy's general approach of taking views +where possible, and performs much faster as well. Deprecations diff --git a/numpy/add_newdocs.py b/numpy/add_newdocs.py index 95379d611..d992a7122 100644 --- a/numpy/add_newdocs.py +++ b/numpy/add_newdocs.py @@ -3252,7 +3252,7 @@ add_newdoc('numpy.core.multiarray', 'ndarray', ('diagonal', Return specified diagonals. - Refer to `numpy.diagonal` for full documentation. + Refer to :func:`numpy.diagonal` for full documentation. See Also -------- diff --git a/numpy/core/fromnumeric.py b/numpy/core/fromnumeric.py index 8825dc6b7..2364fbfe8 100644 --- a/numpy/core/fromnumeric.py +++ b/numpy/core/fromnumeric.py @@ -911,6 +911,8 @@ def diagonal(a, offset=0, axis1=0, axis2=1): removing `axis1` and `axis2` and appending an index to the right equal to the size of the resulting diagonals. + As of NumPy 1.7, this function always returns a view into `a`. + Parameters ---------- a : array_like diff --git a/numpy/core/src/multiarray/item_selection.c b/numpy/core/src/multiarray/item_selection.c index b1d5350df..827b8e4d4 100644 --- a/numpy/core/src/multiarray/item_selection.c +++ b/numpy/core/src/multiarray/item_selection.c @@ -1657,146 +1657,153 @@ PyArray_SearchSorted(PyArrayObject *op1, PyObject *op2, NPY_SEARCHSIDE side) /*NUMPY_API * Diagonal + * + * As of NumPy 1.7, this function always returns a view into 'self'. */ NPY_NO_EXPORT PyObject * PyArray_Diagonal(PyArrayObject *self, int offset, int axis1, int axis2) { - int n = PyArray_NDIM(self); - PyObject *new; - PyArray_Dims newaxes; - intp dims[MAX_DIMS]; - int i, pos; - - newaxes.ptr = dims; - if (n < 2) { + int i, idim, ndim = PyArray_NDIM(self); + npy_intp *strides, *maskna_strides = NULL; + npy_intp stride1, stride2, maskna_stride1 = 0, maskna_stride2 = 0; + npy_intp *shape, dim1, dim2; + int self_has_maskna = PyArray_HASMASKNA(self); + + char *data, *maskna_data; + npy_intp diag_size; + PyArrayObject *ret; + PyArray_Descr *dtype; + npy_intp ret_shape[NPY_MAXDIMS], ret_strides[NPY_MAXDIMS]; + + if (ndim < 2) { PyErr_SetString(PyExc_ValueError, - "array.ndim must be >= 2"); + "diag requires an array of at least two dimensions"); return NULL; } + + /* Handle negative axes with standard Python indexing rules */ if (axis1 < 0) { - axis1 += n; + axis1 += ndim; } if (axis2 < 0) { - axis2 += n; - } - if ((axis1 == axis2) || (axis1 < 0) || (axis1 >= n) || - (axis2 < 0) || (axis2 >= n)) { - PyErr_Format(PyExc_ValueError, "axis1(=%d) and axis2(=%d) "\ - "must be different and within range (nd=%d)", - axis1, axis2, n); - return NULL; + axis2 += ndim; } - newaxes.len = n; - /* insert at the end */ - newaxes.ptr[n-2] = axis1; - newaxes.ptr[n-1] = axis2; - pos = 0; - for (i = 0; i < n; i++) { - if ((i==axis1) || (i==axis2)) { - continue; - } - newaxes.ptr[pos++] = i; + /* Error check the two axes */ + if (axis1 == axis2) { + PyErr_SetString(PyExc_ValueError, + "axis1 and axis2 cannot be the same"); + return NULL; } - new = PyArray_Transpose(self, &newaxes); - if (new == NULL) { + else if (axis1 < 0 || axis1 >= ndim || axis2 < 0 || axis2 >= ndim) { + PyErr_Format(PyExc_ValueError, + "axis1(=%d) and axis2(=%d) " + "must be within range (ndim=%d)", + axis1, axis2, ndim); return NULL; } - self = (PyArrayObject *)new; - if (n == 2) { - PyObject *a = NULL, *ret = NULL; - PyArrayObject *indices = NULL; - intp n1, n2, start, stop, step, count; - intp *dptr; + /* Get the shape and strides of the two axes */ + shape = PyArray_SHAPE(self); + dim1 = shape[axis1]; + dim2 = shape[axis2]; + strides = PyArray_STRIDES(self); + stride1 = strides[axis1]; + stride2 = strides[axis2]; + if (self_has_maskna) { + maskna_strides = PyArray_MASKNA_STRIDES(self); + maskna_stride1 = maskna_strides[axis1]; + maskna_stride2 = maskna_strides[axis2]; + } - n1 = PyArray_DIMS(self)[0]; - n2 = PyArray_DIMS(self)[1]; - step = n2 + 1; - if (offset < 0) { - start = -n2 * offset; - stop = MIN(n2, n1+offset)*(n2+1) - n2*offset; + /* Compute the data pointers and diag_size for the view */ + data = PyArray_DATA(self); + maskna_data = PyArray_MASKNA_DATA(self); + if (offset > 0) { + if (offset >= dim2) { + diag_size = 0; } else { - start = offset; - stop = MIN(n1, n2-offset)*(n2+1) + offset; - } + data += offset * stride2; + maskna_data += offset * maskna_stride2; - /* count = ceil((stop-start)/step) */ - count = ((stop-start) / step) + (((stop-start) % step) != 0); - indices = (PyArrayObject *)PyArray_New(&PyArray_Type, 1, &count, - PyArray_INTP, NULL, NULL, 0, 0, NULL); - if (indices == NULL) { - Py_DECREF(self); - return NULL; + diag_size = dim2 - offset; + if (dim1 < diag_size) { + diag_size = dim1; + } } - dptr = (intp *)PyArray_DATA(indices); - for (n1 = start; n1 < stop; n1 += step) { - *dptr++ = n1; + } + else if (offset < 0) { + offset = -offset; + if (offset >= dim1) { + diag_size = 0; } - a = PyArray_IterNew((PyObject *)self); - Py_DECREF(self); - if (a == NULL) { - Py_DECREF(indices); - return NULL; + else { + data += offset * stride1; + maskna_data += offset * maskna_stride1; + + diag_size = dim1 - offset; + if (dim2 < diag_size) { + diag_size = dim2; + } } - ret = PyObject_GetItem(a, (PyObject *)indices); - Py_DECREF(a); - Py_DECREF(indices); - return ret; } - else { - /* - * my_diagonal = [] - * for i in range (s [0]) : - * my_diagonal.append (diagonal (a [i], offset)) - * return array (my_diagonal) - */ - PyObject *mydiagonal = NULL, *ret = NULL, *sel = NULL; - intp n1; - int res; - PyArray_Descr *typecode; - - new = NULL; - - typecode = PyArray_DESCR(self); - mydiagonal = PyList_New(0); - if (mydiagonal == NULL) { - Py_DECREF(self); - return NULL; + diag_size = dim1 < dim2 ? dim1 : dim2; + } + + /* Build the new shape and strides for the main data */ + i = 0; + for (idim = 0; idim < ndim; ++idim) { + if (idim != axis1 && idim != axis2) { + ret_shape[i] = shape[idim]; + ret_strides[i] = strides[idim]; + ++i; } - n1 = PyArray_DIMS(self)[0]; - for (i = 0; i < n1; i++) { - new = PyInt_FromLong((long) i); - sel = PyArray_EnsureAnyArray(PyObject_GetItem((PyObject *)self, new)); - Py_DECREF(new); - if (sel == NULL) { - Py_DECREF(self); - Py_DECREF(mydiagonal); - return NULL; - } - new = PyArray_Diagonal((PyArrayObject *)sel, offset, n-3, n-2); - Py_DECREF(sel); - if (new == NULL) { - Py_DECREF(self); - Py_DECREF(mydiagonal); - return NULL; - } - res = PyList_Append(mydiagonal, new); - Py_DECREF(new); - if (res < 0) { - Py_DECREF(self); - Py_DECREF(mydiagonal); - return NULL; + } + ret_shape[ndim-2] = diag_size; + ret_strides[ndim-2] = stride1 + stride2; + + /* Create the diagonal view */ + dtype = PyArray_DTYPE(self); + Py_INCREF(dtype); + ret = (PyArrayObject *)PyArray_NewFromDescr(Py_TYPE(self), + dtype, + ndim-1, ret_shape, + ret_strides, + data, + PyArray_FLAGS(self) & ~(NPY_ARRAY_MASKNA | NPY_ARRAY_OWNMASKNA), + (PyObject *)self); + if (ret == NULL) { + return NULL; + } + Py_INCREF(self); + if (PyArray_SetBaseObject(ret, (PyObject *)self) < 0) { + Py_DECREF(ret); + return NULL; + } + + /* Take a view of the mask if it exists */ + if (self_has_maskna) { + PyArrayObject_fieldaccess *fret = (PyArrayObject_fieldaccess *)ret; + npy_intp *maskna_strides = PyArray_MASKNA_STRIDES(self); + + fret->maskna_dtype = PyArray_MASKNA_DTYPE(self); + Py_INCREF(fret->maskna_dtype); + fret->maskna_data = maskna_data; + /* Build the strides for the mask */ + i = 0; + for (idim = 0; idim < ndim; ++idim) { + if (idim != axis1 && idim != axis2) { + fret->maskna_strides[i] = maskna_strides[idim]; + ++i; } } - Py_DECREF(self); - Py_INCREF(typecode); - ret = PyArray_FromAny(mydiagonal, typecode, 0, 0, 0, NULL); - Py_DECREF(mydiagonal); - return ret; + fret->maskna_strides[ndim-2] = maskna_stride1 + maskna_stride2; + fret->flags |= NPY_ARRAY_MASKNA; } + + return (PyObject *)ret; } /*NUMPY_API diff --git a/numpy/core/src/multiarray/methods.c b/numpy/core/src/multiarray/methods.c index 5fdeea7f1..649ff734d 100644 --- a/numpy/core/src/multiarray/methods.c +++ b/numpy/core/src/multiarray/methods.c @@ -2076,6 +2076,7 @@ array_diagonal(PyArrayObject *self, PyObject *args, PyObject *kwds) { int axis1 = 0, axis2 = 1, offset = 0; static char *kwlist[] = {"offset", "axis1", "axis2", NULL}; + PyArrayObject *ret; if (!PyArg_ParseTupleAndKeywords(args, kwds, "|iii", kwlist, &offset, @@ -2083,7 +2084,9 @@ array_diagonal(PyArrayObject *self, PyObject *args, PyObject *kwds) &axis2)) { return NULL; } - return PyArray_Return((PyArrayObject *)PyArray_Diagonal(self, offset, axis1, axis2)); + + ret = (PyArrayObject *)PyArray_Diagonal(self, offset, axis1, axis2); + return PyArray_Return(ret); } diff --git a/numpy/core/tests/test_maskna.py b/numpy/core/tests/test_maskna.py index 0e5cc33bc..85bcde99c 100644 --- a/numpy/core/tests/test_maskna.py +++ b/numpy/core/tests/test_maskna.py @@ -890,5 +890,35 @@ def test_array_maskna_conjugate_method(): assert_equal(np.isna(b), [0,0,1,0,1]) assert_equal(b[~np.isna(b)], [-1j, 2-4j, 2+1.5j]) +def test_array_maskna_diagonal(): + # ndarray.diagonal + a = np.arange(6, maskna=True) + a.shape = (2,3) + a[0,1] = np.NA + + # Should produce a view into a + res = a.diagonal() + assert_(res.base is a) + assert_(res.flags.maskna) + assert_(not res.flags.ownmaskna) + assert_equal(res, [0, 4]) + + res = a.diagonal(-1) + assert_equal(res, [3]) + + res = a.diagonal(-2) + assert_equal(res, []) + + # This diagonal has the NA + res = a.diagonal(1) + assert_equal(np.isna(res), [1,0]) + assert_equal(res[~np.isna(res)], [5]) + + res = a.diagonal(2) + assert_equal(res, [2]) + + res = a.diagonal(3) + assert_equal(res, []) + if __name__ == "__main__": run_module_suite() diff --git a/numpy/lib/twodim_base.py b/numpy/lib/twodim_base.py index d95a59e3f..9b3b50d04 100644 --- a/numpy/lib/twodim_base.py +++ b/numpy/lib/twodim_base.py @@ -221,10 +221,12 @@ def diag(v, k=0): """ Extract a diagonal or construct a diagonal array. + As of NumPy 1.7, extracting a diagonal always returns a view into `v`. + Parameters ---------- v : array_like - If `v` is a 2-D array, return a copy of its `k`-th diagonal. + If `v` is a 2-D array, return a view of its `k`-th diagonal. If `v` is a 1-D array, return a 2-D array with `v` on the `k`-th diagonal. k : int, optional @@ -278,16 +280,7 @@ def diag(v, k=0): res[:n-k].flat[i::n+1] = v return res elif len(s) == 2: - if k >= s[1]: - return empty(0, dtype=v.dtype) - if v.flags.f_contiguous: - # faster slicing - v, k, s = v.T, -k, s[::-1] - if k >= 0: - i = k - else: - i = (-k) * s[1] - return v[:s[1]-k].flat[i::s[1]+1] + return v.diagonal(k) else: raise ValueError("Input must be 1- or 2-d.") |