summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--doc/release/2.0.0-notes.rst22
-rw-r--r--numpy/add_newdocs.py2
-rw-r--r--numpy/core/fromnumeric.py2
-rw-r--r--numpy/core/src/multiarray/item_selection.c229
-rw-r--r--numpy/core/src/multiarray/methods.c5
-rw-r--r--numpy/core/tests/test_maskna.py30
-rw-r--r--numpy/lib/twodim_base.py15
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.")