summaryrefslogtreecommitdiff
path: root/numpy
diff options
context:
space:
mode:
Diffstat (limited to 'numpy')
-rw-r--r--numpy/core/src/umath/ufunc_object.c32
-rw-r--r--numpy/lib/histograms.py32
-rw-r--r--numpy/lib/tests/test_histograms.py41
-rw-r--r--numpy/lib/twodim_base.py2
-rw-r--r--numpy/linalg/linalg.py117
-rw-r--r--numpy/linalg/tests/test_linalg.py33
-rw-r--r--numpy/matrixlib/defmatrix.py115
-rw-r--r--numpy/matrixlib/tests/test_defmatrix.py2
8 files changed, 201 insertions, 173 deletions
diff --git a/numpy/core/src/umath/ufunc_object.c b/numpy/core/src/umath/ufunc_object.c
index 84a329475..c1e8e5a77 100644
--- a/numpy/core/src/umath/ufunc_object.c
+++ b/numpy/core/src/umath/ufunc_object.c
@@ -4282,11 +4282,9 @@ static PyObject *
ufunc_generic_call(PyUFuncObject *ufunc, PyObject *args, PyObject *kwds)
{
int i;
- PyTupleObject *ret;
PyArrayObject *mps[NPY_MAXARGS];
PyObject *retobj[NPY_MAXARGS];
PyObject *wraparr[NPY_MAXARGS];
- PyObject *res;
PyObject *override = NULL;
ufunc_full_args full_args = {NULL, NULL};
int errval;
@@ -4363,13 +4361,17 @@ ufunc_generic_call(PyUFuncObject *ufunc, PyObject *args, PyObject *kwds)
int j = ufunc->nin+i;
PyObject *wrap = wraparr[i];
- if (wrap != NULL) {
+ if (wrap == NULL) {
+ /* default behavior */
+ retobj[i] = PyArray_Return(mps[j]);
+ }
+ else if (wrap == Py_None) {
+ Py_DECREF(wrap);
+ retobj[i] = (PyObject *)mps[j];
+ }
+ else {
+ PyObject *res;
PyObject *args_tup;
- if (wrap == Py_None) {
- Py_DECREF(wrap);
- retobj[i] = (PyObject *)mps[j];
- continue;
- }
/* Call the method with appropriate context */
args_tup = _get_wrap_prepare_args(full_args);
@@ -4389,15 +4391,9 @@ ufunc_generic_call(PyUFuncObject *ufunc, PyObject *args, PyObject *kwds)
if (res == NULL) {
goto fail;
}
- else {
- Py_DECREF(mps[j]);
- retobj[i] = res;
- continue;
- }
- }
- else {
- /* default behavior */
- retobj[i] = PyArray_Return(mps[j]);
+
+ Py_DECREF(mps[j]);
+ retobj[i] = res;
}
}
@@ -4408,6 +4404,8 @@ ufunc_generic_call(PyUFuncObject *ufunc, PyObject *args, PyObject *kwds)
return retobj[0];
}
else {
+ PyTupleObject *ret;
+
ret = (PyTupleObject *)PyTuple_New(ufunc->nout);
for (i = 0; i < ufunc->nout; i++) {
PyTuple_SET_ITEM(ret, i, retobj[i]);
diff --git a/numpy/lib/histograms.py b/numpy/lib/histograms.py
index 90e19769e..2922b3a86 100644
--- a/numpy/lib/histograms.py
+++ b/numpy/lib/histograms.py
@@ -877,12 +877,6 @@ def histogramdd(sample, bins=10, range=None, normed=False, weights=None):
# bins is an integer
bins = D*[bins]
- # avoid rounding issues for comparisons when dealing with inexact types
- if np.issubdtype(sample.dtype, np.inexact):
- edge_dt = sample.dtype
- else:
- edge_dt = float
-
# normalize the range argument
if range is None:
range = (None,) * D
@@ -896,13 +890,12 @@ def histogramdd(sample, bins=10, range=None, normed=False, weights=None):
raise ValueError(
'`bins[{}]` must be positive, when an integer'.format(i))
smin, smax = _get_outer_edges(sample[:,i], range[i])
- edges[i] = np.linspace(smin, smax, bins[i] + 1, dtype=edge_dt)
+ edges[i] = np.linspace(smin, smax, bins[i] + 1)
elif np.ndim(bins[i]) == 1:
- edges[i] = np.asarray(bins[i], edge_dt)
- # not just monotonic, due to the use of mindiff below
- if np.any(edges[i][:-1] >= edges[i][1:]):
+ edges[i] = np.asarray(bins[i])
+ if np.any(edges[i][:-1] > edges[i][1:]):
raise ValueError(
- '`bins[{}]` must be strictly increasing, when an array'
+ '`bins[{}]` must be monotonically increasing, when an array'
.format(i))
else:
raise ValueError(
@@ -913,7 +906,8 @@ def histogramdd(sample, bins=10, range=None, normed=False, weights=None):
# Compute the bin number each sample falls into.
Ncount = tuple(
- np.digitize(sample[:, i], edges[i])
+ # avoid np.digitize to work around gh-11022
+ np.searchsorted(edges[i], sample[:, i], side='right')
for i in _range(D)
)
@@ -921,16 +915,10 @@ def histogramdd(sample, bins=10, range=None, normed=False, weights=None):
# For the rightmost bin, we want values equal to the right edge to be
# counted in the last bin, and not as an outlier.
for i in _range(D):
- # Rounding precision
- mindiff = dedges[i].min()
- if not np.isinf(mindiff):
- decimal = int(-np.log10(mindiff)) + 6
- # Find which points are on the rightmost edge.
- not_smaller_than_edge = (sample[:, i] >= edges[i][-1])
- on_edge = (np.around(sample[:, i], decimal) ==
- np.around(edges[i][-1], decimal))
- # Shift these points one bin to the left.
- Ncount[i][on_edge & not_smaller_than_edge] -= 1
+ # Find which points are on the rightmost edge.
+ on_edge = (sample[:, i] == edges[i][-1])
+ # Shift these points one bin to the left.
+ Ncount[i][on_edge] -= 1
# Compute the sample indices in the flattened histogram matrix.
# This raises an error if the array is too large.
diff --git a/numpy/lib/tests/test_histograms.py b/numpy/lib/tests/test_histograms.py
index 597b5b376..e16ae12c2 100644
--- a/numpy/lib/tests/test_histograms.py
+++ b/numpy/lib/tests/test_histograms.py
@@ -613,8 +613,6 @@ class TestHistogramdd(object):
assert_raises(ValueError, np.histogramdd, x, bins=[-1, 2, 4, 5])
assert_raises(ValueError, np.histogramdd, x, bins=[1, 0.99, 1, 1])
assert_raises(
- ValueError, np.histogramdd, x, bins=[1, 1, 1, [1, 2, 2, 3]])
- assert_raises(
ValueError, np.histogramdd, x, bins=[1, 1, 1, [1, 2, 3, -3]])
assert_(np.histogramdd(x, bins=[1, 1, 1, [1, 2, 3, 4]]))
@@ -646,7 +644,7 @@ class TestHistogramdd(object):
bins = [[0., 0.5, 1.0]]
hist, _ = histogramdd(x, bins=bins)
assert_(hist[0] == 0.0)
- assert_(hist[1] == 1.)
+ assert_(hist[1] == 0.0)
x = [1.0001]
bins = [[0., 0.5, 1.0]]
hist, _ = histogramdd(x, bins=bins)
@@ -660,3 +658,40 @@ class TestHistogramdd(object):
range=[[0.0, 1.0], [0.25, 0.75], [0.25, np.inf]])
assert_raises(ValueError, histogramdd, vals,
range=[[0.0, 1.0], [np.nan, 0.75], [0.25, 0.5]])
+
+ def test_equal_edges(self):
+ """ Test that adjacent entries in an edge array can be equal """
+ x = np.array([0, 1, 2])
+ y = np.array([0, 1, 2])
+ x_edges = np.array([0, 2, 2])
+ y_edges = 1
+ hist, edges = histogramdd((x, y), bins=(x_edges, y_edges))
+
+ hist_expected = np.array([
+ [2.],
+ [1.], # x == 2 falls in the final bin
+ ])
+ assert_equal(hist, hist_expected)
+
+ def test_edge_dtype(self):
+ """ Test that if an edge array is input, its type is preserved """
+ x = np.array([0, 10, 20])
+ y = x / 10
+ x_edges = np.array([0, 5, 15, 20])
+ y_edges = x_edges / 10
+ hist, edges = histogramdd((x, y), bins=(x_edges, y_edges))
+
+ assert_equal(edges[0].dtype, x_edges.dtype)
+ assert_equal(edges[1].dtype, y_edges.dtype)
+
+ def test_large_integers(self):
+ big = 2**60 # Too large to represent with a full precision float
+
+ x = np.array([0], np.int64)
+ x_edges = np.array([-1, +1], np.int64)
+ y = big + x
+ y_edges = big + x_edges
+
+ hist, edges = histogramdd((x, y), bins=(x_edges, y_edges))
+
+ assert_equal(hist[0, 0], 1)
diff --git a/numpy/lib/twodim_base.py b/numpy/lib/twodim_base.py
index 402c18850..cca316e9a 100644
--- a/numpy/lib/twodim_base.py
+++ b/numpy/lib/twodim_base.py
@@ -650,7 +650,7 @@ def histogram2d(x, y, bins=10, range=None, normed=False, weights=None):
N = 1
if N != 1 and N != 2:
- xedges = yedges = asarray(bins, float)
+ xedges = yedges = asarray(bins)
bins = [xedges, yedges]
hist, edges = histogramdd([x, y], bins, range, normed, weights)
return hist, edges[0], edges[1]
diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py
index 5ee230f92..5757b1827 100644
--- a/numpy/linalg/linalg.py
+++ b/numpy/linalg/linalg.py
@@ -16,20 +16,20 @@ __all__ = ['matrix_power', 'solve', 'tensorsolve', 'tensorinv', 'inv',
'svd', 'eig', 'eigh', 'lstsq', 'norm', 'qr', 'cond', 'matrix_rank',
'LinAlgError', 'multi_dot']
+import operator
import warnings
from numpy.core import (
array, asarray, zeros, empty, empty_like, intc, single, double,
- csingle, cdouble, inexact, complexfloating, newaxis, ravel, all, Inf, dot,
- add, multiply, sqrt, maximum, fastCopyAndTranspose, sum, isfinite, size,
- finfo, errstate, geterrobj, longdouble, moveaxis, amin, amax, product, abs,
- broadcast, atleast_2d, intp, asanyarray, object_, ones, matmul,
- swapaxes, divide, count_nonzero, ndarray, isnan
+ csingle, cdouble, inexact, complexfloating, newaxis, all, Inf, dot,
+ add, multiply, sqrt, fastCopyAndTranspose, sum, isfinite,
+ finfo, errstate, geterrobj, moveaxis, amin, amax, product, abs,
+ atleast_2d, intp, asanyarray, object_, matmul,
+ swapaxes, divide, count_nonzero, isnan
)
from numpy.core.multiarray import normalize_axis_index
-from numpy.lib import triu, asfarray
+from numpy.lib.twodim_base import triu, eye
from numpy.linalg import lapack_lite, _umath_linalg
-from numpy.matrixlib.defmatrix import matrix_power
# For Python2/3 compatibility
_N = b'N'
@@ -532,6 +532,109 @@ def inv(a):
return wrap(ainv.astype(result_t, copy=False))
+def matrix_power(a, n):
+ """
+ Raise a square matrix to the (integer) power `n`.
+
+ For positive integers `n`, the power is computed by repeated matrix
+ squarings and matrix multiplications. If ``n == 0``, the identity matrix
+ of the same shape as M is returned. If ``n < 0``, the inverse
+ is computed and then raised to the ``abs(n)``.
+
+ Parameters
+ ----------
+ a : (..., M, M) array_like
+ Matrix to be "powered."
+ n : int
+ The exponent can be any integer or long integer, positive,
+ negative, or zero.
+
+ Returns
+ -------
+ a**n : (..., M, M) ndarray or matrix object
+ The return value is the same shape and type as `M`;
+ if the exponent is positive or zero then the type of the
+ elements is the same as those of `M`. If the exponent is
+ negative the elements are floating-point.
+
+ Raises
+ ------
+ LinAlgError
+ For matrices that are not square or that (for negative powers) cannot
+ be inverted numerically.
+
+ Examples
+ --------
+ >>> from numpy.linalg import matrix_power
+ >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit
+ >>> matrix_power(i, 3) # should = -i
+ array([[ 0, -1],
+ [ 1, 0]])
+ >>> matrix_power(i, 0)
+ array([[1, 0],
+ [0, 1]])
+ >>> matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements
+ array([[ 0., 1.],
+ [-1., 0.]])
+
+ Somewhat more sophisticated example
+
+ >>> q = np.zeros((4, 4))
+ >>> q[0:2, 0:2] = -i
+ >>> q[2:4, 2:4] = i
+ >>> q # one of the three quaternion units not equal to 1
+ array([[ 0., -1., 0., 0.],
+ [ 1., 0., 0., 0.],
+ [ 0., 0., 0., 1.],
+ [ 0., 0., -1., 0.]])
+ >>> matrix_power(q, 2) # = -np.eye(4)
+ array([[-1., 0., 0., 0.],
+ [ 0., -1., 0., 0.],
+ [ 0., 0., -1., 0.],
+ [ 0., 0., 0., -1.]])
+
+ """
+ a = asanyarray(a)
+ _assertRankAtLeast2(a)
+ _assertNdSquareness(a)
+
+ try:
+ n = operator.index(n)
+ except TypeError:
+ raise TypeError("exponent must be an integer")
+
+ if n == 0:
+ a = empty_like(a)
+ a[...] = eye(a.shape[-2], dtype=a.dtype)
+ return a
+
+ elif n < 0:
+ a = inv(a)
+ n = abs(n)
+
+ # short-cuts.
+ if n == 1:
+ return a
+
+ elif n == 2:
+ return matmul(a, a)
+
+ elif n == 3:
+ return matmul(matmul(a, a), a)
+
+ # Use binary decomposition to reduce the number of matrix multiplications.
+ # Here, we iterate over the bits of n, from LSB to MSB, raise `a` to
+ # increasing powers of 2, and multiply into the result as needed.
+ z = result = None
+ while n > 0:
+ z = a if z is None else matmul(z, z)
+ n, bit = divmod(n, 2)
+ if bit:
+ result = z if result is None else matmul(result, z)
+
+ return result
+
+
# Cholesky decomposition
def cholesky(a):
diff --git a/numpy/linalg/tests/test_linalg.py b/numpy/linalg/tests/test_linalg.py
index b3dd2e4ae..5ed1ff1c0 100644
--- a/numpy/linalg/tests/test_linalg.py
+++ b/numpy/linalg/tests/test_linalg.py
@@ -10,8 +10,8 @@ import traceback
import pytest
import numpy as np
-from numpy import array, single, double, csingle, cdouble, dot, identity
-from numpy import multiply, atleast_2d, inf, asarray
+from numpy import array, single, double, csingle, cdouble, dot, identity, matmul
+from numpy import multiply, atleast_2d, inf, asarray, matrix
from numpy import linalg
from numpy.linalg import matrix_power, norm, matrix_rank, multi_dot, LinAlgError
from numpy.linalg.linalg import _multi_dot_matrix_chain_order
@@ -445,8 +445,7 @@ def identity_like_generalized(a):
a = asarray(a)
if a.ndim >= 3:
r = np.empty(a.shape, dtype=a.dtype)
- for c in itertools.product(*map(range, a.shape[:-2])):
- r[c] = identity(a.shape[-2])
+ r[...] = identity(a.shape[-2])
return r
else:
return identity(a.shape[0])
@@ -927,16 +926,21 @@ class TestMatrixPower(object):
R90 = array([[0, 1], [-1, 0]])
Arb22 = array([[4, -7], [-2, 10]])
noninv = array([[1, 0], [0, 0]])
- arbfloat = array([[0.1, 3.2], [1.2, 0.7]])
+ arbfloat = array([[[0.1, 3.2], [1.2, 0.7]],
+ [[0.2, 6.4], [2.4, 1.4]]])
large = identity(10)
t = large[1, :].copy()
- large[1, :] = large[0,:]
+ large[1, :] = large[0, :]
large[0, :] = t
def test_large_power(self):
assert_equal(
matrix_power(self.R90, 2 ** 100 + 2 ** 10 + 2 ** 5 + 1), self.R90)
+ assert_equal(
+ matrix_power(self.R90, 2 ** 100 + 2 ** 10 + 1), self.R90)
+ assert_equal(
+ matrix_power(self.R90, 2 ** 100 + 2 + 1), -self.R90)
def test_large_power_trailing_zero(self):
assert_equal(
@@ -945,7 +949,7 @@ class TestMatrixPower(object):
def testip_zero(self):
def tz(M):
mz = matrix_power(M, 0)
- assert_equal(mz, identity(M.shape[0]))
+ assert_equal(mz, identity_like_generalized(M))
assert_equal(mz.dtype, M.dtype)
for M in [self.Arb22, self.arbfloat, self.large]:
tz(M)
@@ -961,7 +965,7 @@ class TestMatrixPower(object):
def testip_two(self):
def tz(M):
mz = matrix_power(M, 2)
- assert_equal(mz, dot(M, M))
+ assert_equal(mz, matmul(M, M))
assert_equal(mz.dtype, M.dtype)
for M in [self.Arb22, self.arbfloat, self.large]:
tz(M)
@@ -969,14 +973,19 @@ class TestMatrixPower(object):
def testip_invert(self):
def tz(M):
mz = matrix_power(M, -1)
- assert_almost_equal(identity(M.shape[0]), dot(mz, M))
+ assert_almost_equal(matmul(mz, M), identity_like_generalized(M))
for M in [self.R90, self.Arb22, self.arbfloat, self.large]:
tz(M)
def test_invert_noninvertible(self):
- import numpy.linalg
- assert_raises(numpy.linalg.linalg.LinAlgError,
- lambda: matrix_power(self.noninv, -1))
+ assert_raises(LinAlgError, matrix_power, self.noninv, -1)
+
+ def test_invalid(self):
+ assert_raises(TypeError, matrix_power, self.R90, 1.5)
+ assert_raises(TypeError, matrix_power, self.R90, [1])
+ assert_raises(LinAlgError, matrix_power, np.array([1]), 1)
+ assert_raises(LinAlgError, matrix_power, np.array([[1], [2]]), 1)
+ assert_raises(LinAlgError, matrix_power, np.ones((4, 3, 2)), 1)
class TestBoolPower(object):
diff --git a/numpy/matrixlib/defmatrix.py b/numpy/matrixlib/defmatrix.py
index 1f5c94921..9909fec8d 100644
--- a/numpy/matrixlib/defmatrix.py
+++ b/numpy/matrixlib/defmatrix.py
@@ -5,8 +5,11 @@ __all__ = ['matrix', 'bmat', 'mat', 'asmatrix']
import sys
import ast
import numpy.core.numeric as N
-from numpy.core.numeric import concatenate, isscalar, binary_repr, identity, asanyarray
-from numpy.core.numerictypes import issubdtype
+from numpy.core.numeric import concatenate, isscalar
+# While not in __all__, matrix_power used to be defined here, so we import
+# it for backward compatibility.
+from numpy.linalg import matrix_power
+
def _convert_from_string(data):
for char in '[]':
@@ -63,114 +66,6 @@ def asmatrix(data, dtype=None):
"""
return matrix(data, dtype=dtype, copy=False)
-def matrix_power(M, n):
- """
- Raise a square matrix to the (integer) power `n`.
-
- For positive integers `n`, the power is computed by repeated matrix
- squarings and matrix multiplications. If ``n == 0``, the identity matrix
- of the same shape as M is returned. If ``n < 0``, the inverse
- is computed and then raised to the ``abs(n)``.
-
- Parameters
- ----------
- M : ndarray or matrix object
- Matrix to be "powered." Must be square, i.e. ``M.shape == (m, m)``,
- with `m` a positive integer.
- n : int
- The exponent can be any integer or long integer, positive,
- negative, or zero.
-
- Returns
- -------
- M**n : ndarray or matrix object
- The return value is the same shape and type as `M`;
- if the exponent is positive or zero then the type of the
- elements is the same as those of `M`. If the exponent is
- negative the elements are floating-point.
-
- Raises
- ------
- LinAlgError
- If the matrix is not numerically invertible.
-
- See Also
- --------
- matrix
- Provides an equivalent function as the exponentiation operator
- (``**``, not ``^``).
-
- Examples
- --------
- >>> from numpy import linalg as LA
- >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit
- >>> LA.matrix_power(i, 3) # should = -i
- array([[ 0, -1],
- [ 1, 0]])
- >>> LA.matrix_power(np.matrix(i), 3) # matrix arg returns matrix
- matrix([[ 0, -1],
- [ 1, 0]])
- >>> LA.matrix_power(i, 0)
- array([[1, 0],
- [0, 1]])
- >>> LA.matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements
- array([[ 0., 1.],
- [-1., 0.]])
-
- Somewhat more sophisticated example
-
- >>> q = np.zeros((4, 4))
- >>> q[0:2, 0:2] = -i
- >>> q[2:4, 2:4] = i
- >>> q # one of the three quaternion units not equal to 1
- array([[ 0., -1., 0., 0.],
- [ 1., 0., 0., 0.],
- [ 0., 0., 0., 1.],
- [ 0., 0., -1., 0.]])
- >>> LA.matrix_power(q, 2) # = -np.eye(4)
- array([[-1., 0., 0., 0.],
- [ 0., -1., 0., 0.],
- [ 0., 0., -1., 0.],
- [ 0., 0., 0., -1.]])
-
- """
- M = asanyarray(M)
- if M.ndim != 2 or M.shape[0] != M.shape[1]:
- raise ValueError("input must be a square array")
- if not issubdtype(type(n), N.integer):
- raise TypeError("exponent must be an integer")
-
- from numpy.linalg import inv
-
- if n==0:
- M = M.copy()
- M[:] = identity(M.shape[0])
- return M
- elif n<0:
- M = inv(M)
- n *= -1
-
- result = M
- if n <= 3:
- for _ in range(n-1):
- result=N.dot(result, M)
- return result
-
- # binary decomposition to reduce the number of Matrix
- # multiplications for n > 3.
- beta = binary_repr(n)
- Z, q, t = M, 0, len(beta)
- while beta[t-q-1] == '0':
- Z = N.dot(Z, Z)
- q += 1
- result = Z
- for k in range(q+1, t):
- Z = N.dot(Z, Z)
- if beta[t-k-1] == '1':
- result = N.dot(result, Z)
- return result
-
-
class matrix(N.ndarray):
"""
matrix(data, dtype=None, copy=True)
diff --git a/numpy/matrixlib/tests/test_defmatrix.py b/numpy/matrixlib/tests/test_defmatrix.py
index a02a05c09..d160490b3 100644
--- a/numpy/matrixlib/tests/test_defmatrix.py
+++ b/numpy/matrixlib/tests/test_defmatrix.py
@@ -13,7 +13,7 @@ from numpy.testing import (
assert_, assert_equal, assert_almost_equal, assert_array_equal,
assert_array_almost_equal, assert_raises
)
-from numpy.matrixlib.defmatrix import matrix_power
+from numpy.linalg import matrix_power
from numpy.matrixlib import mat
class TestCtor(object):
View Lorry Controller Status