diff options
36 files changed, 1073 insertions, 314 deletions
diff --git a/doc/release/1.12.0-notes.rst b/doc/release/1.12.0-notes.rst index 028d77663..ffab3ccd4 100644 --- a/doc/release/1.12.0-notes.rst +++ b/doc/release/1.12.0-notes.rst @@ -18,6 +18,9 @@ Future Changes * In 1.13 NAT will always compare False except for ``NAT != NAT``, which will be True. In short, NAT will behave like NaN +* In 1.13 np.average will preserve subclasses, to match the behavior of most + other numpy functions such as np.mean. In particular, this means calls which + returned a scalar may return a 0-d subclass object instead. Compatibility notes @@ -86,6 +89,21 @@ FutureWarning to changed behavior * ``np.full`` now returns an array of the fill-value's dtype if no dtype is given, instead of defaulting to float. +* np.average will emit a warning if the argument is a subclass of ndarray, + as the subclass will be preserved starting in 1.13. (see Future Changes) + +Greater consistancy in ``assert_almost_equal`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +The precision check for scalars has been changed to match that for arrays. It +is now + + abs(actual - desired) < 1.5 * 10**(-decimal) + +Note that this is looser than previously documented, but agrees with the +previous implementation used in ``assert_array_almost_equal``. Due to the +change in implementation some very delicate tests may fail that did not +fail before. + C API ~~~~~ @@ -101,6 +119,12 @@ keyword argument. It can be set to False when no write operation to the returned array is expected to avoid accidental unpredictable writes. + +``axes`` keyword argument for ``rot90`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +The ``axes`` keyword argument in ``rot90`` determines the plane in which the +array is rotated. It defaults to ``axes=(0,1)`` as in the originial function. + Generalized ``flip`` ~~~~~~~~~~~~~~~~~~~~ ``flipud`` and ``fliplr`` reverse the elements of an array along axis=0 and @@ -128,7 +152,7 @@ file that will remain empty (bar a docstring) in the standard numpy source, but that can be overwritten by people making binary distributions of numpy. New nanfunctions ``nancumsum`` and ``nancumprod`` added -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Nanfunctions ``nancumsum`` and ``nancumprod`` have been added to compute ``cumsum`` and ``cumprod`` by ignoring nans. @@ -137,6 +161,20 @@ compute ``cumsum`` and ``cumprod`` by ignoring nans. ``np.lib.interp(x, xp, fp)`` now allows the interpolated array ``fp`` to be complex and will interpolate at ``complex128`` precision. +New polynomial evaluation function ``polyvalfromroots`` added +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +The new function ``polyvalfromroots`` evaluates a polynomial at given points +from the roots of the polynomial. This is useful for higher order polynomials, +where expansion into polynomial coefficients is inaccurate at machine +precision. + +New array creation function ``geomspace`` added +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +The new function ``geomspace`` generates a geometric sequence. It is similar +to ``logspace``, but with start and stop specified directly: +``geomspace(start, stop)`` behaves the same as +``logspace(log10(start), log10(stop))``. + Improvements ============ @@ -192,6 +230,14 @@ longer grow without bounds. They have been replaced with LRU (least recently used) caches that automatically evict no longer needed items if either the memory size or item count limit has been reached. +Improved handling of zero-width string/unicode dtypes +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Fixed several interfaces that explicitly disallowed arrays with zero-width +string dtypes (i.e. ``dtype('S0')`` or ``dtype('U0')``, and fixed several +bugs where such dtypes were not handled properly. In particular, changed +``ndarray.__new__`` to not implicitly convert ``dtype('S0')`` to +``dtype('S1')`` (and likewise for unicode) when creating new arrays. + Changes ======= @@ -240,4 +286,5 @@ If a 'width' parameter is passed into ``binary_repr`` that is insufficient to represent the number in base 2 (positive) or 2's complement (negative) form, the function used to silently ignore the parameter and return a representation using the minimal number of bits needed for the form in question. Such behavior -is now considered unsafe from a user perspective and will raise an error in the future. +is now considered unsafe from a user perspective and will raise an error in the +future. diff --git a/doc/source/reference/arrays.datetime.rst b/doc/source/reference/arrays.datetime.rst index f5b454875..cbc696ae8 100644 --- a/doc/source/reference/arrays.datetime.rst +++ b/doc/source/reference/arrays.datetime.rst @@ -177,9 +177,9 @@ And here are the time units: ======== ================ ======================= ========================== h hour +/- 1.0e15 years [1.0e15 BC, 1.0e15 AD] m minute +/- 1.7e13 years [1.7e13 BC, 1.7e13 AD] - s second +/- 2.9e12 years [ 2.9e9 BC, 2.9e9 AD] - ms millisecond +/- 2.9e9 years [ 2.9e6 BC, 2.9e6 AD] - us microsecond +/- 2.9e6 years [290301 BC, 294241 AD] + s second +/- 2.9e11 years [2.9e11 BC, 2.9e11 AD] + ms millisecond +/- 2.9e8 years [ 2.9e8 BC, 2.9e8 AD] + us microsecond +/- 2.9e5 years [290301 BC, 294241 AD] ns nanosecond +/- 292 years [ 1678 AD, 2262 AD] ps picosecond +/- 106 days [ 1969 AD, 1970 AD] fs femtosecond +/- 2.6 hours [ 1969 AD, 1970 AD] diff --git a/doc/source/reference/routines.array-creation.rst b/doc/source/reference/routines.array-creation.rst index c7c6ab815..e718f0052 100644 --- a/doc/source/reference/routines.array-creation.rst +++ b/doc/source/reference/routines.array-creation.rst @@ -80,6 +80,7 @@ Numerical ranges arange linspace logspace + geomspace meshgrid mgrid ogrid diff --git a/doc/source/reference/routines.polynomials.polynomial.rst b/doc/source/reference/routines.polynomials.polynomial.rst index 431856622..8194ca867 100644 --- a/doc/source/reference/routines.polynomials.polynomial.rst +++ b/doc/source/reference/routines.polynomials.polynomial.rst @@ -32,6 +32,7 @@ Basics polygrid3d polyroots polyfromroots + polyvalfromroots Fitting ------- diff --git a/doc/source/reference/ufuncs.rst b/doc/source/reference/ufuncs.rst index 46ca1f128..8faf97c4f 100644 --- a/doc/source/reference/ufuncs.rst +++ b/doc/source/reference/ufuncs.rst @@ -521,7 +521,6 @@ Math operations sqrt square reciprocal - ones_like .. tip:: diff --git a/numpy/add_newdocs.py b/numpy/add_newdocs.py index 16799887f..4fc9bca7d 100644 --- a/numpy/add_newdocs.py +++ b/numpy/add_newdocs.py @@ -552,9 +552,26 @@ add_newdoc('numpy.core', 'broadcast', ('iters', """)) +add_newdoc('numpy.core', 'broadcast', ('ndim', + """ + Number of dimensions of broadcasted result. Alias for `nd`. + + .. versionadded:: 1.12.0 + + Examples + -------- + >>> x = np.array([1, 2, 3]) + >>> y = np.array([[4], [5], [6]]) + >>> b = np.broadcast(x, y) + >>> b.ndim + 2 + + """)) + add_newdoc('numpy.core', 'broadcast', ('nd', """ - Number of dimensions of broadcasted result. + Number of dimensions of broadcasted result. For code intended for Numpy + 1.12.0 and later the more consistent `ndim` is preferred. Examples -------- @@ -647,35 +664,43 @@ add_newdoc('numpy.core', 'broadcast', ('reset', add_newdoc('numpy.core.multiarray', 'array', """ - array(object, dtype=None, copy=True, order=None, subok=False, ndmin=0) + array(object, dtype=None, copy=True, order='K', subok=False, ndmin=0) Create an array. Parameters ---------- object : array_like - An array, any object exposing the array interface, an - object whose __array__ method returns an array, or any - (nested) sequence. + An array, any object exposing the array interface, an object whose + __array__ method returns an array, or any (nested) sequence. dtype : data-type, optional - The desired data-type for the array. If not given, then - the type will be determined as the minimum type required - to hold the objects in the sequence. This argument can only - be used to 'upcast' the array. For downcasting, use the - .astype(t) method. + The desired data-type for the array. If not given, then the type will + be determined as the minimum type required to hold the objects in the + sequence. This argument can only be used to 'upcast' the array. For + downcasting, use the .astype(t) method. copy : bool, optional - If true (default), then the object is copied. Otherwise, a copy - will only be made if __array__ returns a copy, if obj is a - nested sequence, or if a copy is needed to satisfy any of the other - requirements (`dtype`, `order`, etc.). - order : {'C', 'F', 'A'}, optional - Specify the order of the array. If order is 'C', then the array - will be in C-contiguous order (last-index varies the fastest). - If order is 'F', then the returned array will be in - Fortran-contiguous order (first-index varies the fastest). - If order is 'A' (default), then the returned array may be - in any order (either C-, Fortran-contiguous, or even discontiguous), - unless a copy is required, in which case it will be C-contiguous. + If true (default), then the object is copied. Otherwise, a copy will + only be made if __array__ returns a copy, if obj is a nested sequence, + or if a copy is needed to satisfy any of the other requirements + (`dtype`, `order`, etc.). + order : {'K', 'A', 'C', 'F'}, optional + Specify the memory layout of the array. If object is not an array, the + newly created array will be in C order (row major) unless 'F' is + specified, in which case it will be in Fortran order (column major). + If object is an array the following holds. + + ===== ========= =================================================== + order no copy copy=True + ===== ========= =================================================== + 'K' unchanged F & C order preserved, otherwise most similar order + 'A' unchanged F order if input is F and not C, otherwise C order + 'C' C order C order + 'F' F order F order + ===== ========= =================================================== + + When ``copy=False`` and a copy is made for other reasons, the result is + the same as if ``copy=True``, with some exceptions for `A`, see the + Notes section. The default order is 'K'. subok : bool, optional If True, then sub-classes will be passed-through, otherwise the returned array will be forced to be a base-class array (default). @@ -693,6 +718,12 @@ add_newdoc('numpy.core.multiarray', 'array', -------- empty, empty_like, zeros, zeros_like, ones, ones_like, full, full_like + Notes + ----- + When order is 'A' and `object` is an array in neither 'C' nor 'F' order, + and a copy is forced by a change in dtype, then the order of the result is + not necessarily 'C' as expected. This is likely a bug. + Examples -------- >>> np.array([1, 2, 3]) @@ -2411,7 +2442,7 @@ add_newdoc('numpy.core.multiarray', 'ndarray', a low-level method (`ndarray(...)`) for instantiating an array. For more information, refer to the `numpy` module and examine the - the methods and attributes of an array. + methods and attributes of an array. Parameters ---------- diff --git a/numpy/core/fromnumeric.py b/numpy/core/fromnumeric.py index 261a0e2fe..28ee4fffa 100644 --- a/numpy/core/fromnumeric.py +++ b/numpy/core/fromnumeric.py @@ -2855,7 +2855,7 @@ def mean(a, axis=None, dtype=None, out=None, keepdims=np._NoValue): >>> a[0, :] = 1.0 >>> a[1, :] = 0.1 >>> np.mean(a) - 0.546875 + 0.54999924 Computing the mean in float64 is more accurate: diff --git a/numpy/core/function_base.py b/numpy/core/function_base.py index 5b1237a32..17a36eb4c 100644 --- a/numpy/core/function_base.py +++ b/numpy/core/function_base.py @@ -3,10 +3,12 @@ from __future__ import division, absolute_import, print_function import warnings import operator -__all__ = ['logspace', 'linspace'] - from . import numeric as _nx -from .numeric import result_type, NaN, shares_memory, MAY_SHARE_BOUNDS, TooHardError +from .numeric import (result_type, NaN, shares_memory, MAY_SHARE_BOUNDS, + TooHardError) + +__all__ = ['logspace', 'linspace', 'geomspace'] + def _index_deprecate(i, stacklevel=2): try: @@ -73,11 +75,11 @@ def linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None): Examples -------- >>> np.linspace(2.0, 3.0, num=5) - array([ 2. , 2.25, 2.5 , 2.75, 3. ]) + array([ 2. , 2.25, 2.5 , 2.75, 3. ]) >>> np.linspace(2.0, 3.0, num=5, endpoint=False) - array([ 2. , 2.2, 2.4, 2.6, 2.8]) + array([ 2. , 2.2, 2.4, 2.6, 2.8]) >>> np.linspace(2.0, 3.0, num=5, retstep=True) - (array([ 2. , 2.25, 2.5 , 2.75, 3. ]), 0.25) + (array([ 2. , 2.25, 2.5 , 2.75, 3. ]), 0.25) Graphical illustration: @@ -102,8 +104,8 @@ def linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None): div = (num - 1) if endpoint else num # Convert float/complex array scalars to float, gh-3504 - start = start * 1. - stop = stop * 1. + start = start * 1.0 + stop = stop * 1.0 dt = result_type(start, stop, float(num)) if dtype is None: @@ -156,7 +158,7 @@ def logspace(start, stop, num=50, endpoint=True, base=10.0, dtype=None): ``base ** stop`` is the final value of the sequence, unless `endpoint` is False. In that case, ``num + 1`` values are spaced over the interval in log-space, of which all but the last (a sequence of - length ``num``) are returned. + length `num`) are returned. num : integer, optional Number of samples to generate. Default is 50. endpoint : boolean, optional @@ -182,6 +184,7 @@ def logspace(start, stop, num=50, endpoint=True, base=10.0, dtype=None): endpoint may or may not be included. linspace : Similar to logspace, but with the samples uniformly distributed in linear space, instead of log space. + geomspace : Similar to logspace, but with endpoints specified directly. Notes ----- @@ -195,11 +198,11 @@ def logspace(start, stop, num=50, endpoint=True, base=10.0, dtype=None): Examples -------- >>> np.logspace(2.0, 3.0, num=4) - array([ 100. , 215.443469 , 464.15888336, 1000. ]) + array([ 100. , 215.443469 , 464.15888336, 1000. ]) >>> np.logspace(2.0, 3.0, num=4, endpoint=False) - array([ 100. , 177.827941 , 316.22776602, 562.34132519]) + array([ 100. , 177.827941 , 316.22776602, 562.34132519]) >>> np.logspace(2.0, 3.0, num=4, base=2.0) - array([ 4. , 5.0396842 , 6.34960421, 8. ]) + array([ 4. , 5.0396842 , 6.34960421, 8. ]) Graphical illustration: @@ -221,3 +224,127 @@ def logspace(start, stop, num=50, endpoint=True, base=10.0, dtype=None): if dtype is None: return _nx.power(base, y) return _nx.power(base, y).astype(dtype) + + +def geomspace(start, stop, num=50, endpoint=True, dtype=None): + """ + Return numbers spaced evenly on a log scale (a geometric progression). + + This is similar to `logspace`, but with endpoints specified directly. + Each output sample is a constant multiple of the previous. + + Parameters + ---------- + start : scalar + The starting value of the sequence. + stop : scalar + The final value of the sequence, unless `endpoint` is False. + In that case, ``num + 1`` values are spaced over the + interval in log-space, of which all but the last (a sequence of + length `num`) are returned. + num : integer, optional + Number of samples to generate. Default is 50. + endpoint : boolean, optional + If true, `stop` is the last sample. Otherwise, it is not included. + Default is True. + dtype : dtype + The type of the output array. If `dtype` is not given, infer the data + type from the other input arguments. + + Returns + ------- + samples : ndarray + `num` samples, equally spaced on a log scale. + + See Also + -------- + logspace : Similar to geomspace, but with endpoints specified using log + and base. + linspace : Similar to geomspace, but with arithmetic instead of geometric + progression. + arange : Similar to linspace, with the step size specified instead of the + number of samples. + + Notes + ----- + If the inputs or dtype are complex, the output will follow a logarithmic + spiral in the complex plane. (There are an infinite number of spirals + passing through two points; the output will follow the shortest such path.) + + Examples + -------- + >>> np.geomspace(1, 1000, num=4) + array([ 1., 10., 100., 1000.]) + >>> np.geomspace(1, 1000, num=3, endpoint=False) + array([ 1., 10., 100.]) + >>> np.geomspace(1, 1000, num=4, endpoint=False) + array([ 1. , 5.62341325, 31.6227766 , 177.827941 ]) + >>> np.geomspace(1, 256, num=9) + array([ 1., 2., 4., 8., 16., 32., 64., 128., 256.]) + + Note that the above may not produce exact integers: + + >>> np.geomspace(1, 256, num=9, dtype=int) + array([ 1, 2, 4, 7, 16, 32, 63, 127, 256]) + >>> np.around(np.geomspace(1, 256, num=9)).astype(int) + array([ 1, 2, 4, 8, 16, 32, 64, 128, 256]) + + Negative, decreasing, and complex inputs are allowed: + + >>> geomspace(1000, 1, num=4) + array([ 1000., 100., 10., 1.]) + >>> geomspace(-1000, -1, num=4) + array([-1000., -100., -10., -1.]) + >>> geomspace(1j, 1000j, num=4) # Straight line + array([ 0. +1.j, 0. +10.j, 0. +100.j, 0.+1000.j]) + >>> geomspace(-1+0j, 1+0j, num=5) # Circle + array([-1.00000000+0.j , -0.70710678+0.70710678j, + 0.00000000+1.j , 0.70710678+0.70710678j, + 1.00000000+0.j ]) + + Graphical illustration of ``endpoint`` parameter: + + >>> import matplotlib.pyplot as plt + >>> N = 10 + >>> y = np.zeros(N) + >>> plt.semilogx(np.geomspace(1, 1000, N, endpoint=True), y + 1, 'o') + >>> plt.semilogx(np.geomspace(1, 1000, N, endpoint=False), y + 2, 'o') + >>> plt.axis([0.5, 2000, 0, 3]) + >>> plt.grid(True, color='0.7', linestyle='-', which='both', axis='both') + >>> plt.show() + + """ + if start == 0 or stop == 0: + raise ValueError('Geometric sequence cannot include zero') + + dt = result_type(start, stop, float(num)) + if dtype is None: + dtype = dt + else: + # complex to dtype('complex128'), for instance + dtype = _nx.dtype(dtype) + + # Avoid negligible real or imaginary parts in output by rotating to + # positive real, calculating, then undoing rotation + out_sign = 1 + if start.real == stop.real == 0: + start, stop = start.imag, stop.imag + out_sign = 1j * out_sign + if _nx.sign(start) == _nx.sign(stop) == -1: + start, stop = -start, -stop + out_sign = -out_sign + + # Promote both arguments to the same dtype in case, for instance, one is + # complex and another is negative and log would produce NaN otherwise + start = start + (stop - stop) + stop = stop + (start - start) + if _nx.issubdtype(dtype, complex): + start = start + 0j + stop = stop + 0j + + log_start = _nx.log10(start) + log_stop = _nx.log10(stop) + result = out_sign * logspace(log_start, log_stop, num=num, + endpoint=endpoint, base=10.0, dtype=dtype) + + return result.astype(dtype) diff --git a/numpy/core/src/multiarray/ctors.c b/numpy/core/src/multiarray/ctors.c index bf060e8a5..a03bacceb 100644 --- a/numpy/core/src/multiarray/ctors.c +++ b/numpy/core/src/multiarray/ctors.c @@ -892,7 +892,8 @@ discover_dimensions(PyObject *obj, int *maxndim, npy_intp *d, int check_it, * Generic new array creation routine. * Internal variant with calloc argument for PyArray_Zeros. * - * steals a reference to descr (even on failure) + * steals a reference to descr. On failure or descr->subarray, descr will + * be decrefed. */ NPY_NO_EXPORT PyObject * PyArray_NewFromDescr_int(PyTypeObject *subtype, PyArray_Descr *descr, int nd, @@ -1128,7 +1129,8 @@ PyArray_NewFromDescr_int(PyTypeObject *subtype, PyArray_Descr *descr, int nd, /*NUMPY_API * Generic new array creation routine. * - * steals a reference to descr (even on failure) + * steals a reference to descr. On failure or when dtype->subarray is + * true, dtype will be decrefed. */ NPY_NO_EXPORT PyObject * PyArray_NewFromDescr(PyTypeObject *subtype, PyArray_Descr *descr, int nd, @@ -1153,7 +1155,8 @@ PyArray_NewFromDescr(PyTypeObject *subtype, PyArray_Descr *descr, int nd, * subok - If 1, use the prototype's array subtype, otherwise * always create a base-class array. * - * NOTE: If dtype is not NULL, steals the dtype reference. + * NOTE: If dtype is not NULL, steals the dtype reference. On failure or when + * dtype->subarray is true, dtype will be decrefed. */ NPY_NO_EXPORT PyObject * PyArray_NewLikeArray(PyArrayObject *prototype, NPY_ORDER order, @@ -2842,7 +2845,8 @@ PyArray_CheckAxis(PyArrayObject *arr, int *axis, int flags) /*NUMPY_API * Zeros * - * steal a reference + * steals a reference to type. On failure or when dtype->subarray is + * true, dtype will be decrefed. * accepts NULL type */ NPY_NO_EXPORT PyObject * @@ -3260,17 +3264,21 @@ array_fromfile_binary(FILE *fp, PyArray_Descr *dtype, npy_intp num, size_t *nrea } num = numbytes / dtype->elsize; } - r = (PyArrayObject *)PyArray_NewFromDescr(&PyArray_Type, - dtype, - 1, &num, - NULL, NULL, - 0, NULL); + /* + * When dtype->subarray is true, PyArray_NewFromDescr will decref dtype + * even on success, so make sure it stays around until exit. + */ + Py_INCREF(dtype); + r = (PyArrayObject *)PyArray_NewFromDescr(&PyArray_Type, dtype, 1, &num, + NULL, NULL, 0, NULL); if (r == NULL) { + Py_DECREF(dtype); return NULL; } NPY_BEGIN_ALLOW_THREADS; *nread = fread(PyArray_DATA(r), dtype->elsize, num, fp); NPY_END_ALLOW_THREADS; + Py_DECREF(dtype); return r; } @@ -3293,13 +3301,17 @@ array_from_text(PyArray_Descr *dtype, npy_intp num, char *sep, size_t *nread, npy_intp bytes, totalbytes; size = (num >= 0) ? num : FROM_BUFFER_SIZE; + + /* + * When dtype->subarray is true, PyArray_NewFromDescr will decref dtype + * even on success, so make sure it stays around until exit. + */ + Py_INCREF(dtype); r = (PyArrayObject *) - PyArray_NewFromDescr(&PyArray_Type, - dtype, - 1, &size, - NULL, NULL, - 0, NULL); + PyArray_NewFromDescr(&PyArray_Type, dtype, 1, &size, + NULL, NULL, 0, NULL); if (r == NULL) { + Py_DECREF(dtype); return NULL; } clean_sep = swab_separator(sep); @@ -3348,6 +3360,7 @@ array_from_text(PyArray_Descr *dtype, npy_intp num, char *sep, size_t *nread, free(clean_sep); fail: + Py_DECREF(dtype); if (err == 1) { PyErr_NoMemory(); } @@ -3365,7 +3378,8 @@ fail: * array corresponding to the data encoded in that file. * * If the dtype is NULL, the default array type is used (double). - * If non-null, the reference is stolen. + * If non-null, the reference is stolen and if dtype->subarray is true dtype + * will be decrefed even on success. * * The number of elements to read is given as ``num``; if it is < 0, then * then as many as possible are read. diff --git a/numpy/core/src/multiarray/descriptor.c b/numpy/core/src/multiarray/descriptor.c index 7f8be5356..fbfda72d7 100644 --- a/numpy/core/src/multiarray/descriptor.c +++ b/numpy/core/src/multiarray/descriptor.c @@ -1539,6 +1539,31 @@ finish: } #endif if (item) { + /* Check for a deprecated Numeric-style typecode */ + if (PyBytes_Check(obj)) { + char *type = NULL; + Py_ssize_t len = 0; + char *dep_tps[] = {"Bool", "Complex", "Float", "Int", + "Object0", "String0", "Timedelta64", + "Unicode0", "UInt", "Void0"}; + int ndep_tps = sizeof(dep_tps) / sizeof(dep_tps[0]); + int i; + + if (PyBytes_AsStringAndSize(obj, &type, &len) < 0) { + goto error; + } + for (i = 0; i < ndep_tps; ++i) { + char *dep_tp = dep_tps[i]; + + if (strncmp(type, dep_tp, strlen(dep_tp)) == 0) { + if (DEPRECATE("Numeric-style type codes are " + "deprecated and will result in " + "an error in the future.") < 0) { + goto fail; + } + } + } + } return PyArray_DescrConverter(item, at); } } diff --git a/numpy/core/src/multiarray/iterators.c b/numpy/core/src/multiarray/iterators.c index 3dd7b0ebb..50f1cb1f4 100644 --- a/numpy/core/src/multiarray/iterators.c +++ b/numpy/core/src/multiarray/iterators.c @@ -1696,6 +1696,10 @@ static PyMemberDef arraymultiter_members[] = { T_INT, offsetof(PyArrayMultiIterObject, nd), READONLY, NULL}, + {"ndim", + T_INT, + offsetof(PyArrayMultiIterObject, nd), + READONLY, NULL}, {NULL, 0, 0, 0, NULL}, }; diff --git a/numpy/core/tests/test_deprecations.py b/numpy/core/tests/test_deprecations.py index e6d3cd261..547280b23 100644 --- a/numpy/core/tests/test_deprecations.py +++ b/numpy/core/tests/test_deprecations.py @@ -198,7 +198,7 @@ class _DeprecationTestCase(object): (warning.category,)) if num is not None and num_found != num: msg = "%i warnings found but %i expected." % (len(self.log), num) - lst = [w.category for w in self.log] + lst = [str(w.category) for w in self.log] raise AssertionError("\n".join([msg] + lst)) with warnings.catch_warnings(): @@ -609,6 +609,29 @@ class TestBinaryReprInsufficientWidthParameterForRepresentation(_DeprecationTest self.assert_deprecated(np.binary_repr, args=args, kwargs=kwargs) +class TestNumericStyleTypecodes(_DeprecationTestCase): + """ + Deprecate the old numeric-style dtypes, which are especially + confusing for complex types, e.g. Complex32 -> complex64. When the + deprecation cycle is complete, the check for the strings should be + removed from PyArray_DescrConverter in descriptor.c, and the + deprecated keys should not be added as capitalized aliases in + _add_aliases in numerictypes.py. + """ + def test_all_dtypes(self): + deprecated_types = [ + 'Bool', 'Complex32', 'Complex64', 'Float16', 'Float32', 'Float64', + 'Int8', 'Int16', 'Int32', 'Int64', 'Object0', 'Timedelta64', + 'UInt8', 'UInt16', 'UInt32', 'UInt64', 'Void0' + ] + if sys.version_info[0] < 3: + deprecated_types.extend(['Unicode0', 'String0']) + + for dt in deprecated_types: + self.assert_deprecated(np.dtype, exceptions=(TypeError,), + args=(dt,)) + + class TestTestDeprecated(object): def test_assert_deprecated(self): test_case_instance = _DeprecationTestCase() diff --git a/numpy/core/tests/test_function_base.py b/numpy/core/tests/test_function_base.py index 6b5430611..0fabb2588 100644 --- a/numpy/core/tests/test_function_base.py +++ b/numpy/core/tests/test_function_base.py @@ -1,13 +1,45 @@ from __future__ import division, absolute_import, print_function -from numpy import (logspace, linspace, dtype, array, finfo, typecodes, arange, - isnan, ndarray) +from numpy import (logspace, linspace, geomspace, dtype, array, finfo, + typecodes, arange, isnan, ndarray, sqrt) from numpy.testing import ( TestCase, run_module_suite, assert_, assert_equal, assert_raises, - assert_array_equal + assert_array_equal, assert_allclose ) +class PhysicalQuantity(float): + def __new__(cls, value): + return float.__new__(cls, value) + + def __add__(self, x): + assert_(isinstance(x, PhysicalQuantity)) + return PhysicalQuantity(float(x) + float(self)) + __radd__ = __add__ + + def __sub__(self, x): + assert_(isinstance(x, PhysicalQuantity)) + return PhysicalQuantity(float(self) - float(x)) + + def __rsub__(self, x): + assert_(isinstance(x, PhysicalQuantity)) + return PhysicalQuantity(float(x) - float(self)) + + def __mul__(self, x): + return PhysicalQuantity(float(x) * float(self)) + __rmul__ = __mul__ + + def __div__(self, x): + return PhysicalQuantity(float(self) / float(x)) + + def __rdiv__(self, x): + return PhysicalQuantity(float(x) / float(self)) + + +class PhysicalQuantity2(ndarray): + __array_priority__ = 10 + + class TestLogspace(TestCase): def test_basic(self): @@ -28,6 +60,136 @@ class TestLogspace(TestCase): y = logspace(0, 6, dtype='int32') assert_equal(y.dtype, dtype('int32')) + def test_physical_quantities(self): + a = PhysicalQuantity(1.0) + b = PhysicalQuantity(5.0) + assert_equal(logspace(a, b), logspace(1.0, 5.0)) + + def test_subclass(self): + a = array(1).view(PhysicalQuantity2) + b = array(7).view(PhysicalQuantity2) + ls = logspace(a, b) + assert type(ls) is PhysicalQuantity2 + assert_equal(ls, logspace(1.0, 7.0)) + ls = logspace(a, b, 1) + assert type(ls) is PhysicalQuantity2 + assert_equal(ls, logspace(1.0, 7.0, 1)) + + +class TestGeomspace(TestCase): + + def test_basic(self): + y = geomspace(1, 1e6) + assert_(len(y) == 50) + y = geomspace(1, 1e6, num=100) + assert_(y[-1] == 10 ** 6) + y = geomspace(1, 1e6, endpoint=False) + assert_(y[-1] < 10 ** 6) + y = geomspace(1, 1e6, num=7) + assert_array_equal(y, [1, 10, 100, 1e3, 1e4, 1e5, 1e6]) + + y = geomspace(8, 2, num=3) + assert_allclose(y, [8, 4, 2]) + assert_array_equal(y.imag, 0) + + y = geomspace(-1, -100, num=3) + assert_array_equal(y, [-1, -10, -100]) + assert_array_equal(y.imag, 0) + + y = geomspace(-100, -1, num=3) + assert_array_equal(y, [-100, -10, -1]) + assert_array_equal(y.imag, 0) + + def test_complex(self): + # Purely imaginary + y = geomspace(1j, 16j, num=5) + assert_allclose(y, [1j, 2j, 4j, 8j, 16j]) + assert_array_equal(y.real, 0) + + y = geomspace(-4j, -324j, num=5) + assert_allclose(y, [-4j, -12j, -36j, -108j, -324j]) + assert_array_equal(y.real, 0) + + y = geomspace(1+1j, 1000+1000j, num=4) + assert_allclose(y, [1+1j, 10+10j, 100+100j, 1000+1000j]) + + y = geomspace(-1+1j, -1000+1000j, num=4) + assert_allclose(y, [-1+1j, -10+10j, -100+100j, -1000+1000j]) + + # Logarithmic spirals + y = geomspace(-1, 1, num=3, dtype=complex) + assert_allclose(y, [-1, 1j, +1]) + + y = geomspace(0+3j, -3+0j, 3) + assert_allclose(y, [0+3j, -3/sqrt(2)+3j/sqrt(2), -3+0j]) + y = geomspace(0+3j, 3+0j, 3) + assert_allclose(y, [0+3j, 3/sqrt(2)+3j/sqrt(2), 3+0j]) + y = geomspace(-3+0j, 0-3j, 3) + assert_allclose(y, [-3+0j, -3/sqrt(2)-3j/sqrt(2), 0-3j]) + y = geomspace(0+3j, -3+0j, 3) + assert_allclose(y, [0+3j, -3/sqrt(2)+3j/sqrt(2), -3+0j]) + y = geomspace(-2-3j, 5+7j, 7) + assert_allclose(y, [-2-3j, -0.29058977-4.15771027j, + 2.08885354-4.34146838j, 4.58345529-3.16355218j, + 6.41401745-0.55233457j, 6.75707386+3.11795092j, + 5+7j]) + + # Type promotion should prevent the -5 from becoming a NaN + y = geomspace(3j, -5, 2) + assert_allclose(y, [3j, -5]) + y = geomspace(-5, 3j, 2) + assert_allclose(y, [-5, 3j]) + + def test_dtype(self): + y = geomspace(1, 1e6, dtype='float32') + assert_equal(y.dtype, dtype('float32')) + y = geomspace(1, 1e6, dtype='float64') + assert_equal(y.dtype, dtype('float64')) + y = geomspace(1, 1e6, dtype='int32') + assert_equal(y.dtype, dtype('int32')) + + # Native types + y = geomspace(1, 1e6, dtype=float) + assert_equal(y.dtype, dtype('float_')) + y = geomspace(1, 1e6, dtype=complex) + assert_equal(y.dtype, dtype('complex')) + + def test_array_scalar(self): + lim1 = array([120, 100], dtype="int8") + lim2 = array([-120, -100], dtype="int8") + lim3 = array([1200, 1000], dtype="uint16") + t1 = geomspace(lim1[0], lim1[1], 5) + t2 = geomspace(lim2[0], lim2[1], 5) + t3 = geomspace(lim3[0], lim3[1], 5) + t4 = geomspace(120.0, 100.0, 5) + t5 = geomspace(-120.0, -100.0, 5) + t6 = geomspace(1200.0, 1000.0, 5) + + # t3 uses float32, t6 uses float64 + assert_allclose(t1, t4, rtol=1e-2) + assert_allclose(t2, t5, rtol=1e-2) + assert_allclose(t3, t6, rtol=1e-5) + + def test_physical_quantities(self): + a = PhysicalQuantity(1.0) + b = PhysicalQuantity(5.0) + assert_equal(geomspace(a, b), geomspace(1.0, 5.0)) + + def test_subclass(self): + a = array(1).view(PhysicalQuantity2) + b = array(7).view(PhysicalQuantity2) + gs = geomspace(a, b) + assert type(gs) is PhysicalQuantity2 + assert_equal(gs, geomspace(1.0, 7.0)) + gs = geomspace(a, b, 1) + assert type(gs) is PhysicalQuantity2 + assert_equal(gs, geomspace(1.0, 7.0, 1)) + + def test_bounds(self): + assert_raises(ValueError, geomspace, 0, 10) + assert_raises(ValueError, geomspace, 10, 0) + assert_raises(ValueError, geomspace, 0, 0) + class TestLinspace(TestCase): @@ -77,48 +239,18 @@ class TestLinspace(TestCase): def test_complex(self): lim1 = linspace(1 + 2j, 3 + 4j, 5) - t1 = array([ 1.0+2.j, 1.5+2.5j, 2.0+3.j, 2.5+3.5j, 3.0+4.j]) + t1 = array([1.0+2.j, 1.5+2.5j, 2.0+3j, 2.5+3.5j, 3.0+4j]) lim2 = linspace(1j, 10, 5) - t2 = array([ 0.0+1.j, 2.5+0.75j, 5.0+0.5j, 7.5+0.25j, 10.0+0.j]) + t2 = array([0.0+1.j, 2.5+0.75j, 5.0+0.5j, 7.5+0.25j, 10.0+0j]) assert_equal(lim1, t1) assert_equal(lim2, t2) def test_physical_quantities(self): - class PhysicalQuantity(float): - def __new__(cls, value): - return float.__new__(cls, value) - - def __add__(self, x): - assert_(isinstance(x, PhysicalQuantity)) - return PhysicalQuantity(float(x) + float(self)) - __radd__ = __add__ - - def __sub__(self, x): - assert_(isinstance(x, PhysicalQuantity)) - return PhysicalQuantity(float(self) - float(x)) - - def __rsub__(self, x): - assert_(isinstance(x, PhysicalQuantity)) - return PhysicalQuantity(float(x) - float(self)) - - def __mul__(self, x): - return PhysicalQuantity(float(x) * float(self)) - __rmul__ = __mul__ - - def __div__(self, x): - return PhysicalQuantity(float(self) / float(x)) - - def __rdiv__(self, x): - return PhysicalQuantity(float(x) / float(self)) - a = PhysicalQuantity(0.0) b = PhysicalQuantity(1.0) assert_equal(linspace(a, b), linspace(0.0, 1.0)) def test_subclass(self): - class PhysicalQuantity2(ndarray): - __array_priority__ = 10 - a = array(0).view(PhysicalQuantity2) b = array(1).view(PhysicalQuantity2) ls = linspace(a, b) diff --git a/numpy/core/tests/test_numeric.py b/numpy/core/tests/test_numeric.py index dd9c83b25..7fdbbc930 100644 --- a/numpy/core/tests/test_numeric.py +++ b/numpy/core/tests/test_numeric.py @@ -2495,6 +2495,7 @@ class TestBroadcast(TestCase): np.broadcast(arrs[0], np.broadcast(*arrs[1:-1]), arrs[-1])] for mit in mits: assert_equal(mit.shape, (5, 6, 7)) + assert_equal(mit.ndim, 3) assert_equal(mit.nd, 3) assert_equal(mit.numiter, 4) for a, ia in zip(arrs, mit.iters): @@ -2505,6 +2506,7 @@ class TestBroadcast(TestCase): arrs = [np.empty((5, 6, 7))] mit = np.broadcast(*arrs) assert_equal(mit.shape, (5, 6, 7)) + assert_equal(mit.ndim, 3) assert_equal(mit.nd, 3) assert_equal(mit.numiter, 1) assert_(arrs[0] is mit.iters[0].base) diff --git a/numpy/ctypeslib.py b/numpy/ctypeslib.py index 38173fba4..36bcf2764 100644 --- a/numpy/ctypeslib.py +++ b/numpy/ctypeslib.py @@ -178,7 +178,7 @@ class _ndptr(_ndptr_base): def _check_retval_(self): """This method is called when this class is used as the .restype - asttribute for a shared-library function. It constructs a numpy + attribute for a shared-library function. It constructs a numpy array from a void pointer.""" return array(self) diff --git a/numpy/distutils/misc_util.py b/numpy/distutils/misc_util.py index 79b6f25f3..8136f8f4f 100644 --- a/numpy/distutils/misc_util.py +++ b/numpy/distutils/misc_util.py @@ -126,13 +126,13 @@ def allpath(name): return os.path.join(*splitted) def rel_path(path, parent_path): - """Return path relative to parent_path. - """ - pd = os.path.abspath(parent_path) - apath = os.path.abspath(path) - if len(apath)<len(pd): + """Return path relative to parent_path.""" + # Use realpath to avoid issues with symlinked dirs (see gh-7707) + pd = os.path.realpath(os.path.abspath(parent_path)) + apath = os.path.realpath(os.path.abspath(path)) + if len(apath) < len(pd): return path - if apath==pd: + if apath == pd: return '' if pd == apath[:len(pd)]: assert apath[len(pd)] in [os.sep], repr((path, apath[len(pd)])) diff --git a/numpy/distutils/system_info.py b/numpy/distutils/system_info.py index 9f5f4e6c0..5ff5041a6 100644 --- a/numpy/distutils/system_info.py +++ b/numpy/distutils/system_info.py @@ -971,7 +971,7 @@ class djbfft_info(system_info): class mkl_info(system_info): section = 'mkl' - dir_env_var = 'MKL' + dir_env_var = 'MKLROOT' _lib_mkl = ['mkl', 'vml', 'guide'] def get_mkl_rootdir(self): diff --git a/numpy/distutils/tests/test_misc_util.py b/numpy/distutils/tests/test_misc_util.py index c50b9480b..3e97b6fe2 100644 --- a/numpy/distutils/tests/test_misc_util.py +++ b/numpy/distutils/tests/test_misc_util.py @@ -4,7 +4,7 @@ from __future__ import division, absolute_import, print_function from os.path import join, sep, dirname from numpy.distutils.misc_util import ( - appendpath, minrelpath, gpaths, get_shared_lib_extension + appendpath, minrelpath, gpaths, get_shared_lib_extension, get_info ) from numpy.testing import ( TestCase, run_module_suite, assert_, assert_equal @@ -75,5 +75,12 @@ class TestSharedExtension(TestCase): # just check for no crash assert_(get_shared_lib_extension(is_python_ext=True)) + +def test_installed_npymath_ini(): + # Regression test for gh-7707. If npymath.ini wasn't installed, then this + # will give an error. + info = get_info('npymath') + + if __name__ == "__main__": run_module_suite() diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py index 93f4f2634..1e44345b0 100644 --- a/numpy/lib/function_base.py +++ b/numpy/lib/function_base.py @@ -7,11 +7,11 @@ import operator import numpy as np import numpy.core.numeric as _nx -from numpy.core import linspace, atleast_1d, atleast_2d +from numpy.core import linspace, atleast_1d, atleast_2d, transpose from numpy.core.numeric import ( ones, zeros, arange, concatenate, array, asarray, asanyarray, empty, empty_like, ndarray, around, floor, ceil, take, dot, where, intp, - integer, isscalar + integer, isscalar, absolute ) from numpy.core.umath import ( pi, multiply, add, arctan2, frompyfunc, cos, less_equal, sqrt, sin, @@ -39,7 +39,7 @@ if sys.version_info[0] < 3: __all__ = [ 'select', 'piecewise', 'trim_zeros', 'copy', 'iterable', 'percentile', 'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'disp', 'flip', - 'extract', 'place', 'vectorize', 'asarray_chkfinite', 'average', + 'rot90', 'extract', 'place', 'vectorize', 'asarray_chkfinite', 'average', 'histogram', 'histogramdd', 'bincount', 'digitize', 'cov', 'corrcoef', 'msort', 'median', 'sinc', 'hamming', 'hanning', 'bartlett', 'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc', 'add_docstring', @@ -47,6 +47,92 @@ __all__ = [ ] +def rot90(m, k=1, axes=(0,1)): + """ + Rotate an array by 90 degrees in the plane specified by axes. + + Rotation direction is from the first towards the second axis. + + .. versionadded:: 1.12.0 + + Parameters + ---------- + m : array_like + Array of two or more dimensions. + k : integer + Number of times the array is rotated by 90 degrees. + axes: (2,) array_like + The array is rotated in the plane defined by the axes. + Axes must be different. + + Returns + ------- + y : ndarray + A rotated view of `m`. + + See Also + -------- + flip : Reverse the order of elements in an array along the given axis. + fliplr : Flip an array horizontally. + flipud : Flip an array vertically. + + Notes + ----- + rot90(m, k=1, axes=(1,0)) is the reverse of rot90(m, k=1, axes=(0,1)) + rot90(m, k=1, axes=(1,0)) is equivalent to rot90(m, k=-1, axes=(0,1)) + + Examples + -------- + >>> m = np.array([[1,2],[3,4]], int) + >>> m + array([[1, 2], + [3, 4]]) + >>> np.rot90(m) + array([[2, 4], + [1, 3]]) + >>> np.rot90(m, 2) + array([[4, 3], + [2, 1]]) + >>> m = np.arange(8).reshape((2,2,2)) + >>> np.rot90(m, 1, (1,2)) + array([[[1, 3], + [0, 2]], + + [[5, 7], + [4, 6]]]) + + """ + axes = tuple(axes) + if len(axes) != 2: + raise ValueError("len(axes) must be 2.") + + m = asanyarray(m) + + if axes[0] == axes[1] or absolute(axes[0] - axes[1]) == m.ndim: + raise ValueError("Axes must be different.") + + if (axes[0] >= m.ndim or axes[0] < -m.ndim + or axes[1] >= m.ndim or axes[1] < -m.ndim): + raise ValueError("Axes={} out of range for array of ndim={}." + .format(axes, m.ndim)) + + k %= 4 + + if k == 0: + return m[:] + if k == 2: + return flip(flip(m, axes[0]), axes[1]) + + axes_list = arange(0, m.ndim) + axes_list[axes[0]], axes_list[axes[1]] = axes_list[axes[1]], axes_list[axes[0]] + + if k == 1: + return transpose(flip(m,axes[1]), axes_list) + else: + # k == 3 + return flip(transpose(m, axes_list), axes[1]) + + def flip(m, axis): """ Reverse the order of elements in an array along the given axis. @@ -1002,7 +1088,19 @@ def average(a, axis=None, weights=None, returned=False): TypeError: Axis must be specified when shapes of a and weights differ. """ - a = np.asanyarray(a) + # 3/19/2016 1.12.0: + # replace the next few lines with "a = np.asanyarray(a)" + if (type(a) not in (np.ndarray, np.matrix) and + issubclass(type(a), np.ndarray)): + warnings.warn("np.average currently does not preserve subclasses, but " + "will do so in the future to match the behavior of most " + "other numpy functions such as np.mean. In particular, " + "this means calls which returned a scalar may return a " + "0-d subclass object instead.", + FutureWarning, stacklevel=2) + + if not isinstance(a, np.matrix): + a = np.asarray(a) if weights is None: avg = a.mean(axis) diff --git a/numpy/lib/npyio.py b/numpy/lib/npyio.py index 426f6e823..9e3127c89 100644 --- a/numpy/lib/npyio.py +++ b/numpy/lib/npyio.py @@ -1041,7 +1041,7 @@ def savetxt(fname, X, fmt='%.18e', delimiter=' ', newline='\n', header='', a) a single specifier, `fmt='%.4e'`, resulting in numbers formatted like `' (%s+%sj)' % (fmt, fmt)` b) a full string specifying every real and imaginary part, e.g. - `' %.4e %+.4j %.4e %+.4j %.4e %+.4j'` for 3 columns + `' %.4e %+.4ej %.4e %+.4ej %.4e %+.4ej'` for 3 columns c) a list of specifiers, one per column - in this case, the real and imaginary part must have separate specifiers, e.g. `['%.3e + %.3ej', '(%.15e%+.15ej)']` for 2 columns diff --git a/numpy/lib/polynomial.py b/numpy/lib/polynomial.py index c0ab80ee8..d96b8969f 100644 --- a/numpy/lib/polynomial.py +++ b/numpy/lib/polynomial.py @@ -15,7 +15,7 @@ import numpy.core.numeric as NX from numpy.core import (isscalar, abs, finfo, atleast_1d, hstack, dot, array, ones) from numpy.lib.twodim_base import diag, vander -from numpy.lib.function_base import trim_zeros, sort_complex +from numpy.lib.function_base import trim_zeros from numpy.lib.type_check import iscomplex, real, imag, mintypecode from numpy.linalg import eigvals, lstsq, inv @@ -145,11 +145,7 @@ def poly(seq_of_zeros): if issubclass(a.dtype.type, NX.complexfloating): # if complex roots are all complex conjugates, the roots are real. roots = NX.asarray(seq_of_zeros, complex) - pos_roots = sort_complex(NX.compress(roots.imag > 0, roots)) - neg_roots = NX.conjugate(sort_complex( - NX.compress(roots.imag < 0, roots))) - if (len(pos_roots) == len(neg_roots) and - NX.alltrue(neg_roots == pos_roots)): + if NX.all(NX.sort(roots) == NX.sort(roots.conjugate())): a = a.real.copy() return a @@ -603,6 +599,9 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False): # it is included here because the covariance of Multivariate Student-T # (which is implied by a Bayesian uncertainty analysis) includes it. # Plus, it gives a slightly more conservative estimate of uncertainty. + if len(x) <= order + 2: + raise ValueError("the number of data points must exceed order + 2 " + "for Bayesian estimate the covariance matrix") fac = resids / (len(x) - order - 2.0) if y.ndim == 1: return c, Vbase * fac diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py index 15fbbfbd7..c2bcc62ba 100644 --- a/numpy/lib/tests/test_function_base.py +++ b/numpy/lib/tests/test_function_base.py @@ -16,7 +16,7 @@ from numpy.lib import ( add_newdoc_ufunc, angle, average, bartlett, blackman, corrcoef, cov, delete, diff, digitize, extract, flipud, gradient, hamming, hanning, histogram, histogramdd, i0, insert, interp, kaiser, meshgrid, msort, - piecewise, place, select, setxor1d, sinc, split, trapz, trim_zeros, + piecewise, place, rot90, select, setxor1d, sinc, split, trapz, trim_zeros, unwrap, unique, vectorize, ) @@ -29,7 +29,76 @@ def get_mat(n): return data +class TestRot90(TestCase): + def test_basic(self): + self.assertRaises(ValueError, rot90, np.ones(4)) + assert_raises(ValueError, rot90, np.ones((2,2,2)), axes=(0,1,2)) + assert_raises(ValueError, rot90, np.ones((2,2)), axes=(0,2)) + assert_raises(ValueError, rot90, np.ones((2,2)), axes=(1,1)) + assert_raises(ValueError, rot90, np.ones((2,2,2)), axes=(-2,1)) + + a = [[0, 1, 2], + [3, 4, 5]] + b1 = [[2, 5], + [1, 4], + [0, 3]] + b2 = [[5, 4, 3], + [2, 1, 0]] + b3 = [[3, 0], + [4, 1], + [5, 2]] + b4 = [[0, 1, 2], + [3, 4, 5]] + + for k in range(-3, 13, 4): + assert_equal(rot90(a, k=k), b1) + for k in range(-2, 13, 4): + assert_equal(rot90(a, k=k), b2) + for k in range(-1, 13, 4): + assert_equal(rot90(a, k=k), b3) + for k in range(0, 13, 4): + assert_equal(rot90(a, k=k), b4) + + assert_equal(rot90(rot90(a, axes=(0,1)), axes=(1,0)), a) + assert_equal(rot90(a, k=1, axes=(1,0)), rot90(a, k=-1, axes=(0,1))) + + def test_axes(self): + a = np.ones((50, 40, 3)) + assert_equal(rot90(a).shape, (40, 50, 3)) + assert_equal(rot90(a, axes=(0,2)), rot90(a, axes=(0,-1))) + assert_equal(rot90(a, axes=(1,2)), rot90(a, axes=(-2,-1))) + + def test_rotation_axes(self): + a = np.arange(8).reshape((2,2,2)) + + a_rot90_01 = [[[2, 3], + [6, 7]], + [[0, 1], + [4, 5]]] + a_rot90_12 = [[[1, 3], + [0, 2]], + [[5, 7], + [4, 6]]] + a_rot90_20 = [[[4, 0], + [6, 2]], + [[5, 1], + [7, 3]]] + a_rot90_10 = [[[4, 5], + [0, 1]], + [[6, 7], + [2, 3]]] + + assert_equal(rot90(a, axes=(0, 1)), a_rot90_01) + assert_equal(rot90(a, axes=(1, 0)), a_rot90_10) + assert_equal(rot90(a, axes=(1, 2)), a_rot90_12) + + for k in range(1,5): + assert_equal(rot90(a, k=k, axes=(2, 0)), + rot90(a_rot90_20, k=k-1, axes=(2, 0))) + + class TestFlip(TestCase): + def test_axes(self): self.assertRaises(ValueError, np.flip, np.ones(4), axis=1) self.assertRaises(ValueError, np.flip, np.ones((4, 4)), axis=2) @@ -58,13 +127,11 @@ class TestFlip(TestCase): def test_3d_swap_axis0(self): a = np.array([[[0, 1], [2, 3]], - [[4, 5], [6, 7]]]) b = np.array([[[4, 5], [6, 7]], - [[0, 1], [2, 3]]]) @@ -73,13 +140,11 @@ class TestFlip(TestCase): def test_3d_swap_axis1(self): a = np.array([[[0, 1], [2, 3]], - [[4, 5], [6, 7]]]) b = np.array([[[2, 3], [0, 1]], - [[6, 7], [4, 5]]]) @@ -87,16 +152,14 @@ class TestFlip(TestCase): def test_3d_swap_axis2(self): a = np.array([[[0, 1], - [2, 3]], - - [[4, 5], - [6, 7]]]) + [2, 3]], + [[4, 5], + [6, 7]]]) b = np.array([[[1, 0], - [3, 2]], - - [[5, 4], - [7, 6]]]) + [3, 2]], + [[5, 4], + [7, 6]]]) assert_equal(np.flip(a, 2), b) diff --git a/numpy/lib/tests/test_polynomial.py b/numpy/lib/tests/test_polynomial.py index 5c15941e6..00dffd3d3 100644 --- a/numpy/lib/tests/test_polynomial.py +++ b/numpy/lib/tests/test_polynomial.py @@ -81,7 +81,7 @@ poly1d([ 2.]) import numpy as np from numpy.testing import ( run_module_suite, TestCase, assert_, assert_equal, assert_array_equal, - assert_almost_equal, rundocs + assert_almost_equal, assert_array_almost_equal, assert_raises, rundocs ) @@ -89,6 +89,30 @@ class TestDocs(TestCase): def test_doctests(self): return rundocs() + def test_poly(self): + assert_array_almost_equal(np.poly([3, -np.sqrt(2), np.sqrt(2)]), + [1, -3, -2, 6]) + + # From matlab docs + A = [[1, 2, 3], [4, 5, 6], [7, 8, 0]] + assert_array_almost_equal(np.poly(A), [1, -6, -72, -27]) + + # Should produce real output for perfect conjugates + assert_(np.isrealobj(np.poly([+1.082j, +2.613j, -2.613j, -1.082j]))) + assert_(np.isrealobj(np.poly([0+1j, -0+-1j, 1+2j, + 1-2j, 1.+3.5j, 1-3.5j]))) + assert_(np.isrealobj(np.poly([1j, -1j, 1+2j, 1-2j, 1+3j, 1-3.j]))) + assert_(np.isrealobj(np.poly([1j, -1j, 1+2j, 1-2j]))) + assert_(np.isrealobj(np.poly([1j, -1j, 2j, -2j]))) + assert_(np.isrealobj(np.poly([1j, -1j]))) + assert_(np.isrealobj(np.poly([1, -1]))) + + assert_(np.iscomplexobj(np.poly([1j, -1.0000001j]))) + + np.random.seed(42) + a = np.random.randn(100) + 1j*np.random.randn(100) + assert_(np.isrealobj(np.poly(np.concatenate((a, np.conjugate(a)))))) + def test_roots(self): assert_array_equal(np.roots([1, 0, 0]), [0, 0]) @@ -111,6 +135,12 @@ class TestDocs(TestCase): err = [1, -1, 1, -1, 1, -1, 1] weights = np.arange(8, 1, -1)**2/7.0 + # Check exception when too few points for variance estimate. Note that + # the Bayesian estimate requires the number of data points to exceed + # degree + 3. + assert_raises(ValueError, np.polyfit, + [0, 1, 3], [0, 1, 3], deg=0, cov=True) + # check 1D case m, cov = np.polyfit(x, y+err, 2, cov=True) est = [3.8571, 0.2857, 1.619] diff --git a/numpy/lib/tests/test_twodim_base.py b/numpy/lib/tests/test_twodim_base.py index b8e7301d7..31925d5fe 100644 --- a/numpy/lib/tests/test_twodim_base.py +++ b/numpy/lib/tests/test_twodim_base.py @@ -9,7 +9,7 @@ from numpy.testing import ( ) from numpy import ( - arange, rot90, add, fliplr, flipud, zeros, ones, eye, array, diag, + arange, add, fliplr, flipud, zeros, ones, eye, array, diag, histogram2d, tri, mask_indices, triu_indices, triu_indices_from, tril_indices, tril_indices_from, vander, ) @@ -169,37 +169,6 @@ class TestFlipud(TestCase): assert_equal(flipud(a), b) -class TestRot90(TestCase): - def test_basic(self): - self.assertRaises(ValueError, rot90, ones(4)) - - a = [[0, 1, 2], - [3, 4, 5]] - b1 = [[2, 5], - [1, 4], - [0, 3]] - b2 = [[5, 4, 3], - [2, 1, 0]] - b3 = [[3, 0], - [4, 1], - [5, 2]] - b4 = [[0, 1, 2], - [3, 4, 5]] - - for k in range(-3, 13, 4): - assert_equal(rot90(a, k=k), b1) - for k in range(-2, 13, 4): - assert_equal(rot90(a, k=k), b2) - for k in range(-1, 13, 4): - assert_equal(rot90(a, k=k), b3) - for k in range(0, 13, 4): - assert_equal(rot90(a, k=k), b4) - - def test_axes(self): - a = ones((50, 40, 3)) - assert_equal(rot90(a).shape, (40, 50, 3)) - - class TestHistogram2d(TestCase): def test_simple(self): x = array( diff --git a/numpy/lib/twodim_base.py b/numpy/lib/twodim_base.py index aefe8d64b..8858f5bad 100644 --- a/numpy/lib/twodim_base.py +++ b/numpy/lib/twodim_base.py @@ -4,14 +4,14 @@ from __future__ import division, absolute_import, print_function from numpy.core.numeric import ( - asanyarray, arange, zeros, greater_equal, multiply, ones, asarray, - where, int8, int16, int32, int64, empty, promote_types, diagonal, + absolute, asanyarray, arange, zeros, greater_equal, multiply, ones, + asarray, where, int8, int16, int32, int64, empty, promote_types, diagonal, ) -from numpy.core import iinfo +from numpy.core import iinfo, transpose __all__ = [ - 'diag', 'diagflat', 'eye', 'fliplr', 'flipud', 'rot90', 'tri', 'triu', + 'diag', 'diagflat', 'eye', 'fliplr', 'flipud', 'tri', 'triu', 'tril', 'vander', 'histogram2d', 'mask_indices', 'tril_indices', 'tril_indices_from', 'triu_indices', 'triu_indices_from', ] @@ -136,59 +136,6 @@ def flipud(m): return m[::-1, ...] -def rot90(m, k=1): - """ - Rotate an array by 90 degrees in the counter-clockwise direction. - - The first two dimensions are rotated; therefore, the array must be at - least 2-D. - - Parameters - ---------- - m : array_like - Array of two or more dimensions. - k : integer - Number of times the array is rotated by 90 degrees. - - Returns - ------- - y : ndarray - Rotated array. - - See Also - -------- - fliplr : Flip an array horizontally. - flipud : Flip an array vertically. - - Examples - -------- - >>> m = np.array([[1,2],[3,4]], int) - >>> m - array([[1, 2], - [3, 4]]) - >>> np.rot90(m) - array([[2, 4], - [1, 3]]) - >>> np.rot90(m, 2) - array([[4, 3], - [2, 1]]) - - """ - m = asanyarray(m) - if m.ndim < 2: - raise ValueError("Input must >= 2-d.") - k = k % 4 - if k == 0: - return m - elif k == 1: - return fliplr(m).swapaxes(0, 1) - elif k == 2: - return fliplr(flipud(m)) - else: - # k == 3 - return fliplr(m.swapaxes(0, 1)) - - def eye(N, M=None, k=0, dtype=float): """ Return a 2-D array with ones on the diagonal and zeros elsewhere. diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py index 1f3ba7c54..c0a84fc4f 100644 --- a/numpy/linalg/linalg.py +++ b/numpy/linalg/linalg.py @@ -1804,8 +1804,9 @@ def lstsq(a, b, rcond=-1): of `b`. rcond : float, optional Cut-off ratio for small singular values of `a`. - Singular values are set to zero if they are smaller than `rcond` - times the largest singular value of `a`. + For the purposes of rank determination, singular values are treated + as zero if they are smaller than `rcond` times the largest singular + value of `a`. Returns ------- diff --git a/numpy/ma/core.py b/numpy/ma/core.py index 360ef694e..29b818c06 100644 --- a/numpy/ma/core.py +++ b/numpy/ma/core.py @@ -2970,25 +2970,28 @@ class MaskedArray(ndarray): else: with np.errstate(divide='ignore', invalid='ignore'): d = filled(domain(*args), True) - # Fill the result where the domain is wrong - try: - # Binary domain: take the last value - fill_value = ufunc_fills[func][-1] - except TypeError: - # Unary domain: just use this one - fill_value = ufunc_fills[func] - except KeyError: - # Domain not recognized, use fill_value instead - fill_value = self.fill_value - result = result.copy() - np.copyto(result, fill_value, where=d) - # Update the mask - if m is nomask: - if d is not nomask: + + if d.any(): + # Fill the result where the domain is wrong + try: + # Binary domain: take the last value + fill_value = ufunc_fills[func][-1] + except TypeError: + # Unary domain: just use this one + fill_value = ufunc_fills[func] + except KeyError: + # Domain not recognized, use fill_value instead + fill_value = self.fill_value + result = result.copy() + np.copyto(result, fill_value, where=d) + + # Update the mask + if m is nomask: m = d - else: - # Don't modify inplace, we risk back-propagation - m = (m | d) + else: + # Don't modify inplace, we risk back-propagation + m = (m | d) + # Make sure the mask has the proper size if result.shape == () and m: return masked diff --git a/numpy/ma/tests/test_core.py b/numpy/ma/tests/test_core.py index 0a3226f27..5c7ae4356 100644 --- a/numpy/ma/tests/test_core.py +++ b/numpy/ma/tests/test_core.py @@ -1210,6 +1210,12 @@ class TestMaskedArrayArithmetic(TestCase): a /= 1. assert_equal(a.mask, [0, 0, 0]) + def test_ufunc_nomask(self): + # check the case ufuncs should set the mask to false + m = np.ma.array([1]) + # check we don't get array([False], dtype=bool) + assert_equal(np.true_divide(m, 5).mask.shape, ()) + def test_noshink_on_creation(self): # Check that the mask is not shrunk on array creation when not wanted a = np.ma.masked_values([1., 2.5, 3.1], 1.5, shrink=False) diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py index 5d05f5991..c310d659d 100644 --- a/numpy/polynomial/polynomial.py +++ b/numpy/polynomial/polynomial.py @@ -36,6 +36,7 @@ Misc Functions -------------- - `polyfromroots` -- create a polynomial with specified roots. - `polyroots` -- find the roots of a polynomial. +- `polyvalfromroots` -- evalute a polynomial at given points from roots. - `polyvander` -- Vandermonde-like matrix for powers. - `polyvander2d` -- Vandermonde-like matrix for 2D power series. - `polyvander3d` -- Vandermonde-like matrix for 3D power series. @@ -58,8 +59,8 @@ from __future__ import division, absolute_import, print_function __all__ = [ 'polyzero', 'polyone', 'polyx', 'polydomain', 'polyline', 'polyadd', 'polysub', 'polymulx', 'polymul', 'polydiv', 'polypow', 'polyval', - 'polyder', 'polyint', 'polyfromroots', 'polyvander', 'polyfit', - 'polytrim', 'polyroots', 'Polynomial', 'polyval2d', 'polyval3d', + 'polyvalfromroots', 'polyder', 'polyint', 'polyfromroots', 'polyvander', + 'polyfit', 'polytrim', 'polyroots', 'Polynomial', 'polyval2d', 'polyval3d', 'polygrid2d', 'polygrid3d', 'polyvander2d', 'polyvander3d'] import warnings @@ -780,6 +781,94 @@ def polyval(x, c, tensor=True): return c0 +def polyvalfromroots(x, r, tensor=True): + """ + Evaluate a polynomial specified by its roots at points x. + + If `r` is of length `N`, this function returns the value + + .. math:: p(x) = \prod_{n=1}^{N} (x - r_n) + + The parameter `x` is converted to an array only if it is a tuple or a + list, otherwise it is treated as a scalar. In either case, either `x` + or its elements must support multiplication and addition both with + themselves and with the elements of `r`. + + If `r` is a 1-D array, then `p(x)` will have the same shape as `x`. If `r` + is multidimensional, then the shape of the result depends on the value of + `tensor`. If `tensor is ``True`` the shape will be r.shape[1:] + x.shape; + that is, each polynomial is evaluated at every value of `x`. If `tensor` is + ``False``, the shape will be r.shape[1:]; that is, each polynomial is + evaluated only for the corresponding broadcast value of `x`. Note that + scalars have shape (,). + + .. versionadded:: 1.12 + + Parameters + ---------- + x : array_like, compatible object + If `x` is a list or tuple, it is converted to an ndarray, otherwise + it is left unchanged and treated as a scalar. In either case, `x` + or its elements must support addition and multiplication with + with themselves and with the elements of `r`. + r : array_like + Array of roots. If `r` is multidimensional the first index is the + root index, while the remaining indices enumerate multiple + polynomials. For instance, in the two dimensional case the roots + of each polynomial may be thought of as stored in the columns of `r`. + tensor : boolean, optional + If True, the shape of the roots array is extended with ones on the + right, one for each dimension of `x`. Scalars have dimension 0 for this + action. The result is that every column of coefficients in `r` is + evaluated for every element of `x`. If False, `x` is broadcast over the + columns of `r` for the evaluation. This keyword is useful when `r` is + multidimensional. The default value is True. + + Returns + ------- + values : ndarray, compatible object + The shape of the returned array is described above. + + See Also + -------- + polyroots, polyfromroots, polyval + + Examples + -------- + >>> from numpy.polynomial.polynomial import polyvalfromroots + >>> polyvalfromroots(1, [1,2,3]) + 0.0 + >>> a = np.arange(4).reshape(2,2) + >>> a + array([[0, 1], + [2, 3]]) + >>> polyvalfromroots(a, [-1, 0, 1]) + array([[ -0., 0.], + [ 6., 24.]]) + >>> r = np.arange(-2, 2).reshape(2,2) # multidimensional coefficients + >>> r # each column of r defines one polynomial + array([[-2, -1], + [ 0, 1]]) + >>> b = [-2, 1] + >>> polyvalfromroots(b, r, tensor=True) + array([[-0., 3.], + [ 3., 0.]]) + >>> polyvalfromroots(b, r, tensor=False) + array([-0., 0.]) + """ + r = np.array(r, ndmin=1, copy=0) + if r.dtype.char in '?bBhHiIlLqQpP': + r = r.astype(np.double) + if isinstance(x, (tuple, list)): + x = np.asarray(x) + if isinstance(x, np.ndarray): + if tensor: + r = r.reshape(r.shape + (1,)*x.ndim) + elif x.ndim >= r.ndim: + raise ValueError("x.ndim must be < r.ndim when tensor == False") + return np.prod(x - r, axis=0) + + def polyval2d(x, y, c): """ Evaluate a 2-D polynomial at points (x, y). diff --git a/numpy/polynomial/tests/test_polynomial.py b/numpy/polynomial/tests/test_polynomial.py index 0e6a2e8a0..037be5927 100644 --- a/numpy/polynomial/tests/test_polynomial.py +++ b/numpy/polynomial/tests/test_polynomial.py @@ -136,6 +136,70 @@ class TestEvaluation(TestCase): assert_equal(poly.polyval(x, [1, 0]).shape, dims) assert_equal(poly.polyval(x, [1, 0, 0]).shape, dims) + def test_polyvalfromroots(self): + # check exception for broadcasting x values over root array with + # too few dimensions + assert_raises(ValueError, poly.polyvalfromroots, + [1], [1], tensor=False) + + # check empty input + assert_equal(poly.polyvalfromroots([], [1]).size, 0) + assert_(poly.polyvalfromroots([], [1]).shape == (0,)) + + # check empty input + multidimensional roots + assert_equal(poly.polyvalfromroots([], [[1] * 5]).size, 0) + assert_(poly.polyvalfromroots([], [[1] * 5]).shape == (5, 0)) + + # check scalar input + assert_equal(poly.polyvalfromroots(1, 1), 0) + assert_(poly.polyvalfromroots(1, np.ones((3, 3))).shape == (3,)) + + # check normal input) + x = np.linspace(-1, 1) + y = [x**i for i in range(5)] + for i in range(1, 5): + tgt = y[i] + res = poly.polyvalfromroots(x, [0]*i) + assert_almost_equal(res, tgt) + tgt = x*(x - 1)*(x + 1) + res = poly.polyvalfromroots(x, [-1, 0, 1]) + assert_almost_equal(res, tgt) + + # check that shape is preserved + for i in range(3): + dims = [2]*i + x = np.zeros(dims) + assert_equal(poly.polyvalfromroots(x, [1]).shape, dims) + assert_equal(poly.polyvalfromroots(x, [1, 0]).shape, dims) + assert_equal(poly.polyvalfromroots(x, [1, 0, 0]).shape, dims) + + # check compatibility with factorization + ptest = [15, 2, -16, -2, 1] + r = poly.polyroots(ptest) + x = np.linspace(-1, 1) + assert_almost_equal(poly.polyval(x, ptest), + poly.polyvalfromroots(x, r)) + + # check multidimensional arrays of roots and values + # check tensor=False + rshape = (3, 5) + x = np.arange(-3, 2) + r = np.random.randint(-5, 5, size=rshape) + res = poly.polyvalfromroots(x, r, tensor=False) + tgt = np.empty(r.shape[1:]) + for ii in range(tgt.size): + tgt[ii] = poly.polyvalfromroots(x[ii], r[:, ii]) + assert_equal(res, tgt) + + # check tensor=True + x = np.vstack([x, 2*x]) + res = poly.polyvalfromroots(x, r, tensor=True) + tgt = np.empty(r.shape[1:] + x.shape) + for ii in range(r.shape[1]): + for jj in range(x.shape[0]): + tgt[ii, jj, :] = poly.polyvalfromroots(x[jj], r[:, ii]) + assert_equal(res, tgt) + def test_polyval2d(self): x1, x2, x3 = self.x y1, y2, y3 = self.y diff --git a/numpy/random/mtrand/distributions.c b/numpy/random/mtrand/distributions.c index 7c44088a7..e195700d4 100644 --- a/numpy/random/mtrand/distributions.c +++ b/numpy/random/mtrand/distributions.c @@ -500,6 +500,11 @@ long rk_poisson_mult(rk_state *state, double lam) } } +/* + * The transformed rejection method for generating Poisson random variables + * W. Hoermann + * Insurance: Mathematics and Economics 12, 39-45 (1993) + */ #define LS2PI 0.91893853320467267 #define TWELFTH 0.083333333333333333333333 long rk_poisson_ptrs(rk_state *state, double lam) diff --git a/numpy/random/mtrand/mtrand.pyx b/numpy/random/mtrand/mtrand.pyx index 44ae956c8..abf7a4102 100644 --- a/numpy/random/mtrand/mtrand.pyx +++ b/numpy/random/mtrand/mtrand.pyx @@ -1903,15 +1903,13 @@ cdef class RandomState: if oloc.shape == oscale.shape == (): floc = PyFloat_AsDouble(loc) fscale = PyFloat_AsDouble(scale) - - if fscale <= 0: - raise ValueError("scale <= 0") - + if np.signbit(fscale): + raise ValueError("scale < 0") return cont2_array_sc(self.internal_state, rk_normal, size, floc, fscale, self.lock) - if np.any(np.less_equal(oscale, 0)): - raise ValueError("scale <= 0") + if np.any(np.signbit(oscale)): + raise ValueError("scale < 0") return cont2_array(self.internal_state, rk_normal, size, oloc, oscale, self.lock) @@ -2029,14 +2027,13 @@ cdef class RandomState: if oscale.shape == (): fscale = PyFloat_AsDouble(scale) - - if fscale <= 0: - raise ValueError("scale <= 0") + if np.signbit(fscale): + raise ValueError("scale < 0") return cont1_array_sc(self.internal_state, rk_exponential, size, fscale, self.lock) - if np.any(np.less_equal(oscale, 0.0)): - raise ValueError("scale <= 0") + if np.any(np.signbit(oscale)): + raise ValueError("scale < 0") return cont1_array(self.internal_state, rk_exponential, size, oscale, self.lock) @@ -2147,14 +2144,13 @@ cdef class RandomState: if oshape.shape == (): fshape = PyFloat_AsDouble(shape) - - if fshape <= 0: - raise ValueError("shape <= 0") + if np.signbit(fshape): + raise ValueError("shape < 0") return cont1_array_sc(self.internal_state, rk_standard_gamma, size, fshape, self.lock) - if np.any(np.less_equal(oshape, 0.0)): - raise ValueError("shape <= 0") + if np.any(np.signbit(oshape)): + raise ValueError("shape < 0") return cont1_array(self.internal_state, rk_standard_gamma, size, oshape, self.lock) @@ -2240,18 +2236,17 @@ cdef class RandomState: if oshape.shape == oscale.shape == (): fshape = PyFloat_AsDouble(shape) fscale = PyFloat_AsDouble(scale) - - if fshape <= 0: - raise ValueError("shape <= 0") - if fscale <= 0: - raise ValueError("scale <= 0") + if np.signbit(fshape): + raise ValueError("shape < 0") + if np.signbit(fscale): + raise ValueError("scale < 0") return cont2_array_sc(self.internal_state, rk_gamma, size, fshape, fscale, self.lock) - if np.any(np.less_equal(oshape, 0.0)): - raise ValueError("shape <= 0") - if np.any(np.less_equal(oscale, 0.0)): - raise ValueError("scale <= 0") + if np.any(np.signbit(oshape)): + raise ValueError("shape < 0") + if np.any(np.signbit(oscale)): + raise ValueError("scale < 0") return cont2_array(self.internal_state, rk_gamma, size, oshape, oscale, self.lock) @@ -3122,14 +3117,13 @@ cdef class RandomState: if oa.shape == (): fa = PyFloat_AsDouble(a) - - if fa <= 0: - raise ValueError("a <= 0") + if np.signbit(fa): + raise ValueError("a < 0") return cont1_array_sc(self.internal_state, rk_weibull, size, fa, self.lock) - if np.any(np.less_equal(oa, 0.0)): - raise ValueError("a <= 0") + if np.any(np.signbit(oa)): + raise ValueError("a < 0") return cont1_array(self.internal_state, rk_weibull, size, oa, self.lock) @@ -3235,14 +3229,13 @@ cdef class RandomState: if oa.shape == (): fa = PyFloat_AsDouble(a) - - if fa <= 0: - raise ValueError("a <= 0") + if np.signbit(fa): + raise ValueError("a < 0") return cont1_array_sc(self.internal_state, rk_power, size, fa, self.lock) - if np.any(np.less_equal(oa, 0.0)): - raise ValueError("a <= 0") + if np.any(np.signbit(oa)): + raise ValueError("a < 0") return cont1_array(self.internal_state, rk_power, size, oa, self.lock) def laplace(self, loc=0.0, scale=1.0, size=None): @@ -3333,14 +3326,13 @@ cdef class RandomState: if oloc.shape == oscale.shape == (): floc = PyFloat_AsDouble(loc) fscale = PyFloat_AsDouble(scale) - - if fscale <= 0: - raise ValueError("scale <= 0") + if np.signbit(fscale): + raise ValueError("scale < 0") return cont2_array_sc(self.internal_state, rk_laplace, size, floc, fscale, self.lock) - if np.any(np.less_equal(oscale, 0.0)): - raise ValueError("scale <= 0") + if np.any(np.signbit(oscale)): + raise ValueError("scale < 0") return cont2_array(self.internal_state, rk_laplace, size, oloc, oscale, self.lock) @@ -3465,14 +3457,13 @@ cdef class RandomState: if oloc.shape == oscale.shape == (): floc = PyFloat_AsDouble(loc) fscale = PyFloat_AsDouble(scale) - - if fscale <= 0: - raise ValueError("scale <= 0") + if np.signbit(fscale): + raise ValueError("scale < 0") return cont2_array_sc(self.internal_state, rk_gumbel, size, floc, fscale, self.lock) - if np.any(np.less_equal(oscale, 0.0)): - raise ValueError("scale <= 0") + if np.any(np.signbit(oscale)): + raise ValueError("scale < 0") return cont2_array(self.internal_state, rk_gumbel, size, oloc, oscale, self.lock) @@ -3559,14 +3550,13 @@ cdef class RandomState: if oloc.shape == oscale.shape == (): floc = PyFloat_AsDouble(loc) fscale = PyFloat_AsDouble(scale) - - if fscale <= 0: - raise ValueError("scale <= 0") + if np.signbit(fscale): + raise ValueError("scale < 0") return cont2_array_sc(self.internal_state, rk_logistic, size, floc, fscale, self.lock) - if np.any(np.less_equal(oscale, 0.0)): - raise ValueError("scale <= 0") + if np.any(np.signbit(oscale)): + raise ValueError("scale < 0") return cont2_array(self.internal_state, rk_logistic, size, oloc, oscale, self.lock) @@ -3684,14 +3674,13 @@ cdef class RandomState: if omean.shape == osigma.shape == (): fmean = PyFloat_AsDouble(mean) fsigma = PyFloat_AsDouble(sigma) - - if fsigma <= 0: - raise ValueError("sigma <= 0") + if np.signbit(fsigma): + raise ValueError("sigma < 0") return cont2_array_sc(self.internal_state, rk_lognormal, size, fmean, fsigma, self.lock) - if np.any(np.less_equal(osigma, 0.0)): - raise ValueError("sigma <= 0.0") + if np.any(np.signbit(osigma)): + raise ValueError("sigma < 0.0") return cont2_array(self.internal_state, rk_lognormal, size, omean, osigma, self.lock) @@ -3764,14 +3753,13 @@ cdef class RandomState: if oscale.shape == (): fscale = PyFloat_AsDouble(scale) - - if fscale <= 0: - raise ValueError("scale <= 0") + if np.signbit(fscale): + raise ValueError("scale < 0") return cont1_array_sc(self.internal_state, rk_rayleigh, size, fscale, self.lock) - if np.any(np.less_equal(oscale, 0.0)): - raise ValueError("scale <= 0.0") + if np.any(np.signbit(oscale)): + raise ValueError("scale < 0.0") return cont1_array(self.internal_state, rk_rayleigh, size, oscale, self.lock) diff --git a/numpy/random/tests/test_random.py b/numpy/random/tests/test_random.py index 013205835..a06de58e3 100644 --- a/numpy/random/tests/test_random.py +++ b/numpy/random/tests/test_random.py @@ -485,6 +485,10 @@ class TestRandomDist(TestCase): [0.68717433461363442, 1.69175666993575979]]) assert_array_almost_equal(actual, desired, decimal=15) + def test_exponential_0(self): + assert_equal(np.random.exponential(scale=0), 0) + assert_raises(ValueError, np.random.exponential, scale=-0.) + def test_f(self): np.random.seed(self.seed) actual = np.random.f(12, 77, size=(3, 2)) @@ -501,6 +505,10 @@ class TestRandomDist(TestCase): [31.71863275789960568, 33.30143302795922011]]) assert_array_almost_equal(actual, desired, decimal=14) + def test_gamma_0(self): + assert_equal(np.random.gamma(shape=0, scale=0), 0) + assert_raises(ValueError, np.random.gamma, shape=-0., scale=-0.) + def test_geometric(self): np.random.seed(self.seed) actual = np.random.geometric(.123456789, size=(3, 2)) @@ -517,6 +525,10 @@ class TestRandomDist(TestCase): [1.10651090478803416, -0.69535848626236174]]) assert_array_almost_equal(actual, desired, decimal=15) + def test_gumbel_0(self): + assert_equal(np.random.gumbel(scale=0), 0) + assert_raises(ValueError, np.random.gumbel, scale=-0.) + def test_hypergeometric(self): np.random.seed(self.seed) actual = np.random.hypergeometric(10.1, 5.5, 14, size=(3, 2)) @@ -551,6 +563,10 @@ class TestRandomDist(TestCase): [-0.05391065675859356, 1.74901336242837324]]) assert_array_almost_equal(actual, desired, decimal=15) + def test_laplace_0(self): + assert_equal(np.random.laplace(scale=0), 0) + assert_raises(ValueError, np.random.laplace, scale=-0.) + def test_logistic(self): np.random.seed(self.seed) actual = np.random.logistic(loc=.123456789, scale=2.0, size=(3, 2)) @@ -559,6 +575,10 @@ class TestRandomDist(TestCase): [-0.21682183359214885, 2.63373365386060332]]) assert_array_almost_equal(actual, desired, decimal=15) + def test_laplace_0(self): + assert_(np.random.laplace(scale=0) in [0, 1]) + assert_raises(ValueError, np.random.laplace, scale=-0.) + def test_lognormal(self): np.random.seed(self.seed) actual = np.random.lognormal(mean=.123456789, sigma=2.0, size=(3, 2)) @@ -567,6 +587,10 @@ class TestRandomDist(TestCase): [65.72798501792723869, 86.84341601437161273]]) assert_array_almost_equal(actual, desired, decimal=13) + def test_lognormal_0(self): + assert_equal(np.random.lognormal(sigma=0), 1) + assert_raises(ValueError, np.random.lognormal, sigma=-0.) + def test_logseries(self): np.random.seed(self.seed) actual = np.random.logseries(p=.923456789, size=(3, 2)) @@ -657,6 +681,10 @@ class TestRandomDist(TestCase): [4.18552478636557357, 4.46410668111310471]]) assert_array_almost_equal(actual, desired, decimal=15) + def test_normal_0(self): + assert_equal(np.random.normal(scale=0), 0) + assert_raises(ValueError, np.random.normal, scale=-0.) + def test_pareto(self): np.random.seed(self.seed) actual = np.random.pareto(a=.123456789, size=(3, 2)) @@ -704,6 +732,10 @@ class TestRandomDist(TestCase): [11.06066537006854311, 17.35468505778271009]]) assert_array_almost_equal(actual, desired, decimal=14) + def test_rayleigh_0(self): + assert_equal(np.random.rayleigh(scale=0), 0) + assert_raises(ValueError, np.random.rayleigh, scale=-0.) + def test_standard_cauchy(self): np.random.seed(self.seed) actual = np.random.standard_cauchy(size=(3, 2)) @@ -728,6 +760,10 @@ class TestRandomDist(TestCase): [7.54838614231317084, 8.012756093271868]]) assert_array_almost_equal(actual, desired, decimal=14) + def test_standard_gamma_0(self): + assert_equal(np.random.standard_gamma(shape=0), 0) + assert_raises(ValueError, np.random.standard_gamma, shape=-0.) + def test_standard_normal(self): np.random.seed(self.seed) actual = np.random.standard_normal(size=(3, 2)) @@ -803,6 +839,10 @@ class TestRandomDist(TestCase): [0.67057783752390987, 1.39494046635066793]]) assert_array_almost_equal(actual, desired, decimal=15) + def test_weibull_0(self): + assert_equal(np.random.weibull(a=0), 0) + assert_raises(ValueError, np.random.weibull, a=-0.) + def test_zipf(self): np.random.seed(self.seed) actual = np.random.zipf(a=1.23, size=(3, 2)) diff --git a/numpy/testing/tests/test_utils.py b/numpy/testing/tests/test_utils.py index fe1f411c4..842d55b37 100644 --- a/numpy/testing/tests/test_utils.py +++ b/numpy/testing/tests/test_utils.py @@ -11,7 +11,7 @@ from numpy.testing import ( assert_warns, assert_no_warnings, assert_allclose, assert_approx_equal, assert_array_almost_equal_nulp, assert_array_max_ulp, clear_and_catch_warnings, run_module_suite, - assert_string_equal, assert_, tempdir, temppath, + assert_string_equal, assert_, tempdir, temppath, ) import unittest @@ -246,6 +246,23 @@ class TestArrayAlmostEqual(_GenericTest, unittest.TestCase): def setUp(self): self._assert_func = assert_array_almost_equal + def test_closeness(self): + # Note that in the course of time we ended up with + # `abs(x - y) < 1.5 * 10**(-decimal)` + # instead of the previously documented + # `abs(x - y) < 0.5 * 10**(-decimal)` + # so this check serves to preserve the wrongness. + + # test scalars + self._assert_func(1.499999, 0.0, decimal=0) + self.assertRaises(AssertionError, + lambda: self._assert_func(1.5, 0.0, decimal=0)) + + # test arrays + self._assert_func([1.499999], [0.0], decimal=0) + self.assertRaises(AssertionError, + lambda: self._assert_func([1.5], [0.0], decimal=0)) + def test_simple(self): x = np.array([1234.2222]) y = np.array([1234.2223]) @@ -288,6 +305,23 @@ class TestAlmostEqual(_GenericTest, unittest.TestCase): def setUp(self): self._assert_func = assert_almost_equal + def test_closeness(self): + # Note that in the course of time we ended up with + # `abs(x - y) < 1.5 * 10**(-decimal)` + # instead of the previously documented + # `abs(x - y) < 0.5 * 10**(-decimal)` + # so this check serves to preserve the wrongness. + + # test scalars + self._assert_func(1.499999, 0.0, decimal=0) + self.assertRaises(AssertionError, + lambda: self._assert_func(1.5, 0.0, decimal=0)) + + # test arrays + self._assert_func([1.499999], [0.0], decimal=0) + self.assertRaises(AssertionError, + lambda: self._assert_func([1.5], [0.0], decimal=0)) + def test_nan_item(self): self._assert_func(np.nan, np.nan) self.assertRaises(AssertionError, diff --git a/numpy/testing/utils.py b/numpy/testing/utils.py index 3f4a8568a..dfed5d148 100644 --- a/numpy/testing/utils.py +++ b/numpy/testing/utils.py @@ -424,11 +424,14 @@ def assert_almost_equal(actual,desired,decimal=7,err_msg='',verbose=True): instead of this function for more consistent floating point comparisons. - The test is equivalent to ``abs(desired-actual) < 0.5 * 10**(-decimal)``. + The test verifies that the elements of ``actual`` and ``desired`` satisfy. - Given two objects (numbers or ndarrays), check that all elements of these - objects are almost equal. An exception is raised at conflicting values. - For ndarrays this delegates to assert_array_almost_equal + ``abs(desired-actual) < 1.5 * 10**(-decimal)`` + + That is a looser test than originally documented, but agrees with what the + actual implementation in `assert_array_almost_equal` did up to rounding + vagaries. An exception is raised at conflicting values. For ndarrays this + delegates to assert_array_almost_equal Parameters ---------- @@ -529,7 +532,7 @@ def assert_almost_equal(actual,desired,decimal=7,err_msg='',verbose=True): return except (NotImplementedError, TypeError): pass - if round(abs(desired - actual), decimal) != 0: + if abs(desired - actual) >= 1.5 * 10.0**(-decimal): raise AssertionError(_build_err_msg()) @@ -819,14 +822,16 @@ def assert_array_almost_equal(x, y, decimal=6, err_msg='', verbose=True): instead of this function for more consistent floating point comparisons. - The test verifies identical shapes and verifies values with - ``abs(desired-actual) < 0.5 * 10**(-decimal)``. + The test verifies identical shapes and that the elements of ``actual`` and + ``desired`` satisfy. - Given two array_like objects, check that the shape is equal and all - elements of these objects are almost equal. An exception is raised at - shape mismatch or conflicting values. In contrast to the standard usage - in numpy, NaNs are compared like numbers, no assertion is raised if - both objects have NaNs in the same positions. + ``abs(desired-actual) < 1.5 * 10**(-decimal)`` + + That is a looser test than originally documented, but agrees with what the + actual implementation did up to rounding vagaries. An exception is raised + at shape mismatch or conflicting values. In contrast to the standard usage + in numpy, NaNs are compared like numbers, no assertion is raised if both + objects have NaNs in the same positions. Parameters ---------- @@ -903,12 +908,12 @@ def assert_array_almost_equal(x, y, decimal=6, err_msg='', verbose=True): # casting of x later. dtype = result_type(y, 1.) y = array(y, dtype=dtype, copy=False, subok=True) - z = abs(x-y) + z = abs(x - y) if not issubdtype(z.dtype, number): z = z.astype(float_) # handle object arrays - return around(z, decimal) <= 10.0**(-decimal) + return z < 1.5 * 10.0**(-decimal) assert_array_compare(compare, x, y, err_msg=err_msg, verbose=verbose, header=('Arrays are not almost equal to %d decimals' % decimal), diff --git a/tools/swig/numpy.i b/tools/swig/numpy.i index 67a519e6d..b8fdaeb1f 100644 --- a/tools/swig/numpy.i +++ b/tools/swig/numpy.i @@ -81,6 +81,7 @@ %#define array_descr(a) (((PyArrayObject*)a)->descr) %#define array_flags(a) (((PyArrayObject*)a)->flags) %#define array_enableflags(a,f) (((PyArrayObject*)a)->flags) = f +%#define array_is_fortran(a) (PyArray_ISFORTRAN((PyArrayObject*)a)) %#else %#define is_array(a) ((a) && PyArray_Check(a)) %#define array_type(a) PyArray_TYPE((PyArrayObject*)a) @@ -93,10 +94,10 @@ %#define array_descr(a) PyArray_DESCR((PyArrayObject*)a) %#define array_flags(a) PyArray_FLAGS((PyArrayObject*)a) %#define array_enableflags(a,f) PyArray_ENABLEFLAGS((PyArrayObject*)a,f) +%#define array_is_fortran(a) (PyArray_IS_F_CONTIGUOUS((PyArrayObject*)a)) %#endif %#define array_is_contiguous(a) (PyArray_ISCONTIGUOUS((PyArrayObject*)a)) %#define array_is_native(a) (PyArray_ISNOTSWAPPED((PyArrayObject*)a)) -%#define array_is_fortran(a) (PyArray_IS_F_CONTIGUOUS((PyArrayObject*)a)) } /**********************************************************************/ @@ -295,7 +296,11 @@ Py_INCREF(array_descr(ary)); result = (PyArrayObject*) PyArray_FromArray(ary, array_descr(ary), +%#if NPY_API_VERSION < 0x00000007 + NPY_FORTRANORDER); +%#else NPY_ARRAY_F_CONTIGUOUS); +%#endif *is_new_object = 1; } return result; |