diff options
Diffstat (limited to 'numpy/lib')
| -rw-r--r-- | numpy/lib/__init__.py | 1 | ||||
| -rw-r--r-- | numpy/lib/_version.py | 155 | ||||
| -rw-r--r-- | numpy/lib/arraysetops.py | 73 | ||||
| -rw-r--r-- | numpy/lib/bscript | 6 | ||||
| -rw-r--r-- | numpy/lib/financial.py | 20 | ||||
| -rw-r--r-- | numpy/lib/format.py | 175 | ||||
| -rw-r--r-- | numpy/lib/function_base.py | 474 | ||||
| -rw-r--r-- | numpy/lib/nanfunctions.py | 339 | ||||
| -rw-r--r-- | numpy/lib/npyio.py | 80 | ||||
| -rw-r--r-- | numpy/lib/polynomial.py | 10 | ||||
| -rw-r--r-- | numpy/lib/shape_base.py | 32 | ||||
| -rw-r--r-- | numpy/lib/src/_compiled_base.c | 148 | ||||
| -rw-r--r-- | numpy/lib/stride_tricks.py | 3 | ||||
| -rw-r--r-- | numpy/lib/tests/test__version.py | 57 | ||||
| -rw-r--r-- | numpy/lib/tests/test_arraysetops.py | 64 | ||||
| -rw-r--r-- | numpy/lib/tests/test_financial.py | 15 | ||||
| -rw-r--r-- | numpy/lib/tests/test_format.py | 84 | ||||
| -rw-r--r-- | numpy/lib/tests/test_function_base.py | 309 | ||||
| -rw-r--r-- | numpy/lib/tests/test_io.py | 111 | ||||
| -rw-r--r-- | numpy/lib/tests/test_nanfunctions.py | 266 | ||||
| -rw-r--r-- | numpy/lib/tests/test_stride_tricks.py | 17 | ||||
| -rw-r--r-- | numpy/lib/tests/test_twodim_base.py | 80 | ||||
| -rw-r--r-- | numpy/lib/tests/test_utils.py | 9 | ||||
| -rw-r--r-- | numpy/lib/twodim_base.py | 175 | ||||
| -rw-r--r-- | numpy/lib/utils.py | 7 |
25 files changed, 2207 insertions, 503 deletions
diff --git a/numpy/lib/__init__.py b/numpy/lib/__init__.py index 73e4b2306..8c420b0c3 100644 --- a/numpy/lib/__init__.py +++ b/numpy/lib/__init__.py @@ -23,6 +23,7 @@ from .npyio import * from .financial import * from .arrayterator import * from .arraypad import * +from ._version import * __all__ = ['emath', 'math'] __all__ += type_check.__all__ diff --git a/numpy/lib/_version.py b/numpy/lib/_version.py new file mode 100644 index 000000000..4eaabd0ff --- /dev/null +++ b/numpy/lib/_version.py @@ -0,0 +1,155 @@ +"""Utility to compare (Numpy) version strings. + +The NumpyVersion class allows properly comparing numpy version strings. +The LooseVersion and StrictVersion classes that distutils provides don't +work; they don't recognize anything like alpha/beta/rc/dev versions. + +""" +from __future__ import division, absolute_import, print_function + +import re + +from numpy.compat import basestring + + +__all__ = ['NumpyVersion'] + + +class NumpyVersion(): + """Parse and compare numpy version strings. + + Numpy has the following versioning scheme (numbers given are examples; they + can be > 9) in principle): + + - Released version: '1.8.0', '1.8.1', etc. + - Alpha: '1.8.0a1', '1.8.0a2', etc. + - Beta: '1.8.0b1', '1.8.0b2', etc. + - Release candidates: '1.8.0rc1', '1.8.0rc2', etc. + - Development versions: '1.8.0.dev-f1234afa' (git commit hash appended) + - Development versions after a1: '1.8.0a1.dev-f1234afa', + '1.8.0b2.dev-f1234afa', + '1.8.1rc1.dev-f1234afa', etc. + - Development versions (no git hash available): '1.8.0.dev-Unknown' + + Comparing needs to be done against a valid version string or other + `NumpyVersion` instance. Note that all development versions of the same + (pre-)release compare equal. + + .. versionadded:: 1.9.0 + + Parameters + ---------- + vstring : str + Numpy version string (``np.__version__``). + + Examples + -------- + >>> from numpy.lib import NumpyVersion + >>> if NumpyVersion(np.__version__) < '1.7.0'): + ... print('skip') + skip + + >>> NumpyVersion('1.7') # raises ValueError, add ".0" + + """ + def __init__(self, vstring): + self.vstring = vstring + ver_main = re.match(r'\d[.]\d+[.]\d+', vstring) + if not ver_main: + raise ValueError("Not a valid numpy version string") + + self.version = ver_main.group() + self.major, self.minor, self.bugfix = [int(x) for x in + self.version.split('.')] + if len(vstring) == ver_main.end(): + self.pre_release = 'final' + else: + alpha = re.match(r'a\d', vstring[ver_main.end():]) + beta = re.match(r'b\d', vstring[ver_main.end():]) + rc = re.match(r'rc\d', vstring[ver_main.end():]) + pre_rel = [m for m in [alpha, beta, rc] if m is not None] + if pre_rel: + self.pre_release = pre_rel[0].group() + else: + self.pre_release = '' + + self.is_devversion = bool(re.search(r'.dev', vstring)) + + def _compare_version(self, other): + """Compare major.minor.bugfix""" + if self.major == other.major: + if self.minor == other.minor: + if self.bugfix == other.bugfix: + vercmp = 0 + elif self.bugfix > other.bugfix: + vercmp = 1 + else: + vercmp = -1 + elif self.minor > other.minor: + vercmp = 1 + else: + vercmp = -1 + elif self.major > other.major: + vercmp = 1 + else: + vercmp = -1 + + return vercmp + + def _compare_pre_release(self, other): + """Compare alpha/beta/rc/final.""" + if self.pre_release == other.pre_release: + vercmp = 0 + elif self.pre_release == 'final': + vercmp = 1 + elif other.pre_release == 'final': + vercmp = -1 + elif self.pre_release > other.pre_release: + vercmp = 1 + else: + vercmp = -1 + + return vercmp + + def _compare(self, other): + if not isinstance(other, (basestring, NumpyVersion)): + raise ValueError("Invalid object to compare with NumpyVersion.") + + if isinstance(other, basestring): + other = NumpyVersion(other) + + vercmp = self._compare_version(other) + if vercmp == 0: + # Same x.y.z version, check for alpha/beta/rc + vercmp = self._compare_pre_release(other) + if vercmp == 0: + # Same version and same pre-release, check if dev version + if self.is_devversion is other.is_devversion: + vercmp = 0 + elif self.is_devversion: + vercmp = -1 + else: + vercmp = 1 + + return vercmp + + def __lt__(self, other): + return self._compare(other) < 0 + + def __le__(self, other): + return self._compare(other) <= 0 + + def __eq__(self, other): + return self._compare(other) == 0 + + def __ne__(self, other): + return self._compare(other) != 0 + + def __gt__(self, other): + return self._compare(other) > 0 + + def __ge__(self, other): + return self._compare(other) >= 0 + + def __repr(self): + return "NumpyVersion(%s)" % self.vstring diff --git a/numpy/lib/arraysetops.py b/numpy/lib/arraysetops.py index 5cd535703..005703d16 100644 --- a/numpy/lib/arraysetops.py +++ b/numpy/lib/arraysetops.py @@ -90,7 +90,7 @@ def ediff1d(ary, to_end=None, to_begin=None): return ed -def unique(ar, return_index=False, return_inverse=False): +def unique(ar, return_index=False, return_inverse=False, return_counts=False): """ Find the unique elements of an array. @@ -109,6 +109,10 @@ def unique(ar, return_index=False, return_inverse=False): return_inverse : bool, optional If True, also return the indices of the unique array that can be used to reconstruct `ar`. + return_counts : bool, optional + .. versionadded:: 1.9.0 + If True, also return the number of times each unique value comes up + in `ar`. Returns ------- @@ -120,6 +124,10 @@ def unique(ar, return_index=False, return_inverse=False): unique_inverse : ndarray, optional The indices to reconstruct the (flattened) original array from the unique array. Only provided if `return_inverse` is True. + unique_counts : ndarray, optional + .. versionadded:: 1.9.0 + The number of times each of the unique values comes up in the + original array. Only provided if `return_counts` is True. See Also -------- @@ -159,45 +167,46 @@ def unique(ar, return_index=False, return_inverse=False): array([1, 2, 6, 4, 2, 3, 2]) """ - try: - ar = ar.flatten() - except AttributeError: - if not return_inverse and not return_index: - items = sorted(set(ar)) - return np.asarray(items) - else: - ar = np.asanyarray(ar).flatten() + ar = np.asanyarray(ar).flatten() + + optional_indices = return_index or return_inverse + optional_returns = optional_indices or return_counts if ar.size == 0: - if return_inverse and return_index: - return ar, np.empty(0, np.bool), np.empty(0, np.bool) - elif return_inverse or return_index: - return ar, np.empty(0, np.bool) + if not optional_returns: + ret = ar else: - return ar + ret = (ar,) + if return_index: + ret += (np.empty(0, np.bool),) + if return_inverse: + ret += (np.empty(0, np.bool),) + if return_counts: + ret += (np.empty(0, np.intp),) + return ret + + if optional_indices: + perm = ar.argsort(kind='mergesort' if return_index else 'quicksort') + aux = ar[perm] + else: + ar.sort() + aux = ar + flag = np.concatenate(([True], aux[1:] != aux[:-1])) - if return_inverse or return_index: + if not optional_returns: + ret = aux[flag] + else: + ret = (aux[flag],) if return_index: - perm = ar.argsort(kind='mergesort') - else: - perm = ar.argsort() - aux = ar[perm] - flag = np.concatenate(([True], aux[1:] != aux[:-1])) + ret += (perm[flag],) if return_inverse: iflag = np.cumsum(flag) - 1 iperm = perm.argsort() - if return_index: - return aux[flag], perm[flag], iflag[iperm] - else: - return aux[flag], iflag[iperm] - else: - return aux[flag], perm[flag] - - else: - ar.sort() - flag = np.concatenate(([True], ar[1:] != ar[:-1])) - return ar[flag] - + ret += (np.take(iflag, iperm),) + if return_counts: + idx = np.concatenate(np.nonzero(flag) + ([ar.size],)) + ret += (np.diff(idx),) + return ret def intersect1d(ar1, ar2, assume_unique=False): """ diff --git a/numpy/lib/bscript b/numpy/lib/bscript index a9200d043..61debfe41 100644 --- a/numpy/lib/bscript +++ b/numpy/lib/bscript @@ -4,4 +4,8 @@ from bento.commands import hooks def build(context): context.tweak_extension("_compiled_base", includes=["../core/include", "../core/include/numpy", "../core", - "../core/src/private"]) + "../core/src/private"], + defines=['_FILE_OFFSET_BITS=64', + '_LARGEFILE_SOURCE=1', + '_LARGEFILE64_SOURCE=1'] + ) diff --git a/numpy/lib/financial.py b/numpy/lib/financial.py index ec642afd3..c085a5d53 100644 --- a/numpy/lib/financial.py +++ b/numpy/lib/financial.py @@ -628,21 +628,29 @@ def irr(values): Examples -------- - >>> print round(np.irr([-100, 39, 59, 55, 20]), 5) + >>> round(irr([-100, 39, 59, 55, 20]), 5) 0.28095 + >>> round(irr([-100, 0, 0, 74]), 5) + -0.0955 + >>> round(irr([-100, 100, 0, -7]), 5) + -0.0833 + >>> round(irr([-100, 100, 0, 7]), 5) + 0.06206 + >>> round(irr([-5, 10.5, 1, -8, 1]), 5) + 0.0886 (Compare with the Example given for numpy.lib.financial.npv) """ res = np.roots(values[::-1]) - # Find the root(s) between 0 and 1 - mask = (res.imag == 0) & (res.real > 0) & (res.real <= 1) - res = res[mask].real + mask = (res.imag == 0) & (res.real > 0) if res.size == 0: return np.nan + res = res[mask].real + # NPV(rate) = 0 can have more than one solution so we return + # only the solution closest to zero. rate = 1.0/res - 1 - if rate.size == 1: - rate = rate.item() + rate = rate.item(np.argmin(np.abs(rate))) return rate def npv(rate, values): diff --git a/numpy/lib/format.py b/numpy/lib/format.py index 4cfbbe05d..6083312de 100644 --- a/numpy/lib/format.py +++ b/numpy/lib/format.py @@ -139,6 +139,7 @@ from __future__ import division, absolute_import, print_function import numpy import sys import io +import warnings from numpy.lib.utils import safe_eval from numpy.compat import asbytes, isfileobj, long, basestring @@ -151,6 +152,14 @@ MAGIC_PREFIX = asbytes('\x93NUMPY') MAGIC_LEN = len(MAGIC_PREFIX) + 2 BUFFER_SIZE = 2 ** 18 #size of buffer for reading npz files in bytes +# difference between version 1.0 and 2.0 is a 4 byte (I) header length +# instead of 2 bytes (H) allowing storage of large structured arrays + +def _check_version(version): + if version not in [(1, 0), (2, 0), None]: + msg = "we only support format version (1,0) and (2, 0), not %s" + raise ValueError(msg % (version,)) + def magic(major, minor): """ Return the magic string for the given file format version. @@ -258,8 +267,8 @@ def header_data_from_array_1_0(array): d['descr'] = dtype_to_descr(array.dtype) return d -def write_array_header_1_0(fp, d): - """ Write the header for an array using the 1.0 format. +def _write_array_header(fp, d, version=None): + """ Write the header for an array and returns the version used Parameters ---------- @@ -267,6 +276,14 @@ def write_array_header_1_0(fp, d): d : dict This has the appropriate entries for writing its string representation to the header of the file. + version: tuple or None + None means use oldest that works + explicit version will raise a ValueError if the format does not + allow saving this data. Default: None + Returns + ------- + version : tuple of int + the file version which needs to be used to store the data """ import struct header = ["{"] @@ -282,11 +299,52 @@ def write_array_header_1_0(fp, d): current_header_len = MAGIC_LEN + 2 + len(header) + 1 # 1 for the newline topad = 16 - (current_header_len % 16) header = asbytes(header + ' '*topad + '\n') - if len(header) >= (256*256): - raise ValueError("header does not fit inside %s bytes" % (256*256)) - header_len_str = struct.pack('<H', len(header)) + + if len(header) >= (256*256) and version == (1, 0): + raise ValueError("header does not fit inside %s bytes required by the" + " 1.0 format" % (256*256)) + if len(header) < (256*256): + header_len_str = struct.pack('<H', len(header)) + version = (1, 0) + elif len(header) < (2**32): + header_len_str = struct.pack('<I', len(header)) + version = (2, 0) + else: + raise ValueError("header does not fit inside 4 GiB required by " + "the 2.0 format") + + fp.write(magic(*version)) fp.write(header_len_str) fp.write(header) + return version + +def write_array_header_1_0(fp, d): + """ Write the header for an array using the 1.0 format. + + Parameters + ---------- + fp : filelike object + d : dict + This has the appropriate entries for writing its string representation + to the header of the file. + """ + _write_array_header(fp, d, (1, 0)) + + +def write_array_header_2_0(fp, d): + """ Write the header for an array using the 2.0 format. + The 2.0 format allows storing very large structured arrays. + + .. versionadded:: 1.9.0 + + Parameters + ---------- + fp : filelike object + d : dict + This has the appropriate entries for writing its string representation + to the header of the file. + """ + _write_array_header(fp, d, (2, 0)) def read_array_header_1_0(fp): """ @@ -317,12 +375,58 @@ def read_array_header_1_0(fp): If the data is invalid. """ + _read_array_header(fp, version=(1, 0)) + +def read_array_header_2_0(fp): + """ + Read an array header from a filelike object using the 2.0 file format + version. + + This will leave the file object located just after the header. + + .. versionadded:: 1.9.0 + + Parameters + ---------- + fp : filelike object + A file object or something with a `.read()` method like a file. + + Returns + ------- + shape : tuple of int + The shape of the array. + fortran_order : bool + The array data will be written out directly if it is either C-contiguous + or Fortran-contiguous. Otherwise, it will be made contiguous before + writing it out. + dtype : dtype + The dtype of the file's data. + + Raises + ------ + ValueError + If the data is invalid. + + """ + _read_array_header(fp, version=(2, 0)) + +def _read_array_header(fp, version): + """ + see read_array_header_1_0 + """ # Read an unsigned, little-endian short int which has the length of the # header. import struct - hlength_str = _read_bytes(fp, 2, "array header length") - header_length = struct.unpack('<H', hlength_str)[0] - header = _read_bytes(fp, header_length, "array header") + if version == (1, 0): + hlength_str = _read_bytes(fp, 2, "array header length") + header_length = struct.unpack('<H', hlength_str)[0] + header = _read_bytes(fp, header_length, "array header") + elif version == (2, 0): + hlength_str = _read_bytes(fp, 4, "array header length") + header_length = struct.unpack('<I', hlength_str)[0] + header = _read_bytes(fp, header_length, "array header") + else: + raise ValueError("Invalid version %r" % version) # The header is a pretty-printed string representation of a literal Python # dictionary with trailing newlines padded to a 16-byte boundary. The keys @@ -359,7 +463,7 @@ def read_array_header_1_0(fp): return d['shape'], d['fortran_order'], dtype -def write_array(fp, array, version=(1, 0)): +def write_array(fp, array, version=None): """ Write an array to an NPY file, including a header. @@ -374,8 +478,9 @@ def write_array(fp, array, version=(1, 0)): method. array : ndarray The array to write to disk. - version : (int, int), optional - The version number of the format. Default: (1, 0) + version : (int, int) or None, optional + The version number of the format. None means use the oldest supported + version that is able to store the data. Default: None Raises ------ @@ -387,11 +492,13 @@ def write_array(fp, array, version=(1, 0)): are not picklable. """ - if version != (1, 0): - msg = "we only support format version (1,0), not %s" - raise ValueError(msg % (version,)) - fp.write(magic(*version)) - write_array_header_1_0(fp, header_data_from_array_1_0(array)) + _check_version(version) + used_ver = _write_array_header(fp, header_data_from_array_1_0(array), + version) + # this warning can be removed when 1.9 has aged enough + if version != (2, 0) and used_ver == (2, 0): + warnings.warn("Stored array in format 2.0. It can only be" + "read by NumPy >= 1.9", UserWarning) # Set buffer size to 16 MiB to hide the Python loop overhead. buffersize = max(16 * 1024 ** 2 // array.itemsize, 1) @@ -407,7 +514,7 @@ def write_array(fp, array, version=(1, 0)): for chunk in numpy.nditer( array, flags=['external_loop', 'buffered', 'zerosize_ok'], buffersize=buffersize, order='F'): - fp.write(chunk.tostring('C')) + fp.write(chunk.tobytes('C')) else: if isfileobj(fp): array.tofile(fp) @@ -415,7 +522,7 @@ def write_array(fp, array, version=(1, 0)): for chunk in numpy.nditer( array, flags=['external_loop', 'buffered', 'zerosize_ok'], buffersize=buffersize, order='C'): - fp.write(chunk.tostring('C')) + fp.write(chunk.tobytes('C')) def read_array(fp): @@ -440,10 +547,8 @@ def read_array(fp): """ version = read_magic(fp) - if version != (1, 0): - msg = "only support version (1,0) of file format, not %r" - raise ValueError(msg % (version,)) - shape, fortran_order, dtype = read_array_header_1_0(fp) + _check_version(version) + shape, fortran_order, dtype = _read_array_header(fp, version) if len(shape) == 0: count = 1 else: @@ -486,7 +591,7 @@ def read_array(fp): def open_memmap(filename, mode='r+', dtype=None, shape=None, - fortran_order=False, version=(1, 0)): + fortran_order=False, version=None): """ Open a .npy file as a memory-mapped array. @@ -513,9 +618,11 @@ def open_memmap(filename, mode='r+', dtype=None, shape=None, Whether the array should be Fortran-contiguous (True) or C-contiguous (False, the default) if we are creating a new file in "write" mode. - version : tuple of int (major, minor) + version : tuple of int (major, minor) or None If the mode is a "write" mode, then this is the version of the file - format used to create the file. Default: (1,0) + format used to create the file. + None means use the oldest supported version that is able to store the + data. Default: None Returns ------- @@ -541,9 +648,7 @@ def open_memmap(filename, mode='r+', dtype=None, shape=None, if 'w' in mode: # We are creating the file, not reading it. # Check if we ought to create the file. - if version != (1, 0): - msg = "only support version (1,0) of file format, not %r" - raise ValueError(msg % (version,)) + _check_version(version) # Ensure that the given dtype is an authentic dtype object rather than # just something that can be interpreted as a dtype object. dtype = numpy.dtype(dtype) @@ -558,8 +663,11 @@ def open_memmap(filename, mode='r+', dtype=None, shape=None, # If we got here, then it should be safe to create the file. fp = open(filename, mode+'b') try: - fp.write(magic(*version)) - write_array_header_1_0(fp, d) + used_ver = _write_array_header(fp, d, version) + # this warning can be removed when 1.9 has aged enough + if version != (2, 0) and used_ver == (2, 0): + warnings.warn("Stored array in format 2.0. It can only be" + "read by NumPy >= 1.9", UserWarning) offset = fp.tell() finally: fp.close() @@ -568,10 +676,9 @@ def open_memmap(filename, mode='r+', dtype=None, shape=None, fp = open(filename, 'rb') try: version = read_magic(fp) - if version != (1, 0): - msg = "only support version (1,0) of file format, not %r" - raise ValueError(msg % (version,)) - shape, fortran_order, dtype = read_array_header_1_0(fp) + _check_version(version) + + shape, fortran_order, dtype = _read_array_header(fp, version) if dtype.hasobject: msg = "Array can't be memory-mapped: Python objects in dtype." raise ValueError(msg) diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py index 63b191b07..f625bcb90 100644 --- a/numpy/lib/function_base.py +++ b/numpy/lib/function_base.py @@ -13,6 +13,7 @@ __all__ = [ import warnings import sys import collections +import operator import numpy as np import numpy.core.numeric as _nx @@ -20,7 +21,7 @@ from numpy.core import linspace, atleast_1d, atleast_2d from numpy.core.numeric import ( ones, zeros, arange, concatenate, array, asarray, asanyarray, empty, empty_like, ndarray, around, floor, ceil, take, ScalarType, dot, where, - newaxis, intp, integer, isscalar + intp, integer, isscalar ) from numpy.core.umath import ( pi, multiply, add, arctan2, frompyfunc, cos, less_equal, sqrt, sin, @@ -256,22 +257,22 @@ def histogramdd(sample, bins=10, range=None, normed=False, weights=None): range : sequence, optional A sequence of lower and upper bin edges to be used if the edges are - not given explicitely in `bins`. Defaults to the minimum and maximum + not given explicitly in `bins`. Defaults to the minimum and maximum values along each dimension. normed : bool, optional - If False, returns the number of samples in each bin. If True, returns - the bin density, ie, the bin count divided by the bin hypervolume. + If False, returns the number of samples in each bin. If True, + returns the bin density ``bin_count / sample_count / bin_volume``. weights : array_like (N,), optional An array of values `w_i` weighing each sample `(x_i, y_i, z_i, ...)`. - Weights are normalized to 1 if normed is True. If normed is False, the - values of the returned histogram are equal to the sum of the weights - belonging to the samples falling into each bin. + Weights are normalized to 1 if normed is True. If normed is False, + the values of the returned histogram are equal to the sum of the + weights belonging to the samples falling into each bin. Returns ------- H : ndarray - The multidimensional histogram of sample x. See normed and weights for - the different possible semantics. + The multidimensional histogram of sample x. See normed and weights + for the different possible semantics. edges : list A list of D arrays describing the bin edges for each dimension. @@ -335,6 +336,11 @@ def histogramdd(sample, bins=10, range=None, normed=False, weights=None): smin[i] = smin[i] - .5 smax[i] = smax[i] + .5 + # 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 # Create edge arrays for i in arange(D): if isscalar(bins[i]): @@ -343,9 +349,9 @@ def histogramdd(sample, bins=10, range=None, normed=False, weights=None): "Element at index %s in `bins` should be a positive " "integer." % i) nbin[i] = bins[i] + 2 # +2 for outlier bins - edges[i] = linspace(smin[i], smax[i], nbin[i]-1) + edges[i] = linspace(smin[i], smax[i], nbin[i]-1, dtype=edge_dt) else: - edges[i] = asarray(bins[i], float) + edges[i] = asarray(bins[i], edge_dt) nbin[i] = len(edges[i]) + 1 # +1 for outlier bins dedges[i] = diff(edges[i]) if np.any(np.asarray(dedges[i]) <= 0): @@ -373,10 +379,10 @@ def histogramdd(sample, bins=10, range=None, normed=False, weights=None): if not np.isinf(mindiff): decimal = int(-log10(mindiff)) + 6 # Find which points are on the rightmost edge. - on_edge = where(around(sample[:, i], decimal) == - around(edges[i][-1], decimal))[0] + not_smaller_than_edge = (sample[:, i] >= edges[i][-1]) + on_edge = (around(sample[:, i], decimal) == around(edges[i][-1], decimal)) # Shift these points one bin to the left. - Ncount[i][on_edge] -= 1 + Ncount[i][where(on_edge & not_smaller_than_edge)[0]] -= 1 # Flattened histogram matrix (1D) # Reshape is used so that overlarge arrays @@ -650,7 +656,7 @@ def piecewise(x, condlist, funclist, *args, **kw): The output is the same shape and type as x and is found by calling the functions in `funclist` on the appropriate portions of `x`, as defined by the boolean arrays in `condlist`. Portions not covered - by any condition have undefined values. + by any condition have a default value of 0. See Also @@ -692,32 +698,24 @@ def piecewise(x, condlist, funclist, *args, **kw): if (isscalar(condlist) or not (isinstance(condlist[0], list) or isinstance(condlist[0], ndarray))): condlist = [condlist] - condlist = [asarray(c, dtype=bool) for c in condlist] + condlist = array(condlist, dtype=bool) n = len(condlist) - if n == n2 - 1: # compute the "otherwise" condition. - totlist = condlist[0] - for k in range(1, n): - totlist |= condlist[k] - condlist.append(~totlist) - n += 1 - if (n != n2): - raise ValueError( - "function list and condition list must be the same") - zerod = False # This is a hack to work around problems with NumPy's # handling of 0-d arrays and boolean indexing with # numpy.bool_ scalars + zerod = False if x.ndim == 0: x = x[None] zerod = True - newcondlist = [] - for k in range(n): - if condlist[k].ndim == 0: - condition = condlist[k][None] - else: - condition = condlist[k] - newcondlist.append(condition) - condlist = newcondlist + if condlist.shape[-1] != 1: + condlist = condlist.T + if n == n2 - 1: # compute the "otherwise" condition. + totlist = np.logical_or.reduce(condlist, axis=0) + condlist = np.vstack([condlist, ~totlist]) + n += 1 + if (n != n2): + raise ValueError( + "function list and condition list must be the same") y = zeros(x.shape, x.dtype) for k in range(n): @@ -770,29 +768,68 @@ def select(condlist, choicelist, default=0): array([ 0, 1, 2, 0, 0, 0, 36, 49, 64, 81]) """ - n = len(condlist) - n2 = len(choicelist) - if n2 != n: + # Check the size of condlist and choicelist are the same, or abort. + if len(condlist) != len(choicelist): raise ValueError( - "list of cases must be same length as list of conditions") - choicelist = [default] + choicelist - S = 0 - pfac = 1 - for k in range(1, n+1): - S += k * pfac * asarray(condlist[k-1]) - if k < n: - pfac *= (1-asarray(condlist[k-1])) - # handle special case of a 1-element condition but - # a multi-element choice - if type(S) in ScalarType or max(asarray(S).shape) == 1: - pfac = asarray(1) - for k in range(n2+1): - pfac = pfac + asarray(choicelist[k]) - if type(S) in ScalarType: - S = S*ones(asarray(pfac).shape, type(S)) - else: - S = S*ones(asarray(pfac).shape, S.dtype) - return choose(S, tuple(choicelist)) + 'list of cases must be same length as list of conditions') + + # Now that the dtype is known, handle the deprecated select([], []) case + if len(condlist) == 0: + warnings.warn("select with an empty condition list is not possible" + "and will be deprecated", + DeprecationWarning) + return np.asarray(default)[()] + + choicelist = [np.asarray(choice) for choice in choicelist] + choicelist.append(np.asarray(default)) + + # need to get the result type before broadcasting for correct scalar + # behaviour + dtype = np.result_type(*choicelist) + + # Convert conditions to arrays and broadcast conditions and choices + # as the shape is needed for the result. Doing it seperatly optimizes + # for example when all choices are scalars. + condlist = np.broadcast_arrays(*condlist) + choicelist = np.broadcast_arrays(*choicelist) + + # If cond array is not an ndarray in boolean format or scalar bool, abort. + deprecated_ints = False + for i in range(len(condlist)): + cond = condlist[i] + if cond.dtype.type is not np.bool_: + if np.issubdtype(cond.dtype, np.integer): + # A previous implementation accepted int ndarrays accidentally. + # Supported here deliberately, but deprecated. + condlist[i] = condlist[i].astype(bool) + deprecated_ints = True + else: + raise ValueError( + 'invalid entry in choicelist: should be boolean ndarray') + + if deprecated_ints: + msg = "select condlists containing integer ndarrays is deprecated " \ + "and will be removed in the future. Use `.astype(bool)` to " \ + "convert to bools." + warnings.warn(msg, DeprecationWarning) + + if choicelist[0].ndim == 0: + # This may be common, so avoid the call. + result_shape = condlist[0].shape + else: + result_shape = np.broadcast_arrays(condlist[0], choicelist[0])[0].shape + + result = np.full(result_shape, choicelist[-1], dtype) + + # Use np.copyto to burn each choicelist array onto result, using the + # corresponding condlist as a boolean mask. This is done in reverse + # order since the first choice should take precedence. + choicelist = choicelist[-2::-1] + condlist = condlist[::-1] + for choice, cond in zip(choicelist, condlist): + np.copyto(result, choice, where=cond) + + return result def copy(a, order='K'): @@ -1101,7 +1138,7 @@ def interp(x, xp, fp, left=None, right=None): ----- Does not check that the x-coordinate sequence `xp` is increasing. If `xp` is not increasing, the results are nonsense. - A simple check for increasingness is:: + A simple check for increasing is:: np.all(np.diff(xp) > 0) @@ -1578,20 +1615,22 @@ class vectorize(object): The `vectorize` function is provided primarily for convenience, not for performance. The implementation is essentially a for loop. - If `otypes` is not specified, then a call to the function with the first - argument will be used to determine the number of outputs. The results of - this call will be cached if `cache` is `True` to prevent calling the - function twice. However, to implement the cache, the original function - must be wrapped which will slow down subsequent calls, so only do this if - your function is expensive. + If `otypes` is not specified, then a call to the function with the + first argument will be used to determine the number of outputs. The + results of this call will be cached if `cache` is `True` to prevent + calling the function twice. However, to implement the cache, the + original function must be wrapped which will slow down subsequent + calls, so only do this if your function is expensive. + + The new keyword argument interface and `excluded` argument support + further degrades performance. - The new keyword argument interface and `excluded` argument support further - degrades performance. """ def __init__(self, pyfunc, otypes='', doc=None, excluded=None, cache=False): self.pyfunc = pyfunc self.cache = cache + self._ufunc = None # Caching to improve default performance if doc is None: self.__doc__ = pyfunc.__doc__ @@ -1615,9 +1654,6 @@ class vectorize(object): excluded = set() self.excluded = set(excluded) - if self.otypes and not self.excluded: - self._ufunc = None # Caching to improve default performance - def __call__(self, *args, **kwargs): """ Return arrays with the results of `pyfunc` broadcast (vectorized) over @@ -1651,7 +1687,8 @@ class vectorize(object): def _get_ufunc_and_otypes(self, func, args): """Return (ufunc, otypes).""" # frompyfunc will fail if args is empty - assert args + if not args: + raise ValueError('args can not be empty') if self.otypes: otypes = self.otypes @@ -1810,9 +1847,11 @@ def cov(m, y=None, rowvar=1, bias=0, ddof=None): "ddof must be integer") # Handles complex arrays too + m = np.asarray(m) if y is None: dtype = np.result_type(m, np.float64) else: + y = np.asarray(y) dtype = np.result_type(m, y, np.float64) X = array(m, ndmin=2, dtype=dtype) @@ -1821,11 +1860,9 @@ def cov(m, y=None, rowvar=1, bias=0, ddof=None): if rowvar: N = X.shape[1] axis = 0 - tup = (slice(None), newaxis) else: N = X.shape[0] axis = 1 - tup = (newaxis, slice(None)) # check ddof if ddof is None: @@ -1842,7 +1879,7 @@ def cov(m, y=None, rowvar=1, bias=0, ddof=None): y = array(y, copy=False, ndmin=2, dtype=dtype) X = concatenate((X, y), axis) - X -= X.mean(axis=1-axis)[tup] + X -= X.mean(axis=1-axis, keepdims=True) if not rowvar: return (dot(X.T, X.conj()) / fact).squeeze() else: @@ -1909,7 +1946,7 @@ def blackman(M): """ Return the Blackman window. - The Blackman window is a taper formed by using the the first three + The Blackman window is a taper formed by using the first three terms of a summation of cosines. It was designed to have close to the minimal leakage possible. It is close to optimal, only slightly worse than a Kaiser window. @@ -2139,9 +2176,10 @@ def hanning(M): .. math:: w(n) = 0.5 - 0.5cos\\left(\\frac{2\\pi{n}}{M-1}\\right) \\qquad 0 \\leq n \\leq M-1 - The Hanning was named for Julius van Hann, an Austrian meterologist. It is - also known as the Cosine Bell. Some authors prefer that it be called a - Hann window, to help avoid confusion with the very similar Hamming window. + The Hanning was named for Julius van Hann, an Austrian meteorologist. + It is also known as the Cosine Bell. Some authors prefer that it be + called a Hann window, to help avoid confusion with the very similar + Hamming window. Most references to the Hanning window come from the signal processing literature, where it is used as one of many windowing functions for @@ -2238,9 +2276,9 @@ def hamming(M): .. math:: w(n) = 0.54 - 0.46cos\\left(\\frac{2\\pi{n}}{M-1}\\right) \\qquad 0 \\leq n \\leq M-1 - The Hamming was named for R. W. Hamming, an associate of J. W. Tukey and - is described in Blackman and Tukey. It was recommended for smoothing the - truncated autocovariance function in the time domain. + The Hamming was named for R. W. Hamming, an associate of J. W. Tukey + and is described in Blackman and Tukey. It was recommended for + smoothing the truncated autocovariance function in the time domain. Most references to the Hamming window come from the signal processing literature, where it is used as one of many windowing functions for smoothing values. It is also known as an apodization (which means @@ -2420,11 +2458,11 @@ def i0(x): Notes ----- We use the algorithm published by Clenshaw [1]_ and referenced by - Abramowitz and Stegun [2]_, for which the function domain is partitioned - into the two intervals [0,8] and (8,inf), and Chebyshev polynomial - expansions are employed in each interval. Relative error on the domain - [0,30] using IEEE arithmetic is documented [3]_ as having a peak of 5.8e-16 - with an rms of 1.4e-16 (n = 30000). + Abramowitz and Stegun [2]_, for which the function domain is + partitioned into the two intervals [0,8] and (8,inf), and Chebyshev + polynomial expansions are employed in each interval. Relative error on + the domain [0,30] using IEEE arithmetic is documented [3]_ as having a + peak of 5.8e-16 with an rms of 1.4e-16 (n = 30000). References ---------- @@ -2494,12 +2532,11 @@ def kaiser(M, beta): where :math:`I_0` is the modified zeroth-order Bessel function. - The Kaiser was named for Jim Kaiser, who discovered a simple approximation - to the DPSS window based on Bessel functions. - The Kaiser window is a very good approximation to the Digital Prolate - Spheroidal Sequence, or Slepian window, which is the transform which - maximizes the energy in the main lobe of the window relative to total - energy. + The Kaiser was named for Jim Kaiser, who discovered a simple + approximation to the DPSS window based on Bessel functions. The Kaiser + window is a very good approximation to the Digital Prolate Spheroidal + Sequence, or Slepian window, which is the transform which maximizes the + energy in the main lobe of the window relative to total energy. The Kaiser can approximate many other windows by varying the beta parameter. @@ -2609,8 +2646,8 @@ def sinc(x): The name sinc is short for "sine cardinal" or "sinus cardinalis". The sinc function is used in various signal processing applications, - including in anti-aliasing, in the construction of a - Lanczos resampling filter, and in interpolation. + including in anti-aliasing, in the construction of a Lanczos resampling + filter, and in interpolation. For bandlimited interpolation of discrete-time signals, the ideal interpolation kernel is proportional to the sinc function. @@ -2692,7 +2729,67 @@ def msort(a): return b -def median(a, axis=None, out=None, overwrite_input=False): +def _ureduce(a, func, **kwargs): + """ + Internal Function. + Call `func` with `a` as first argument swapping the axes to use extended + axis on functions that don't support it natively. + + Returns result and a.shape with axis dims set to 1. + + Parameters + ---------- + a : array_like + Input array or object that can be converted to an array. + func : callable + Reduction function Kapable of receiving an axis argument. + It is is called with `a` as first argument followed by `kwargs`. + kwargs : keyword arguments + additional keyword arguments to pass to `func`. + + Returns + ------- + result : tuple + Result of func(a, **kwargs) and a.shape with axis dims set to 1 + which can be used to reshape the result to the same shape a ufunc with + keepdims=True would produce. + + """ + a = np.asanyarray(a) + axis = kwargs.get('axis', None) + if axis is not None: + keepdim = list(a.shape) + nd = a.ndim + try: + axis = operator.index(axis) + if axis >= nd or axis < -nd: + raise IndexError("axis %d out of bounds (%d)" % (axis, a.ndim)) + keepdim[axis] = 1 + except TypeError: + sax = set() + for x in axis: + if x >= nd or x < -nd: + raise IndexError("axis %d out of bounds (%d)" % (x, nd)) + if x in sax: + raise ValueError("duplicate value in axis") + sax.add(x % nd) + keepdim[x] = 1 + keep = sax.symmetric_difference(frozenset(range(nd))) + nkeep = len(keep) + # swap axis that should not be reduced to front + for i, s in enumerate(sorted(keep)): + a = a.swapaxes(i, s) + # merge reduced axis + a = a.reshape(a.shape[:nkeep] + (-1,)) + kwargs['axis'] = -1 + else: + keepdim = [1] * a.ndim + + r = func(a, **kwargs) + return r, keepdim + + +def median(a, axis=None, out=None, overwrite_input=False, keepdims=False): """ Compute the median along the specified axis. @@ -2702,27 +2799,35 @@ def median(a, axis=None, out=None, overwrite_input=False): ---------- a : array_like Input array or object that can be converted to an array. - axis : int, optional + axis : int or sequence of int, optional Axis along which the medians are computed. The default (axis=None) is to compute the median along a flattened version of the array. + A sequence of axes is supported since version 1.9.0. out : ndarray, optional - Alternative output array in which to place the result. It must - have the same shape and buffer length as the expected output, - but the type (of the output) will be cast if necessary. + Alternative output array in which to place the result. It must have + the same shape and buffer length as the expected output, but the + type (of the output) will be cast if necessary. overwrite_input : bool, optional If True, then allow use of memory of input array (a) for calculations. The input array will be modified by the call to - median. This will save memory when you do not need to preserve - the contents of the input array. Treat the input as undefined, - but it will probably be fully or partially sorted. Default is - False. Note that, if `overwrite_input` is True and the input - is not already an ndarray, an error will be raised. + median. This will save memory when you do not need to preserve the + contents of the input array. Treat the input as undefined, but it + will probably be fully or partially sorted. Default is False. Note + that, if `overwrite_input` is True and the input is not already an + ndarray, an error will be raised. + keepdims : bool, optional + If this is set to True, the axes which are reduced are left + in the result as dimensions with size one. With this option, + the result will broadcast correctly against the original `arr`. + + .. versionadded:: 1.9.0 + Returns ------- median : ndarray - A new array holding the result (unless `out` is specified, in - which case that array is returned instead). If the input contains + A new array holding the result (unless `out` is specified, in which + case that array is returned instead). If the input contains integers, or floats of smaller precision than 64, then the output data-type is float64. Otherwise, the output data-type is the same as that of the input. @@ -2766,6 +2871,16 @@ def median(a, axis=None, out=None, overwrite_input=False): >>> assert not np.all(a==b) """ + r, k = _ureduce(a, func=_median, axis=axis, out=out, + overwrite_input=overwrite_input) + if keepdims: + return r.reshape(k) + else: + return r + +def _median(a, axis=None, out=None, overwrite_input=False): + # can't be reasonably be implemented in terms of percentile as we have to + # call mean to not break astropy a = np.asanyarray(a) if axis is not None and axis >= a.ndim: raise IndexError( @@ -2815,7 +2930,7 @@ def median(a, axis=None, out=None, overwrite_input=False): def percentile(a, q, axis=None, out=None, - overwrite_input=False, interpolation='linear'): + overwrite_input=False, interpolation='linear', keepdims=False): """ Compute the qth percentile of the data along the specified axis. @@ -2827,9 +2942,10 @@ def percentile(a, q, axis=None, out=None, Input array or object that can be converted to an array. q : float in range of [0,100] (or sequence of floats) Percentile to compute which must be between 0 and 100 inclusive. - axis : int, optional + axis : int or sequence of int, optional Axis along which the percentiles are computed. The default (None) is to compute the percentiles along a flattened version of the array. + A sequence of axes is supported since version 1.9.0. out : ndarray, optional Alternative output array in which to place the result. It must have the same shape and buffer length as the expected output, @@ -2855,17 +2971,23 @@ def percentile(a, q, axis=None, out=None, * midpoint: (`i` + `j`) / 2. .. versionadded:: 1.9.0 + keepdims : bool, optional + If this is set to True, the axes which are reduced are left + in the result as dimensions with size one. With this option, + the result will broadcast correctly against the original `arr`. + + .. versionadded:: 1.9.0 Returns ------- percentile : scalar or ndarray - If a single percentile `q` is given and axis=None a scalar is returned. - If multiple percentiles `q` are given an array holding the result is - returned. The results are listed in the first axis. - (If `out` is specified, in which case that array is returned instead). - If the input contains integers, or floats of smaller precision than 64, - then the output data-type is float64. Otherwise, the output data-type - is the same as that of the input. + If a single percentile `q` is given and axis=None a scalar is + returned. If multiple percentiles `q` are given an array holding + the result is returned. The results are listed in the first axis. + (If `out` is specified, in which case that array is returned + instead). If the input contains integers, or floats of smaller + precision than 64, then the output data-type is float64. Otherwise, + the output data-type is the same as that of the input. See Also -------- @@ -2873,12 +2995,12 @@ def percentile(a, q, axis=None, out=None, Notes ----- - Given a vector V of length N, the qth percentile of V is the qth ranked - value in a sorted copy of V. The values and distances of the two nearest - neighbors as well as the `interpolation` parameter will determine the - percentile if the normalized ranking does not match q exactly. This - function is the same as the median if ``q=50``, the same as the minimum - if ``q=0``and the same as the maximum if ``q=100``. + Given a vector V of length N, the q-th percentile of V is the q-th ranked + value in a sorted copy of V. The values and distances of the two + nearest neighbors as well as the `interpolation` parameter will + determine the percentile if the normalized ranking does not match q + exactly. This function is the same as the median if ``q=50``, the same + as the minimum if ``q=0``and the same as the maximum if ``q=100``. Examples -------- @@ -2911,8 +3033,22 @@ def percentile(a, q, axis=None, out=None, array([ 3.5]) """ + q = asarray(q, dtype=np.float64) + r, k = _ureduce(a, func=_percentile, q=q, axis=axis, out=out, + overwrite_input=overwrite_input, + interpolation=interpolation) + if keepdims: + if q.ndim == 0: + return r.reshape(k) + else: + return r.reshape([len(q)] + k) + else: + return r + + +def _percentile(a, q, axis=None, out=None, + overwrite_input=False, interpolation='linear', keepdims=False): a = asarray(a) - q = asarray(q) if q.ndim == 0: # Do not allow 0-d arrays because following code fails for scalar zerod = True @@ -2920,10 +3056,17 @@ def percentile(a, q, axis=None, out=None, else: zerod = False - q = q / 100.0 - if (q < 0).any() or (q > 1).any(): - raise ValueError( - "Percentiles must be in the range [0,100]") + # avoid expensive reductions, relevant for arrays with < O(1000) elements + if q.size < 10: + for i in range(q.size): + if q[i] < 0. or q[i] > 100.: + raise ValueError("Percentiles must be in the range [0,100]") + q[i] /= 100. + else: + # faster than any() + if np.count_nonzero(q < 0.) or np.count_nonzero(q > 100.): + raise ValueError("Percentiles must be in the range [0,100]") + q /= 100. # prepare a for partioning if overwrite_input: @@ -3029,10 +3172,11 @@ def trapz(y, x=None, dx=1.0, axis=-1): Notes ----- - Image [2]_ illustrates trapezoidal rule -- y-axis locations of points will - be taken from `y` array, by default x-axis distances between points will be - 1.0, alternatively they can be provided with `x` array or with `dx` scalar. - Return value will be equal to combined area under the red lines. + Image [2]_ illustrates trapezoidal rule -- y-axis locations of points + will be taken from `y` array, by default x-axis distances between + points will be 1.0, alternatively they can be provided with `x` array + or with `dx` scalar. Return value will be equal to combined area under + the red lines. References @@ -3130,7 +3274,7 @@ def meshgrid(*xi, **kwargs): Make N-D coordinate arrays for vectorized evaluations of N-D scalar/vector fields over N-D grids, given one-dimensional coordinate arrays x1, x2,..., xn. - + .. versionchanged:: 1.9 1-D and 0-D cases are allowed. @@ -3141,16 +3285,22 @@ def meshgrid(*xi, **kwargs): indexing : {'xy', 'ij'}, optional Cartesian ('xy', default) or matrix ('ij') indexing of output. See Notes for more details. + + .. versionadded:: 1.7.0 sparse : bool, optional - If True a sparse grid is returned in order to conserve memory. - Default is False. + If True a sparse grid is returned in order to conserve memory. + Default is False. + + .. versionadded:: 1.7.0 copy : bool, optional - If False, a view into the original arrays are returned in - order to conserve memory. Default is True. Please note that - ``sparse=False, copy=False`` will likely return non-contiguous arrays. - Furthermore, more than one element of a broadcast array may refer to - a single memory location. If you need to write to the arrays, make - copies first. + If False, a view into the original arrays are returned in order to + conserve memory. Default is True. Please note that + ``sparse=False, copy=False`` will likely return non-contiguous + arrays. Furthermore, more than one element of a broadcast array + may refer to a single memory location. If you need to write to the + arrays, make copies first. + + .. versionadded:: 1.7.0 Returns ------- @@ -3164,13 +3314,13 @@ def meshgrid(*xi, **kwargs): Notes ----- This function supports both indexing conventions through the indexing - keyword argument. Giving the string 'ij' returns a meshgrid with matrix - indexing, while 'xy' returns a meshgrid with Cartesian indexing. In the - 2-D case with inputs of length M and N, the outputs are of shape (N, M) for - 'xy' indexing and (M, N) for 'ij' indexing. In the 3-D case with inputs of - length M, N and P, outputs are of shape (N, M, P) for 'xy' indexing and - (M, N, P) for 'ij' indexing. The difference is illustrated by the - following code snippet:: + keyword argument. Giving the string 'ij' returns a meshgrid with + matrix indexing, while 'xy' returns a meshgrid with Cartesian indexing. + In the 2-D case with inputs of length M and N, the outputs are of shape + (N, M) for 'xy' indexing and (M, N) for 'ij' indexing. In the 3-D case + with inputs of length M, N and P, outputs are of shape (N, M, P) for + 'xy' indexing and (M, N, P) for 'ij' indexing. The difference is + illustrated by the following code snippet:: xv, yv = meshgrid(x, y, sparse=False, indexing='ij') for i in range(nx): @@ -3181,7 +3331,7 @@ def meshgrid(*xi, **kwargs): for i in range(nx): for j in range(ny): # treat xv[j,i], yv[j,i] - + In the 1-D and 0-D case, the indexing and sparse keywords have no effect. See Also @@ -3221,9 +3371,13 @@ def meshgrid(*xi, **kwargs): """ ndim = len(xi) - copy_ = kwargs.get('copy', True) - sparse = kwargs.get('sparse', False) - indexing = kwargs.get('indexing', 'xy') + copy_ = kwargs.pop('copy', True) + sparse = kwargs.pop('sparse', False) + indexing = kwargs.pop('indexing', 'xy') + + if kwargs: + raise TypeError("meshgrid() got an unexpected keyword argument '%s'" + % (list(kwargs)[0],)) if not indexing in ['xy', 'ij']: raise ValueError( @@ -3258,7 +3412,8 @@ def meshgrid(*xi, **kwargs): def delete(arr, obj, axis=None): """ Return a new array with sub-arrays along an axis deleted. For a one - dimensional array, this returns those entries not returned by `arr[obj]`. + dimensional array, this returns those entries not returned by + `arr[obj]`. Parameters ---------- @@ -3285,9 +3440,11 @@ def delete(arr, obj, axis=None): Notes ----- Often it is preferable to use a boolean mask. For example: + >>> mask = np.ones(len(arr), dtype=bool) >>> mask[[0,2,4]] = False >>> result = arr[mask,...] + Is equivalent to `np.delete(arr, [0,2,4], axis=0)`, but allows further use of `mask`. @@ -3465,7 +3622,8 @@ def insert(arr, obj, values, axis=None): .. versionadded:: 1.8.0 Support for multiple insertions when `obj` is a single scalar or a - sequence with one element (similar to calling insert multiple times). + sequence with one element (similar to calling insert multiple + times). values : array_like Values to insert into `arr`. If the type of `values` is different from that of `arr`, `values` is converted to the type of `arr`. @@ -3664,19 +3822,19 @@ def append(arr, values, axis=None): Values are appended to a copy of this array. values : array_like These values are appended to a copy of `arr`. It must be of the - correct shape (the same shape as `arr`, excluding `axis`). If `axis` - is not specified, `values` can be any shape and will be flattened - before use. + correct shape (the same shape as `arr`, excluding `axis`). If + `axis` is not specified, `values` can be any shape and will be + flattened before use. axis : int, optional - The axis along which `values` are appended. If `axis` is not given, - both `arr` and `values` are flattened before use. + The axis along which `values` are appended. If `axis` is not + given, both `arr` and `values` are flattened before use. Returns ------- append : ndarray - A copy of `arr` with `values` appended to `axis`. Note that `append` - does not occur in-place: a new array is allocated and filled. If - `axis` is None, `out` is a flattened array. + A copy of `arr` with `values` appended to `axis`. Note that + `append` does not occur in-place: a new array is allocated and + filled. If `axis` is None, `out` is a flattened array. See Also -------- diff --git a/numpy/lib/nanfunctions.py b/numpy/lib/nanfunctions.py index 5766084ab..7120760b5 100644 --- a/numpy/lib/nanfunctions.py +++ b/numpy/lib/nanfunctions.py @@ -17,12 +17,14 @@ Functions from __future__ import division, absolute_import, print_function import warnings +import operator import numpy as np - +from numpy.core.fromnumeric import partition +from numpy.lib.function_base import _ureduce as _ureduce __all__ = [ 'nansum', 'nanmax', 'nanmin', 'nanargmax', 'nanargmin', 'nanmean', - 'nanvar', 'nanstd' + 'nanmedian', 'nanpercentile', 'nanvar', 'nanstd' ] @@ -228,7 +230,7 @@ def nanmin(a, axis=None, out=None, keepdims=False): # Check for all-NaN axis mask = np.all(mask, axis=axis, keepdims=keepdims) if np.any(mask): - res = _copyto(res, mask, np.nan) + res = _copyto(res, np.nan, mask) warnings.warn("All-NaN axis encountered", RuntimeWarning) return res @@ -327,7 +329,7 @@ def nanmax(a, axis=None, out=None, keepdims=False): # Check for all-NaN axis mask = np.all(mask, axis=axis, keepdims=keepdims) if np.any(mask): - res = _copyto(res, mask, np.nan) + res = _copyto(res, np.nan, mask) warnings.warn("All-NaN axis encountered", RuntimeWarning) return res @@ -601,6 +603,335 @@ def nanmean(a, axis=None, dtype=None, out=None, keepdims=False): return avg +def _nanmedian1d(arr1d, overwrite_input=False): + """ + Private function for rank 1 arrays. Compute the median ignoring NaNs. + See nanmedian for parameter usage + """ + c = np.isnan(arr1d) + s = np.where(c)[0] + if s.size == arr1d.size: + warnings.warn("All-NaN slice encountered", RuntimeWarning) + return np.nan + elif s.size == 0: + return np.median(arr1d, overwrite_input=overwrite_input) + else: + if overwrite_input: + x = arr1d + else: + x = arr1d.copy() + # select non-nans at end of array + enonan = arr1d[-s.size:][~c[-s.size:]] + # fill nans in beginning of array with non-nans of end + x[s[:enonan.size]] = enonan + # slice nans away + return np.median(x[:-s.size], overwrite_input=True) + + +def _nanmedian(a, axis=None, out=None, overwrite_input=False): + """ + Private function that doesn't support extended axis or keepdims. + These methods are extended to this function using _ureduce + See nanmedian for parameter usage + + """ + if axis is None or a.ndim == 1: + part = a.ravel() + if out is None: + return _nanmedian1d(part, overwrite_input) + else: + out[...] = _nanmedian1d(part, overwrite_input) + return out + else: + # for small medians use sort + indexing which is still faster than + # apply_along_axis + if a.shape[axis] < 400: + return _nanmedian_small(a, axis, out, overwrite_input) + result = np.apply_along_axis(_nanmedian1d, axis, a, overwrite_input) + if out is not None: + out[...] = result + return result + +def _nanmedian_small(a, axis=None, out=None, overwrite_input=False): + """ + sort + indexing median, faster for small medians along multiple dimensions + due to the high overhead of apply_along_axis + see nanmedian for parameter usage + """ + a = np.ma.masked_array(a, np.isnan(a)) + m = np.ma.median(a, axis=axis, overwrite_input=overwrite_input) + for i in range(np.count_nonzero(m.mask.ravel())): + warnings.warn("All-NaN slice encountered", RuntimeWarning) + if out is not None: + out[...] = m.filled(np.nan) + return out + return m.filled(np.nan) + +def nanmedian(a, axis=None, out=None, overwrite_input=False, keepdims=False): + """ + Compute the median along the specified axis, while ignoring NaNs. + + Returns the median of the array elements. + + .. versionadded:: 1.9.0 + + Parameters + ---------- + a : array_like + Input array or object that can be converted to an array. + axis : int, optional + Axis along which the medians are computed. The default (axis=None) + is to compute the median along a flattened version of the array. + A sequence of axes is supported since version 1.9.0. + out : ndarray, optional + Alternative output array in which to place the result. It must have + the same shape and buffer length as the expected output, but the + type (of the output) will be cast if necessary. + overwrite_input : bool, optional + If True, then allow use of memory of input array (a) for + calculations. The input array will be modified by the call to + median. This will save memory when you do not need to preserve + the contents of the input array. Treat the input as undefined, + but it will probably be fully or partially sorted. Default is + False. Note that, if `overwrite_input` is True and the input + is not already an ndarray, an error will be raised. + keepdims : bool, optional + If this is set to True, the axes which are reduced are left + in the result as dimensions with size one. With this option, + the result will broadcast correctly against the original `arr`. + + + + Returns + ------- + median : ndarray + A new array holding the result. If the input contains integers, or + floats of smaller precision than 64, then the output data-type is + float64. Otherwise, the output data-type is the same as that of the + input. + + See Also + -------- + mean, median, percentile + + Notes + ----- + Given a vector V of length N, the median of V is the middle value of + a sorted copy of V, ``V_sorted`` - i.e., ``V_sorted[(N-1)/2]``, when N is + odd. When N is even, it is the average of the two middle values of + ``V_sorted``. + + Examples + -------- + >>> a = np.array([[10.0, 7, 4], [3, 2, 1]]) + >>> a[0, 1] = np.nan + >>> a + array([[ 10., nan, 4.], + [ 3., 2., 1.]]) + >>> np.median(a) + nan + >>> np.nanmedian(a) + 3.0 + >>> np.nanmedian(a, axis=0) + array([ 6.5, 2., 2.5]) + >>> np.median(a, axis=1) + array([ 7., 2.]) + >>> b = a.copy() + >>> np.nanmedian(b, axis=1, overwrite_input=True) + array([ 7., 2.]) + >>> assert not np.all(a==b) + >>> b = a.copy() + >>> np.nanmedian(b, axis=None, overwrite_input=True) + 3.0 + >>> assert not np.all(a==b) + + """ + a = np.asanyarray(a) + # apply_along_axis in _nanmedian doesn't handle empty arrays well, + # so deal them upfront + if a.size == 0: + return np.nanmean(a, axis, out=out, keepdims=keepdims) + + r, k = _ureduce(a, func=_nanmedian, axis=axis, out=out, + overwrite_input=overwrite_input) + if keepdims: + return r.reshape(k) + else: + return r + + +def nanpercentile(a, q, axis=None, out=None, overwrite_input=False, + interpolation='linear', keepdims=False): + """ + Compute the qth percentile of the data along the specified axis, while + ignoring nan values. + + Returns the qth percentile of the array elements. + + Parameters + ---------- + a : array_like + Input array or object that can be converted to an array. + q : float in range of [0,100] (or sequence of floats) + Percentile to compute which must be between 0 and 100 inclusive. + axis : int or sequence of int, optional + Axis along which the percentiles are computed. The default (None) + is to compute the percentiles along a flattened version of the array. + A sequence of axes is supported since version 1.9.0. + out : ndarray, optional + Alternative output array in which to place the result. It must + have the same shape and buffer length as the expected output, + but the type (of the output) will be cast if necessary. + overwrite_input : bool, optional + If True, then allow use of memory of input array `a` for + calculations. The input array will be modified by the call to + percentile. This will save memory when you do not need to preserve + the contents of the input array. In this case you should not make + any assumptions about the content of the passed in array `a` after + this function completes -- treat it as undefined. Default is False. + Note that, if the `a` input is not already an array this parameter + will have no effect, `a` will be converted to an array internally + regardless of the value of this parameter. + interpolation : {'linear', 'lower', 'higher', 'midpoint', 'nearest'} + This optional parameter specifies the interpolation method to use, + when the desired quantile lies between two data points `i` and `j`: + * linear: `i + (j - i) * fraction`, where `fraction` is the + fractional part of the index surrounded by `i` and `j`. + * lower: `i`. + * higher: `j`. + * nearest: `i` or `j` whichever is nearest. + * midpoint: (`i` + `j`) / 2. + + keepdims : bool, optional + If this is set to True, the axes which are reduced are left + in the result as dimensions with size one. With this option, + the result will broadcast correctly against the original `arr`. + + + Returns + ------- + nanpercentile : scalar or ndarray + If a single percentile `q` is given and axis=None a scalar is + returned. If multiple percentiles `q` are given an array holding + the result is returned. The results are listed in the first axis. + (If `out` is specified, in which case that array is returned + instead). If the input contains integers, or floats of smaller + precision than 64, then the output data-type is float64. Otherwise, + the output data-type is the same as that of the input. + + See Also + -------- + nanmean, nanmedian, percentile, median, mean + + Notes + ----- + Given a vector V of length N, the q-th percentile of V is the q-th ranked + value in a sorted copy of V. The values and distances of the two + nearest neighbors as well as the `interpolation` parameter will + determine the percentile if the normalized ranking does not match q + exactly. This function is the same as the median if ``q=50``, the same + as the minimum if ``q=0``and the same as the maximum if ``q=100``. + + Examples + -------- + >>> a = np.array([[10., 7., 4.], [3., 2., 1.]]) + >>> a[0][1] = np.nan + >>> a + array([[ 10., nan, 4.], + [ 3., 2., 1.]]) + >>> np.percentile(a, 50) + nan + >>> np.nanpercentile(a, 50) + 3.5 + >>> np.nanpercentile(a, 50, axis=0) + array([[ 6.5, 4.5, 2.5]]) + >>> np.nanpercentile(a, 50, axis=1) + array([[ 7.], + [ 2.]]) + >>> m = np.nanpercentile(a, 50, axis=0) + >>> out = np.zeros_like(m) + >>> np.nanpercentile(a, 50, axis=0, out=m) + array([[ 6.5, 4.5, 2.5]]) + >>> m + array([[ 6.5, 4.5, 2.5]]) + >>> b = a.copy() + >>> np.nanpercentile(b, 50, axis=1, overwrite_input=True) + array([[ 7.], + [ 2.]]) + >>> assert not np.all(a==b) + >>> b = a.copy() + >>> np.nanpercentile(b, 50, axis=None, overwrite_input=True) + array([ 3.5]) + + """ + + a = np.asanyarray(a) + q = np.asanyarray(q) + # apply_along_axis in _nanpercentile doesn't handle empty arrays well, + # so deal them upfront + if a.size == 0: + return np.nanmean(a, axis, out=out, keepdims=keepdims) + + r, k = _ureduce(a, func=_nanpercentile, q=q, axis=axis, out=out, + overwrite_input=overwrite_input, + interpolation=interpolation) + if keepdims: + if q.ndim == 0: + return r.reshape(k) + else: + return r.reshape([len(q)] + k) + else: + return r + + +def _nanpercentile(a, q, axis=None, out=None, overwrite_input=False, + interpolation='linear', keepdims=False): + """ + Private function that doesn't support extended axis or keepdims. + These methods are extended to this function using _ureduce + See nanpercentile for parameter usage + + """ + if axis is None: + part = a.ravel() + result = _nanpercentile1d(part, q, overwrite_input, interpolation) + else: + result = np.apply_along_axis(_nanpercentile1d, axis, a, q, + overwrite_input, interpolation) + + if out is not None: + out[...] = result + return result + + +def _nanpercentile1d(arr1d, q, overwrite_input=False, interpolation='linear'): + """ + Private function for rank 1 arrays. Compute percentile ignoring NaNs. + See nanpercentile for parameter usage + + """ + c = np.isnan(arr1d) + s = np.where(c)[0] + if s.size == arr1d.size: + warnings.warn("All-NaN slice encountered", RuntimeWarning) + return np.nan + elif s.size == 0: + return np.percentile(arr1d, q, overwrite_input=overwrite_input, + interpolation=interpolation) + else: + if overwrite_input: + x = arr1d + else: + x = arr1d.copy() + # select non-nans at end of array + enonan = arr1d[-s.size:][~c[-s.size:]] + # fill nans in beginning of array with non-nans of end + x[s[:enonan.size]] = enonan + # slice nans away + return np.percentile(x[:-s.size], q, overwrite_input=True, + interpolation=interpolation) + + def nanvar(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False): """ Compute the variance along the specified axis, while ignoring NaNs. diff --git a/numpy/lib/npyio.py b/numpy/lib/npyio.py index af259cfef..42a539f78 100644 --- a/numpy/lib/npyio.py +++ b/numpy/lib/npyio.py @@ -285,21 +285,22 @@ def load(file, mmap_mode=None): Parameters ---------- file : file-like object or string - The file to read. It must support ``seek()`` and ``read()`` methods. - If the filename extension is ``.gz``, the file is first decompressed. + The file to read. File-like objects must support the + ``seek()`` and ``read()`` methods. Pickled files require that the + file-like object support the ``readline()`` method as well. mmap_mode : {None, 'r+', 'r', 'w+', 'c'}, optional - If not None, then memory-map the file, using the given mode - (see `numpy.memmap` for a detailed description of the modes). - A memory-mapped array is kept on disk. However, it can be accessed - and sliced like any ndarray. Memory mapping is especially useful for - accessing small fragments of large files without reading the entire - file into memory. + If not None, then memory-map the file, using the given mode (see + `numpy.memmap` for a detailed description of the modes). A + memory-mapped array is kept on disk. However, it can be accessed + and sliced like any ndarray. Memory mapping is especially useful + for accessing small fragments of large files without reading the + entire file into memory. Returns ------- result : array, tuple, dict, etc. - Data stored in the file. For ``.npz`` files, the returned instance of - NpzFile class must be closed to avoid leaking file descriptors. + Data stored in the file. For ``.npz`` files, the returned instance + of NpzFile class must be closed to avoid leaking file descriptors. Raises ------ @@ -308,7 +309,7 @@ def load(file, mmap_mode=None): See Also -------- - save, savez, loadtxt + save, savez, savez_compressed, loadtxt memmap : Create a memory-map to an array stored in a file on disk. Notes @@ -319,13 +320,14 @@ def load(file, mmap_mode=None): - If the file is a ``.npz`` file, then a dictionary-like object is returned, containing ``{filename: array}`` key-value pairs, one for each file in the archive. - - If the file is a ``.npz`` file, the returned value supports the context - manager protocol in a similar fashion to the open function:: + - If the file is a ``.npz`` file, the returned value supports the + context manager protocol in a similar fashion to the open function:: with load('foo.npz') as data: a = data['a'] - The underlying file descriptor is closed when exiting the 'with' block. + The underlying file descriptor is closed when exiting the 'with' + block. Examples -------- @@ -549,6 +551,7 @@ def savez_compressed(file, *args, **kwds): See Also -------- numpy.savez : Save several arrays into an uncompressed ``.npz`` file format + numpy.load : Load the files created by savez_compressed. """ _savez(file, args, kwds, True) @@ -668,6 +671,7 @@ def loadtxt(fname, dtype=float, comments='#', delimiter=None, The returned array will have at least `ndmin` dimensions. Otherwise mono-dimensional axes will be squeezed. Legal values: 0 (default), 1 or 2. + .. versionadded:: 1.6.0 Returns @@ -840,6 +844,11 @@ def loadtxt(fname, dtype=float, comments='#', delimiter=None, continue if usecols: vals = [vals[i] for i in usecols] + if len(vals) != N: + line_num = i + skiprows + 1 + raise ValueError("Wrong number of columns at line %d" + % line_num) + # Convert each value according to its column and store items = [conv(val) for (conv, val) in zip(converters, vals)] # Then pack it according to the dtype's nesting @@ -906,22 +915,26 @@ def savetxt(fname, X, fmt='%.18e', delimiter=' ', newline='\n', header='', and imaginary part must have separate specifiers, e.g. `['%.3e + %.3ej', '(%.15e%+.15ej)']` for 2 columns delimiter : str, optional - Character separating columns. + String or character separating columns. newline : str, optional + String or character separating lines. + .. versionadded:: 1.5.0 header : str, optional String that will be written at the beginning of the file. + .. versionadded:: 1.7.0 footer : str, optional String that will be written at the end of the file. + .. versionadded:: 1.7.0 comments : str, optional String that will be prepended to the ``header`` and ``footer`` strings, to mark them as comments. Default: '# ', as expected by e.g. ``numpy.loadtxt``. + .. versionadded:: 1.7.0 - Character separating lines. See Also -------- @@ -1191,14 +1204,20 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None, The string used to separate values. By default, any consecutive whitespaces act as delimiter. An integer or sequence of integers can also be provided as width(s) of each field. + skip_rows : int, optional + `skip_rows` was deprecated in numpy 1.5, and will be removed in + numpy 2.0. Please use `skip_header` instead. skip_header : int, optional - The numbers of lines to skip at the beginning of the file. + The number of lines to skip at the beginning of the file. skip_footer : int, optional - The numbers of lines to skip at the end of the file + The number of lines to skip at the end of the file. converters : variable, optional The set of functions that convert the data of a column to a value. The converters can also be used to provide a default value for missing data: ``converters = {3: lambda s: float(s or 0)}``. + missing : variable, optional + `missing` was deprecated in numpy 1.5, and will be removed in + numpy 2.0. Please use `missing_values` instead. missing_values : variable, optional The set of strings corresponding to missing data. filling_values : variable, optional @@ -1236,6 +1255,8 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None, usemask : bool, optional If True, return a masked array. If False, return a regular array. + loose : bool, optional + If True, do not raise errors for invalid values. invalid_raise : bool, optional If True, an exception is raised if an inconsistency is detected in the number of columns. @@ -1840,7 +1861,7 @@ def recfromtxt(fname, **kwargs): array will be determined from the data. """ - kwargs.update(dtype=kwargs.get('dtype', None)) + kwargs.setdefault("dtype", None) usemask = kwargs.get('usemask', False) output = genfromtxt(fname, **kwargs) if usemask: @@ -1867,17 +1888,20 @@ def recfromcsv(fname, **kwargs): -------- numpy.genfromtxt : generic function to load ASCII data. + Notes + ----- + By default, `dtype` is None, which means that the data-type of the output + array will be determined from the data. + """ - case_sensitive = kwargs.get('case_sensitive', "lower") or "lower" - names = kwargs.get('names', True) - if names is None: - names = True - kwargs.update(dtype=kwargs.get('update', None), - delimiter=kwargs.get('delimiter', ",") or ",", - names=names, - case_sensitive=case_sensitive) - usemask = kwargs.get("usemask", False) + # Set default kwargs for genfromtxt as relevant to csv import. + kwargs.setdefault("case_sensitive", "lower") + kwargs.setdefault("names", True) + kwargs.setdefault("delimiter", ",") + kwargs.setdefault("dtype", None) output = genfromtxt(fname, **kwargs) + + usemask = kwargs.get("usemask", False) if usemask: from numpy.ma.mrecords import MaskedRecords output = output.view(MaskedRecords) diff --git a/numpy/lib/polynomial.py b/numpy/lib/polynomial.py index 8545375cd..10ae32a60 100644 --- a/numpy/lib/polynomial.py +++ b/numpy/lib/polynomial.py @@ -127,7 +127,7 @@ def poly(seq_of_zeros): elif len(sh) == 1: pass else: - raise ValueError("input must be 1d or square 2d array.") + raise ValueError("input must be 1d or non-empty square 2d array.") if len(seq_of_zeros) == 0: return 1.0 @@ -806,6 +806,8 @@ def polymul(a1, a2): poly1d : A one-dimensional polynomial class. poly, polyadd, polyder, polydiv, polyfit, polyint, polysub, polyval + convolve : Array convolution. Same output as polymul, but has parameter + for overlap mode. Examples -------- @@ -1191,10 +1193,12 @@ class poly1d(object): __rtruediv__ = __rdiv__ def __eq__(self, other): - return NX.alltrue(self.coeffs == other.coeffs) + if self.coeffs.shape != other.coeffs.shape: + return False + return (self.coeffs == other.coeffs).all() def __ne__(self, other): - return NX.any(self.coeffs != other.coeffs) + return not self.__eq__(other) def __setattr__(self, key, val): raise ValueError("Attributes cannot be changed this way.") diff --git a/numpy/lib/shape_base.py b/numpy/lib/shape_base.py index 4fdaba349..a6d391728 100644 --- a/numpy/lib/shape_base.py +++ b/numpy/lib/shape_base.py @@ -12,7 +12,7 @@ from numpy.core.numeric import asarray, zeros, newaxis, outer, \ from numpy.core.fromnumeric import product, reshape from numpy.core import hstack, vstack, atleast_3d -def apply_along_axis(func1d,axis,arr,*args): +def apply_along_axis(func1d, axis, arr, *args, **kwargs): """ Apply a function to 1-D slices along the given axis. @@ -30,6 +30,11 @@ def apply_along_axis(func1d,axis,arr,*args): Input array. args : any Additional arguments to `func1d`. + kwargs: any + Additional named arguments to `func1d`. + + .. versionadded:: 1.9.0 + Returns ------- @@ -78,7 +83,7 @@ def apply_along_axis(func1d,axis,arr,*args): i[axis] = slice(None, None) outshape = asarray(arr.shape).take(indlist) i.put(indlist, ind) - res = func1d(arr[tuple(i.tolist())],*args) + res = func1d(arr[tuple(i.tolist())], *args, **kwargs) # if res is a number, then we have a smaller output array if isscalar(res): outarr = zeros(outshape, asarray(res).dtype) @@ -94,7 +99,7 @@ def apply_along_axis(func1d,axis,arr,*args): ind[n] = 0 n -= 1 i.put(indlist, ind) - res = func1d(arr[tuple(i.tolist())],*args) + res = func1d(arr[tuple(i.tolist())], *args, **kwargs) outarr[tuple(ind)] = res k += 1 return outarr @@ -115,7 +120,7 @@ def apply_along_axis(func1d,axis,arr,*args): ind[n] = 0 n -= 1 i.put(indlist, ind) - res = func1d(arr[tuple(i.tolist())],*args) + res = func1d(arr[tuple(i.tolist())], *args, **kwargs) outarr[tuple(i.tolist())] = res k += 1 return outarr @@ -153,6 +158,12 @@ def apply_over_axes(func, a, axes): apply_along_axis : Apply a function to 1-D slices of an array along the given axis. + Notes + ------ + This function is equivalent to tuple axis arguments to reorderable ufuncs + with keepdims=True. Tuple axis arguments to ufuncs have been availabe since + version 1.7.0. + Examples -------- >>> a = np.arange(24).reshape(2,3,4) @@ -172,6 +183,13 @@ def apply_over_axes(func, a, axes): [ 92], [124]]]) + Tuple axis arguments to ufuncs are equivalent: + + >>> np.sum(a, axis=(0,2), keepdims=True) + array([[[ 60], + [ 92], + [124]]]) + """ val = asarray(a) N = a.ndim @@ -274,10 +292,6 @@ def column_stack(tup): -------- hstack, vstack, concatenate - Notes - ----- - This function is equivalent to ``np.vstack(tup).T``. - Examples -------- >>> a = np.array((1,2,3)) @@ -638,7 +652,7 @@ def dsplit(ary, indices_or_sections): """ if len(_nx.shape(ary)) < 3: - raise ValueError('vsplit only works on arrays of 3 or more dimensions') + raise ValueError('dsplit only works on arrays of 3 or more dimensions') return split(ary, indices_or_sections, 2) def get_array_prepare(*args): diff --git a/numpy/lib/src/_compiled_base.c b/numpy/lib/src/_compiled_base.c index 328fc2d14..a461613e3 100644 --- a/numpy/lib/src/_compiled_base.c +++ b/numpy/lib/src/_compiled_base.c @@ -60,29 +60,44 @@ decr_slot_right_(double x, double * bins, npy_intp lbins) return 0; } -/** +/* * 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(double * a, int lena) +check_array_monotonic(const double *a, npy_int lena) { - int i; + 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; + } - if (a [0] <= a [1]) { - /* possibly monotonic increasing */ - for (i = 1; i < lena - 1; i ++) { - if (a [i] > a [i + 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 { - /* possibly monotonic decreasing */ - for (i = 1; i < lena - 1; i ++) { - if (a [i] < a [i + 1]) { + /* last > next, possibly monotonic decreasing */ + for (i += 1; i < lena; i++) { + last = next; + next = a[i]; + if (last < next) { return 0; } } @@ -90,39 +105,26 @@ check_array_monotonic(double * a, int lena) } } - - -/* find the index of the maximum element of an integer array */ -static npy_intp -mxx (npy_intp *i , npy_intp len) +/* 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 mx = 0, max = i[0]; - npy_intp j; + npy_intp min = *data; + npy_intp max = *data; - for ( j = 1; j < len; j ++ ) { - if ( i [j] > max ) { - max = i [j]; - mx = j; + while (--data_len) { + const npy_intp val = *(++data); + if (val < min) { + min = val; + } + else if (val > max) { + max = val; } } - return mx; -} -/* find the index of the minimum element of an integer array */ -static npy_intp -mnx (npy_intp *i , npy_intp len) -{ - npy_intp mn = 0, min = i [0]; - npy_intp j; - - for ( j = 1; j < len; j ++ ) - if ( i [j] < min ) - {min = i [j]; - mn = j;} - return mn; + *mn = min; + *mx = max; } - - /* * arr_bincount is registered as bincount. * @@ -142,7 +144,7 @@ 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 , mxi, mni, ans_size, minlength; + npy_intp *numbers, *ians, len , mx, mn, ans_size, minlength; int i; double *weights , *dans; static char *kwlist[] = {"list", "weights", "minlength", NULL}; @@ -159,14 +161,22 @@ arr_bincount(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds) len = PyArray_SIZE(lst); type = PyArray_DescrFromType(NPY_INTP); - /* handle empty list */ - if (len < 1) { - if (mlength == Py_None) { - minlength = 0; - } - else if (!(minlength = PyArray_PyIntAsIntp(mlength))) { + 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; } @@ -175,24 +185,14 @@ arr_bincount(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds) } numbers = (npy_intp *) PyArray_DATA(lst); - mxi = mxx(numbers, len); - mni = mnx(numbers, len); - if (numbers[mni] < 0) { + 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 = numbers [mxi] + 1; + ans_size = mx + 1; if (mlength != Py_None) { - if (!(minlength = PyArray_PyIntAsIntp(mlength))) { - goto fail; - } - if (minlength <= 0) { - /* superfluous, but may catch incorrect usage */ - PyErr_SetString(PyExc_ValueError, - "minlength must be positive"); - goto fail; - } if (ans_size < minlength) { ans_size = minlength; } @@ -680,8 +680,15 @@ arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict) slopes[i] = (dy[i + 1] - dy[i])/(dx[i + 1] - dx[i]); } for (i = 0; i < lenx; i++) { - npy_intp j = binary_search(dz[i], dx, lenxp); + 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; } @@ -692,7 +699,7 @@ arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict) dres[i] = rval; } else { - dres[i] = slopes[j]*(dz[i] - dx[j]) + dy[j]; + dres[i] = slopes[j]*(x - dx[j]) + dy[j]; } } NPY_END_ALLOW_THREADS; @@ -701,8 +708,15 @@ arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict) else { NPY_BEGIN_ALLOW_THREADS; for (i = 0; i < lenx; i++) { - npy_intp j = binary_search(dz[i], dx, lenxp); + 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; } @@ -713,8 +727,8 @@ arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict) dres[i] = rval; } else { - double slope = (dy[j + 1] - dy[j])/(dx[j + 1] - dx[j]); - dres[i] = slope*(dz[i] - dx[j]) + dy[j]; + 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; @@ -1402,6 +1416,9 @@ _packbits( void *In, npy_intp out_Nm1; int maxi, remain, nonzero, j; char *outptr,*inptr; + NPY_BEGIN_THREADS_DEF; + + NPY_BEGIN_THREADS_THRESHOLDED(out_N); outptr = Out; /* pointer to output buffer */ inptr = In; /* pointer to input buffer */ @@ -1436,6 +1453,8 @@ _packbits( void *In, *outptr = build; outptr += out_stride; } + + NPY_END_THREADS; return; } @@ -1453,6 +1472,9 @@ _unpackbits(void *In, unsigned char mask; int i, index; char *inptr, *outptr; + NPY_BEGIN_THREADS_DEF; + + NPY_BEGIN_THREADS_THRESHOLDED(in_N); outptr = Out; inptr = In; @@ -1465,6 +1487,8 @@ _unpackbits(void *In, } inptr += in_stride; } + + NPY_END_THREADS; return; } diff --git a/numpy/lib/stride_tricks.py b/numpy/lib/stride_tricks.py index d092f92a8..1ffaaee36 100644 --- a/numpy/lib/stride_tricks.py +++ b/numpy/lib/stride_tricks.py @@ -29,7 +29,8 @@ def as_strided(x, shape=None, strides=None): interface['strides'] = tuple(strides) array = np.asarray(DummyArray(interface, base=x)) # Make sure dtype is correct in case of custom dtype - array.dtype = x.dtype + if array.dtype.kind == 'V': + array.dtype = x.dtype return array def broadcast_arrays(*args): diff --git a/numpy/lib/tests/test__version.py b/numpy/lib/tests/test__version.py new file mode 100644 index 000000000..bbafe68eb --- /dev/null +++ b/numpy/lib/tests/test__version.py @@ -0,0 +1,57 @@ +"""Tests for the NumpyVersion class. + +""" +from __future__ import division, absolute_import, print_function + +from numpy.testing import assert_, run_module_suite, assert_raises +from numpy.lib import NumpyVersion + + +def test_main_versions(): + assert_(NumpyVersion('1.8.0') == '1.8.0') + for ver in ['1.9.0', '2.0.0', '1.8.1']: + assert_(NumpyVersion('1.8.0') < ver) + + for ver in ['1.7.0', '1.7.1', '0.9.9']: + assert_(NumpyVersion('1.8.0') > ver) + + +def test_version_1_point_10(): + # regression test for gh-2998. + assert_(NumpyVersion('1.9.0') < '1.10.0') + assert_(NumpyVersion('1.11.0') < '1.11.1') + assert_(NumpyVersion('1.11.0') == '1.11.0') + assert_(NumpyVersion('1.99.11') < '1.99.12') + + +def test_alpha_beta_rc(): + assert_(NumpyVersion('1.8.0rc1') == '1.8.0rc1') + for ver in ['1.8.0', '1.8.0rc2']: + assert_(NumpyVersion('1.8.0rc1') < ver) + + for ver in ['1.8.0a2', '1.8.0b3', '1.7.2rc4']: + assert_(NumpyVersion('1.8.0rc1') > ver) + + assert_(NumpyVersion('1.8.0b1') > '1.8.0a2') + + +def test_dev_version(): + assert_(NumpyVersion('1.9.0.dev-Unknown') < '1.9.0') + for ver in ['1.9.0', '1.9.0a1', '1.9.0b2', '1.9.0b2.dev-ffffffff']: + assert_(NumpyVersion('1.9.0.dev-f16acvda') < ver) + + assert_(NumpyVersion('1.9.0.dev-f16acvda') == '1.9.0.dev-11111111') + + +def test_dev_a_b_rc_mixed(): + assert_(NumpyVersion('1.9.0a2.dev-f16acvda') == '1.9.0a2.dev-11111111') + assert_(NumpyVersion('1.9.0a2.dev-6acvda54') < '1.9.0a2') + + +def test_raises(): + for ver in ['1.9', '1,9.0', '1.7.x']: + assert_raises(ValueError, NumpyVersion, ver) + + +if __name__ == "__main__": + run_module_suite() diff --git a/numpy/lib/tests/test_arraysetops.py b/numpy/lib/tests/test_arraysetops.py index 5934ca05a..271943abc 100644 --- a/numpy/lib/tests/test_arraysetops.py +++ b/numpy/lib/tests/test_arraysetops.py @@ -14,31 +14,59 @@ class TestSetOps(TestCase): def test_unique(self): - def check_all(a, b, i1, i2, dt): - msg = "check values failed for type '%s'" % dt + def check_all(a, b, i1, i2, c, dt): + base_msg = 'check {0} failed for type {1}' + + msg = base_msg.format('values', dt) v = unique(a) assert_array_equal(v, b, msg) - msg = "check indexes failed for type '%s'" % dt - v, j = unique(a, 1, 0) + msg = base_msg.format('return_index', dt) + v, j = unique(a, 1, 0, 0) assert_array_equal(v, b, msg) assert_array_equal(j, i1, msg) - msg = "check reverse indexes failed for type '%s'" % dt - v, j = unique(a, 0, 1) + msg = base_msg.format('return_inverse', dt) + v, j = unique(a, 0, 1, 0) assert_array_equal(v, b, msg) assert_array_equal(j, i2, msg) - msg = "check with all indexes failed for type '%s'" % dt - v, j1, j2 = unique(a, 1, 1) + msg = base_msg.format('return_counts', dt) + v, j = unique(a, 0, 0, 1) + assert_array_equal(v, b, msg) + assert_array_equal(j, c, msg) + + msg = base_msg.format('return_index and return_inverse', dt) + v, j1, j2 = unique(a, 1, 1, 0) + assert_array_equal(v, b, msg) + assert_array_equal(j1, i1, msg) + assert_array_equal(j2, i2, msg) + + msg = base_msg.format('return_index and return_counts', dt) + v, j1, j2 = unique(a, 1, 0, 1) + assert_array_equal(v, b, msg) + assert_array_equal(j1, i1, msg) + assert_array_equal(j2, c, msg) + + msg = base_msg.format('return_inverse and return_counts', dt) + v, j1, j2 = unique(a, 0, 1, 1) + assert_array_equal(v, b, msg) + assert_array_equal(j1, i2, msg) + assert_array_equal(j2, c, msg) + + msg = base_msg.format(('return_index, return_inverse ' + 'and return_counts'), dt) + v, j1, j2, j3 = unique(a, 1, 1, 1) assert_array_equal(v, b, msg) assert_array_equal(j1, i1, msg) assert_array_equal(j2, i2, msg) + assert_array_equal(j3, c, msg) a = [5, 7, 1, 2, 1, 5, 7]*10 b = [1, 2, 5, 7] i1 = [2, 3, 0, 1] i2 = [2, 3, 0, 1, 0, 2, 3]*10 + c = np.multiply([2, 1, 2, 2], 10) # test for numeric arrays types = [] @@ -49,7 +77,7 @@ class TestSetOps(TestCase): for dt in types: aa = np.array(a, dt) bb = np.array(b, dt) - check_all(aa, bb, i1, i2, dt) + check_all(aa, bb, i1, i2, c, dt) # test for object arrays dt = 'O' @@ -57,13 +85,27 @@ class TestSetOps(TestCase): aa[:] = a bb = np.empty(len(b), dt) bb[:] = b - check_all(aa, bb, i1, i2, dt) + check_all(aa, bb, i1, i2, c, dt) # test for structured arrays dt = [('', 'i'), ('', 'i')] aa = np.array(list(zip(a, a)), dt) bb = np.array(list(zip(b, b)), dt) - check_all(aa, bb, i1, i2, dt) + check_all(aa, bb, i1, i2, c, dt) + + # test for ticket #2799 + aa = [1.+0.j, 1- 1.j, 1] + assert_array_equal(np.unique(aa), [ 1.-1.j, 1.+0.j]) + + # test for ticket #4785 + a = [(1, 2), (1, 2), (2, 3)] + unq = [1, 2, 3] + inv = [0, 1, 0, 1, 1, 2] + a1 = unique(a) + assert_array_equal(a1, unq) + a2, a2_inv = unique(a, return_inverse=True) + assert_array_equal(a2, unq) + assert_array_equal(a2_inv, inv) def test_intersect1d(self): # unique inputs diff --git a/numpy/lib/tests/test_financial.py b/numpy/lib/tests/test_financial.py index 6b7c6ef53..41a060a3f 100644 --- a/numpy/lib/tests/test_financial.py +++ b/numpy/lib/tests/test_financial.py @@ -13,6 +13,21 @@ class TestFinancial(TestCase): v = [-150000, 15000, 25000, 35000, 45000, 60000] assert_almost_equal(np.irr(v), 0.0524, 2) + v = [-100, 0, 0, 74] + assert_almost_equal(np.irr(v), + -0.0955, 2) + v = [-100, 39, 59, 55, 20] + assert_almost_equal(np.irr(v), + 0.28095, 2) + v = [-100, 100, 0, -7] + assert_almost_equal(np.irr(v), + -0.0833, 2) + v = [-100, 100, 0, 7] + assert_almost_equal(np.irr(v), + 0.06206, 2) + v = [-5, 10.5, 1, -8, 1] + assert_almost_equal(np.irr(v), + 0.0886, 2) def test_pv(self): assert_almost_equal(np.pv(0.07, 20, 12000, 0), diff --git a/numpy/lib/tests/test_format.py b/numpy/lib/tests/test_format.py index b9be643c8..1034b5125 100644 --- a/numpy/lib/tests/test_format.py +++ b/numpy/lib/tests/test_format.py @@ -280,6 +280,7 @@ import sys import os import shutil import tempfile +import warnings from io import BytesIO import numpy as np @@ -521,19 +522,71 @@ def test_compressed_roundtrip(): assert_array_equal(arr, arr1) -def test_write_version_1_0(): +def test_version_2_0(): + f = BytesIO() + # requires more than 2 byte for header + dt = [(("%d" % i) * 100, float) for i in range(500)] + d = np.ones(1000, dtype=dt) + + format.write_array(f, d, version=(2, 0)) + with warnings.catch_warnings(record=True) as w: + warnings.filterwarnings('always', '', UserWarning) + format.write_array(f, d) + assert_(w[0].category is UserWarning) + + f.seek(0) + n = format.read_array(f) + assert_array_equal(d, n) + + # 1.0 requested but data cannot be saved this way + assert_raises(ValueError, format.write_array, f, d, (1, 0)) + + +def test_version_2_0_memmap(): + # requires more than 2 byte for header + dt = [(("%d" % i) * 100, float) for i in range(500)] + d = np.ones(1000, dtype=dt) + tf = tempfile.mktemp('', 'mmap', dir=tempdir) + + # 1.0 requested but data cannot be saved this way + assert_raises(ValueError, format.open_memmap, tf, mode='w+', dtype=d.dtype, + shape=d.shape, version=(1, 0)) + + ma = format.open_memmap(tf, mode='w+', dtype=d.dtype, + shape=d.shape, version=(2, 0)) + ma[...] = d + del ma + + with warnings.catch_warnings(record=True) as w: + warnings.filterwarnings('always', '', UserWarning) + ma = format.open_memmap(tf, mode='w+', dtype=d.dtype, + shape=d.shape, version=None) + assert_(w[0].category is UserWarning) + ma[...] = d + del ma + + ma = format.open_memmap(tf, mode='r') + assert_array_equal(ma, d) + + +def test_write_version(): f = BytesIO() arr = np.arange(1) # These should pass. format.write_array(f, arr, version=(1, 0)) format.write_array(f, arr) + format.write_array(f, arr, version=None) + format.write_array(f, arr) + + format.write_array(f, arr, version=(2, 0)) + format.write_array(f, arr) + # These should all fail. bad_versions = [ (1, 1), (0, 0), (0, 1), - (2, 0), (2, 2), (255, 255), ] @@ -620,5 +673,32 @@ def test_bad_header(): format.write_array_header_1_0(s, d) assert_raises(ValueError, format.read_array_header_1_0, s) + +def test_large_file_support(): + from nose import SkipTest + # try creating a large sparse file + with tempfile.NamedTemporaryFile() as tf: + try: + # seek past end would work too, but linux truncate somewhat + # increases the chances that we have a sparse filesystem and can + # avoid actually writing 5GB + import subprocess as sp + sp.check_call(["truncate", "-s", "5368709120", tf.name]) + except: + raise SkipTest("Could not create 5GB large file") + # write a small array to the end + f = open(tf.name, "wb") + f.seek(5368709120) + d = np.arange(5) + np.save(f, d) + f.close() + # read it back + f = open(tf.name, "rb") + f.seek(5368709120) + r = np.load(f) + f.close() + assert_array_equal(r, d) + + if __name__ == "__main__": run_module_suite() diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py index 494b512f7..ac677a308 100644 --- a/numpy/lib/tests/test_function_base.py +++ b/numpy/lib/tests/test_function_base.py @@ -6,7 +6,7 @@ import numpy as np from numpy.testing import ( run_module_suite, TestCase, assert_, assert_equal, assert_array_equal, assert_almost_equal, assert_array_almost_equal, assert_raises, - assert_allclose, assert_array_max_ulp, assert_warns + assert_allclose, assert_array_max_ulp, assert_warns, assert_raises_regex ) from numpy.random import rand from numpy.lib import * @@ -150,6 +150,13 @@ class TestAverage(TestCase): class TestSelect(TestCase): + choices = [np.array([1, 2, 3]), + np.array([4, 5, 6]), + np.array([7, 8, 9])] + conditions = [np.array([False, False, False]), + np.array([False, True, False]), + np.array([False, False, True])] + def _select(self, cond, values, default=0): output = [] for m in range(len(cond)): @@ -157,18 +164,62 @@ class TestSelect(TestCase): return output def test_basic(self): - choices = [np.array([1, 2, 3]), - np.array([4, 5, 6]), - np.array([7, 8, 9])] - conditions = [np.array([0, 0, 0]), - np.array([0, 1, 0]), - np.array([0, 0, 1])] + choices = self.choices + conditions = self.conditions assert_array_equal(select(conditions, choices, default=15), self._select(conditions, choices, default=15)) assert_equal(len(choices), 3) assert_equal(len(conditions), 3) + def test_broadcasting(self): + conditions = [np.array(True), np.array([False, True, False])] + choices = [1, np.arange(12).reshape(4, 3)] + assert_array_equal(select(conditions, choices), np.ones((4, 3))) + # default can broadcast too: + assert_equal(select([True], [0], default=[0]).shape, (1,)) + + def test_return_dtype(self): + assert_equal(select(self.conditions, self.choices, 1j).dtype, + np.complex_) + # But the conditions need to be stronger then the scalar default + # if it is scalar. + choices = [choice.astype(np.int8) for choice in self.choices] + assert_equal(select(self.conditions, choices).dtype, np.int8) + + d = np.array([1, 2, 3, np.nan, 5, 7]) + m = np.isnan(d) + assert_equal(select([m], [d]), [0, 0, 0, np.nan, 0, 0]) + + def test_deprecated_empty(self): + with warnings.catch_warnings(record=True): + warnings.simplefilter("always") + assert_equal(select([], [], 3j), 3j) + + with warnings.catch_warnings(): + warnings.simplefilter("always") + assert_warns(DeprecationWarning, select, [], []) + warnings.simplefilter("error") + assert_raises(DeprecationWarning, select, [], []) + + def test_non_bool_deprecation(self): + choices = self.choices + conditions = self.conditions[:] + with warnings.catch_warnings(): + warnings.filterwarnings("always") + conditions[0] = conditions[0].astype(np.int_) + assert_warns(DeprecationWarning, select, conditions, choices) + conditions[0] = conditions[0].astype(np.uint8) + assert_warns(DeprecationWarning, select, conditions, choices) + warnings.filterwarnings("error") + assert_raises(DeprecationWarning, select, conditions, choices) + + def test_many_arguments(self): + # This used to be limited by NPY_MAXARGS == 32 + conditions = [np.array([False])] * 100 + choices = [np.array([1])] * 100 + select(conditions, choices) + class TestInsert(TestCase): def test_basic(self): @@ -249,7 +300,7 @@ class TestInsert(TestCase): assert_(isinstance(np.insert(a, [], []), SubClass)) assert_(isinstance(np.insert(a, [0, 1], [1, 2]), SubClass)) assert_(isinstance(np.insert(a, slice(1, 2), [1, 2]), SubClass)) - assert_(isinstance(np.insert(a, slice(1, -2), []), SubClass)) + assert_(isinstance(np.insert(a, slice(1, -2, -1), []), SubClass)) # This is an error in the future: a = np.array(1).view(SubClass) assert_(isinstance(np.insert(a, 0, [0]), SubClass)) @@ -721,6 +772,12 @@ class TestVectorize(TestCase): assert_array_equal(f(x), x*x) assert_equal(_calls[0], len(x)) + def test_otypes(self): + f = np.vectorize(lambda x: x) + f.otypes = 'i' + x = np.arange(5) + assert_array_equal(f(x), x) + class TestDigitize(TestCase): def test_forward(self): @@ -761,6 +818,22 @@ class TestDigitize(TestCase): bins = np.linspace(x.min(), x.max(), 10) assert_(np.all(digitize(x, bins, True) != 10)) + def test_monotonic(self): + x = [-1, 0, 1, 2] + bins = [0, 0, 1] + assert_array_equal(digitize(x, bins, False), [0, 2, 3, 3]) + assert_array_equal(digitize(x, bins, True), [0, 0, 2, 3]) + bins = [1, 1, 0] + assert_array_equal(digitize(x, bins, False), [3, 2, 0, 0]) + assert_array_equal(digitize(x, bins, True), [3, 3, 2, 0]) + bins = [1, 1, 1, 1] + assert_array_equal(digitize(x, bins, False), [0, 0, 4, 4]) + assert_array_equal(digitize(x, bins, True), [0, 0, 0, 4]) + bins = [0, 0, 1, 0] + assert_raises(ValueError, digitize, x, bins) + bins = [1, 1, 0, 1] + assert_raises(ValueError, digitize, x, bins) + class TestUnwrap(TestCase): def test_simple(self): @@ -997,6 +1070,13 @@ class TestHistogram(TestCase): h, b = histogram(a, weights=np.ones(10, float)) assert_(issubdtype(h.dtype, float)) + def test_f32_rounding(self): + # gh-4799, check that the rounding of the edges works with float32 + x = np.array([276.318359 , -69.593948 , 21.329449], dtype=np.float32) + y = np.array([5005.689453, 4481.327637, 6010.369629], dtype=np.float32) + counts_hist, xedges, yedges = np.histogram2d(x, y, bins=100) + assert_equal(counts_hist.sum(), 3.) + def test_weights(self): v = rand(100) w = np.ones(100) * 5 @@ -1139,6 +1219,30 @@ class TestHistogramdd(TestCase): h, e = np.histogramdd(x, bins=[3, [-np.inf, 3, np.inf]]) assert_allclose(h, expected) + def test_rightmost_binedge(self): + """Test event very close to rightmost binedge. + See Github issue #4266""" + x = [0.9999999995] + bins = [[0.,0.5,1.0]] + hist, _ = histogramdd(x, bins=bins) + assert_(hist[0] == 0.0) + assert_(hist[1] == 1.) + x = [1.0] + bins = [[0.,0.5,1.0]] + hist, _ = histogramdd(x, bins=bins) + assert_(hist[0] == 0.0) + assert_(hist[1] == 1.) + x = [1.0000000001] + bins = [[0.,0.5,1.0]] + hist, _ = histogramdd(x, bins=bins) + assert_(hist[0] == 0.0) + assert_(hist[1] == 1.) + x = [1.0001] + bins = [[0.,0.5,1.0]] + hist, _ = histogramdd(x, bins=bins) + assert_(hist[0] == 0.0) + assert_(hist[1] == 0.0) + class TestUnique(TestCase): def test_simple(self): @@ -1188,6 +1292,10 @@ class TestCorrCoef(TestCase): [0.66318558, 0.88157256, 0.71483595, -0.51366032, 1., 0.98317823], [0.51532523, 0.78052386, 0.83053601, -0.66173113, 0.98317823, 1.]]) + def test_non_array(self): + assert_almost_equal(np.corrcoef([0, 1, 0], [1, 0, 1]), + [[1., -1.], [-1., 1.]]) + def test_simple(self): assert_almost_equal(corrcoef(self.A), self.res1) assert_almost_equal(corrcoef(self.A, self.B), self.res2) @@ -1206,8 +1314,8 @@ class TestCorrCoef(TestCase): assert_allclose(np.corrcoef(x, y), np.array([[1., -1.j], [1.j, 1.]])) def test_empty(self): - with warnings.catch_warnings(): - warnings.simplefilter('ignore', RuntimeWarning) + with warnings.catch_warnings(record=True): + warnings.simplefilter('always', RuntimeWarning) assert_array_equal(corrcoef(np.array([])), np.nan) assert_array_equal(corrcoef(np.array([]).reshape(0, 2)), np.array([]).reshape(0, 0)) @@ -1216,8 +1324,8 @@ class TestCorrCoef(TestCase): def test_wrong_ddof(self): x = np.array([[0, 2], [1, 1], [2, 0]]).T - with warnings.catch_warnings(): - warnings.simplefilter('ignore', RuntimeWarning) + with warnings.catch_warnings(record=True): + warnings.simplefilter('always', RuntimeWarning) assert_array_equal(corrcoef(x, ddof=5), np.array([[np.nan, np.nan], [np.nan, np.nan]])) @@ -1237,8 +1345,8 @@ class TestCov(TestCase): assert_allclose(cov(x, y), np.array([[1., -1.j], [1.j, 1.]])) def test_empty(self): - with warnings.catch_warnings(): - warnings.simplefilter('ignore', RuntimeWarning) + with warnings.catch_warnings(record=True): + warnings.simplefilter('always', RuntimeWarning) assert_array_equal(cov(np.array([])), np.nan) assert_array_equal(cov(np.array([]).reshape(0, 2)), np.array([]).reshape(0, 0)) @@ -1247,8 +1355,8 @@ class TestCov(TestCase): def test_wrong_ddof(self): x = np.array([[0, 2], [1, 1], [2, 0]]).T - with warnings.catch_warnings(): - warnings.simplefilter('ignore', RuntimeWarning) + with warnings.catch_warnings(record=True): + warnings.simplefilter('always', RuntimeWarning) assert_array_equal(cov(x, ddof=5), np.array([[np.inf, -np.inf], [-np.inf, np.inf]])) @@ -1353,6 +1461,13 @@ class TestMeshgrid(TestCase): assert_array_equal(X, np.array([[1, 2, 3]])) assert_array_equal(Y, np.array([[4], [5], [6], [7]])) + def test_invalid_arguments(self): + # Test that meshgrid complains about invalid arguments + # Regression test for issue #4755: + # https://github.com/numpy/numpy/issues/4755 + assert_raises(TypeError, meshgrid, + [1, 2, 3], [4, 5, 6, 7], indices='ij') + class TestPiecewise(TestCase): def test_simple(self): @@ -1379,6 +1494,7 @@ class TestPiecewise(TestCase): x = piecewise([0, 0], [[False, True]], [lambda x:-1]) assert_array_equal(x, [0, -1]) + def test_two_conditions(self): x = piecewise([1, 2], [[True, False], [False, True]], [3, 4]) assert_array_equal(x, [3, 4]) @@ -1397,6 +1513,15 @@ class TestPiecewise(TestCase): assert_(y.ndim == 0) assert_(y == 0) + x = 5 + y = piecewise(x, [[True], [False]], [1, 0]) + assert_(y.ndim == 0) + assert_(y == 1) + + def test_0d_comparison(self): + x = 3 + y = piecewise(x, [x <= 3, x > 3], [4, 0]) + class TestBincount(TestCase): def test_simple(self): @@ -1445,6 +1570,23 @@ class TestBincount(TestCase): y = np.bincount(x, minlength=5) assert_array_equal(y, np.zeros(5, dtype=int)) + def test_with_incorrect_minlength(self): + x = np.array([], dtype=int) + assert_raises_regex(TypeError, "an integer is required", + lambda: np.bincount(x, minlength="foobar")) + assert_raises_regex(ValueError, "must be positive", + lambda: np.bincount(x, minlength=-1)) + assert_raises_regex(ValueError, "must be positive", + lambda: np.bincount(x, minlength=0)) + + x = np.arange(5) + assert_raises_regex(TypeError, "an integer is required", + lambda: np.bincount(x, minlength="foobar")) + assert_raises_regex(ValueError, "minlength must be positive", + lambda: np.bincount(x, minlength=-1)) + assert_raises_regex(ValueError, "minlength must be positive", + lambda: np.bincount(x, minlength=0)) + class TestInterp(TestCase): def test_exceptions(self): @@ -1474,6 +1616,8 @@ class TestInterp(TestCase): assert_almost_equal(np.interp(x0, x, y), x0) x0 = np.float64(.3) assert_almost_equal(np.interp(x0, x, y), x0) + x0 = np.nan + assert_almost_equal(np.interp(x0, x, y), x0) def test_zero_dimensional_interpolation_point(self): x = np.linspace(0, 1, 5) @@ -1636,6 +1780,8 @@ class TestScoreatpercentile(TestCase): interpolation='foobar') assert_raises(ValueError, np.percentile, [1], 101) assert_raises(ValueError, np.percentile, [1], -1) + assert_raises(ValueError, np.percentile, [1], list(range(50)) + [101]) + assert_raises(ValueError, np.percentile, [1], list(range(50)) + [-0.1]) def test_percentile_list(self): assert_equal(np.percentile([1, 2, 3], 0), 1) @@ -1727,6 +1873,65 @@ class TestScoreatpercentile(TestCase): b = np.percentile([2, 3, 4, 1], [50], overwrite_input=True) assert_equal(b, np.array([2.5])) + def test_extended_axis(self): + o = np.random.normal(size=(71, 23)) + x = np.dstack([o] * 10) + assert_equal(np.percentile(x, 30, axis=(0, 1)), np.percentile(o, 30)) + x = np.rollaxis(x, -1, 0) + assert_equal(np.percentile(x, 30, axis=(-2, -1)), np.percentile(o, 30)) + x = x.swapaxes(0, 1).copy() + assert_equal(np.percentile(x, 30, axis=(0, -1)), np.percentile(o, 30)) + x = x.swapaxes(0, 1).copy() + + assert_equal(np.percentile(x, [25, 60], axis=(0, 1, 2)), + np.percentile(x, [25, 60], axis=None)) + assert_equal(np.percentile(x, [25, 60], axis=(0,)), + np.percentile(x, [25, 60], axis=0)) + + d = np.arange(3 * 5 * 7 * 11).reshape(3, 5, 7, 11) + np.random.shuffle(d) + assert_equal(np.percentile(d, 25, axis=(0, 1, 2))[0], + np.percentile(d[:, :, :, 0].flatten(), 25)) + assert_equal(np.percentile(d, [10, 90], axis=(0, 1, 3))[:, 1], + np.percentile(d[:, :, 1, :].flatten(), [10, 90])) + assert_equal(np.percentile(d, 25, axis=(3, 1, -4))[2], + np.percentile(d[:, :, 2, :].flatten(), 25)) + assert_equal(np.percentile(d, 25, axis=(3, 1, 2))[2], + np.percentile(d[2, :, :, :].flatten(), 25)) + assert_equal(np.percentile(d, 25, axis=(3, 2))[2, 1], + np.percentile(d[2, 1, :, :].flatten(), 25)) + assert_equal(np.percentile(d, 25, axis=(1, -2))[2, 1], + np.percentile(d[2, :, :, 1].flatten(), 25)) + assert_equal(np.percentile(d, 25, axis=(1, 3))[2, 2], + np.percentile(d[2, :, 2, :].flatten(), 25)) + + def test_extended_axis_invalid(self): + d = np.ones((3, 5, 7, 11)) + assert_raises(IndexError, np.percentile, d, axis=-5, q=25) + assert_raises(IndexError, np.percentile, d, axis=(0, -5), q=25) + assert_raises(IndexError, np.percentile, d, axis=4, q=25) + assert_raises(IndexError, np.percentile, d, axis=(0, 4), q=25) + assert_raises(ValueError, np.percentile, d, axis=(1, 1), q=25) + + def test_keepdims(self): + d = np.ones((3, 5, 7, 11)) + assert_equal(np.percentile(d, 7, axis=None, keepdims=True).shape, + (1, 1, 1, 1)) + assert_equal(np.percentile(d, 7, axis=(0, 1), keepdims=True).shape, + (1, 1, 7, 11)) + assert_equal(np.percentile(d, 7, axis=(0, 3), keepdims=True).shape, + (1, 5, 7, 1)) + assert_equal(np.percentile(d, 7, axis=(1,), keepdims=True).shape, + (3, 1, 7, 11)) + assert_equal(np.percentile(d, 7, (0, 1, 2, 3), keepdims=True).shape, + (1, 1, 1, 1)) + assert_equal(np.percentile(d, 7, axis=(0, 1, 3), keepdims=True).shape, + (1, 1, 7, 1)) + + assert_equal(np.percentile(d, [1, 7], axis=(0, 1, 3), + keepdims=True).shape, (2, 1, 1, 7, 1)) + assert_equal(np.percentile(d, [1, 7], axis=(0, 3), + keepdims=True).shape, (2, 1, 5, 7, 1)) class TestMedian(TestCase): @@ -1734,19 +1939,23 @@ class TestMedian(TestCase): a0 = np.array(1) a1 = np.arange(2) a2 = np.arange(6).reshape(2, 3) - assert_allclose(np.median(a0), 1) + assert_equal(np.median(a0), 1) assert_allclose(np.median(a1), 0.5) assert_allclose(np.median(a2), 2.5) assert_allclose(np.median(a2, axis=0), [1.5, 2.5, 3.5]) - assert_allclose(np.median(a2, axis=1), [1, 4]) + assert_equal(np.median(a2, axis=1), [1, 4]) assert_allclose(np.median(a2, axis=None), 2.5) a = np.array([0.0444502, 0.0463301, 0.141249, 0.0606775]) assert_almost_equal((a[1] + a[3]) / 2., np.median(a)) a = np.array([0.0463301, 0.0444502, 0.141249]) - assert_almost_equal(a[0], np.median(a)) + assert_equal(a[0], np.median(a)) a = np.array([0.0444502, 0.141249, 0.0463301]) - assert_almost_equal(a[-1], np.median(a)) + assert_equal(a[-1], np.median(a)) + # check array scalar result + assert_equal(np.median(a).ndim, 0) + a[1] = np.nan + assert_equal(np.median(a).ndim, 0) def test_axis_keyword(self): a3 = np.array([[2, 3], @@ -1820,6 +2029,66 @@ class TestMedian(TestCase): a = MySubClass([1,2,3]) assert_equal(np.median(a), -7) + def test_object(self): + o = np.arange(7.); + assert_(type(np.median(o.astype(object))), float) + o[2] = np.nan + assert_(type(np.median(o.astype(object))), float) + + def test_extended_axis(self): + o = np.random.normal(size=(71, 23)) + x = np.dstack([o] * 10) + assert_equal(np.median(x, axis=(0, 1)), np.median(o)) + x = np.rollaxis(x, -1, 0) + assert_equal(np.median(x, axis=(-2, -1)), np.median(o)) + x = x.swapaxes(0, 1).copy() + assert_equal(np.median(x, axis=(0, -1)), np.median(o)) + + assert_equal(np.median(x, axis=(0, 1, 2)), np.median(x, axis=None)) + assert_equal(np.median(x, axis=(0, )), np.median(x, axis=0)) + assert_equal(np.median(x, axis=(-1, )), np.median(x, axis=-1)) + + d = np.arange(3 * 5 * 7 * 11).reshape(3, 5, 7, 11) + np.random.shuffle(d) + assert_equal(np.median(d, axis=(0, 1, 2))[0], + np.median(d[:, :, :, 0].flatten())) + assert_equal(np.median(d, axis=(0, 1, 3))[1], + np.median(d[:, :, 1, :].flatten())) + assert_equal(np.median(d, axis=(3, 1, -4))[2], + np.median(d[:, :, 2, :].flatten())) + assert_equal(np.median(d, axis=(3, 1, 2))[2], + np.median(d[2, :, :, :].flatten())) + assert_equal(np.median(d, axis=(3, 2))[2, 1], + np.median(d[2, 1, :, :].flatten())) + assert_equal(np.median(d, axis=(1, -2))[2, 1], + np.median(d[2, :, :, 1].flatten())) + assert_equal(np.median(d, axis=(1, 3))[2, 2], + np.median(d[2, :, 2, :].flatten())) + + def test_extended_axis_invalid(self): + d = np.ones((3, 5, 7, 11)) + assert_raises(IndexError, np.median, d, axis=-5) + assert_raises(IndexError, np.median, d, axis=(0, -5)) + assert_raises(IndexError, np.median, d, axis=4) + assert_raises(IndexError, np.median, d, axis=(0, 4)) + assert_raises(ValueError, np.median, d, axis=(1, 1)) + + def test_keepdims(self): + d = np.ones((3, 5, 7, 11)) + assert_equal(np.median(d, axis=None, keepdims=True).shape, + (1, 1, 1, 1)) + assert_equal(np.median(d, axis=(0, 1), keepdims=True).shape, + (1, 1, 7, 11)) + assert_equal(np.median(d, axis=(0, 3), keepdims=True).shape, + (1, 5, 7, 1)) + assert_equal(np.median(d, axis=(1,), keepdims=True).shape, + (3, 1, 7, 11)) + assert_equal(np.median(d, axis=(0, 1, 2, 3), keepdims=True).shape, + (1, 1, 1, 1)) + assert_equal(np.median(d, axis=(0, 1, 3), keepdims=True).shape, + (1, 1, 7, 1)) + + class TestAdd_newdoc_ufunc(TestCase): diff --git a/numpy/lib/tests/test_io.py b/numpy/lib/tests/test_io.py index e3ccb391c..03e238261 100644 --- a/numpy/lib/tests/test_io.py +++ b/numpy/lib/tests/test_io.py @@ -4,7 +4,7 @@ import sys import gzip import os import threading -from tempfile import mkstemp, mktemp, NamedTemporaryFile +from tempfile import mkstemp, NamedTemporaryFile import time import warnings import gc @@ -17,9 +17,12 @@ from numpy.lib._iotools import (ConverterError, ConverterLockError, ConversionWarning) from numpy.compat import asbytes, asbytes_nested, bytes, asstr from nose import SkipTest -from numpy.ma.testutils import (TestCase, assert_equal, assert_array_equal, - assert_raises, run_module_suite) +from numpy.ma.testutils import ( + TestCase, assert_equal, assert_array_equal, + assert_raises, assert_raises_regex, run_module_suite +) from numpy.testing import assert_warns, assert_, build_err_msg +from numpy.testing.utils import tempdir class TextIO(BytesIO): @@ -77,32 +80,32 @@ class RoundtripTest(object): file_on_disk = kwargs.get('file_on_disk', False) if file_on_disk: - # Do not delete the file on windows, because we can't - # reopen an already opened file on that platform, so we - # need to close the file and reopen it, implying no - # automatic deletion. - if sys.platform == 'win32' and MAJVER >= 2 and MINVER >= 6: - target_file = NamedTemporaryFile(delete=False) - else: - target_file = NamedTemporaryFile() + target_file = NamedTemporaryFile(delete=False) load_file = target_file.name else: target_file = BytesIO() load_file = target_file - arr = args + try: + arr = args - save_func(target_file, *arr, **save_kwds) - target_file.flush() - target_file.seek(0) + save_func(target_file, *arr, **save_kwds) + target_file.flush() + target_file.seek(0) - if sys.platform == 'win32' and not isinstance(target_file, BytesIO): - target_file.close() + if sys.platform == 'win32' and not isinstance(target_file, BytesIO): + target_file.close() - arr_reloaded = np.load(load_file, **load_kwds) + arr_reloaded = np.load(load_file, **load_kwds) - self.arr = arr - self.arr_reloaded = arr_reloaded + self.arr = arr + self.arr_reloaded = arr_reloaded + finally: + if not isinstance(target_file, BytesIO): + target_file.close() + # holds an open file descriptor so it can't be deleted on win + if not isinstance(arr_reloaded, np.lib.npyio.NpzFile): + os.remove(target_file.name) def check_roundtrips(self, a): self.roundtrip(a) @@ -155,6 +158,13 @@ class RoundtripTest(object): a = np.array([(1, 2), (3, 4)], dtype=[('x', 'i4'), ('y', 'i4')]) self.check_roundtrips(a) + def test_format_2_0(self): + dt = [(("%d" % i) * 100, float) for i in range(500)] + a = np.ones(1000, dtype=dt) + with warnings.catch_warnings(record=True): + warnings.filterwarnings('always', '', UserWarning) + self.check_roundtrips(a) + class TestSaveLoad(RoundtripTest, TestCase): def roundtrip(self, *args, **kwargs): @@ -167,24 +177,30 @@ class TestSaveLoad(RoundtripTest, TestCase): class TestSavezLoad(RoundtripTest, TestCase): def roundtrip(self, *args, **kwargs): RoundtripTest.roundtrip(self, np.savez, *args, **kwargs) - for n, arr in enumerate(self.arr): - reloaded = self.arr_reloaded['arr_%d' % n] - assert_equal(arr, reloaded) - assert_equal(arr.dtype, reloaded.dtype) - assert_equal(arr.flags.fnc, reloaded.flags.fnc) + try: + for n, arr in enumerate(self.arr): + reloaded = self.arr_reloaded['arr_%d' % n] + assert_equal(arr, reloaded) + assert_equal(arr.dtype, reloaded.dtype) + assert_equal(arr.flags.fnc, reloaded.flags.fnc) + finally: + # delete tempfile, must be done here on windows + if self.arr_reloaded.fid: + self.arr_reloaded.fid.close() + os.remove(self.arr_reloaded.fid.name) @np.testing.dec.skipif(not IS_64BIT, "Works only with 64bit systems") @np.testing.dec.slow def test_big_arrays(self): L = (1 << 31) + 100000 - tmp = mktemp(suffix='.npz') a = np.empty(L, dtype=np.uint8) - np.savez(tmp, a=a) - del a - npfile = np.load(tmp) - a = npfile['a'] - npfile.close() - os.remove(tmp) + with tempdir(prefix="numpy_test_big_arrays_") as tmpdir: + tmp = os.path.join(tmpdir, "file.npz") + np.savez(tmp, a=a) + del a + npfile = np.load(tmp) + a = npfile['a'] + npfile.close() def test_multiple_arrays(self): a = np.array([[1, 2], [3, 4]], float) @@ -287,13 +303,14 @@ class TestSavezLoad(RoundtripTest, TestCase): # Check that zipfile owns file and can close it. # This needs to pass a file name to load for the # test. - fd, tmp = mkstemp(suffix='.npz') - os.close(fd) - np.savez(tmp, lab='place holder') - data = np.load(tmp) - fp = data.zip.fp - data.close() - assert_(fp.closed) + with tempdir(prefix="numpy_test_closing_zipfile_after_load_") as tmpdir: + fd, tmp = mkstemp(suffix='.npz', dir=tmpdir) + os.close(fd) + np.savez(tmp, lab='place holder') + data = np.load(tmp) + fp = data.zip.fp + data.close() + assert_(fp.closed) class TestSaveTxt(TestCase): @@ -750,6 +767,14 @@ class TestLoadTxt(TestCase): res = np.loadtxt(count()) assert_array_equal(res, np.arange(10)) + def test_bad_line(self): + c = TextIO() + c.write('1 2 3\n4 5 6\n2 3') + c.seek(0) + + # Check for exception and that exception contains line number + assert_raises_regex(ValueError, "3", np.loadtxt, c) + class Testfromregex(TestCase): # np.fromregex expects files opened in binary mode. @@ -1563,6 +1588,14 @@ M 33 21.99 dtype=[('a', np.int), ('b', np.int)]) self.assertTrue(isinstance(test, np.recarray)) assert_equal(test, control) + # + data = TextIO('A,B\n0,1\n2,3') + dtype = [('a', np.int), ('b', np.float)] + test = np.recfromcsv(data, missing_values='N/A', dtype=dtype) + control = np.array([(0, 1), (2, 3)], + dtype=dtype) + self.assertTrue(isinstance(test, np.recarray)) + assert_equal(test, control) def test_gft_using_filename(self): # Test that we can load data from a filename as well as a file object diff --git a/numpy/lib/tests/test_nanfunctions.py b/numpy/lib/tests/test_nanfunctions.py index af01a7167..c5af61434 100644 --- a/numpy/lib/tests/test_nanfunctions.py +++ b/numpy/lib/tests/test_nanfunctions.py @@ -5,21 +5,21 @@ import warnings import numpy as np from numpy.testing import ( run_module_suite, TestCase, assert_, assert_equal, assert_almost_equal, - assert_raises + assert_raises, assert_array_equal ) # Test data _ndat = np.array([[0.6244, np.nan, 0.2692, 0.0116, np.nan, 0.1170], - [0.5351, 0.9403, np.nan, 0.2100, 0.4759, 0.2833], - [np.nan, np.nan, np.nan, 0.1042, np.nan, 0.5954], + [0.5351, -0.9403, np.nan, 0.2100, 0.4759, 0.2833], + [np.nan, np.nan, np.nan, 0.1042, np.nan, -0.5954], [0.1610, np.nan, np.nan, 0.1859, 0.3146, np.nan]]) # Rows of _ndat with nans removed _rdat = [np.array([ 0.6244, 0.2692, 0.0116, 0.1170]), - np.array([ 0.5351, 0.9403, 0.2100, 0.4759, 0.2833]), - np.array([ 0.1042, 0.5954]), + np.array([ 0.5351, -0.9403, 0.2100, 0.4759, 0.2833]), + np.array([ 0.1042, -0.5954]), np.array([ 0.1610, 0.1859, 0.3146])] @@ -114,6 +114,31 @@ class TestNanFunctions_MinMax(TestCase): assert_(res.shape == (3, 1)) res = f(mat) assert_(np.isscalar(res)) + # check that rows of nan are dealt with for subclasses (#4628) + mat[1] = np.nan + for f in self.nanfuncs: + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter('always') + res = f(mat, axis=0) + assert_(isinstance(res, np.matrix)) + assert_(not np.any(np.isnan(res))) + assert_(len(w) == 0) + + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter('always') + res = f(mat, axis=1) + assert_(isinstance(res, np.matrix)) + assert_(np.isnan(res[1, 0]) and not np.isnan(res[0, 0]) + and not np.isnan(res[2, 0])) + assert_(len(w) == 1, 'no warning raised') + assert_(issubclass(w[0].category, RuntimeWarning)) + + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter('always') + res = f(mat) + assert_(np.isscalar(res)) + assert_(res != np.nan) + assert_(len(w) == 0) class TestNanFunctions_ArgminArgmax(TestCase): @@ -130,8 +155,8 @@ class TestNanFunctions_ArgminArgmax(TestCase): def test_result_values(self): for f, fcmp in zip(self.nanfuncs, [np.greater, np.less]): for row in _ndat: - with warnings.catch_warnings(): - warnings.simplefilter('ignore') + with warnings.catch_warnings(record=True): + warnings.simplefilter('always') ind = f(row) val = row[ind] # comparing with NaN is tricky as the result @@ -502,5 +527,232 @@ class TestNanFunctions_MeanVarStd(TestCase): assert_(np.isscalar(res)) +class TestNanFunctions_Median(TestCase): + + def test_mutation(self): + # Check that passed array is not modified. + ndat = _ndat.copy() + np.nanmedian(ndat) + assert_equal(ndat, _ndat) + + def test_keepdims(self): + mat = np.eye(3) + for axis in [None, 0, 1]: + tgt = np.median(mat, axis=axis, out=None, overwrite_input=False) + res = np.nanmedian(mat, axis=axis, out=None, overwrite_input=False) + assert_(res.ndim == tgt.ndim) + + d = np.ones((3, 5, 7, 11)) + # Randomly set some elements to NaN: + w = np.random.random((4, 200)) * np.array(d.shape)[:, None] + w = w.astype(np.intp) + d[tuple(w)] = np.nan + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter('always', RuntimeWarning) + res = np.nanmedian(d, axis=None, keepdims=True) + assert_equal(res.shape, (1, 1, 1, 1)) + res = np.nanmedian(d, axis=(0, 1), keepdims=True) + assert_equal(res.shape, (1, 1, 7, 11)) + res = np.nanmedian(d, axis=(0, 3), keepdims=True) + assert_equal(res.shape, (1, 5, 7, 1)) + res = np.nanmedian(d, axis=(1,), keepdims=True) + assert_equal(res.shape, (3, 1, 7, 11)) + res = np.nanmedian(d, axis=(0, 1, 2, 3), keepdims=True) + assert_equal(res.shape, (1, 1, 1, 1)) + res = np.nanmedian(d, axis=(0, 1, 3), keepdims=True) + assert_equal(res.shape, (1, 1, 7, 1)) + + def test_out(self): + mat = np.random.rand(3, 3) + nan_mat = np.insert(mat, [0, 2], np.nan, axis=1) + resout = np.zeros(3) + tgt = np.median(mat, axis=1) + res = np.nanmedian(nan_mat, axis=1, out=resout) + assert_almost_equal(res, resout) + assert_almost_equal(res, tgt) + # 0-d output: + resout = np.zeros(()) + tgt = np.median(mat, axis=None) + res = np.nanmedian(nan_mat, axis=None, out=resout) + assert_almost_equal(res, resout) + assert_almost_equal(res, tgt) + res = np.nanmedian(nan_mat, axis=(0, 1), out=resout) + assert_almost_equal(res, resout) + assert_almost_equal(res, tgt) + + def test_small_large(self): + # test the small and large code paths, current cutoff 400 elements + for s in [5, 20, 51, 200, 1000]: + d = np.random.randn(4, s) + # Randomly set some elements to NaN: + w = np.random.randint(0, d.size, size=d.size // 5) + d.ravel()[w] = np.nan + d[:,0] = 1. # ensure at least one good value + # use normal median without nans to compare + tgt = [] + for x in d: + nonan = np.compress(~np.isnan(x), x) + tgt.append(np.median(nonan, overwrite_input=True)) + + assert_array_equal(np.nanmedian(d, axis=-1), tgt) + + def test_result_values(self): + tgt = [np.median(d) for d in _rdat] + res = np.nanmedian(_ndat, axis=1) + assert_almost_equal(res, tgt) + + def test_allnans(self): + mat = np.array([np.nan]*9).reshape(3, 3) + for axis in [None, 0, 1]: + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter('always') + assert_(np.isnan(np.nanmedian(mat, axis=axis)).all()) + if axis is None: + assert_(len(w) == 1) + else: + assert_(len(w) == 3) + assert_(issubclass(w[0].category, RuntimeWarning)) + # Check scalar + assert_(np.isnan(np.nanmedian(np.nan))) + if axis is None: + assert_(len(w) == 2) + else: + assert_(len(w) == 4) + assert_(issubclass(w[0].category, RuntimeWarning)) + + def test_empty(self): + mat = np.zeros((0, 3)) + for axis in [0, None]: + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter('always') + assert_(np.isnan(np.nanmedian(mat, axis=axis)).all()) + assert_(len(w) == 1) + assert_(issubclass(w[0].category, RuntimeWarning)) + for axis in [1]: + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter('always') + assert_equal(np.nanmedian(mat, axis=axis), np.zeros([])) + assert_(len(w) == 0) + + def test_scalar(self): + assert_(np.nanmedian(0.) == 0.) + + def test_extended_axis_invalid(self): + d = np.ones((3, 5, 7, 11)) + assert_raises(IndexError, np.nanmedian, d, axis=-5) + assert_raises(IndexError, np.nanmedian, d, axis=(0, -5)) + assert_raises(IndexError, np.nanmedian, d, axis=4) + assert_raises(IndexError, np.nanmedian, d, axis=(0, 4)) + assert_raises(ValueError, np.nanmedian, d, axis=(1, 1)) + + +class TestNanFunctions_Percentile(TestCase): + + def test_mutation(self): + # Check that passed array is not modified. + ndat = _ndat.copy() + np.nanpercentile(ndat, 30) + assert_equal(ndat, _ndat) + + def test_keepdims(self): + mat = np.eye(3) + for axis in [None, 0, 1]: + tgt = np.percentile(mat, 70, axis=axis, out=None, + overwrite_input=False) + res = np.nanpercentile(mat, 70, axis=axis, out=None, + overwrite_input=False) + assert_(res.ndim == tgt.ndim) + + d = np.ones((3, 5, 7, 11)) + # Randomly set some elements to NaN: + w = np.random.random((4, 200)) * np.array(d.shape)[:, None] + w = w.astype(np.intp) + d[tuple(w)] = np.nan + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter('always', RuntimeWarning) + res = np.nanpercentile(d, 90, axis=None, keepdims=True) + assert_equal(res.shape, (1, 1, 1, 1)) + res = np.nanpercentile(d, 90, axis=(0, 1), keepdims=True) + assert_equal(res.shape, (1, 1, 7, 11)) + res = np.nanpercentile(d, 90, axis=(0, 3), keepdims=True) + assert_equal(res.shape, (1, 5, 7, 1)) + res = np.nanpercentile(d, 90, axis=(1,), keepdims=True) + assert_equal(res.shape, (3, 1, 7, 11)) + res = np.nanpercentile(d, 90, axis=(0, 1, 2, 3), keepdims=True) + assert_equal(res.shape, (1, 1, 1, 1)) + res = np.nanpercentile(d, 90, axis=(0, 1, 3), keepdims=True) + assert_equal(res.shape, (1, 1, 7, 1)) + + def test_out(self): + mat = np.random.rand(3, 3) + nan_mat = np.insert(mat, [0, 2], np.nan, axis=1) + resout = np.zeros(3) + tgt = np.percentile(mat, 42, axis=1) + res = np.nanpercentile(nan_mat, 42, axis=1, out=resout) + assert_almost_equal(res, resout) + assert_almost_equal(res, tgt) + # 0-d output: + resout = np.zeros(()) + tgt = np.percentile(mat, 42, axis=None) + res = np.nanpercentile(nan_mat, 42, axis=None, out=resout) + assert_almost_equal(res, resout) + assert_almost_equal(res, tgt) + res = np.nanpercentile(nan_mat, 42, axis=(0, 1), out=resout) + assert_almost_equal(res, resout) + assert_almost_equal(res, tgt) + + def test_result_values(self): + tgt = [np.percentile(d, 28) for d in _rdat] + res = np.nanpercentile(_ndat, 28, axis=1) + assert_almost_equal(res, tgt) + tgt = [np.percentile(d, (28, 98)) for d in _rdat] + res = np.nanpercentile(_ndat, (28, 98), axis=1) + assert_almost_equal(res, tgt) + + def test_allnans(self): + mat = np.array([np.nan]*9).reshape(3, 3) + for axis in [None, 0, 1]: + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter('always') + assert_(np.isnan(np.nanpercentile(mat, 60, axis=axis)).all()) + if axis is None: + assert_(len(w) == 1) + else: + assert_(len(w) == 3) + assert_(issubclass(w[0].category, RuntimeWarning)) + # Check scalar + assert_(np.isnan(np.nanpercentile(np.nan, 60))) + if axis is None: + assert_(len(w) == 2) + else: + assert_(len(w) == 4) + assert_(issubclass(w[0].category, RuntimeWarning)) + + def test_empty(self): + mat = np.zeros((0, 3)) + for axis in [0, None]: + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter('always') + assert_(np.isnan(np.nanpercentile(mat, 40, axis=axis)).all()) + assert_(len(w) == 1) + assert_(issubclass(w[0].category, RuntimeWarning)) + for axis in [1]: + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter('always') + assert_equal(np.nanpercentile(mat, 40, axis=axis), np.zeros([])) + assert_(len(w) == 0) + + def test_scalar(self): + assert_(np.nanpercentile(0., 100) == 0.) + + def test_extended_axis_invalid(self): + d = np.ones((3, 5, 7, 11)) + assert_raises(IndexError, np.nanpercentile, d, q=5, axis=-5) + assert_raises(IndexError, np.nanpercentile, d, q=5, axis=(0, -5)) + assert_raises(IndexError, np.nanpercentile, d, q=5, axis=4) + assert_raises(IndexError, np.nanpercentile, d, q=5, axis=(0, 4)) + assert_raises(ValueError, np.nanpercentile, d, q=5, axis=(1, 1)) + + if __name__ == "__main__": run_module_suite() diff --git a/numpy/lib/tests/test_stride_tricks.py b/numpy/lib/tests/test_stride_tricks.py index 4d57d2ca7..3adcf0b41 100644 --- a/numpy/lib/tests/test_stride_tricks.py +++ b/numpy/lib/tests/test_stride_tricks.py @@ -3,6 +3,7 @@ from __future__ import division, absolute_import, print_function import numpy as np from numpy.testing import * from numpy.lib.stride_tricks import broadcast_arrays +from numpy.lib.stride_tricks import as_strided def assert_shapes_correct(input_shapes, expected_shape): @@ -214,6 +215,22 @@ def test_same_as_ufunc(): assert_same_as_ufunc(input_shapes[0], input_shapes[1], False, True) assert_same_as_ufunc(input_shapes[0], input_shapes[1], True, True) +def test_as_strided(): + a = np.array([None]) + a_view = as_strided(a) + expected = np.array([None]) + assert_array_equal(a_view, np.array([None])) + + a = np.array([1, 2, 3, 4]) + a_view = as_strided(a, shape=(2,), strides=(2 * a.itemsize,)) + expected = np.array([1, 3]) + assert_array_equal(a_view, expected) + + a = np.array([1, 2, 3, 4]) + a_view = as_strided(a, shape=(3, 4), strides=(0, 1 * a.itemsize)) + expected = np.array([[1, 2, 3, 4], [1, 2, 3, 4], [1, 2, 3, 4]]) + assert_array_equal(a_view, expected) + if __name__ == "__main__": run_module_suite() diff --git a/numpy/lib/tests/test_twodim_base.py b/numpy/lib/tests/test_twodim_base.py index 8d0275a25..f5b8fab4a 100644 --- a/numpy/lib/tests/test_twodim_base.py +++ b/numpy/lib/tests/test_twodim_base.py @@ -254,7 +254,7 @@ class TestHistogram2d(TestCase): assert_array_almost_equal(H, answer, 3) def test_all_outliers(self): - r = rand(100)+1. + r = rand(100) + 1. + 1e6 # histogramdd rounds by decimal=6 H, xed, yed = histogram2d(r, r, (4, 5), range=([0, 1], [0, 1])) assert_array_equal(H, 0) @@ -275,16 +275,41 @@ class TestTri(TestCase): assert_array_equal(tri(3, dtype=bool), out.astype(bool)) -def test_tril_triu(): +def test_tril_triu_ndim2(): for dtype in np.typecodes['AllFloat'] + np.typecodes['AllInteger']: a = np.ones((2, 2), dtype=dtype) b = np.tril(a) c = np.triu(a) - assert_array_equal(b, [[1, 0], [1, 1]]) - assert_array_equal(c, b.T) + yield assert_array_equal, b, [[1, 0], [1, 1]] + yield assert_array_equal, c, b.T # should return the same dtype as the original array - assert_equal(b.dtype, a.dtype) - assert_equal(c.dtype, a.dtype) + yield assert_equal, b.dtype, a.dtype + yield assert_equal, c.dtype, a.dtype + + +def test_tril_triu_ndim3(): + for dtype in np.typecodes['AllFloat'] + np.typecodes['AllInteger']: + a = np.array([ + [[1, 1], [1, 1]], + [[1, 1], [1, 0]], + [[1, 1], [0, 0]], + ], dtype=dtype) + a_tril_desired = np.array([ + [[1, 0], [1, 1]], + [[1, 0], [1, 0]], + [[1, 0], [0, 0]], + ], dtype=dtype) + a_triu_desired = np.array([ + [[1, 1], [0, 1]], + [[1, 1], [0, 0]], + [[1, 1], [0, 0]], + ], dtype=dtype) + a_triu_observed = np.triu(a) + a_tril_observed = np.tril(a) + yield assert_array_equal, a_triu_observed, a_triu_desired + yield assert_array_equal, a_tril_observed, a_tril_desired + yield assert_equal, a_triu_observed.dtype, a.dtype + yield assert_equal, a_tril_observed.dtype, a.dtype def test_mask_indices(): @@ -300,16 +325,21 @@ def test_mask_indices(): def test_tril_indices(): # indices without and with offset il1 = tril_indices(4) - il2 = tril_indices(4, 2) + il2 = tril_indices(4, k=2) + il3 = tril_indices(4, m=5) + il4 = tril_indices(4, k=2, m=5) a = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]]) + b = np.arange(1, 21).reshape(4, 5) # indexing: yield (assert_array_equal, a[il1], array([1, 5, 6, 9, 10, 11, 13, 14, 15, 16])) + yield (assert_array_equal, b[il3], + array([1, 6, 7, 11, 12, 13, 16, 17, 18, 19])) # And for assigning values: a[il1] = -1 @@ -318,7 +348,12 @@ def test_tril_indices(): [-1, -1, 7, 8], [-1, -1, -1, 12], [-1, -1, -1, -1]])) - + b[il3] = -1 + yield (assert_array_equal, b, + array([[-1, 2, 3, 4, 5], + [-1, -1, 8, 9, 10], + [-1, -1, -1, 14, 15], + [-1, -1, -1, -1, 20]])) # These cover almost the whole array (two diagonals right of the main one): a[il2] = -10 yield (assert_array_equal, a, @@ -326,21 +361,32 @@ def test_tril_indices(): [-10, -10, -10, -10], [-10, -10, -10, -10], [-10, -10, -10, -10]])) + b[il4] = -10 + yield (assert_array_equal, b, + array([[-10, -10, -10, 4, 5], + [-10, -10, -10, -10, 10], + [-10, -10, -10, -10, -10], + [-10, -10, -10, -10, -10]])) class TestTriuIndices(object): def test_triu_indices(self): iu1 = triu_indices(4) - iu2 = triu_indices(4, 2) + iu2 = triu_indices(4, k=2) + iu3 = triu_indices(4, m=5) + iu4 = triu_indices(4, k=2, m=5) a = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]]) + b = np.arange(1, 21).reshape(4, 5) # Both for indexing: yield (assert_array_equal, a[iu1], array([1, 2, 3, 4, 6, 7, 8, 11, 12, 16])) + yield (assert_array_equal, b[iu3], + array([1, 2, 3, 4, 5, 7, 8, 9, 10, 13, 14, 15, 19, 20])) # And for assigning values: a[iu1] = -1 @@ -349,6 +395,12 @@ class TestTriuIndices(object): [5, -1, -1, -1], [9, 10, -1, -1], [13, 14, 15, -1]])) + b[iu3] = -1 + yield (assert_array_equal, b, + array([[-1, -1, -1, -1, -1], + [ 6, -1, -1, -1, -1], + [11, 12, -1, -1, -1], + [16, 17, 18, -1, -1]])) # These cover almost the whole array (two diagonals right of the # main one): @@ -358,20 +410,26 @@ class TestTriuIndices(object): [5, -1, -1, -10], [9, 10, -1, -1], [13, 14, 15, -1]])) + b[iu4] = -10 + yield (assert_array_equal, b, + array([[-1, -1, -10, -10, -10], + [6, -1, -1, -10, -10], + [11, 12, -1, -1, -10], + [16, 17, 18, -1, -1]])) class TestTrilIndicesFrom(object): def test_exceptions(self): assert_raises(ValueError, tril_indices_from, np.ones((2,))) assert_raises(ValueError, tril_indices_from, np.ones((2, 2, 2))) - assert_raises(ValueError, tril_indices_from, np.ones((2, 3))) + # assert_raises(ValueError, tril_indices_from, np.ones((2, 3))) class TestTriuIndicesFrom(object): def test_exceptions(self): assert_raises(ValueError, triu_indices_from, np.ones((2,))) assert_raises(ValueError, triu_indices_from, np.ones((2, 2, 2))) - assert_raises(ValueError, triu_indices_from, np.ones((2, 3))) + # assert_raises(ValueError, triu_indices_from, np.ones((2, 3))) class TestVander(object): diff --git a/numpy/lib/tests/test_utils.py b/numpy/lib/tests/test_utils.py index 93c674766..36d5d6428 100644 --- a/numpy/lib/tests/test_utils.py +++ b/numpy/lib/tests/test_utils.py @@ -1,6 +1,7 @@ from __future__ import division, absolute_import, print_function import sys +from numpy.core import arange from numpy.testing import * import numpy.lib.utils as utils from numpy.lib import deprecate @@ -46,9 +47,17 @@ def test_deprecate_fn(): assert_('old_func3' in new_func3.__doc__) assert_('new_func3' in new_func3.__doc__) + def test_safe_eval_nameconstant(): # Test if safe_eval supports Python 3.4 _ast.NameConstant utils.safe_eval('None') + +def test_byte_bounds(): + a = arange(12).reshape(3, 4) + low, high = utils.byte_bounds(a) + assert_equal(high - low, a.size * a.itemsize) + + if __name__ == "__main__": run_module_suite() diff --git a/numpy/lib/twodim_base.py b/numpy/lib/twodim_base.py index 12c0f9bb3..a8925592a 100644 --- a/numpy/lib/twodim_base.py +++ b/numpy/lib/twodim_base.py @@ -11,8 +11,23 @@ __all__ = ['diag', 'diagflat', 'eye', 'fliplr', 'flipud', 'rot90', 'tri', from numpy.core.numeric import ( asanyarray, subtract, arange, zeros, greater_equal, multiply, ones, - asarray, where, + asarray, where, less, int8, int16, int32, int64, empty, promote_types ) +from numpy.core import iinfo + + +i1 = iinfo(int8) +i2 = iinfo(int16) +i4 = iinfo(int32) +def _min_int(low, high): + """ get small int that fits the range """ + if high <= i1.max and low >= i1.min: + return int8 + if high <= i2.max and low >= i2.min: + return int16 + if high <= i4.max and low >= i4.min: + return int32 + return int64 def fliplr(m): @@ -25,7 +40,7 @@ def fliplr(m): Parameters ---------- m : array_like - Input array. + Input array, must be at least 2-D. Returns ------- @@ -40,8 +55,7 @@ def fliplr(m): Notes ----- - Equivalent to A[:,::-1]. Does not require the array to be - two-dimensional. + Equivalent to A[:,::-1]. Requires the array to be at least 2-D. Examples -------- @@ -373,6 +387,7 @@ def tri(N, M=None, k=0, dtype=float): dtype : dtype, optional Data type of the returned array. The default is float. + Returns ------- tri : ndarray of shape (N, M) @@ -394,8 +409,14 @@ def tri(N, M=None, k=0, dtype=float): """ if M is None: M = N - m = greater_equal(subtract.outer(arange(N), arange(M)), -k) - return m.astype(dtype) + + m = greater_equal.outer(arange(N, dtype=_min_int(0, N)), + arange(-k, M-k, dtype=_min_int(-k, M - k))) + + # Avoid making a copy if the requested type is already bool + m = m.astype(dtype, copy=False) + + return m def tril(m, k=0): @@ -431,8 +452,7 @@ def tril(m, k=0): """ m = asanyarray(m) - out = multiply(tri(m.shape[0], m.shape[1], k=k, dtype=m.dtype), m) - return out + return multiply(tri(*m.shape[-2:], k=k, dtype=bool), m, dtype=m.dtype) def triu(m, k=0): @@ -458,22 +478,20 @@ def triu(m, k=0): """ m = asanyarray(m) - out = multiply((1 - tri(m.shape[0], m.shape[1], k - 1, dtype=m.dtype)), m) - return out + return multiply(~tri(*m.shape[-2:], k=k-1, dtype=bool), m, dtype=m.dtype) # Originally borrowed from John Hunter and matplotlib -def vander(x, N=None, order='decreasing'): +def vander(x, N=None, increasing=False): """ Generate a Vandermonde matrix. - The columns of the output matrix are powers of the input vector. The - order of the powers is determined by the `order` argument, either - "decreasing" (the default) or "increasing". Specifically, when - `order` is "decreasing", the `i`-th output column is the input vector - raised element-wise to the power of ``N - i - 1``. Such a matrix with - a geometric progression in each row is named for Alexandre-Theophile - Vandermonde. + The columns of the output matrix are powers of the input vector. The + order of the powers is determined by the `increasing` boolean argument. + Specifically, when `increasing` is False, the `i`-th output column is + the input vector raised element-wise to the power of ``N - i - 1``. Such + a matrix with a geometric progression in each row is named for Alexandre- + Theophile Vandermonde. Parameters ---------- @@ -482,16 +500,18 @@ def vander(x, N=None, order='decreasing'): N : int, optional Number of columns in the output. If `N` is not specified, a square array is returned (``N = len(x)``). - order : str, optional - Order of the powers of the columns. Must be either 'decreasing' - (the default) or 'increasing'. + increasing : bool, optional + Order of the powers of the columns. If True, the powers increase + from left to right, if False (the default) they are reversed. + + .. versionadded:: 1.9.0 Returns ------- out : ndarray - Vandermonde matrix. If `order` is "decreasing", the first column is - ``x^(N-1)``, the second ``x^(N-2)`` and so forth. If `order` is - "increasing", the columns are ``x^0, x^1, ..., x^(N-1)``. + Vandermonde matrix. If `increasing` is False, the first column is + ``x^(N-1)``, the second ``x^(N-2)`` and so forth. If `increasing` is + True, the columns are ``x^0, x^1, ..., x^(N-1)``. See Also -------- @@ -519,7 +539,7 @@ def vander(x, N=None, order='decreasing'): [ 8, 4, 2, 1], [ 27, 9, 3, 1], [125, 25, 5, 1]]) - >>> np.vander(x, order='increasing') + >>> np.vander(x, increasing=True) array([[ 1, 1, 1, 1], [ 1, 2, 4, 8], [ 1, 3, 9, 27], @@ -534,22 +554,22 @@ def vander(x, N=None, order='decreasing'): 48 """ - if order not in ['decreasing', 'increasing']: - raise ValueError("Invalid order %r; order must be either " - "'decreasing' or 'increasing'." % (order,)) x = asarray(x) if x.ndim != 1: raise ValueError("x must be a one-dimensional array or sequence.") if N is None: N = len(x) - if order == "decreasing": - powers = arange(N - 1, -1, -1) - else: - powers = arange(N) - V = x.reshape(-1, 1) ** powers + v = empty((len(x), N), dtype=promote_types(x.dtype, int)) + tmp = v[:, ::-1] if not increasing else v + + if N > 0: + tmp[:, 0] = 1 + if N > 1: + tmp[:, 1:] = x[:, None] + multiply.accumulate(tmp[:, 1:], out=tmp[:, 1:], axis=1) - return V + return v def histogram2d(x, y, bins=10, range=None, normed=False, weights=None): @@ -559,9 +579,11 @@ def histogram2d(x, y, bins=10, range=None, normed=False, weights=None): Parameters ---------- x : array_like, shape (N,) - An array containing the x coordinates of the points to be histogrammed. + An array containing the x coordinates of the points to be + histogrammed. y : array_like, shape (N,) - An array containing the y coordinates of the points to be histogrammed. + An array containing the y coordinates of the points to be + histogrammed. bins : int or [int, int] or array_like or [array, array], optional The bin specification: @@ -579,13 +601,13 @@ def histogram2d(x, y, bins=10, range=None, normed=False, weights=None): ``[[xmin, xmax], [ymin, ymax]]``. All values outside of this range will be considered outliers and not tallied in the histogram. normed : bool, optional - If False, returns the number of samples in each bin. If True, returns - the bin density, i.e. the bin count divided by the bin area. + If False, returns the number of samples in each bin. If True, + returns the bin density ``bin_count / sample_count / bin_area``. weights : array_like, shape(N,), optional - An array of values ``w_i`` weighing each sample ``(x_i, y_i)``. Weights - are normalized to 1 if `normed` is True. If `normed` is False, the - values of the returned histogram are equal to the sum of the weights - belonging to the samples falling into each bin. + An array of values ``w_i`` weighing each sample ``(x_i, y_i)``. + Weights are normalized to 1 if `normed` is True. If `normed` is + False, the values of the returned histogram are equal to the sum of + the weights belonging to the samples falling into each bin. Returns ------- @@ -605,20 +627,15 @@ def histogram2d(x, y, bins=10, range=None, normed=False, weights=None): Notes ----- - When `normed` is True, then the returned histogram is the sample density, - defined such that: - - .. math:: - \\sum_{i=0}^{nx-1} \\sum_{j=0}^{ny-1} H_{i,j} \\Delta x_i \\Delta y_j = 1 - - where `H` is the histogram array and :math:`\\Delta x_i \\Delta y_i` - the area of bin ``{i,j}``. + When `normed` is True, then the returned histogram is the sample + density, defined such that the sum over bins of the product + ``bin_value * bin_area`` is 1. Please note that the histogram does not follow the Cartesian convention - where `x` values are on the abcissa and `y` values on the ordinate axis. - Rather, `x` is histogrammed along the first dimension of the array - (vertical), and `y` along the second dimension of the array (horizontal). - This ensures compatibility with `histogramdd`. + where `x` values are on the abscissa and `y` values on the ordinate + axis. Rather, `x` is histogrammed along the first dimension of the + array (vertical), and `y` along the second dimension of the array + (horizontal). This ensures compatibility with `histogramdd`. Examples -------- @@ -650,14 +667,14 @@ def histogram2d(x, y, bins=10, range=None, normed=False, weights=None): >>> fig = plt.figure(figsize=(7, 3)) >>> ax = fig.add_subplot(131) - >>> ax.set_title('imshow:\nequidistant') + >>> ax.set_title('imshow: equidistant') >>> im = plt.imshow(H, interpolation='nearest', origin='low', - extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]]) + extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]]) pcolormesh can display exact bin edges: >>> ax = fig.add_subplot(132) - >>> ax.set_title('pcolormesh:\nexact bin edges') + >>> ax.set_title('pcolormesh: exact bin edges') >>> X, Y = np.meshgrid(xedges, yedges) >>> ax.pcolormesh(X, Y, H) >>> ax.set_aspect('equal') @@ -665,7 +682,7 @@ def histogram2d(x, y, bins=10, range=None, normed=False, weights=None): NonUniformImage displays exact bin edges with interpolation: >>> ax = fig.add_subplot(133) - >>> ax.set_title('NonUniformImage:\ninterpolated') + >>> ax.set_title('NonUniformImage: interpolated') >>> im = mpl.image.NonUniformImage(ax, interpolation='bilinear') >>> xcenters = xedges[:-1] + 0.5 * (xedges[1:] - xedges[:-1]) >>> ycenters = yedges[:-1] + 0.5 * (yedges[1:] - yedges[:-1]) @@ -761,17 +778,24 @@ def mask_indices(n, mask_func, k=0): return where(a != 0) -def tril_indices(n, k=0): +def tril_indices(n, k=0, m=None): """ - Return the indices for the lower-triangle of an (n, n) array. + Return the indices for the lower-triangle of an (n, m) array. Parameters ---------- n : int - The row dimension of the square arrays for which the returned + The row dimension of the arrays for which the returned indices will be valid. k : int, optional Diagonal offset (see `tril` for details). + m : int, optional + .. versionadded:: 1.9.0 + + The column dimension of the arrays for which the returned + arrays will be valid. + By default `m` is taken equal to `n`. + Returns ------- @@ -831,7 +855,7 @@ def tril_indices(n, k=0): [-10, -10, -10, -10]]) """ - return mask_indices(n, tril, k) + return where(tri(n, m, k=k, dtype=bool)) def tril_indices_from(arr, k=0): @@ -857,14 +881,14 @@ def tril_indices_from(arr, k=0): .. versionadded:: 1.4.0 """ - if not (arr.ndim == 2 and arr.shape[0] == arr.shape[1]): - raise ValueError("input array must be 2-d and square") - return tril_indices(arr.shape[0], k) + if arr.ndim != 2: + raise ValueError("input array must be 2-d") + return tril_indices(arr.shape[-2], k=k, m=arr.shape[-1]) -def triu_indices(n, k=0): +def triu_indices(n, k=0, m=None): """ - Return the indices for the upper-triangle of an (n, n) array. + Return the indices for the upper-triangle of an (n, m) array. Parameters ---------- @@ -873,6 +897,13 @@ def triu_indices(n, k=0): be valid. k : int, optional Diagonal offset (see `triu` for details). + m : int, optional + .. versionadded:: 1.9.0 + + The column dimension of the arrays for which the returned + arrays will be valid. + By default `m` is taken equal to `n`. + Returns ------- @@ -934,12 +965,12 @@ def triu_indices(n, k=0): [ 12, 13, 14, -1]]) """ - return mask_indices(n, triu, k) + return where(~tri(n, m, k=k-1, dtype=bool)) def triu_indices_from(arr, k=0): """ - Return the indices for the upper-triangle of a (N, N) array. + Return the indices for the upper-triangle of arr. See `triu_indices` for full details. @@ -964,6 +995,6 @@ def triu_indices_from(arr, k=0): .. versionadded:: 1.4.0 """ - if not (arr.ndim == 2 and arr.shape[0] == arr.shape[1]): - raise ValueError("input array must be 2-d and square") - return triu_indices(arr.shape[0], k) + if arr.ndim != 2: + raise ValueError("input array must be 2-d") + return triu_indices(arr.shape[-2], k=k, m=arr.shape[-1]) diff --git a/numpy/lib/utils.py b/numpy/lib/utils.py index 1f1cdfc8a..6d41e8eb4 100644 --- a/numpy/lib/utils.py +++ b/numpy/lib/utils.py @@ -6,7 +6,7 @@ import types import re from numpy.core.numerictypes import issubclass_, issubsctype, issubdtype -from numpy.core import product, ndarray, ufunc +from numpy.core import product, ndarray, ufunc, asarray __all__ = [ 'issubclass_', 'issubsctype', 'issubdtype', 'deprecate', @@ -210,8 +210,9 @@ def byte_bounds(a): a_data = ai['data'][0] astrides = ai['strides'] ashape = ai['shape'] - bytes_a = int(ai['typestr'][2:]) - + + bytes_a = asarray(a).dtype.itemsize + a_low = a_high = a_data if astrides is None: # contiguous case a_high += a.size * bytes_a |