You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
4816 lines
153 KiB
Python
4816 lines
153 KiB
Python
from __future__ import division, absolute_import, print_function
|
|
|
|
try:
|
|
# Accessing collections abstract classes from collections
|
|
# has been deprecated since Python 3.3
|
|
import collections.abc as collections_abc
|
|
except ImportError:
|
|
import collections as collections_abc
|
|
import functools
|
|
import re
|
|
import sys
|
|
import warnings
|
|
|
|
import numpy as np
|
|
import numpy.core.numeric as _nx
|
|
from numpy.core import atleast_1d, 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, absolute
|
|
)
|
|
from numpy.core.umath import (
|
|
pi, add, arctan2, frompyfunc, cos, less_equal, sqrt, sin,
|
|
mod, exp, not_equal, subtract
|
|
)
|
|
from numpy.core.fromnumeric import (
|
|
ravel, nonzero, partition, mean, any, sum
|
|
)
|
|
from numpy.core.numerictypes import typecodes
|
|
from numpy.core.overrides import set_module
|
|
from numpy.core import overrides
|
|
from numpy.core.function_base import add_newdoc
|
|
from numpy.lib.twodim_base import diag
|
|
from numpy.core.multiarray import (
|
|
_insert, add_docstring, bincount, normalize_axis_index, _monotonicity,
|
|
interp as compiled_interp, interp_complex as compiled_interp_complex
|
|
)
|
|
from numpy.core.umath import _add_newdoc_ufunc as add_newdoc_ufunc
|
|
from numpy.compat import long
|
|
|
|
if sys.version_info[0] < 3:
|
|
# Force range to be a generator, for np.delete's usage.
|
|
range = xrange
|
|
import __builtin__ as builtins
|
|
else:
|
|
import builtins
|
|
|
|
|
|
array_function_dispatch = functools.partial(
|
|
overrides.array_function_dispatch, module='numpy')
|
|
|
|
|
|
# needed in this module for compatibility
|
|
from numpy.lib.histograms import histogram, histogramdd
|
|
|
|
__all__ = [
|
|
'select', 'piecewise', 'trim_zeros', 'copy', 'iterable', 'percentile',
|
|
'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'disp', 'flip',
|
|
'rot90', 'extract', 'place', 'vectorize', 'asarray_chkfinite', 'average',
|
|
'bincount', 'digitize', 'cov', 'corrcoef',
|
|
'msort', 'median', 'sinc', 'hamming', 'hanning', 'bartlett',
|
|
'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc', 'add_docstring',
|
|
'meshgrid', 'delete', 'insert', 'append', 'interp', 'add_newdoc_ufunc',
|
|
'quantile'
|
|
]
|
|
|
|
|
|
def _rot90_dispatcher(m, k=None, axes=None):
|
|
return (m,)
|
|
|
|
|
|
@array_function_dispatch(_rot90_dispatcher)
|
|
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.
|
|
|
|
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.
|
|
|
|
.. versionadded:: 1.12.0
|
|
|
|
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_dispatcher(m, axis=None):
|
|
return (m,)
|
|
|
|
|
|
@array_function_dispatch(_flip_dispatcher)
|
|
def flip(m, axis=None):
|
|
"""
|
|
Reverse the order of elements in an array along the given axis.
|
|
|
|
The shape of the array is preserved, but the elements are reordered.
|
|
|
|
.. versionadded:: 1.12.0
|
|
|
|
Parameters
|
|
----------
|
|
m : array_like
|
|
Input array.
|
|
axis : None or int or tuple of ints, optional
|
|
Axis or axes along which to flip over. The default,
|
|
axis=None, will flip over all of the axes of the input array.
|
|
If axis is negative it counts from the last to the first axis.
|
|
|
|
If axis is a tuple of ints, flipping is performed on all of the axes
|
|
specified in the tuple.
|
|
|
|
.. versionchanged:: 1.15.0
|
|
None and tuples of axes are supported
|
|
|
|
Returns
|
|
-------
|
|
out : array_like
|
|
A view of `m` with the entries of axis reversed. Since a view is
|
|
returned, this operation is done in constant time.
|
|
|
|
See Also
|
|
--------
|
|
flipud : Flip an array vertically (axis=0).
|
|
fliplr : Flip an array horizontally (axis=1).
|
|
|
|
Notes
|
|
-----
|
|
flip(m, 0) is equivalent to flipud(m).
|
|
|
|
flip(m, 1) is equivalent to fliplr(m).
|
|
|
|
flip(m, n) corresponds to ``m[...,::-1,...]`` with ``::-1`` at position n.
|
|
|
|
flip(m) corresponds to ``m[::-1,::-1,...,::-1]`` with ``::-1`` at all
|
|
positions.
|
|
|
|
flip(m, (0, 1)) corresponds to ``m[::-1,::-1,...]`` with ``::-1`` at
|
|
position 0 and position 1.
|
|
|
|
Examples
|
|
--------
|
|
>>> A = np.arange(8).reshape((2,2,2))
|
|
>>> A
|
|
array([[[0, 1],
|
|
[2, 3]],
|
|
[[4, 5],
|
|
[6, 7]]])
|
|
>>> np.flip(A, 0)
|
|
array([[[4, 5],
|
|
[6, 7]],
|
|
[[0, 1],
|
|
[2, 3]]])
|
|
>>> np.flip(A, 1)
|
|
array([[[2, 3],
|
|
[0, 1]],
|
|
[[6, 7],
|
|
[4, 5]]])
|
|
>>> np.flip(A)
|
|
array([[[7, 6],
|
|
[5, 4]],
|
|
[[3, 2],
|
|
[1, 0]]])
|
|
>>> np.flip(A, (0, 2))
|
|
array([[[5, 4],
|
|
[7, 6]],
|
|
[[1, 0],
|
|
[3, 2]]])
|
|
>>> A = np.random.randn(3,4,5)
|
|
>>> np.all(np.flip(A,2) == A[:,:,::-1,...])
|
|
True
|
|
"""
|
|
if not hasattr(m, 'ndim'):
|
|
m = asarray(m)
|
|
if axis is None:
|
|
indexer = (np.s_[::-1],) * m.ndim
|
|
else:
|
|
axis = _nx.normalize_axis_tuple(axis, m.ndim)
|
|
indexer = [np.s_[:]] * m.ndim
|
|
for ax in axis:
|
|
indexer[ax] = np.s_[::-1]
|
|
indexer = tuple(indexer)
|
|
return m[indexer]
|
|
|
|
|
|
@set_module('numpy')
|
|
def iterable(y):
|
|
"""
|
|
Check whether or not an object can be iterated over.
|
|
|
|
Parameters
|
|
----------
|
|
y : object
|
|
Input object.
|
|
|
|
Returns
|
|
-------
|
|
b : bool
|
|
Return ``True`` if the object has an iterator method or is a
|
|
sequence and ``False`` otherwise.
|
|
|
|
|
|
Examples
|
|
--------
|
|
>>> np.iterable([1, 2, 3])
|
|
True
|
|
>>> np.iterable(2)
|
|
False
|
|
|
|
"""
|
|
try:
|
|
iter(y)
|
|
except TypeError:
|
|
return False
|
|
return True
|
|
|
|
|
|
def _average_dispatcher(a, axis=None, weights=None, returned=None):
|
|
return (a, weights)
|
|
|
|
|
|
@array_function_dispatch(_average_dispatcher)
|
|
def average(a, axis=None, weights=None, returned=False):
|
|
"""
|
|
Compute the weighted average along the specified axis.
|
|
|
|
Parameters
|
|
----------
|
|
a : array_like
|
|
Array containing data to be averaged. If `a` is not an array, a
|
|
conversion is attempted.
|
|
axis : None or int or tuple of ints, optional
|
|
Axis or axes along which to average `a`. The default,
|
|
axis=None, will average over all of the elements of the input array.
|
|
If axis is negative it counts from the last to the first axis.
|
|
|
|
.. versionadded:: 1.7.0
|
|
|
|
If axis is a tuple of ints, averaging is performed on all of the axes
|
|
specified in the tuple instead of a single axis or all the axes as
|
|
before.
|
|
weights : array_like, optional
|
|
An array of weights associated with the values in `a`. Each value in
|
|
`a` contributes to the average according to its associated weight.
|
|
The weights array can either be 1-D (in which case its length must be
|
|
the size of `a` along the given axis) or of the same shape as `a`.
|
|
If `weights=None`, then all data in `a` are assumed to have a
|
|
weight equal to one.
|
|
returned : bool, optional
|
|
Default is `False`. If `True`, the tuple (`average`, `sum_of_weights`)
|
|
is returned, otherwise only the average is returned.
|
|
If `weights=None`, `sum_of_weights` is equivalent to the number of
|
|
elements over which the average is taken.
|
|
|
|
|
|
Returns
|
|
-------
|
|
retval, [sum_of_weights] : array_type or double
|
|
Return the average along the specified axis. When `returned` is `True`,
|
|
return a tuple with the average as the first element and the sum
|
|
of the weights as the second element. `sum_of_weights` is of the
|
|
same type as `retval`. The result dtype follows a genereal pattern.
|
|
If `weights` is None, the result dtype will be that of `a` , or ``float64``
|
|
if `a` is integral. Otherwise, if `weights` is not None and `a` is non-
|
|
integral, the result type will be the type of lowest precision capable of
|
|
representing values of both `a` and `weights`. If `a` happens to be
|
|
integral, the previous rules still applies but the result dtype will
|
|
at least be ``float64``.
|
|
|
|
Raises
|
|
------
|
|
ZeroDivisionError
|
|
When all weights along axis are zero. See `numpy.ma.average` for a
|
|
version robust to this type of error.
|
|
TypeError
|
|
When the length of 1D `weights` is not the same as the shape of `a`
|
|
along axis.
|
|
|
|
See Also
|
|
--------
|
|
mean
|
|
|
|
ma.average : average for masked arrays -- useful if your data contains
|
|
"missing" values
|
|
numpy.result_type : Returns the type that results from applying the
|
|
numpy type promotion rules to the arguments.
|
|
|
|
Examples
|
|
--------
|
|
>>> data = list(range(1,5))
|
|
>>> data
|
|
[1, 2, 3, 4]
|
|
>>> np.average(data)
|
|
2.5
|
|
>>> np.average(range(1,11), weights=range(10,0,-1))
|
|
4.0
|
|
|
|
>>> data = np.arange(6).reshape((3,2))
|
|
>>> data
|
|
array([[0, 1],
|
|
[2, 3],
|
|
[4, 5]])
|
|
>>> np.average(data, axis=1, weights=[1./4, 3./4])
|
|
array([0.75, 2.75, 4.75])
|
|
>>> np.average(data, weights=[1./4, 3./4])
|
|
Traceback (most recent call last):
|
|
...
|
|
TypeError: Axis must be specified when shapes of a and weights differ.
|
|
|
|
>>> a = np.ones(5, dtype=np.float128)
|
|
>>> w = np.ones(5, dtype=np.complex64)
|
|
>>> avg = np.average(a, weights=w)
|
|
>>> print(avg.dtype)
|
|
complex256
|
|
"""
|
|
a = np.asanyarray(a)
|
|
|
|
if weights is None:
|
|
avg = a.mean(axis)
|
|
scl = avg.dtype.type(a.size/avg.size)
|
|
else:
|
|
wgt = np.asanyarray(weights)
|
|
|
|
if issubclass(a.dtype.type, (np.integer, np.bool_)):
|
|
result_dtype = np.result_type(a.dtype, wgt.dtype, 'f8')
|
|
else:
|
|
result_dtype = np.result_type(a.dtype, wgt.dtype)
|
|
|
|
# Sanity checks
|
|
if a.shape != wgt.shape:
|
|
if axis is None:
|
|
raise TypeError(
|
|
"Axis must be specified when shapes of a and weights "
|
|
"differ.")
|
|
if wgt.ndim != 1:
|
|
raise TypeError(
|
|
"1D weights expected when shapes of a and weights differ.")
|
|
if wgt.shape[0] != a.shape[axis]:
|
|
raise ValueError(
|
|
"Length of weights not compatible with specified axis.")
|
|
|
|
# setup wgt to broadcast along axis
|
|
wgt = np.broadcast_to(wgt, (a.ndim-1)*(1,) + wgt.shape)
|
|
wgt = wgt.swapaxes(-1, axis)
|
|
|
|
scl = wgt.sum(axis=axis, dtype=result_dtype)
|
|
if np.any(scl == 0.0):
|
|
raise ZeroDivisionError(
|
|
"Weights sum to zero, can't be normalized")
|
|
|
|
avg = np.multiply(a, wgt, dtype=result_dtype).sum(axis)/scl
|
|
|
|
if returned:
|
|
if scl.shape != avg.shape:
|
|
scl = np.broadcast_to(scl, avg.shape).copy()
|
|
return avg, scl
|
|
else:
|
|
return avg
|
|
|
|
|
|
@set_module('numpy')
|
|
def asarray_chkfinite(a, dtype=None, order=None):
|
|
"""Convert the input to an array, checking for NaNs or Infs.
|
|
|
|
Parameters
|
|
----------
|
|
a : array_like
|
|
Input data, in any form that can be converted to an array. This
|
|
includes lists, lists of tuples, tuples, tuples of tuples, tuples
|
|
of lists and ndarrays. Success requires no NaNs or Infs.
|
|
dtype : data-type, optional
|
|
By default, the data-type is inferred from the input data.
|
|
order : {'C', 'F'}, optional
|
|
Whether to use row-major (C-style) or
|
|
column-major (Fortran-style) memory representation.
|
|
Defaults to 'C'.
|
|
|
|
Returns
|
|
-------
|
|
out : ndarray
|
|
Array interpretation of `a`. No copy is performed if the input
|
|
is already an ndarray. If `a` is a subclass of ndarray, a base
|
|
class ndarray is returned.
|
|
|
|
Raises
|
|
------
|
|
ValueError
|
|
Raises ValueError if `a` contains NaN (Not a Number) or Inf (Infinity).
|
|
|
|
See Also
|
|
--------
|
|
asarray : Create and array.
|
|
asanyarray : Similar function which passes through subclasses.
|
|
ascontiguousarray : Convert input to a contiguous array.
|
|
asfarray : Convert input to a floating point ndarray.
|
|
asfortranarray : Convert input to an ndarray with column-major
|
|
memory order.
|
|
fromiter : Create an array from an iterator.
|
|
fromfunction : Construct an array by executing a function on grid
|
|
positions.
|
|
|
|
Examples
|
|
--------
|
|
Convert a list into an array. If all elements are finite
|
|
``asarray_chkfinite`` is identical to ``asarray``.
|
|
|
|
>>> a = [1, 2]
|
|
>>> np.asarray_chkfinite(a, dtype=float)
|
|
array([1., 2.])
|
|
|
|
Raises ValueError if array_like contains Nans or Infs.
|
|
|
|
>>> a = [1, 2, np.inf]
|
|
>>> try:
|
|
... np.asarray_chkfinite(a)
|
|
... except ValueError:
|
|
... print('ValueError')
|
|
...
|
|
ValueError
|
|
|
|
"""
|
|
a = asarray(a, dtype=dtype, order=order)
|
|
if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
|
|
raise ValueError(
|
|
"array must not contain infs or NaNs")
|
|
return a
|
|
|
|
|
|
def _piecewise_dispatcher(x, condlist, funclist, *args, **kw):
|
|
yield x
|
|
# support the undocumented behavior of allowing scalars
|
|
if np.iterable(condlist):
|
|
for c in condlist:
|
|
yield c
|
|
|
|
|
|
@array_function_dispatch(_piecewise_dispatcher)
|
|
def piecewise(x, condlist, funclist, *args, **kw):
|
|
"""
|
|
Evaluate a piecewise-defined function.
|
|
|
|
Given a set of conditions and corresponding functions, evaluate each
|
|
function on the input data wherever its condition is true.
|
|
|
|
Parameters
|
|
----------
|
|
x : ndarray or scalar
|
|
The input domain.
|
|
condlist : list of bool arrays or bool scalars
|
|
Each boolean array corresponds to a function in `funclist`. Wherever
|
|
`condlist[i]` is True, `funclist[i](x)` is used as the output value.
|
|
|
|
Each boolean array in `condlist` selects a piece of `x`,
|
|
and should therefore be of the same shape as `x`.
|
|
|
|
The length of `condlist` must correspond to that of `funclist`.
|
|
If one extra function is given, i.e. if
|
|
``len(funclist) == len(condlist) + 1``, then that extra function
|
|
is the default value, used wherever all conditions are false.
|
|
funclist : list of callables, f(x,*args,**kw), or scalars
|
|
Each function is evaluated over `x` wherever its corresponding
|
|
condition is True. It should take a 1d array as input and give an 1d
|
|
array or a scalar value as output. If, instead of a callable,
|
|
a scalar is provided then a constant function (``lambda x: scalar``) is
|
|
assumed.
|
|
args : tuple, optional
|
|
Any further arguments given to `piecewise` are passed to the functions
|
|
upon execution, i.e., if called ``piecewise(..., ..., 1, 'a')``, then
|
|
each function is called as ``f(x, 1, 'a')``.
|
|
kw : dict, optional
|
|
Keyword arguments used in calling `piecewise` are passed to the
|
|
functions upon execution, i.e., if called
|
|
``piecewise(..., ..., alpha=1)``, then each function is called as
|
|
``f(x, alpha=1)``.
|
|
|
|
Returns
|
|
-------
|
|
out : ndarray
|
|
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 a default value of 0.
|
|
|
|
|
|
See Also
|
|
--------
|
|
choose, select, where
|
|
|
|
Notes
|
|
-----
|
|
This is similar to choose or select, except that functions are
|
|
evaluated on elements of `x` that satisfy the corresponding condition from
|
|
`condlist`.
|
|
|
|
The result is::
|
|
|
|
|--
|
|
|funclist[0](x[condlist[0]])
|
|
out = |funclist[1](x[condlist[1]])
|
|
|...
|
|
|funclist[n2](x[condlist[n2]])
|
|
|--
|
|
|
|
Examples
|
|
--------
|
|
Define the sigma function, which is -1 for ``x < 0`` and +1 for ``x >= 0``.
|
|
|
|
>>> x = np.linspace(-2.5, 2.5, 6)
|
|
>>> np.piecewise(x, [x < 0, x >= 0], [-1, 1])
|
|
array([-1., -1., -1., 1., 1., 1.])
|
|
|
|
Define the absolute value, which is ``-x`` for ``x <0`` and ``x`` for
|
|
``x >= 0``.
|
|
|
|
>>> np.piecewise(x, [x < 0, x >= 0], [lambda x: -x, lambda x: x])
|
|
array([2.5, 1.5, 0.5, 0.5, 1.5, 2.5])
|
|
|
|
Apply the same function to a scalar value.
|
|
|
|
>>> y = -2
|
|
>>> np.piecewise(y, [y < 0, y >= 0], [lambda x: -x, lambda x: x])
|
|
array(2)
|
|
|
|
"""
|
|
x = asanyarray(x)
|
|
n2 = len(funclist)
|
|
|
|
# undocumented: single condition is promoted to a list of one condition
|
|
if isscalar(condlist) or (
|
|
not isinstance(condlist[0], (list, ndarray)) and x.ndim != 0):
|
|
condlist = [condlist]
|
|
|
|
condlist = array(condlist, dtype=bool)
|
|
n = len(condlist)
|
|
|
|
if n == n2 - 1: # compute the "otherwise" condition.
|
|
condelse = ~np.any(condlist, axis=0, keepdims=True)
|
|
condlist = np.concatenate([condlist, condelse], axis=0)
|
|
n += 1
|
|
elif n != n2:
|
|
raise ValueError(
|
|
"with {} condition(s), either {} or {} functions are expected"
|
|
.format(n, n, n+1)
|
|
)
|
|
|
|
y = zeros(x.shape, x.dtype)
|
|
for k in range(n):
|
|
item = funclist[k]
|
|
if not isinstance(item, collections_abc.Callable):
|
|
y[condlist[k]] = item
|
|
else:
|
|
vals = x[condlist[k]]
|
|
if vals.size > 0:
|
|
y[condlist[k]] = item(vals, *args, **kw)
|
|
|
|
return y
|
|
|
|
|
|
def _select_dispatcher(condlist, choicelist, default=None):
|
|
for c in condlist:
|
|
yield c
|
|
for c in choicelist:
|
|
yield c
|
|
|
|
|
|
@array_function_dispatch(_select_dispatcher)
|
|
def select(condlist, choicelist, default=0):
|
|
"""
|
|
Return an array drawn from elements in choicelist, depending on conditions.
|
|
|
|
Parameters
|
|
----------
|
|
condlist : list of bool ndarrays
|
|
The list of conditions which determine from which array in `choicelist`
|
|
the output elements are taken. When multiple conditions are satisfied,
|
|
the first one encountered in `condlist` is used.
|
|
choicelist : list of ndarrays
|
|
The list of arrays from which the output elements are taken. It has
|
|
to be of the same length as `condlist`.
|
|
default : scalar, optional
|
|
The element inserted in `output` when all conditions evaluate to False.
|
|
|
|
Returns
|
|
-------
|
|
output : ndarray
|
|
The output at position m is the m-th element of the array in
|
|
`choicelist` where the m-th element of the corresponding array in
|
|
`condlist` is True.
|
|
|
|
See Also
|
|
--------
|
|
where : Return elements from one of two arrays depending on condition.
|
|
take, choose, compress, diag, diagonal
|
|
|
|
Examples
|
|
--------
|
|
>>> x = np.arange(10)
|
|
>>> condlist = [x<3, x>5]
|
|
>>> choicelist = [x, x**2]
|
|
>>> np.select(condlist, choicelist)
|
|
array([ 0, 1, 2, ..., 49, 64, 81])
|
|
|
|
"""
|
|
# 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')
|
|
|
|
# Now that the dtype is known, handle the deprecated select([], []) case
|
|
if len(condlist) == 0:
|
|
# 2014-02-24, 1.9
|
|
warnings.warn("select with an empty condition list is not possible"
|
|
"and will be deprecated",
|
|
DeprecationWarning, stacklevel=3)
|
|
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 separately 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 condlist: should be boolean ndarray'.format(i))
|
|
|
|
if deprecated_ints:
|
|
# 2014-02-24, 1.9
|
|
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, stacklevel=3)
|
|
|
|
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_dispatcher(a, order=None):
|
|
return (a,)
|
|
|
|
|
|
@array_function_dispatch(_copy_dispatcher)
|
|
def copy(a, order='K'):
|
|
"""
|
|
Return an array copy of the given object.
|
|
|
|
Parameters
|
|
----------
|
|
a : array_like
|
|
Input data.
|
|
order : {'C', 'F', 'A', 'K'}, optional
|
|
Controls the memory layout of the copy. 'C' means C-order,
|
|
'F' means F-order, 'A' means 'F' if `a` is Fortran contiguous,
|
|
'C' otherwise. 'K' means match the layout of `a` as closely
|
|
as possible. (Note that this function and :meth:`ndarray.copy` are very
|
|
similar, but have different default values for their order=
|
|
arguments.)
|
|
|
|
Returns
|
|
-------
|
|
arr : ndarray
|
|
Array interpretation of `a`.
|
|
|
|
Notes
|
|
-----
|
|
This is equivalent to:
|
|
|
|
>>> np.array(a, copy=True) #doctest: +SKIP
|
|
|
|
Examples
|
|
--------
|
|
Create an array x, with a reference y and a copy z:
|
|
|
|
>>> x = np.array([1, 2, 3])
|
|
>>> y = x
|
|
>>> z = np.copy(x)
|
|
|
|
Note that, when we modify x, y changes, but not z:
|
|
|
|
>>> x[0] = 10
|
|
>>> x[0] == y[0]
|
|
True
|
|
>>> x[0] == z[0]
|
|
False
|
|
|
|
"""
|
|
return array(a, order=order, copy=True)
|
|
|
|
# Basic operations
|
|
|
|
|
|
def _gradient_dispatcher(f, *varargs, **kwargs):
|
|
yield f
|
|
for v in varargs:
|
|
yield v
|
|
|
|
|
|
@array_function_dispatch(_gradient_dispatcher)
|
|
def gradient(f, *varargs, **kwargs):
|
|
"""
|
|
Return the gradient of an N-dimensional array.
|
|
|
|
The gradient is computed using second order accurate central differences
|
|
in the interior points and either first or second order accurate one-sides
|
|
(forward or backwards) differences at the boundaries.
|
|
The returned gradient hence has the same shape as the input array.
|
|
|
|
Parameters
|
|
----------
|
|
f : array_like
|
|
An N-dimensional array containing samples of a scalar function.
|
|
varargs : list of scalar or array, optional
|
|
Spacing between f values. Default unitary spacing for all dimensions.
|
|
Spacing can be specified using:
|
|
|
|
1. single scalar to specify a sample distance for all dimensions.
|
|
2. N scalars to specify a constant sample distance for each dimension.
|
|
i.e. `dx`, `dy`, `dz`, ...
|
|
3. N arrays to specify the coordinates of the values along each
|
|
dimension of F. The length of the array must match the size of
|
|
the corresponding dimension
|
|
4. Any combination of N scalars/arrays with the meaning of 2. and 3.
|
|
|
|
If `axis` is given, the number of varargs must equal the number of axes.
|
|
Default: 1.
|
|
|
|
edge_order : {1, 2}, optional
|
|
Gradient is calculated using N-th order accurate differences
|
|
at the boundaries. Default: 1.
|
|
|
|
.. versionadded:: 1.9.1
|
|
|
|
axis : None or int or tuple of ints, optional
|
|
Gradient is calculated only along the given axis or axes
|
|
The default (axis = None) is to calculate the gradient for all the axes
|
|
of the input array. axis may be negative, in which case it counts from
|
|
the last to the first axis.
|
|
|
|
.. versionadded:: 1.11.0
|
|
|
|
Returns
|
|
-------
|
|
gradient : ndarray or list of ndarray
|
|
A set of ndarrays (or a single ndarray if there is only one dimension)
|
|
corresponding to the derivatives of f with respect to each dimension.
|
|
Each derivative has the same shape as f.
|
|
|
|
Examples
|
|
--------
|
|
>>> f = np.array([1, 2, 4, 7, 11, 16], dtype=float)
|
|
>>> np.gradient(f)
|
|
array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
|
|
>>> np.gradient(f, 2)
|
|
array([0.5 , 0.75, 1.25, 1.75, 2.25, 2.5 ])
|
|
|
|
Spacing can be also specified with an array that represents the coordinates
|
|
of the values F along the dimensions.
|
|
For instance a uniform spacing:
|
|
|
|
>>> x = np.arange(f.size)
|
|
>>> np.gradient(f, x)
|
|
array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
|
|
|
|
Or a non uniform one:
|
|
|
|
>>> x = np.array([0., 1., 1.5, 3.5, 4., 6.], dtype=float)
|
|
>>> np.gradient(f, x)
|
|
array([1. , 3. , 3.5, 6.7, 6.9, 2.5])
|
|
|
|
For two dimensional arrays, the return will be two arrays ordered by
|
|
axis. In this example the first array stands for the gradient in
|
|
rows and the second one in columns direction:
|
|
|
|
>>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float))
|
|
[array([[ 2., 2., -1.],
|
|
[ 2., 2., -1.]]), array([[1. , 2.5, 4. ],
|
|
[1. , 1. , 1. ]])]
|
|
|
|
In this example the spacing is also specified:
|
|
uniform for axis=0 and non uniform for axis=1
|
|
|
|
>>> dx = 2.
|
|
>>> y = [1., 1.5, 3.5]
|
|
>>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), dx, y)
|
|
[array([[ 1. , 1. , -0.5],
|
|
[ 1. , 1. , -0.5]]), array([[2. , 2. , 2. ],
|
|
[2. , 1.7, 0.5]])]
|
|
|
|
It is possible to specify how boundaries are treated using `edge_order`
|
|
|
|
>>> x = np.array([0, 1, 2, 3, 4])
|
|
>>> f = x**2
|
|
>>> np.gradient(f, edge_order=1)
|
|
array([1., 2., 4., 6., 7.])
|
|
>>> np.gradient(f, edge_order=2)
|
|
array([0., 2., 4., 6., 8.])
|
|
|
|
The `axis` keyword can be used to specify a subset of axes of which the
|
|
gradient is calculated
|
|
|
|
>>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), axis=0)
|
|
array([[ 2., 2., -1.],
|
|
[ 2., 2., -1.]])
|
|
|
|
Notes
|
|
-----
|
|
Assuming that :math:`f\\in C^{3}` (i.e., :math:`f` has at least 3 continuous
|
|
derivatives) and let :math:`h_{*}` be a non-homogeneous stepsize, we
|
|
minimize the "consistency error" :math:`\\eta_{i}` between the true gradient
|
|
and its estimate from a linear combination of the neighboring grid-points:
|
|
|
|
.. math::
|
|
|
|
\\eta_{i} = f_{i}^{\\left(1\\right)} -
|
|
\\left[ \\alpha f\\left(x_{i}\\right) +
|
|
\\beta f\\left(x_{i} + h_{d}\\right) +
|
|
\\gamma f\\left(x_{i}-h_{s}\\right)
|
|
\\right]
|
|
|
|
By substituting :math:`f(x_{i} + h_{d})` and :math:`f(x_{i} - h_{s})`
|
|
with their Taylor series expansion, this translates into solving
|
|
the following the linear system:
|
|
|
|
.. math::
|
|
|
|
\\left\\{
|
|
\\begin{array}{r}
|
|
\\alpha+\\beta+\\gamma=0 \\\\
|
|
\\beta h_{d}-\\gamma h_{s}=1 \\\\
|
|
\\beta h_{d}^{2}+\\gamma h_{s}^{2}=0
|
|
\\end{array}
|
|
\\right.
|
|
|
|
The resulting approximation of :math:`f_{i}^{(1)}` is the following:
|
|
|
|
.. math::
|
|
|
|
\\hat f_{i}^{(1)} =
|
|
\\frac{
|
|
h_{s}^{2}f\\left(x_{i} + h_{d}\\right)
|
|
+ \\left(h_{d}^{2} - h_{s}^{2}\\right)f\\left(x_{i}\\right)
|
|
- h_{d}^{2}f\\left(x_{i}-h_{s}\\right)}
|
|
{ h_{s}h_{d}\\left(h_{d} + h_{s}\\right)}
|
|
+ \\mathcal{O}\\left(\\frac{h_{d}h_{s}^{2}
|
|
+ h_{s}h_{d}^{2}}{h_{d}
|
|
+ h_{s}}\\right)
|
|
|
|
It is worth noting that if :math:`h_{s}=h_{d}`
|
|
(i.e., data are evenly spaced)
|
|
we find the standard second order approximation:
|
|
|
|
.. math::
|
|
|
|
\\hat f_{i}^{(1)}=
|
|
\\frac{f\\left(x_{i+1}\\right) - f\\left(x_{i-1}\\right)}{2h}
|
|
+ \\mathcal{O}\\left(h^{2}\\right)
|
|
|
|
With a similar procedure the forward/backward approximations used for
|
|
boundaries can be derived.
|
|
|
|
References
|
|
----------
|
|
.. [1] Quarteroni A., Sacco R., Saleri F. (2007) Numerical Mathematics
|
|
(Texts in Applied Mathematics). New York: Springer.
|
|
.. [2] Durran D. R. (1999) Numerical Methods for Wave Equations
|
|
in Geophysical Fluid Dynamics. New York: Springer.
|
|
.. [3] Fornberg B. (1988) Generation of Finite Difference Formulas on
|
|
Arbitrarily Spaced Grids,
|
|
Mathematics of Computation 51, no. 184 : 699-706.
|
|
`PDF <http://www.ams.org/journals/mcom/1988-51-184/
|
|
S0025-5718-1988-0935077-0/S0025-5718-1988-0935077-0.pdf>`_.
|
|
"""
|
|
f = np.asanyarray(f)
|
|
N = f.ndim # number of dimensions
|
|
|
|
axes = kwargs.pop('axis', None)
|
|
if axes is None:
|
|
axes = tuple(range(N))
|
|
else:
|
|
axes = _nx.normalize_axis_tuple(axes, N)
|
|
|
|
len_axes = len(axes)
|
|
n = len(varargs)
|
|
if n == 0:
|
|
# no spacing argument - use 1 in all axes
|
|
dx = [1.0] * len_axes
|
|
elif n == 1 and np.ndim(varargs[0]) == 0:
|
|
# single scalar for all axes
|
|
dx = varargs * len_axes
|
|
elif n == len_axes:
|
|
# scalar or 1d array for each axis
|
|
dx = list(varargs)
|
|
for i, distances in enumerate(dx):
|
|
if np.ndim(distances) == 0:
|
|
continue
|
|
elif np.ndim(distances) != 1:
|
|
raise ValueError("distances must be either scalars or 1d")
|
|
if len(distances) != f.shape[axes[i]]:
|
|
raise ValueError("when 1d, distances must match "
|
|
"the length of the corresponding dimension")
|
|
diffx = np.diff(distances)
|
|
# if distances are constant reduce to the scalar case
|
|
# since it brings a consistent speedup
|
|
if (diffx == diffx[0]).all():
|
|
diffx = diffx[0]
|
|
dx[i] = diffx
|
|
else:
|
|
raise TypeError("invalid number of arguments")
|
|
|
|
edge_order = kwargs.pop('edge_order', 1)
|
|
if kwargs:
|
|
raise TypeError('"{}" are not valid keyword arguments.'.format(
|
|
'", "'.join(kwargs.keys())))
|
|
if edge_order > 2:
|
|
raise ValueError("'edge_order' greater than 2 not supported")
|
|
|
|
# use central differences on interior and one-sided differences on the
|
|
# endpoints. This preserves second order-accuracy over the full domain.
|
|
|
|
outvals = []
|
|
|
|
# create slice objects --- initially all are [:, :, ..., :]
|
|
slice1 = [slice(None)]*N
|
|
slice2 = [slice(None)]*N
|
|
slice3 = [slice(None)]*N
|
|
slice4 = [slice(None)]*N
|
|
|
|
otype = f.dtype
|
|
if otype.type is np.datetime64:
|
|
# the timedelta dtype with the same unit information
|
|
otype = np.dtype(otype.name.replace('datetime', 'timedelta'))
|
|
# view as timedelta to allow addition
|
|
f = f.view(otype)
|
|
elif otype.type is np.timedelta64:
|
|
pass
|
|
elif np.issubdtype(otype, np.inexact):
|
|
pass
|
|
else:
|
|
# all other types convert to floating point
|
|
otype = np.double
|
|
|
|
for axis, ax_dx in zip(axes, dx):
|
|
if f.shape[axis] < edge_order + 1:
|
|
raise ValueError(
|
|
"Shape of array too small to calculate a numerical gradient, "
|
|
"at least (edge_order + 1) elements are required.")
|
|
# result allocation
|
|
out = np.empty_like(f, dtype=otype)
|
|
|
|
# spacing for the current axis
|
|
uniform_spacing = np.ndim(ax_dx) == 0
|
|
|
|
# Numerical differentiation: 2nd order interior
|
|
slice1[axis] = slice(1, -1)
|
|
slice2[axis] = slice(None, -2)
|
|
slice3[axis] = slice(1, -1)
|
|
slice4[axis] = slice(2, None)
|
|
|
|
if uniform_spacing:
|
|
out[tuple(slice1)] = (f[tuple(slice4)] - f[tuple(slice2)]) / (2. * ax_dx)
|
|
else:
|
|
dx1 = ax_dx[0:-1]
|
|
dx2 = ax_dx[1:]
|
|
a = -(dx2)/(dx1 * (dx1 + dx2))
|
|
b = (dx2 - dx1) / (dx1 * dx2)
|
|
c = dx1 / (dx2 * (dx1 + dx2))
|
|
# fix the shape for broadcasting
|
|
shape = np.ones(N, dtype=int)
|
|
shape[axis] = -1
|
|
a.shape = b.shape = c.shape = shape
|
|
# 1D equivalent -- out[1:-1] = a * f[:-2] + b * f[1:-1] + c * f[2:]
|
|
out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
|
|
|
|
# Numerical differentiation: 1st order edges
|
|
if edge_order == 1:
|
|
slice1[axis] = 0
|
|
slice2[axis] = 1
|
|
slice3[axis] = 0
|
|
dx_0 = ax_dx if uniform_spacing else ax_dx[0]
|
|
# 1D equivalent -- out[0] = (f[1] - f[0]) / (x[1] - x[0])
|
|
out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0
|
|
|
|
slice1[axis] = -1
|
|
slice2[axis] = -1
|
|
slice3[axis] = -2
|
|
dx_n = ax_dx if uniform_spacing else ax_dx[-1]
|
|
# 1D equivalent -- out[-1] = (f[-1] - f[-2]) / (x[-1] - x[-2])
|
|
out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_n
|
|
|
|
# Numerical differentiation: 2nd order edges
|
|
else:
|
|
slice1[axis] = 0
|
|
slice2[axis] = 0
|
|
slice3[axis] = 1
|
|
slice4[axis] = 2
|
|
if uniform_spacing:
|
|
a = -1.5 / ax_dx
|
|
b = 2. / ax_dx
|
|
c = -0.5 / ax_dx
|
|
else:
|
|
dx1 = ax_dx[0]
|
|
dx2 = ax_dx[1]
|
|
a = -(2. * dx1 + dx2)/(dx1 * (dx1 + dx2))
|
|
b = (dx1 + dx2) / (dx1 * dx2)
|
|
c = - dx1 / (dx2 * (dx1 + dx2))
|
|
# 1D equivalent -- out[0] = a * f[0] + b * f[1] + c * f[2]
|
|
out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
|
|
|
|
slice1[axis] = -1
|
|
slice2[axis] = -3
|
|
slice3[axis] = -2
|
|
slice4[axis] = -1
|
|
if uniform_spacing:
|
|
a = 0.5 / ax_dx
|
|
b = -2. / ax_dx
|
|
c = 1.5 / ax_dx
|
|
else:
|
|
dx1 = ax_dx[-2]
|
|
dx2 = ax_dx[-1]
|
|
a = (dx2) / (dx1 * (dx1 + dx2))
|
|
b = - (dx2 + dx1) / (dx1 * dx2)
|
|
c = (2. * dx2 + dx1) / (dx2 * (dx1 + dx2))
|
|
# 1D equivalent -- out[-1] = a * f[-3] + b * f[-2] + c * f[-1]
|
|
out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
|
|
|
|
outvals.append(out)
|
|
|
|
# reset the slice object in this dimension to ":"
|
|
slice1[axis] = slice(None)
|
|
slice2[axis] = slice(None)
|
|
slice3[axis] = slice(None)
|
|
slice4[axis] = slice(None)
|
|
|
|
if len_axes == 1:
|
|
return outvals[0]
|
|
else:
|
|
return outvals
|
|
|
|
|
|
def _diff_dispatcher(a, n=None, axis=None, prepend=None, append=None):
|
|
return (a, prepend, append)
|
|
|
|
|
|
@array_function_dispatch(_diff_dispatcher)
|
|
def diff(a, n=1, axis=-1, prepend=np._NoValue, append=np._NoValue):
|
|
"""
|
|
Calculate the n-th discrete difference along the given axis.
|
|
|
|
The first difference is given by ``out[i] = a[i+1] - a[i]`` along
|
|
the given axis, higher differences are calculated by using `diff`
|
|
recursively.
|
|
|
|
Parameters
|
|
----------
|
|
a : array_like
|
|
Input array
|
|
n : int, optional
|
|
The number of times values are differenced. If zero, the input
|
|
is returned as-is.
|
|
axis : int, optional
|
|
The axis along which the difference is taken, default is the
|
|
last axis.
|
|
prepend, append : array_like, optional
|
|
Values to prepend or append to "a" along axis prior to
|
|
performing the difference. Scalar values are expanded to
|
|
arrays with length 1 in the direction of axis and the shape
|
|
of the input array in along all other axes. Otherwise the
|
|
dimension and shape must match "a" except along axis.
|
|
|
|
Returns
|
|
-------
|
|
diff : ndarray
|
|
The n-th differences. The shape of the output is the same as `a`
|
|
except along `axis` where the dimension is smaller by `n`. The
|
|
type of the output is the same as the type of the difference
|
|
between any two elements of `a`. This is the same as the type of
|
|
`a` in most cases. A notable exception is `datetime64`, which
|
|
results in a `timedelta64` output array.
|
|
|
|
See Also
|
|
--------
|
|
gradient, ediff1d, cumsum
|
|
|
|
Notes
|
|
-----
|
|
Type is preserved for boolean arrays, so the result will contain
|
|
`False` when consecutive elements are the same and `True` when they
|
|
differ.
|
|
|
|
For unsigned integer arrays, the results will also be unsigned. This
|
|
should not be surprising, as the result is consistent with
|
|
calculating the difference directly:
|
|
|
|
>>> u8_arr = np.array([1, 0], dtype=np.uint8)
|
|
>>> np.diff(u8_arr)
|
|
array([255], dtype=uint8)
|
|
>>> u8_arr[1,...] - u8_arr[0,...]
|
|
255
|
|
|
|
If this is not desirable, then the array should be cast to a larger
|
|
integer type first:
|
|
|
|
>>> i16_arr = u8_arr.astype(np.int16)
|
|
>>> np.diff(i16_arr)
|
|
array([-1], dtype=int16)
|
|
|
|
Examples
|
|
--------
|
|
>>> x = np.array([1, 2, 4, 7, 0])
|
|
>>> np.diff(x)
|
|
array([ 1, 2, 3, -7])
|
|
>>> np.diff(x, n=2)
|
|
array([ 1, 1, -10])
|
|
|
|
>>> x = np.array([[1, 3, 6, 10], [0, 5, 6, 8]])
|
|
>>> np.diff(x)
|
|
array([[2, 3, 4],
|
|
[5, 1, 2]])
|
|
>>> np.diff(x, axis=0)
|
|
array([[-1, 2, 0, -2]])
|
|
|
|
>>> x = np.arange('1066-10-13', '1066-10-16', dtype=np.datetime64)
|
|
>>> np.diff(x)
|
|
array([1, 1], dtype='timedelta64[D]')
|
|
|
|
"""
|
|
if n == 0:
|
|
return a
|
|
if n < 0:
|
|
raise ValueError(
|
|
"order must be non-negative but got " + repr(n))
|
|
|
|
a = asanyarray(a)
|
|
nd = a.ndim
|
|
if nd == 0:
|
|
raise ValueError("diff requires input that is at least one dimensional")
|
|
axis = normalize_axis_index(axis, nd)
|
|
|
|
combined = []
|
|
if prepend is not np._NoValue:
|
|
prepend = np.asanyarray(prepend)
|
|
if prepend.ndim == 0:
|
|
shape = list(a.shape)
|
|
shape[axis] = 1
|
|
prepend = np.broadcast_to(prepend, tuple(shape))
|
|
combined.append(prepend)
|
|
|
|
combined.append(a)
|
|
|
|
if append is not np._NoValue:
|
|
append = np.asanyarray(append)
|
|
if append.ndim == 0:
|
|
shape = list(a.shape)
|
|
shape[axis] = 1
|
|
append = np.broadcast_to(append, tuple(shape))
|
|
combined.append(append)
|
|
|
|
if len(combined) > 1:
|
|
a = np.concatenate(combined, axis)
|
|
|
|
slice1 = [slice(None)] * nd
|
|
slice2 = [slice(None)] * nd
|
|
slice1[axis] = slice(1, None)
|
|
slice2[axis] = slice(None, -1)
|
|
slice1 = tuple(slice1)
|
|
slice2 = tuple(slice2)
|
|
|
|
op = not_equal if a.dtype == np.bool_ else subtract
|
|
for _ in range(n):
|
|
a = op(a[slice1], a[slice2])
|
|
|
|
return a
|
|
|
|
|
|
def _interp_dispatcher(x, xp, fp, left=None, right=None, period=None):
|
|
return (x, xp, fp)
|
|
|
|
|
|
@array_function_dispatch(_interp_dispatcher)
|
|
def interp(x, xp, fp, left=None, right=None, period=None):
|
|
"""
|
|
One-dimensional linear interpolation.
|
|
|
|
Returns the one-dimensional piecewise linear interpolant to a function
|
|
with given discrete data points (`xp`, `fp`), evaluated at `x`.
|
|
|
|
Parameters
|
|
----------
|
|
x : array_like
|
|
The x-coordinates at which to evaluate the interpolated values.
|
|
|
|
xp : 1-D sequence of floats
|
|
The x-coordinates of the data points, must be increasing if argument
|
|
`period` is not specified. Otherwise, `xp` is internally sorted after
|
|
normalizing the periodic boundaries with ``xp = xp % period``.
|
|
|
|
fp : 1-D sequence of float or complex
|
|
The y-coordinates of the data points, same length as `xp`.
|
|
|
|
left : optional float or complex corresponding to fp
|
|
Value to return for `x < xp[0]`, default is `fp[0]`.
|
|
|
|
right : optional float or complex corresponding to fp
|
|
Value to return for `x > xp[-1]`, default is `fp[-1]`.
|
|
|
|
period : None or float, optional
|
|
A period for the x-coordinates. This parameter allows the proper
|
|
interpolation of angular x-coordinates. Parameters `left` and `right`
|
|
are ignored if `period` is specified.
|
|
|
|
.. versionadded:: 1.10.0
|
|
|
|
Returns
|
|
-------
|
|
y : float or complex (corresponding to fp) or ndarray
|
|
The interpolated values, same shape as `x`.
|
|
|
|
Raises
|
|
------
|
|
ValueError
|
|
If `xp` and `fp` have different length
|
|
If `xp` or `fp` are not 1-D sequences
|
|
If `period == 0`
|
|
|
|
Notes
|
|
-----
|
|
Does not check that the x-coordinate sequence `xp` is increasing.
|
|
If `xp` is not increasing, the results are nonsense.
|
|
A simple check for increasing is::
|
|
|
|
np.all(np.diff(xp) > 0)
|
|
|
|
Examples
|
|
--------
|
|
>>> xp = [1, 2, 3]
|
|
>>> fp = [3, 2, 0]
|
|
>>> np.interp(2.5, xp, fp)
|
|
1.0
|
|
>>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
|
|
array([3. , 3. , 2.5 , 0.56, 0. ])
|
|
>>> UNDEF = -99.0
|
|
>>> np.interp(3.14, xp, fp, right=UNDEF)
|
|
-99.0
|
|
|
|
Plot an interpolant to the sine function:
|
|
|
|
>>> x = np.linspace(0, 2*np.pi, 10)
|
|
>>> y = np.sin(x)
|
|
>>> xvals = np.linspace(0, 2*np.pi, 50)
|
|
>>> yinterp = np.interp(xvals, x, y)
|
|
>>> import matplotlib.pyplot as plt
|
|
>>> plt.plot(x, y, 'o')
|
|
[<matplotlib.lines.Line2D object at 0x...>]
|
|
>>> plt.plot(xvals, yinterp, '-x')
|
|
[<matplotlib.lines.Line2D object at 0x...>]
|
|
>>> plt.show()
|
|
|
|
Interpolation with periodic x-coordinates:
|
|
|
|
>>> x = [-180, -170, -185, 185, -10, -5, 0, 365]
|
|
>>> xp = [190, -190, 350, -350]
|
|
>>> fp = [5, 10, 3, 4]
|
|
>>> np.interp(x, xp, fp, period=360)
|
|
array([7.5 , 5. , 8.75, 6.25, 3. , 3.25, 3.5 , 3.75])
|
|
|
|
Complex interpolation:
|
|
|
|
>>> x = [1.5, 4.0]
|
|
>>> xp = [2,3,5]
|
|
>>> fp = [1.0j, 0, 2+3j]
|
|
>>> np.interp(x, xp, fp)
|
|
array([0.+1.j , 1.+1.5j])
|
|
|
|
"""
|
|
|
|
fp = np.asarray(fp)
|
|
|
|
if np.iscomplexobj(fp):
|
|
interp_func = compiled_interp_complex
|
|
input_dtype = np.complex128
|
|
else:
|
|
interp_func = compiled_interp
|
|
input_dtype = np.float64
|
|
|
|
if period is not None:
|
|
if period == 0:
|
|
raise ValueError("period must be a non-zero value")
|
|
period = abs(period)
|
|
left = None
|
|
right = None
|
|
|
|
x = np.asarray(x, dtype=np.float64)
|
|
xp = np.asarray(xp, dtype=np.float64)
|
|
fp = np.asarray(fp, dtype=input_dtype)
|
|
|
|
if xp.ndim != 1 or fp.ndim != 1:
|
|
raise ValueError("Data points must be 1-D sequences")
|
|
if xp.shape[0] != fp.shape[0]:
|
|
raise ValueError("fp and xp are not of the same length")
|
|
# normalizing periodic boundaries
|
|
x = x % period
|
|
xp = xp % period
|
|
asort_xp = np.argsort(xp)
|
|
xp = xp[asort_xp]
|
|
fp = fp[asort_xp]
|
|
xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period))
|
|
fp = np.concatenate((fp[-1:], fp, fp[0:1]))
|
|
|
|
return interp_func(x, xp, fp, left, right)
|
|
|
|
|
|
def _angle_dispatcher(z, deg=None):
|
|
return (z,)
|
|
|
|
|
|
@array_function_dispatch(_angle_dispatcher)
|
|
def angle(z, deg=False):
|
|
"""
|
|
Return the angle of the complex argument.
|
|
|
|
Parameters
|
|
----------
|
|
z : array_like
|
|
A complex number or sequence of complex numbers.
|
|
deg : bool, optional
|
|
Return angle in degrees if True, radians if False (default).
|
|
|
|
Returns
|
|
-------
|
|
angle : ndarray or scalar
|
|
The counterclockwise angle from the positive real axis on the complex
|
|
plane in the range ``(-pi, pi]``, with dtype as numpy.float64.
|
|
|
|
..versionchanged:: 1.16.0
|
|
This function works on subclasses of ndarray like `ma.array`.
|
|
|
|
See Also
|
|
--------
|
|
arctan2
|
|
absolute
|
|
|
|
Examples
|
|
--------
|
|
>>> np.angle([1.0, 1.0j, 1+1j]) # in radians
|
|
array([ 0. , 1.57079633, 0.78539816]) # may vary
|
|
>>> np.angle(1+1j, deg=True) # in degrees
|
|
45.0
|
|
|
|
"""
|
|
z = asanyarray(z)
|
|
if issubclass(z.dtype.type, _nx.complexfloating):
|
|
zimag = z.imag
|
|
zreal = z.real
|
|
else:
|
|
zimag = 0
|
|
zreal = z
|
|
|
|
a = arctan2(zimag, zreal)
|
|
if deg:
|
|
a *= 180/pi
|
|
return a
|
|
|
|
|
|
def _unwrap_dispatcher(p, discont=None, axis=None):
|
|
return (p,)
|
|
|
|
|
|
@array_function_dispatch(_unwrap_dispatcher)
|
|
def unwrap(p, discont=pi, axis=-1):
|
|
"""
|
|
Unwrap by changing deltas between values to 2*pi complement.
|
|
|
|
Unwrap radian phase `p` by changing absolute jumps greater than
|
|
`discont` to their 2*pi complement along the given axis.
|
|
|
|
Parameters
|
|
----------
|
|
p : array_like
|
|
Input array.
|
|
discont : float, optional
|
|
Maximum discontinuity between values, default is ``pi``.
|
|
axis : int, optional
|
|
Axis along which unwrap will operate, default is the last axis.
|
|
|
|
Returns
|
|
-------
|
|
out : ndarray
|
|
Output array.
|
|
|
|
See Also
|
|
--------
|
|
rad2deg, deg2rad
|
|
|
|
Notes
|
|
-----
|
|
If the discontinuity in `p` is smaller than ``pi``, but larger than
|
|
`discont`, no unwrapping is done because taking the 2*pi complement
|
|
would only make the discontinuity larger.
|
|
|
|
Examples
|
|
--------
|
|
>>> phase = np.linspace(0, np.pi, num=5)
|
|
>>> phase[3:] += np.pi
|
|
>>> phase
|
|
array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531]) # may vary
|
|
>>> np.unwrap(phase)
|
|
array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ]) # may vary
|
|
|
|
"""
|
|
p = asarray(p)
|
|
nd = p.ndim
|
|
dd = diff(p, axis=axis)
|
|
slice1 = [slice(None, None)]*nd # full slices
|
|
slice1[axis] = slice(1, None)
|
|
slice1 = tuple(slice1)
|
|
ddmod = mod(dd + pi, 2*pi) - pi
|
|
_nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
|
|
ph_correct = ddmod - dd
|
|
_nx.copyto(ph_correct, 0, where=abs(dd) < discont)
|
|
up = array(p, copy=True, dtype='d')
|
|
up[slice1] = p[slice1] + ph_correct.cumsum(axis)
|
|
return up
|
|
|
|
|
|
def _sort_complex(a):
|
|
return (a,)
|
|
|
|
|
|
@array_function_dispatch(_sort_complex)
|
|
def sort_complex(a):
|
|
"""
|
|
Sort a complex array using the real part first, then the imaginary part.
|
|
|
|
Parameters
|
|
----------
|
|
a : array_like
|
|
Input array
|
|
|
|
Returns
|
|
-------
|
|
out : complex ndarray
|
|
Always returns a sorted complex array.
|
|
|
|
Examples
|
|
--------
|
|
>>> np.sort_complex([5, 3, 6, 2, 1])
|
|
array([1.+0.j, 2.+0.j, 3.+0.j, 5.+0.j, 6.+0.j])
|
|
|
|
>>> np.sort_complex([1 + 2j, 2 - 1j, 3 - 2j, 3 - 3j, 3 + 5j])
|
|
array([1.+2.j, 2.-1.j, 3.-3.j, 3.-2.j, 3.+5.j])
|
|
|
|
"""
|
|
b = array(a, copy=True)
|
|
b.sort()
|
|
if not issubclass(b.dtype.type, _nx.complexfloating):
|
|
if b.dtype.char in 'bhBH':
|
|
return b.astype('F')
|
|
elif b.dtype.char == 'g':
|
|
return b.astype('G')
|
|
else:
|
|
return b.astype('D')
|
|
else:
|
|
return b
|
|
|
|
|
|
def _trim_zeros(filt, trim=None):
|
|
return (filt,)
|
|
|
|
|
|
@array_function_dispatch(_trim_zeros)
|
|
def trim_zeros(filt, trim='fb'):
|
|
"""
|
|
Trim the leading and/or trailing zeros from a 1-D array or sequence.
|
|
|
|
Parameters
|
|
----------
|
|
filt : 1-D array or sequence
|
|
Input array.
|
|
trim : str, optional
|
|
A string with 'f' representing trim from front and 'b' to trim from
|
|
back. Default is 'fb', trim zeros from both front and back of the
|
|
array.
|
|
|
|
Returns
|
|
-------
|
|
trimmed : 1-D array or sequence
|
|
The result of trimming the input. The input data type is preserved.
|
|
|
|
Examples
|
|
--------
|
|
>>> a = np.array((0, 0, 0, 1, 2, 3, 0, 2, 1, 0))
|
|
>>> np.trim_zeros(a)
|
|
array([1, 2, 3, 0, 2, 1])
|
|
|
|
>>> np.trim_zeros(a, 'b')
|
|
array([0, 0, 0, ..., 0, 2, 1])
|
|
|
|
The input data type is preserved, list/tuple in means list/tuple out.
|
|
|
|
>>> np.trim_zeros([0, 1, 2, 0])
|
|
[1, 2]
|
|
|
|
"""
|
|
first = 0
|
|
trim = trim.upper()
|
|
if 'F' in trim:
|
|
for i in filt:
|
|
if i != 0.:
|
|
break
|
|
else:
|
|
first = first + 1
|
|
last = len(filt)
|
|
if 'B' in trim:
|
|
for i in filt[::-1]:
|
|
if i != 0.:
|
|
break
|
|
else:
|
|
last = last - 1
|
|
return filt[first:last]
|
|
|
|
def _extract_dispatcher(condition, arr):
|
|
return (condition, arr)
|
|
|
|
|
|
@array_function_dispatch(_extract_dispatcher)
|
|
def extract(condition, arr):
|
|
"""
|
|
Return the elements of an array that satisfy some condition.
|
|
|
|
This is equivalent to ``np.compress(ravel(condition), ravel(arr))``. If
|
|
`condition` is boolean ``np.extract`` is equivalent to ``arr[condition]``.
|
|
|
|
Note that `place` does the exact opposite of `extract`.
|
|
|
|
Parameters
|
|
----------
|
|
condition : array_like
|
|
An array whose nonzero or True entries indicate the elements of `arr`
|
|
to extract.
|
|
arr : array_like
|
|
Input array of the same size as `condition`.
|
|
|
|
Returns
|
|
-------
|
|
extract : ndarray
|
|
Rank 1 array of values from `arr` where `condition` is True.
|
|
|
|
See Also
|
|
--------
|
|
take, put, copyto, compress, place
|
|
|
|
Examples
|
|
--------
|
|
>>> arr = np.arange(12).reshape((3, 4))
|
|
>>> arr
|
|
array([[ 0, 1, 2, 3],
|
|
[ 4, 5, 6, 7],
|
|
[ 8, 9, 10, 11]])
|
|
>>> condition = np.mod(arr, 3)==0
|
|
>>> condition
|
|
array([[ True, False, False, True],
|
|
[False, False, True, False],
|
|
[False, True, False, False]])
|
|
>>> np.extract(condition, arr)
|
|
array([0, 3, 6, 9])
|
|
|
|
|
|
If `condition` is boolean:
|
|
|
|
>>> arr[condition]
|
|
array([0, 3, 6, 9])
|
|
|
|
"""
|
|
return _nx.take(ravel(arr), nonzero(ravel(condition))[0])
|
|
|
|
|
|
def _place_dispatcher(arr, mask, vals):
|
|
return (arr, mask, vals)
|
|
|
|
|
|
@array_function_dispatch(_place_dispatcher)
|
|
def place(arr, mask, vals):
|
|
"""
|
|
Change elements of an array based on conditional and input values.
|
|
|
|
Similar to ``np.copyto(arr, vals, where=mask)``, the difference is that
|
|
`place` uses the first N elements of `vals`, where N is the number of
|
|
True values in `mask`, while `copyto` uses the elements where `mask`
|
|
is True.
|
|
|
|
Note that `extract` does the exact opposite of `place`.
|
|
|
|
Parameters
|
|
----------
|
|
arr : ndarray
|
|
Array to put data into.
|
|
mask : array_like
|
|
Boolean mask array. Must have the same size as `a`.
|
|
vals : 1-D sequence
|
|
Values to put into `a`. Only the first N elements are used, where
|
|
N is the number of True values in `mask`. If `vals` is smaller
|
|
than N, it will be repeated, and if elements of `a` are to be masked,
|
|
this sequence must be non-empty.
|
|
|
|
See Also
|
|
--------
|
|
copyto, put, take, extract
|
|
|
|
Examples
|
|
--------
|
|
>>> arr = np.arange(6).reshape(2, 3)
|
|
>>> np.place(arr, arr>2, [44, 55])
|
|
>>> arr
|
|
array([[ 0, 1, 2],
|
|
[44, 55, 44]])
|
|
|
|
"""
|
|
if not isinstance(arr, np.ndarray):
|
|
raise TypeError("argument 1 must be numpy.ndarray, "
|
|
"not {name}".format(name=type(arr).__name__))
|
|
|
|
return _insert(arr, mask, vals)
|
|
|
|
|
|
def disp(mesg, device=None, linefeed=True):
|
|
"""
|
|
Display a message on a device.
|
|
|
|
Parameters
|
|
----------
|
|
mesg : str
|
|
Message to display.
|
|
device : object
|
|
Device to write message. If None, defaults to ``sys.stdout`` which is
|
|
very similar to ``print``. `device` needs to have ``write()`` and
|
|
``flush()`` methods.
|
|
linefeed : bool, optional
|
|
Option whether to print a line feed or not. Defaults to True.
|
|
|
|
Raises
|
|
------
|
|
AttributeError
|
|
If `device` does not have a ``write()`` or ``flush()`` method.
|
|
|
|
Examples
|
|
--------
|
|
Besides ``sys.stdout``, a file-like object can also be used as it has
|
|
both required methods:
|
|
|
|
>>> from io import StringIO
|
|
>>> buf = StringIO()
|
|
>>> np.disp(u'"Display" in a file', device=buf)
|
|
>>> buf.getvalue()
|
|
'"Display" in a file\\n'
|
|
|
|
"""
|
|
if device is None:
|
|
device = sys.stdout
|
|
if linefeed:
|
|
device.write('%s\n' % mesg)
|
|
else:
|
|
device.write('%s' % mesg)
|
|
device.flush()
|
|
return
|
|
|
|
|
|
# See https://docs.scipy.org/doc/numpy/reference/c-api.generalized-ufuncs.html
|
|
_DIMENSION_NAME = r'\w+'
|
|
_CORE_DIMENSION_LIST = '(?:{0:}(?:,{0:})*)?'.format(_DIMENSION_NAME)
|
|
_ARGUMENT = r'\({}\)'.format(_CORE_DIMENSION_LIST)
|
|
_ARGUMENT_LIST = '{0:}(?:,{0:})*'.format(_ARGUMENT)
|
|
_SIGNATURE = '^{0:}->{0:}$'.format(_ARGUMENT_LIST)
|
|
|
|
|
|
def _parse_gufunc_signature(signature):
|
|
"""
|
|
Parse string signatures for a generalized universal function.
|
|
|
|
Arguments
|
|
---------
|
|
signature : string
|
|
Generalized universal function signature, e.g., ``(m,n),(n,p)->(m,p)``
|
|
for ``np.matmul``.
|
|
|
|
Returns
|
|
-------
|
|
Tuple of input and output core dimensions parsed from the signature, each
|
|
of the form List[Tuple[str, ...]].
|
|
"""
|
|
if not re.match(_SIGNATURE, signature):
|
|
raise ValueError(
|
|
'not a valid gufunc signature: {}'.format(signature))
|
|
return tuple([tuple(re.findall(_DIMENSION_NAME, arg))
|
|
for arg in re.findall(_ARGUMENT, arg_list)]
|
|
for arg_list in signature.split('->'))
|
|
|
|
|
|
def _update_dim_sizes(dim_sizes, arg, core_dims):
|
|
"""
|
|
Incrementally check and update core dimension sizes for a single argument.
|
|
|
|
Arguments
|
|
---------
|
|
dim_sizes : Dict[str, int]
|
|
Sizes of existing core dimensions. Will be updated in-place.
|
|
arg : ndarray
|
|
Argument to examine.
|
|
core_dims : Tuple[str, ...]
|
|
Core dimensions for this argument.
|
|
"""
|
|
if not core_dims:
|
|
return
|
|
|
|
num_core_dims = len(core_dims)
|
|
if arg.ndim < num_core_dims:
|
|
raise ValueError(
|
|
'%d-dimensional argument does not have enough '
|
|
'dimensions for all core dimensions %r'
|
|
% (arg.ndim, core_dims))
|
|
|
|
core_shape = arg.shape[-num_core_dims:]
|
|
for dim, size in zip(core_dims, core_shape):
|
|
if dim in dim_sizes:
|
|
if size != dim_sizes[dim]:
|
|
raise ValueError(
|
|
'inconsistent size for core dimension %r: %r vs %r'
|
|
% (dim, size, dim_sizes[dim]))
|
|
else:
|
|
dim_sizes[dim] = size
|
|
|
|
|
|
def _parse_input_dimensions(args, input_core_dims):
|
|
"""
|
|
Parse broadcast and core dimensions for vectorize with a signature.
|
|
|
|
Arguments
|
|
---------
|
|
args : Tuple[ndarray, ...]
|
|
Tuple of input arguments to examine.
|
|
input_core_dims : List[Tuple[str, ...]]
|
|
List of core dimensions corresponding to each input.
|
|
|
|
Returns
|
|
-------
|
|
broadcast_shape : Tuple[int, ...]
|
|
Common shape to broadcast all non-core dimensions to.
|
|
dim_sizes : Dict[str, int]
|
|
Common sizes for named core dimensions.
|
|
"""
|
|
broadcast_args = []
|
|
dim_sizes = {}
|
|
for arg, core_dims in zip(args, input_core_dims):
|
|
_update_dim_sizes(dim_sizes, arg, core_dims)
|
|
ndim = arg.ndim - len(core_dims)
|
|
dummy_array = np.lib.stride_tricks.as_strided(0, arg.shape[:ndim])
|
|
broadcast_args.append(dummy_array)
|
|
broadcast_shape = np.lib.stride_tricks._broadcast_shape(*broadcast_args)
|
|
return broadcast_shape, dim_sizes
|
|
|
|
|
|
def _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims):
|
|
"""Helper for calculating broadcast shapes with core dimensions."""
|
|
return [broadcast_shape + tuple(dim_sizes[dim] for dim in core_dims)
|
|
for core_dims in list_of_core_dims]
|
|
|
|
|
|
def _create_arrays(broadcast_shape, dim_sizes, list_of_core_dims, dtypes):
|
|
"""Helper for creating output arrays in vectorize."""
|
|
shapes = _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims)
|
|
arrays = tuple(np.empty(shape, dtype=dtype)
|
|
for shape, dtype in zip(shapes, dtypes))
|
|
return arrays
|
|
|
|
|
|
@set_module('numpy')
|
|
class vectorize(object):
|
|
"""
|
|
vectorize(pyfunc, otypes=None, doc=None, excluded=None, cache=False,
|
|
signature=None)
|
|
|
|
Generalized function class.
|
|
|
|
Define a vectorized function which takes a nested sequence of objects or
|
|
numpy arrays as inputs and returns a single numpy array or a tuple of numpy
|
|
arrays. The vectorized function evaluates `pyfunc` over successive tuples
|
|
of the input arrays like the python map function, except it uses the
|
|
broadcasting rules of numpy.
|
|
|
|
The data type of the output of `vectorized` is determined by calling
|
|
the function with the first element of the input. This can be avoided
|
|
by specifying the `otypes` argument.
|
|
|
|
Parameters
|
|
----------
|
|
pyfunc : callable
|
|
A python function or method.
|
|
otypes : str or list of dtypes, optional
|
|
The output data type. It must be specified as either a string of
|
|
typecode characters or a list of data type specifiers. There should
|
|
be one data type specifier for each output.
|
|
doc : str, optional
|
|
The docstring for the function. If `None`, the docstring will be the
|
|
``pyfunc.__doc__``.
|
|
excluded : set, optional
|
|
Set of strings or integers representing the positional or keyword
|
|
arguments for which the function will not be vectorized. These will be
|
|
passed directly to `pyfunc` unmodified.
|
|
|
|
.. versionadded:: 1.7.0
|
|
|
|
cache : bool, optional
|
|
If `True`, then cache the first function call that determines the number
|
|
of outputs if `otypes` is not provided.
|
|
|
|
.. versionadded:: 1.7.0
|
|
|
|
signature : string, optional
|
|
Generalized universal function signature, e.g., ``(m,n),(n)->(m)`` for
|
|
vectorized matrix-vector multiplication. If provided, ``pyfunc`` will
|
|
be called with (and expected to return) arrays with shapes given by the
|
|
size of corresponding core dimensions. By default, ``pyfunc`` is
|
|
assumed to take scalars as input and output.
|
|
|
|
.. versionadded:: 1.12.0
|
|
|
|
Returns
|
|
-------
|
|
vectorized : callable
|
|
Vectorized function.
|
|
|
|
See Also
|
|
--------
|
|
frompyfunc : Takes an arbitrary Python function and returns a ufunc
|
|
|
|
Notes
|
|
-----
|
|
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.
|
|
|
|
The new keyword argument interface and `excluded` argument support
|
|
further degrades performance.
|
|
|
|
References
|
|
----------
|
|
.. [1] NumPy Reference, section `Generalized Universal Function API
|
|
<https://docs.scipy.org/doc/numpy/reference/c-api.generalized-ufuncs.html>`_.
|
|
|
|
Examples
|
|
--------
|
|
>>> def myfunc(a, b):
|
|
... "Return a-b if a>b, otherwise return a+b"
|
|
... if a > b:
|
|
... return a - b
|
|
... else:
|
|
... return a + b
|
|
|
|
>>> vfunc = np.vectorize(myfunc)
|
|
>>> vfunc([1, 2, 3, 4], 2)
|
|
array([3, 4, 1, 2])
|
|
|
|
The docstring is taken from the input function to `vectorize` unless it
|
|
is specified:
|
|
|
|
>>> vfunc.__doc__
|
|
'Return a-b if a>b, otherwise return a+b'
|
|
>>> vfunc = np.vectorize(myfunc, doc='Vectorized `myfunc`')
|
|
>>> vfunc.__doc__
|
|
'Vectorized `myfunc`'
|
|
|
|
The output type is determined by evaluating the first element of the input,
|
|
unless it is specified:
|
|
|
|
>>> out = vfunc([1, 2, 3, 4], 2)
|
|
>>> type(out[0])
|
|
<class 'numpy.int64'>
|
|
>>> vfunc = np.vectorize(myfunc, otypes=[float])
|
|
>>> out = vfunc([1, 2, 3, 4], 2)
|
|
>>> type(out[0])
|
|
<class 'numpy.float64'>
|
|
|
|
The `excluded` argument can be used to prevent vectorizing over certain
|
|
arguments. This can be useful for array-like arguments of a fixed length
|
|
such as the coefficients for a polynomial as in `polyval`:
|
|
|
|
>>> def mypolyval(p, x):
|
|
... _p = list(p)
|
|
... res = _p.pop(0)
|
|
... while _p:
|
|
... res = res*x + _p.pop(0)
|
|
... return res
|
|
>>> vpolyval = np.vectorize(mypolyval, excluded=['p'])
|
|
>>> vpolyval(p=[1, 2, 3], x=[0, 1])
|
|
array([3, 6])
|
|
|
|
Positional arguments may also be excluded by specifying their position:
|
|
|
|
>>> vpolyval.excluded.add(0)
|
|
>>> vpolyval([1, 2, 3], x=[0, 1])
|
|
array([3, 6])
|
|
|
|
The `signature` argument allows for vectorizing functions that act on
|
|
non-scalar arrays of fixed length. For example, you can use it for a
|
|
vectorized calculation of Pearson correlation coefficient and its p-value:
|
|
|
|
>>> import scipy.stats
|
|
>>> pearsonr = np.vectorize(scipy.stats.pearsonr,
|
|
... signature='(n),(n)->(),()')
|
|
>>> pearsonr([[0, 1, 2, 3]], [[1, 2, 3, 4], [4, 3, 2, 1]])
|
|
(array([ 1., -1.]), array([ 0., 0.]))
|
|
|
|
Or for a vectorized convolution:
|
|
|
|
>>> convolve = np.vectorize(np.convolve, signature='(n),(m)->(k)')
|
|
>>> convolve(np.eye(4), [1, 2, 1])
|
|
array([[1., 2., 1., 0., 0., 0.],
|
|
[0., 1., 2., 1., 0., 0.],
|
|
[0., 0., 1., 2., 1., 0.],
|
|
[0., 0., 0., 1., 2., 1.]])
|
|
|
|
"""
|
|
def __init__(self, pyfunc, otypes=None, doc=None, excluded=None,
|
|
cache=False, signature=None):
|
|
self.pyfunc = pyfunc
|
|
self.cache = cache
|
|
self.signature = signature
|
|
self._ufunc = None # Caching to improve default performance
|
|
|
|
if doc is None:
|
|
self.__doc__ = pyfunc.__doc__
|
|
else:
|
|
self.__doc__ = doc
|
|
|
|
if isinstance(otypes, str):
|
|
for char in otypes:
|
|
if char not in typecodes['All']:
|
|
raise ValueError("Invalid otype specified: %s" % (char,))
|
|
elif iterable(otypes):
|
|
otypes = ''.join([_nx.dtype(x).char for x in otypes])
|
|
elif otypes is not None:
|
|
raise ValueError("Invalid otype specification")
|
|
self.otypes = otypes
|
|
|
|
# Excluded variable support
|
|
if excluded is None:
|
|
excluded = set()
|
|
self.excluded = set(excluded)
|
|
|
|
if signature is not None:
|
|
self._in_and_out_core_dims = _parse_gufunc_signature(signature)
|
|
else:
|
|
self._in_and_out_core_dims = None
|
|
|
|
def __call__(self, *args, **kwargs):
|
|
"""
|
|
Return arrays with the results of `pyfunc` broadcast (vectorized) over
|
|
`args` and `kwargs` not in `excluded`.
|
|
"""
|
|
excluded = self.excluded
|
|
if not kwargs and not excluded:
|
|
func = self.pyfunc
|
|
vargs = args
|
|
else:
|
|
# The wrapper accepts only positional arguments: we use `names` and
|
|
# `inds` to mutate `the_args` and `kwargs` to pass to the original
|
|
# function.
|
|
nargs = len(args)
|
|
|
|
names = [_n for _n in kwargs if _n not in excluded]
|
|
inds = [_i for _i in range(nargs) if _i not in excluded]
|
|
the_args = list(args)
|
|
|
|
def func(*vargs):
|
|
for _n, _i in enumerate(inds):
|
|
the_args[_i] = vargs[_n]
|
|
kwargs.update(zip(names, vargs[len(inds):]))
|
|
return self.pyfunc(*the_args, **kwargs)
|
|
|
|
vargs = [args[_i] for _i in inds]
|
|
vargs.extend([kwargs[_n] for _n in names])
|
|
|
|
return self._vectorize_call(func=func, args=vargs)
|
|
|
|
def _get_ufunc_and_otypes(self, func, args):
|
|
"""Return (ufunc, otypes)."""
|
|
# frompyfunc will fail if args is empty
|
|
if not args:
|
|
raise ValueError('args can not be empty')
|
|
|
|
if self.otypes is not None:
|
|
otypes = self.otypes
|
|
nout = len(otypes)
|
|
|
|
# Note logic here: We only *use* self._ufunc if func is self.pyfunc
|
|
# even though we set self._ufunc regardless.
|
|
if func is self.pyfunc and self._ufunc is not None:
|
|
ufunc = self._ufunc
|
|
else:
|
|
ufunc = self._ufunc = frompyfunc(func, len(args), nout)
|
|
else:
|
|
# Get number of outputs and output types by calling the function on
|
|
# the first entries of args. We also cache the result to prevent
|
|
# the subsequent call when the ufunc is evaluated.
|
|
# Assumes that ufunc first evaluates the 0th elements in the input
|
|
# arrays (the input values are not checked to ensure this)
|
|
args = [asarray(arg) for arg in args]
|
|
if builtins.any(arg.size == 0 for arg in args):
|
|
raise ValueError('cannot call `vectorize` on size 0 inputs '
|
|
'unless `otypes` is set')
|
|
|
|
inputs = [arg.flat[0] for arg in args]
|
|
outputs = func(*inputs)
|
|
|
|
# Performance note: profiling indicates that -- for simple
|
|
# functions at least -- this wrapping can almost double the
|
|
# execution time.
|
|
# Hence we make it optional.
|
|
if self.cache:
|
|
_cache = [outputs]
|
|
|
|
def _func(*vargs):
|
|
if _cache:
|
|
return _cache.pop()
|
|
else:
|
|
return func(*vargs)
|
|
else:
|
|
_func = func
|
|
|
|
if isinstance(outputs, tuple):
|
|
nout = len(outputs)
|
|
else:
|
|
nout = 1
|
|
outputs = (outputs,)
|
|
|
|
otypes = ''.join([asarray(outputs[_k]).dtype.char
|
|
for _k in range(nout)])
|
|
|
|
# Performance note: profiling indicates that creating the ufunc is
|
|
# not a significant cost compared with wrapping so it seems not
|
|
# worth trying to cache this.
|
|
ufunc = frompyfunc(_func, len(args), nout)
|
|
|
|
return ufunc, otypes
|
|
|
|
def _vectorize_call(self, func, args):
|
|
"""Vectorized call to `func` over positional `args`."""
|
|
if self.signature is not None:
|
|
res = self._vectorize_call_with_signature(func, args)
|
|
elif not args:
|
|
res = func()
|
|
else:
|
|
ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args)
|
|
|
|
# Convert args to object arrays first
|
|
inputs = [array(a, copy=False, subok=True, dtype=object)
|
|
for a in args]
|
|
|
|
outputs = ufunc(*inputs)
|
|
|
|
if ufunc.nout == 1:
|
|
res = array(outputs, copy=False, subok=True, dtype=otypes[0])
|
|
else:
|
|
res = tuple([array(x, copy=False, subok=True, dtype=t)
|
|
for x, t in zip(outputs, otypes)])
|
|
return res
|
|
|
|
def _vectorize_call_with_signature(self, func, args):
|
|
"""Vectorized call over positional arguments with a signature."""
|
|
input_core_dims, output_core_dims = self._in_and_out_core_dims
|
|
|
|
if len(args) != len(input_core_dims):
|
|
raise TypeError('wrong number of positional arguments: '
|
|
'expected %r, got %r'
|
|
% (len(input_core_dims), len(args)))
|
|
args = tuple(asanyarray(arg) for arg in args)
|
|
|
|
broadcast_shape, dim_sizes = _parse_input_dimensions(
|
|
args, input_core_dims)
|
|
input_shapes = _calculate_shapes(broadcast_shape, dim_sizes,
|
|
input_core_dims)
|
|
args = [np.broadcast_to(arg, shape, subok=True)
|
|
for arg, shape in zip(args, input_shapes)]
|
|
|
|
outputs = None
|
|
otypes = self.otypes
|
|
nout = len(output_core_dims)
|
|
|
|
for index in np.ndindex(*broadcast_shape):
|
|
results = func(*(arg[index] for arg in args))
|
|
|
|
n_results = len(results) if isinstance(results, tuple) else 1
|
|
|
|
if nout != n_results:
|
|
raise ValueError(
|
|
'wrong number of outputs from pyfunc: expected %r, got %r'
|
|
% (nout, n_results))
|
|
|
|
if nout == 1:
|
|
results = (results,)
|
|
|
|
if outputs is None:
|
|
for result, core_dims in zip(results, output_core_dims):
|
|
_update_dim_sizes(dim_sizes, result, core_dims)
|
|
|
|
if otypes is None:
|
|
otypes = [asarray(result).dtype for result in results]
|
|
|
|
outputs = _create_arrays(broadcast_shape, dim_sizes,
|
|
output_core_dims, otypes)
|
|
|
|
for output, result in zip(outputs, results):
|
|
output[index] = result
|
|
|
|
if outputs is None:
|
|
# did not call the function even once
|
|
if otypes is None:
|
|
raise ValueError('cannot call `vectorize` on size 0 inputs '
|
|
'unless `otypes` is set')
|
|
if builtins.any(dim not in dim_sizes
|
|
for dims in output_core_dims
|
|
for dim in dims):
|
|
raise ValueError('cannot call `vectorize` with a signature '
|
|
'including new output dimensions on size 0 '
|
|
'inputs')
|
|
outputs = _create_arrays(broadcast_shape, dim_sizes,
|
|
output_core_dims, otypes)
|
|
|
|
return outputs[0] if nout == 1 else outputs
|
|
|
|
|
|
def _cov_dispatcher(m, y=None, rowvar=None, bias=None, ddof=None,
|
|
fweights=None, aweights=None):
|
|
return (m, y, fweights, aweights)
|
|
|
|
|
|
@array_function_dispatch(_cov_dispatcher)
|
|
def cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None,
|
|
aweights=None):
|
|
"""
|
|
Estimate a covariance matrix, given data and weights.
|
|
|
|
Covariance indicates the level to which two variables vary together.
|
|
If we examine N-dimensional samples, :math:`X = [x_1, x_2, ... x_N]^T`,
|
|
then the covariance matrix element :math:`C_{ij}` is the covariance of
|
|
:math:`x_i` and :math:`x_j`. The element :math:`C_{ii}` is the variance
|
|
of :math:`x_i`.
|
|
|
|
See the notes for an outline of the algorithm.
|
|
|
|
Parameters
|
|
----------
|
|
m : array_like
|
|
A 1-D or 2-D array containing multiple variables and observations.
|
|
Each row of `m` represents a variable, and each column a single
|
|
observation of all those variables. Also see `rowvar` below.
|
|
y : array_like, optional
|
|
An additional set of variables and observations. `y` has the same form
|
|
as that of `m`.
|
|
rowvar : bool, optional
|
|
If `rowvar` is True (default), then each row represents a
|
|
variable, with observations in the columns. Otherwise, the relationship
|
|
is transposed: each column represents a variable, while the rows
|
|
contain observations.
|
|
bias : bool, optional
|
|
Default normalization (False) is by ``(N - 1)``, where ``N`` is the
|
|
number of observations given (unbiased estimate). If `bias` is True,
|
|
then normalization is by ``N``. These values can be overridden by using
|
|
the keyword ``ddof`` in numpy versions >= 1.5.
|
|
ddof : int, optional
|
|
If not ``None`` the default value implied by `bias` is overridden.
|
|
Note that ``ddof=1`` will return the unbiased estimate, even if both
|
|
`fweights` and `aweights` are specified, and ``ddof=0`` will return
|
|
the simple average. See the notes for the details. The default value
|
|
is ``None``.
|
|
|
|
.. versionadded:: 1.5
|
|
fweights : array_like, int, optional
|
|
1-D array of integer frequency weights; the number of times each
|
|
observation vector should be repeated.
|
|
|
|
.. versionadded:: 1.10
|
|
aweights : array_like, optional
|
|
1-D array of observation vector weights. These relative weights are
|
|
typically large for observations considered "important" and smaller for
|
|
observations considered less "important". If ``ddof=0`` the array of
|
|
weights can be used to assign probabilities to observation vectors.
|
|
|
|
.. versionadded:: 1.10
|
|
|
|
Returns
|
|
-------
|
|
out : ndarray
|
|
The covariance matrix of the variables.
|
|
|
|
See Also
|
|
--------
|
|
corrcoef : Normalized covariance matrix
|
|
|
|
Notes
|
|
-----
|
|
Assume that the observations are in the columns of the observation
|
|
array `m` and let ``f = fweights`` and ``a = aweights`` for brevity. The
|
|
steps to compute the weighted covariance are as follows::
|
|
|
|
>>> m = np.arange(10, dtype=np.float64)
|
|
>>> f = np.arange(10) * 2
|
|
>>> a = np.arange(10) ** 2.
|
|
>>> ddof = 9 # N - 1
|
|
>>> w = f * a
|
|
>>> v1 = np.sum(w)
|
|
>>> v2 = np.sum(w * a)
|
|
>>> m -= np.sum(m * w, axis=None, keepdims=True) / v1
|
|
>>> cov = np.dot(m * w, m.T) * v1 / (v1**2 - ddof * v2)
|
|
|
|
Note that when ``a == 1``, the normalization factor
|
|
``v1 / (v1**2 - ddof * v2)`` goes over to ``1 / (np.sum(f) - ddof)``
|
|
as it should.
|
|
|
|
Examples
|
|
--------
|
|
Consider two variables, :math:`x_0` and :math:`x_1`, which
|
|
correlate perfectly, but in opposite directions:
|
|
|
|
>>> x = np.array([[0, 2], [1, 1], [2, 0]]).T
|
|
>>> x
|
|
array([[0, 1, 2],
|
|
[2, 1, 0]])
|
|
|
|
Note how :math:`x_0` increases while :math:`x_1` decreases. The covariance
|
|
matrix shows this clearly:
|
|
|
|
>>> np.cov(x)
|
|
array([[ 1., -1.],
|
|
[-1., 1.]])
|
|
|
|
Note that element :math:`C_{0,1}`, which shows the correlation between
|
|
:math:`x_0` and :math:`x_1`, is negative.
|
|
|
|
Further, note how `x` and `y` are combined:
|
|
|
|
>>> x = [-2.1, -1, 4.3]
|
|
>>> y = [3, 1.1, 0.12]
|
|
>>> X = np.stack((x, y), axis=0)
|
|
>>> np.cov(X)
|
|
array([[11.71 , -4.286 ], # may vary
|
|
[-4.286 , 2.144133]])
|
|
>>> np.cov(x, y)
|
|
array([[11.71 , -4.286 ], # may vary
|
|
[-4.286 , 2.144133]])
|
|
>>> np.cov(x)
|
|
array(11.71)
|
|
|
|
"""
|
|
# Check inputs
|
|
if ddof is not None and ddof != int(ddof):
|
|
raise ValueError(
|
|
"ddof must be integer")
|
|
|
|
# Handles complex arrays too
|
|
m = np.asarray(m)
|
|
if m.ndim > 2:
|
|
raise ValueError("m has more than 2 dimensions")
|
|
|
|
if y is None:
|
|
dtype = np.result_type(m, np.float64)
|
|
else:
|
|
y = np.asarray(y)
|
|
if y.ndim > 2:
|
|
raise ValueError("y has more than 2 dimensions")
|
|
dtype = np.result_type(m, y, np.float64)
|
|
|
|
X = array(m, ndmin=2, dtype=dtype)
|
|
if not rowvar and X.shape[0] != 1:
|
|
X = X.T
|
|
if X.shape[0] == 0:
|
|
return np.array([]).reshape(0, 0)
|
|
if y is not None:
|
|
y = array(y, copy=False, ndmin=2, dtype=dtype)
|
|
if not rowvar and y.shape[0] != 1:
|
|
y = y.T
|
|
X = np.concatenate((X, y), axis=0)
|
|
|
|
if ddof is None:
|
|
if bias == 0:
|
|
ddof = 1
|
|
else:
|
|
ddof = 0
|
|
|
|
# Get the product of frequencies and weights
|
|
w = None
|
|
if fweights is not None:
|
|
fweights = np.asarray(fweights, dtype=float)
|
|
if not np.all(fweights == np.around(fweights)):
|
|
raise TypeError(
|
|
"fweights must be integer")
|
|
if fweights.ndim > 1:
|
|
raise RuntimeError(
|
|
"cannot handle multidimensional fweights")
|
|
if fweights.shape[0] != X.shape[1]:
|
|
raise RuntimeError(
|
|
"incompatible numbers of samples and fweights")
|
|
if any(fweights < 0):
|
|
raise ValueError(
|
|
"fweights cannot be negative")
|
|
w = fweights
|
|
if aweights is not None:
|
|
aweights = np.asarray(aweights, dtype=float)
|
|
if aweights.ndim > 1:
|
|
raise RuntimeError(
|
|
"cannot handle multidimensional aweights")
|
|
if aweights.shape[0] != X.shape[1]:
|
|
raise RuntimeError(
|
|
"incompatible numbers of samples and aweights")
|
|
if any(aweights < 0):
|
|
raise ValueError(
|
|
"aweights cannot be negative")
|
|
if w is None:
|
|
w = aweights
|
|
else:
|
|
w *= aweights
|
|
|
|
avg, w_sum = average(X, axis=1, weights=w, returned=True)
|
|
w_sum = w_sum[0]
|
|
|
|
# Determine the normalization
|
|
if w is None:
|
|
fact = X.shape[1] - ddof
|
|
elif ddof == 0:
|
|
fact = w_sum
|
|
elif aweights is None:
|
|
fact = w_sum - ddof
|
|
else:
|
|
fact = w_sum - ddof*sum(w*aweights)/w_sum
|
|
|
|
if fact <= 0:
|
|
warnings.warn("Degrees of freedom <= 0 for slice",
|
|
RuntimeWarning, stacklevel=3)
|
|
fact = 0.0
|
|
|
|
X -= avg[:, None]
|
|
if w is None:
|
|
X_T = X.T
|
|
else:
|
|
X_T = (X*w).T
|
|
c = dot(X, X_T.conj())
|
|
c *= np.true_divide(1, fact)
|
|
return c.squeeze()
|
|
|
|
|
|
def _corrcoef_dispatcher(x, y=None, rowvar=None, bias=None, ddof=None):
|
|
return (x, y)
|
|
|
|
|
|
@array_function_dispatch(_corrcoef_dispatcher)
|
|
def corrcoef(x, y=None, rowvar=True, bias=np._NoValue, ddof=np._NoValue):
|
|
"""
|
|
Return Pearson product-moment correlation coefficients.
|
|
|
|
Please refer to the documentation for `cov` for more detail. The
|
|
relationship between the correlation coefficient matrix, `R`, and the
|
|
covariance matrix, `C`, is
|
|
|
|
.. math:: R_{ij} = \\frac{ C_{ij} } { \\sqrt{ C_{ii} * C_{jj} } }
|
|
|
|
The values of `R` are between -1 and 1, inclusive.
|
|
|
|
Parameters
|
|
----------
|
|
x : array_like
|
|
A 1-D or 2-D array containing multiple variables and observations.
|
|
Each row of `x` represents a variable, and each column a single
|
|
observation of all those variables. Also see `rowvar` below.
|
|
y : array_like, optional
|
|
An additional set of variables and observations. `y` has the same
|
|
shape as `x`.
|
|
rowvar : bool, optional
|
|
If `rowvar` is True (default), then each row represents a
|
|
variable, with observations in the columns. Otherwise, the relationship
|
|
is transposed: each column represents a variable, while the rows
|
|
contain observations.
|
|
bias : _NoValue, optional
|
|
Has no effect, do not use.
|
|
|
|
.. deprecated:: 1.10.0
|
|
ddof : _NoValue, optional
|
|
Has no effect, do not use.
|
|
|
|
.. deprecated:: 1.10.0
|
|
|
|
Returns
|
|
-------
|
|
R : ndarray
|
|
The correlation coefficient matrix of the variables.
|
|
|
|
See Also
|
|
--------
|
|
cov : Covariance matrix
|
|
|
|
Notes
|
|
-----
|
|
Due to floating point rounding the resulting array may not be Hermitian,
|
|
the diagonal elements may not be 1, and the elements may not satisfy the
|
|
inequality abs(a) <= 1. The real and imaginary parts are clipped to the
|
|
interval [-1, 1] in an attempt to improve on that situation but is not
|
|
much help in the complex case.
|
|
|
|
This function accepts but discards arguments `bias` and `ddof`. This is
|
|
for backwards compatibility with previous versions of this function. These
|
|
arguments had no effect on the return values of the function and can be
|
|
safely ignored in this and previous versions of numpy.
|
|
|
|
"""
|
|
if bias is not np._NoValue or ddof is not np._NoValue:
|
|
# 2015-03-15, 1.10
|
|
warnings.warn('bias and ddof have no effect and are deprecated',
|
|
DeprecationWarning, stacklevel=3)
|
|
c = cov(x, y, rowvar)
|
|
try:
|
|
d = diag(c)
|
|
except ValueError:
|
|
# scalar covariance
|
|
# nan if incorrect value (nan, inf, 0), 1 otherwise
|
|
return c / c
|
|
stddev = sqrt(d.real)
|
|
c /= stddev[:, None]
|
|
c /= stddev[None, :]
|
|
|
|
# Clip real and imaginary parts to [-1, 1]. This does not guarantee
|
|
# abs(a[i,j]) <= 1 for complex arrays, but is the best we can do without
|
|
# excessive work.
|
|
np.clip(c.real, -1, 1, out=c.real)
|
|
if np.iscomplexobj(c):
|
|
np.clip(c.imag, -1, 1, out=c.imag)
|
|
|
|
return c
|
|
|
|
|
|
@set_module('numpy')
|
|
def blackman(M):
|
|
"""
|
|
Return the Blackman window.
|
|
|
|
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.
|
|
|
|
Parameters
|
|
----------
|
|
M : int
|
|
Number of points in the output window. If zero or less, an empty
|
|
array is returned.
|
|
|
|
Returns
|
|
-------
|
|
out : ndarray
|
|
The window, with the maximum value normalized to one (the value one
|
|
appears only if the number of samples is odd).
|
|
|
|
See Also
|
|
--------
|
|
bartlett, hamming, hanning, kaiser
|
|
|
|
Notes
|
|
-----
|
|
The Blackman window is defined as
|
|
|
|
.. math:: w(n) = 0.42 - 0.5 \\cos(2\\pi n/M) + 0.08 \\cos(4\\pi n/M)
|
|
|
|
Most references to the Blackman 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
|
|
"removing the foot", i.e. smoothing discontinuities at the beginning
|
|
and end of the sampled signal) or tapering function. It is known as a
|
|
"near optimal" tapering function, almost as good (by some measures)
|
|
as the kaiser window.
|
|
|
|
References
|
|
----------
|
|
Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra,
|
|
Dover Publications, New York.
|
|
|
|
Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing.
|
|
Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471.
|
|
|
|
Examples
|
|
--------
|
|
>>> import matplotlib.pyplot as plt
|
|
>>> np.blackman(12)
|
|
array([-1.38777878e-17, 3.26064346e-02, 1.59903635e-01, # may vary
|
|
4.14397981e-01, 7.36045180e-01, 9.67046769e-01,
|
|
9.67046769e-01, 7.36045180e-01, 4.14397981e-01,
|
|
1.59903635e-01, 3.26064346e-02, -1.38777878e-17])
|
|
|
|
Plot the window and the frequency response:
|
|
|
|
>>> from numpy.fft import fft, fftshift
|
|
>>> window = np.blackman(51)
|
|
>>> plt.plot(window)
|
|
[<matplotlib.lines.Line2D object at 0x...>]
|
|
>>> plt.title("Blackman window")
|
|
Text(0.5, 1.0, 'Blackman window')
|
|
>>> plt.ylabel("Amplitude")
|
|
Text(0, 0.5, 'Amplitude')
|
|
>>> plt.xlabel("Sample")
|
|
Text(0.5, 0, 'Sample')
|
|
>>> plt.show()
|
|
|
|
>>> plt.figure()
|
|
<Figure size 640x480 with 0 Axes>
|
|
>>> A = fft(window, 2048) / 25.5
|
|
>>> mag = np.abs(fftshift(A))
|
|
>>> freq = np.linspace(-0.5, 0.5, len(A))
|
|
>>> with np.errstate(divide='ignore', invalid='ignore'):
|
|
... response = 20 * np.log10(mag)
|
|
...
|
|
>>> response = np.clip(response, -100, 100)
|
|
>>> plt.plot(freq, response)
|
|
[<matplotlib.lines.Line2D object at 0x...>]
|
|
>>> plt.title("Frequency response of Blackman window")
|
|
Text(0.5, 1.0, 'Frequency response of Blackman window')
|
|
>>> plt.ylabel("Magnitude [dB]")
|
|
Text(0, 0.5, 'Magnitude [dB]')
|
|
>>> plt.xlabel("Normalized frequency [cycles per sample]")
|
|
Text(0.5, 0, 'Normalized frequency [cycles per sample]')
|
|
>>> _ = plt.axis('tight')
|
|
>>> plt.show()
|
|
|
|
"""
|
|
if M < 1:
|
|
return array([])
|
|
if M == 1:
|
|
return ones(1, float)
|
|
n = arange(0, M)
|
|
return 0.42 - 0.5*cos(2.0*pi*n/(M-1)) + 0.08*cos(4.0*pi*n/(M-1))
|
|
|
|
|
|
@set_module('numpy')
|
|
def bartlett(M):
|
|
"""
|
|
Return the Bartlett window.
|
|
|
|
The Bartlett window is very similar to a triangular window, except
|
|
that the end points are at zero. It is often used in signal
|
|
processing for tapering a signal, without generating too much
|
|
ripple in the frequency domain.
|
|
|
|
Parameters
|
|
----------
|
|
M : int
|
|
Number of points in the output window. If zero or less, an
|
|
empty array is returned.
|
|
|
|
Returns
|
|
-------
|
|
out : array
|
|
The triangular window, with the maximum value normalized to one
|
|
(the value one appears only if the number of samples is odd), with
|
|
the first and last samples equal to zero.
|
|
|
|
See Also
|
|
--------
|
|
blackman, hamming, hanning, kaiser
|
|
|
|
Notes
|
|
-----
|
|
The Bartlett window is defined as
|
|
|
|
.. math:: w(n) = \\frac{2}{M-1} \\left(
|
|
\\frac{M-1}{2} - \\left|n - \\frac{M-1}{2}\\right|
|
|
\\right)
|
|
|
|
Most references to the Bartlett window come from the signal
|
|
processing literature, where it is used as one of many windowing
|
|
functions for smoothing values. Note that convolution with this
|
|
window produces linear interpolation. It is also known as an
|
|
apodization (which means"removing the foot", i.e. smoothing
|
|
discontinuities at the beginning and end of the sampled signal) or
|
|
tapering function. The fourier transform of the Bartlett is the product
|
|
of two sinc functions.
|
|
Note the excellent discussion in Kanasewich.
|
|
|
|
References
|
|
----------
|
|
.. [1] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra",
|
|
Biometrika 37, 1-16, 1950.
|
|
.. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
|
|
The University of Alberta Press, 1975, pp. 109-110.
|
|
.. [3] A.V. Oppenheim and R.W. Schafer, "Discrete-Time Signal
|
|
Processing", Prentice-Hall, 1999, pp. 468-471.
|
|
.. [4] Wikipedia, "Window function",
|
|
https://en.wikipedia.org/wiki/Window_function
|
|
.. [5] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
|
|
"Numerical Recipes", Cambridge University Press, 1986, page 429.
|
|
|
|
Examples
|
|
--------
|
|
>>> import matplotlib.pyplot as plt
|
|
>>> np.bartlett(12)
|
|
array([ 0. , 0.18181818, 0.36363636, 0.54545455, 0.72727273, # may vary
|
|
0.90909091, 0.90909091, 0.72727273, 0.54545455, 0.36363636,
|
|
0.18181818, 0. ])
|
|
|
|
Plot the window and its frequency response (requires SciPy and matplotlib):
|
|
|
|
>>> from numpy.fft import fft, fftshift
|
|
>>> window = np.bartlett(51)
|
|
>>> plt.plot(window)
|
|
[<matplotlib.lines.Line2D object at 0x...>]
|
|
>>> plt.title("Bartlett window")
|
|
Text(0.5, 1.0, 'Bartlett window')
|
|
>>> plt.ylabel("Amplitude")
|
|
Text(0, 0.5, 'Amplitude')
|
|
>>> plt.xlabel("Sample")
|
|
Text(0.5, 0, 'Sample')
|
|
>>> plt.show()
|
|
|
|
>>> plt.figure()
|
|
<Figure size 640x480 with 0 Axes>
|
|
>>> A = fft(window, 2048) / 25.5
|
|
>>> mag = np.abs(fftshift(A))
|
|
>>> freq = np.linspace(-0.5, 0.5, len(A))
|
|
>>> with np.errstate(divide='ignore', invalid='ignore'):
|
|
... response = 20 * np.log10(mag)
|
|
...
|
|
>>> response = np.clip(response, -100, 100)
|
|
>>> plt.plot(freq, response)
|
|
[<matplotlib.lines.Line2D object at 0x...>]
|
|
>>> plt.title("Frequency response of Bartlett window")
|
|
Text(0.5, 1.0, 'Frequency response of Bartlett window')
|
|
>>> plt.ylabel("Magnitude [dB]")
|
|
Text(0, 0.5, 'Magnitude [dB]')
|
|
>>> plt.xlabel("Normalized frequency [cycles per sample]")
|
|
Text(0.5, 0, 'Normalized frequency [cycles per sample]')
|
|
>>> _ = plt.axis('tight')
|
|
>>> plt.show()
|
|
|
|
"""
|
|
if M < 1:
|
|
return array([])
|
|
if M == 1:
|
|
return ones(1, float)
|
|
n = arange(0, M)
|
|
return where(less_equal(n, (M-1)/2.0), 2.0*n/(M-1), 2.0 - 2.0*n/(M-1))
|
|
|
|
|
|
@set_module('numpy')
|
|
def hanning(M):
|
|
"""
|
|
Return the Hanning window.
|
|
|
|
The Hanning window is a taper formed by using a weighted cosine.
|
|
|
|
Parameters
|
|
----------
|
|
M : int
|
|
Number of points in the output window. If zero or less, an
|
|
empty array is returned.
|
|
|
|
Returns
|
|
-------
|
|
out : ndarray, shape(M,)
|
|
The window, with the maximum value normalized to one (the value
|
|
one appears only if `M` is odd).
|
|
|
|
See Also
|
|
--------
|
|
bartlett, blackman, hamming, kaiser
|
|
|
|
Notes
|
|
-----
|
|
The Hanning window is defined as
|
|
|
|
.. 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 von 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
|
|
smoothing values. It is also known as an apodization (which means
|
|
"removing the foot", i.e. smoothing discontinuities at the beginning
|
|
and end of the sampled signal) or tapering function.
|
|
|
|
References
|
|
----------
|
|
.. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
|
|
spectra, Dover Publications, New York.
|
|
.. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
|
|
The University of Alberta Press, 1975, pp. 106-108.
|
|
.. [3] Wikipedia, "Window function",
|
|
https://en.wikipedia.org/wiki/Window_function
|
|
.. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
|
|
"Numerical Recipes", Cambridge University Press, 1986, page 425.
|
|
|
|
Examples
|
|
--------
|
|
>>> np.hanning(12)
|
|
array([0. , 0.07937323, 0.29229249, 0.57115742, 0.82743037,
|
|
0.97974649, 0.97974649, 0.82743037, 0.57115742, 0.29229249,
|
|
0.07937323, 0. ])
|
|
|
|
Plot the window and its frequency response:
|
|
|
|
>>> import matplotlib.pyplot as plt
|
|
>>> from numpy.fft import fft, fftshift
|
|
>>> window = np.hanning(51)
|
|
>>> plt.plot(window)
|
|
[<matplotlib.lines.Line2D object at 0x...>]
|
|
>>> plt.title("Hann window")
|
|
Text(0.5, 1.0, 'Hann window')
|
|
>>> plt.ylabel("Amplitude")
|
|
Text(0, 0.5, 'Amplitude')
|
|
>>> plt.xlabel("Sample")
|
|
Text(0.5, 0, 'Sample')
|
|
>>> plt.show()
|
|
|
|
>>> plt.figure()
|
|
<Figure size 640x480 with 0 Axes>
|
|
>>> A = fft(window, 2048) / 25.5
|
|
>>> mag = np.abs(fftshift(A))
|
|
>>> freq = np.linspace(-0.5, 0.5, len(A))
|
|
>>> with np.errstate(divide='ignore', invalid='ignore'):
|
|
... response = 20 * np.log10(mag)
|
|
...
|
|
>>> response = np.clip(response, -100, 100)
|
|
>>> plt.plot(freq, response)
|
|
[<matplotlib.lines.Line2D object at 0x...>]
|
|
>>> plt.title("Frequency response of the Hann window")
|
|
Text(0.5, 1.0, 'Frequency response of the Hann window')
|
|
>>> plt.ylabel("Magnitude [dB]")
|
|
Text(0, 0.5, 'Magnitude [dB]')
|
|
>>> plt.xlabel("Normalized frequency [cycles per sample]")
|
|
Text(0.5, 0, 'Normalized frequency [cycles per sample]')
|
|
>>> plt.axis('tight')
|
|
...
|
|
>>> plt.show()
|
|
|
|
"""
|
|
if M < 1:
|
|
return array([])
|
|
if M == 1:
|
|
return ones(1, float)
|
|
n = arange(0, M)
|
|
return 0.5 - 0.5*cos(2.0*pi*n/(M-1))
|
|
|
|
|
|
@set_module('numpy')
|
|
def hamming(M):
|
|
"""
|
|
Return the Hamming window.
|
|
|
|
The Hamming window is a taper formed by using a weighted cosine.
|
|
|
|
Parameters
|
|
----------
|
|
M : int
|
|
Number of points in the output window. If zero or less, an
|
|
empty array is returned.
|
|
|
|
Returns
|
|
-------
|
|
out : ndarray
|
|
The window, with the maximum value normalized to one (the value
|
|
one appears only if the number of samples is odd).
|
|
|
|
See Also
|
|
--------
|
|
bartlett, blackman, hanning, kaiser
|
|
|
|
Notes
|
|
-----
|
|
The Hamming window is defined as
|
|
|
|
.. 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.
|
|
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
|
|
"removing the foot", i.e. smoothing discontinuities at the beginning
|
|
and end of the sampled signal) or tapering function.
|
|
|
|
References
|
|
----------
|
|
.. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
|
|
spectra, Dover Publications, New York.
|
|
.. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
|
|
University of Alberta Press, 1975, pp. 109-110.
|
|
.. [3] Wikipedia, "Window function",
|
|
https://en.wikipedia.org/wiki/Window_function
|
|
.. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
|
|
"Numerical Recipes", Cambridge University Press, 1986, page 425.
|
|
|
|
Examples
|
|
--------
|
|
>>> np.hamming(12)
|
|
array([ 0.08 , 0.15302337, 0.34890909, 0.60546483, 0.84123594, # may vary
|
|
0.98136677, 0.98136677, 0.84123594, 0.60546483, 0.34890909,
|
|
0.15302337, 0.08 ])
|
|
|
|
Plot the window and the frequency response:
|
|
|
|
>>> import matplotlib.pyplot as plt
|
|
>>> from numpy.fft import fft, fftshift
|
|
>>> window = np.hamming(51)
|
|
>>> plt.plot(window)
|
|
[<matplotlib.lines.Line2D object at 0x...>]
|
|
>>> plt.title("Hamming window")
|
|
Text(0.5, 1.0, 'Hamming window')
|
|
>>> plt.ylabel("Amplitude")
|
|
Text(0, 0.5, 'Amplitude')
|
|
>>> plt.xlabel("Sample")
|
|
Text(0.5, 0, 'Sample')
|
|
>>> plt.show()
|
|
|
|
>>> plt.figure()
|
|
<Figure size 640x480 with 0 Axes>
|
|
>>> A = fft(window, 2048) / 25.5
|
|
>>> mag = np.abs(fftshift(A))
|
|
>>> freq = np.linspace(-0.5, 0.5, len(A))
|
|
>>> response = 20 * np.log10(mag)
|
|
>>> response = np.clip(response, -100, 100)
|
|
>>> plt.plot(freq, response)
|
|
[<matplotlib.lines.Line2D object at 0x...>]
|
|
>>> plt.title("Frequency response of Hamming window")
|
|
Text(0.5, 1.0, 'Frequency response of Hamming window')
|
|
>>> plt.ylabel("Magnitude [dB]")
|
|
Text(0, 0.5, 'Magnitude [dB]')
|
|
>>> plt.xlabel("Normalized frequency [cycles per sample]")
|
|
Text(0.5, 0, 'Normalized frequency [cycles per sample]')
|
|
>>> plt.axis('tight')
|
|
...
|
|
>>> plt.show()
|
|
|
|
"""
|
|
if M < 1:
|
|
return array([])
|
|
if M == 1:
|
|
return ones(1, float)
|
|
n = arange(0, M)
|
|
return 0.54 - 0.46*cos(2.0*pi*n/(M-1))
|
|
|
|
## Code from cephes for i0
|
|
|
|
_i0A = [
|
|
-4.41534164647933937950E-18,
|
|
3.33079451882223809783E-17,
|
|
-2.43127984654795469359E-16,
|
|
1.71539128555513303061E-15,
|
|
-1.16853328779934516808E-14,
|
|
7.67618549860493561688E-14,
|
|
-4.85644678311192946090E-13,
|
|
2.95505266312963983461E-12,
|
|
-1.72682629144155570723E-11,
|
|
9.67580903537323691224E-11,
|
|
-5.18979560163526290666E-10,
|
|
2.65982372468238665035E-9,
|
|
-1.30002500998624804212E-8,
|
|
6.04699502254191894932E-8,
|
|
-2.67079385394061173391E-7,
|
|
1.11738753912010371815E-6,
|
|
-4.41673835845875056359E-6,
|
|
1.64484480707288970893E-5,
|
|
-5.75419501008210370398E-5,
|
|
1.88502885095841655729E-4,
|
|
-5.76375574538582365885E-4,
|
|
1.63947561694133579842E-3,
|
|
-4.32430999505057594430E-3,
|
|
1.05464603945949983183E-2,
|
|
-2.37374148058994688156E-2,
|
|
4.93052842396707084878E-2,
|
|
-9.49010970480476444210E-2,
|
|
1.71620901522208775349E-1,
|
|
-3.04682672343198398683E-1,
|
|
6.76795274409476084995E-1
|
|
]
|
|
|
|
_i0B = [
|
|
-7.23318048787475395456E-18,
|
|
-4.83050448594418207126E-18,
|
|
4.46562142029675999901E-17,
|
|
3.46122286769746109310E-17,
|
|
-2.82762398051658348494E-16,
|
|
-3.42548561967721913462E-16,
|
|
1.77256013305652638360E-15,
|
|
3.81168066935262242075E-15,
|
|
-9.55484669882830764870E-15,
|
|
-4.15056934728722208663E-14,
|
|
1.54008621752140982691E-14,
|
|
3.85277838274214270114E-13,
|
|
7.18012445138366623367E-13,
|
|
-1.79417853150680611778E-12,
|
|
-1.32158118404477131188E-11,
|
|
-3.14991652796324136454E-11,
|
|
1.18891471078464383424E-11,
|
|
4.94060238822496958910E-10,
|
|
3.39623202570838634515E-9,
|
|
2.26666899049817806459E-8,
|
|
2.04891858946906374183E-7,
|
|
2.89137052083475648297E-6,
|
|
6.88975834691682398426E-5,
|
|
3.36911647825569408990E-3,
|
|
8.04490411014108831608E-1
|
|
]
|
|
|
|
|
|
def _chbevl(x, vals):
|
|
b0 = vals[0]
|
|
b1 = 0.0
|
|
|
|
for i in range(1, len(vals)):
|
|
b2 = b1
|
|
b1 = b0
|
|
b0 = x*b1 - b2 + vals[i]
|
|
|
|
return 0.5*(b0 - b2)
|
|
|
|
|
|
def _i0_1(x):
|
|
return exp(x) * _chbevl(x/2.0-2, _i0A)
|
|
|
|
|
|
def _i0_2(x):
|
|
return exp(x) * _chbevl(32.0/x - 2.0, _i0B) / sqrt(x)
|
|
|
|
|
|
def _i0_dispatcher(x):
|
|
return (x,)
|
|
|
|
|
|
@array_function_dispatch(_i0_dispatcher)
|
|
def i0(x):
|
|
"""
|
|
Modified Bessel function of the first kind, order 0.
|
|
|
|
Usually denoted :math:`I_0`. This function does broadcast, but will *not*
|
|
"up-cast" int dtype arguments unless accompanied by at least one float or
|
|
complex dtype argument (see Raises below).
|
|
|
|
Parameters
|
|
----------
|
|
x : array_like, dtype float or complex
|
|
Argument of the Bessel function.
|
|
|
|
Returns
|
|
-------
|
|
out : ndarray, shape = x.shape, dtype = x.dtype
|
|
The modified Bessel function evaluated at each of the elements of `x`.
|
|
|
|
Raises
|
|
------
|
|
TypeError: array cannot be safely cast to required type
|
|
If argument consists exclusively of int dtypes.
|
|
|
|
See Also
|
|
--------
|
|
scipy.special.i0, scipy.special.iv, scipy.special.ive
|
|
|
|
Notes
|
|
-----
|
|
The scipy implementation is recommended over this function: it is a
|
|
proper ufunc written in C, and more than an order of magnitude faster.
|
|
|
|
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).
|
|
|
|
References
|
|
----------
|
|
.. [1] C. W. Clenshaw, "Chebyshev series for mathematical functions", in
|
|
*National Physical Laboratory Mathematical Tables*, vol. 5, London:
|
|
Her Majesty's Stationery Office, 1962.
|
|
.. [2] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical
|
|
Functions*, 10th printing, New York: Dover, 1964, pp. 379.
|
|
http://www.math.sfu.ca/~cbm/aands/page_379.htm
|
|
.. [3] http://kobesearch.cpan.org/htdocs/Math-Cephes/Math/Cephes.html
|
|
|
|
Examples
|
|
--------
|
|
>>> np.i0(0.)
|
|
array(1.0) # may vary
|
|
>>> np.i0([0., 1. + 2j])
|
|
array([ 1.00000000+0.j , 0.18785373+0.64616944j]) # may vary
|
|
|
|
"""
|
|
x = np.asanyarray(x)
|
|
x = np.abs(x)
|
|
return piecewise(x, [x <= 8.0], [_i0_1, _i0_2])
|
|
|
|
## End of cephes code for i0
|
|
|
|
|
|
@set_module('numpy')
|
|
def kaiser(M, beta):
|
|
"""
|
|
Return the Kaiser window.
|
|
|
|
The Kaiser window is a taper formed by using a Bessel function.
|
|
|
|
Parameters
|
|
----------
|
|
M : int
|
|
Number of points in the output window. If zero or less, an
|
|
empty array is returned.
|
|
beta : float
|
|
Shape parameter for window.
|
|
|
|
Returns
|
|
-------
|
|
out : array
|
|
The window, with the maximum value normalized to one (the value
|
|
one appears only if the number of samples is odd).
|
|
|
|
See Also
|
|
--------
|
|
bartlett, blackman, hamming, hanning
|
|
|
|
Notes
|
|
-----
|
|
The Kaiser window is defined as
|
|
|
|
.. math:: w(n) = I_0\\left( \\beta \\sqrt{1-\\frac{4n^2}{(M-1)^2}}
|
|
\\right)/I_0(\\beta)
|
|
|
|
with
|
|
|
|
.. math:: \\quad -\\frac{M-1}{2} \\leq n \\leq \\frac{M-1}{2},
|
|
|
|
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 can approximate many other windows by varying the beta
|
|
parameter.
|
|
|
|
==== =======================
|
|
beta Window shape
|
|
==== =======================
|
|
0 Rectangular
|
|
5 Similar to a Hamming
|
|
6 Similar to a Hanning
|
|
8.6 Similar to a Blackman
|
|
==== =======================
|
|
|
|
A beta value of 14 is probably a good starting point. Note that as beta
|
|
gets large, the window narrows, and so the number of samples needs to be
|
|
large enough to sample the increasingly narrow spike, otherwise NaNs will
|
|
get returned.
|
|
|
|
Most references to the Kaiser 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
|
|
"removing the foot", i.e. smoothing discontinuities at the beginning
|
|
and end of the sampled signal) or tapering function.
|
|
|
|
References
|
|
----------
|
|
.. [1] J. F. Kaiser, "Digital Filters" - Ch 7 in "Systems analysis by
|
|
digital computer", Editors: F.F. Kuo and J.F. Kaiser, p 218-285.
|
|
John Wiley and Sons, New York, (1966).
|
|
.. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
|
|
University of Alberta Press, 1975, pp. 177-178.
|
|
.. [3] Wikipedia, "Window function",
|
|
https://en.wikipedia.org/wiki/Window_function
|
|
|
|
Examples
|
|
--------
|
|
>>> import matplotlib.pyplot as plt
|
|
>>> np.kaiser(12, 14)
|
|
array([7.72686684e-06, 3.46009194e-03, 4.65200189e-02, # may vary
|
|
2.29737120e-01, 5.99885316e-01, 9.45674898e-01,
|
|
9.45674898e-01, 5.99885316e-01, 2.29737120e-01,
|
|
4.65200189e-02, 3.46009194e-03, 7.72686684e-06])
|
|
|
|
|
|
Plot the window and the frequency response:
|
|
|
|
>>> from numpy.fft import fft, fftshift
|
|
>>> window = np.kaiser(51, 14)
|
|
>>> plt.plot(window)
|
|
[<matplotlib.lines.Line2D object at 0x...>]
|
|
>>> plt.title("Kaiser window")
|
|
Text(0.5, 1.0, 'Kaiser window')
|
|
>>> plt.ylabel("Amplitude")
|
|
Text(0, 0.5, 'Amplitude')
|
|
>>> plt.xlabel("Sample")
|
|
Text(0.5, 0, 'Sample')
|
|
>>> plt.show()
|
|
|
|
>>> plt.figure()
|
|
<Figure size 640x480 with 0 Axes>
|
|
>>> A = fft(window, 2048) / 25.5
|
|
>>> mag = np.abs(fftshift(A))
|
|
>>> freq = np.linspace(-0.5, 0.5, len(A))
|
|
>>> response = 20 * np.log10(mag)
|
|
>>> response = np.clip(response, -100, 100)
|
|
>>> plt.plot(freq, response)
|
|
[<matplotlib.lines.Line2D object at 0x...>]
|
|
>>> plt.title("Frequency response of Kaiser window")
|
|
Text(0.5, 1.0, 'Frequency response of Kaiser window')
|
|
>>> plt.ylabel("Magnitude [dB]")
|
|
Text(0, 0.5, 'Magnitude [dB]')
|
|
>>> plt.xlabel("Normalized frequency [cycles per sample]")
|
|
Text(0.5, 0, 'Normalized frequency [cycles per sample]')
|
|
>>> plt.axis('tight')
|
|
(-0.5, 0.5, -100.0, ...) # may vary
|
|
>>> plt.show()
|
|
|
|
"""
|
|
from numpy.dual import i0
|
|
if M == 1:
|
|
return np.array([1.])
|
|
n = arange(0, M)
|
|
alpha = (M-1)/2.0
|
|
return i0(beta * sqrt(1-((n-alpha)/alpha)**2.0))/i0(float(beta))
|
|
|
|
|
|
def _sinc_dispatcher(x):
|
|
return (x,)
|
|
|
|
|
|
@array_function_dispatch(_sinc_dispatcher)
|
|
def sinc(x):
|
|
"""
|
|
Return the sinc function.
|
|
|
|
The sinc function is :math:`\\sin(\\pi x)/(\\pi x)`.
|
|
|
|
Parameters
|
|
----------
|
|
x : ndarray
|
|
Array (possibly multi-dimensional) of values for which to to
|
|
calculate ``sinc(x)``.
|
|
|
|
Returns
|
|
-------
|
|
out : ndarray
|
|
``sinc(x)``, which has the same shape as the input.
|
|
|
|
Notes
|
|
-----
|
|
``sinc(0)`` is the limit value 1.
|
|
|
|
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.
|
|
|
|
For bandlimited interpolation of discrete-time signals, the ideal
|
|
interpolation kernel is proportional to the sinc function.
|
|
|
|
References
|
|
----------
|
|
.. [1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram Web
|
|
Resource. http://mathworld.wolfram.com/SincFunction.html
|
|
.. [2] Wikipedia, "Sinc function",
|
|
https://en.wikipedia.org/wiki/Sinc_function
|
|
|
|
Examples
|
|
--------
|
|
>>> import matplotlib.pyplot as plt
|
|
>>> x = np.linspace(-4, 4, 41)
|
|
>>> np.sinc(x)
|
|
array([-3.89804309e-17, -4.92362781e-02, -8.40918587e-02, # may vary
|
|
-8.90384387e-02, -5.84680802e-02, 3.89804309e-17,
|
|
6.68206631e-02, 1.16434881e-01, 1.26137788e-01,
|
|
8.50444803e-02, -3.89804309e-17, -1.03943254e-01,
|
|
-1.89206682e-01, -2.16236208e-01, -1.55914881e-01,
|
|
3.89804309e-17, 2.33872321e-01, 5.04551152e-01,
|
|
7.56826729e-01, 9.35489284e-01, 1.00000000e+00,
|
|
9.35489284e-01, 7.56826729e-01, 5.04551152e-01,
|
|
2.33872321e-01, 3.89804309e-17, -1.55914881e-01,
|
|
-2.16236208e-01, -1.89206682e-01, -1.03943254e-01,
|
|
-3.89804309e-17, 8.50444803e-02, 1.26137788e-01,
|
|
1.16434881e-01, 6.68206631e-02, 3.89804309e-17,
|
|
-5.84680802e-02, -8.90384387e-02, -8.40918587e-02,
|
|
-4.92362781e-02, -3.89804309e-17])
|
|
|
|
>>> plt.plot(x, np.sinc(x))
|
|
[<matplotlib.lines.Line2D object at 0x...>]
|
|
>>> plt.title("Sinc Function")
|
|
Text(0.5, 1.0, 'Sinc Function')
|
|
>>> plt.ylabel("Amplitude")
|
|
Text(0, 0.5, 'Amplitude')
|
|
>>> plt.xlabel("X")
|
|
Text(0.5, 0, 'X')
|
|
>>> plt.show()
|
|
|
|
It works in 2-D as well:
|
|
|
|
>>> x = np.linspace(-4, 4, 401)
|
|
>>> xx = np.outer(x, x)
|
|
>>> plt.imshow(np.sinc(xx))
|
|
<matplotlib.image.AxesImage object at 0x...>
|
|
|
|
"""
|
|
x = np.asanyarray(x)
|
|
y = pi * where(x == 0, 1.0e-20, x)
|
|
return sin(y)/y
|
|
|
|
|
|
def _msort_dispatcher(a):
|
|
return (a,)
|
|
|
|
|
|
@array_function_dispatch(_msort_dispatcher)
|
|
def msort(a):
|
|
"""
|
|
Return a copy of an array sorted along the first axis.
|
|
|
|
Parameters
|
|
----------
|
|
a : array_like
|
|
Array to be sorted.
|
|
|
|
Returns
|
|
-------
|
|
sorted_array : ndarray
|
|
Array of the same type and shape as `a`.
|
|
|
|
See Also
|
|
--------
|
|
sort
|
|
|
|
Notes
|
|
-----
|
|
``np.msort(a)`` is equivalent to ``np.sort(a, axis=0)``.
|
|
|
|
"""
|
|
b = array(a, subok=True, copy=True)
|
|
b.sort(0)
|
|
return b
|
|
|
|
|
|
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 capable of receiving a single axis argument.
|
|
It 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
|
|
axis = _nx.normalize_axis_tuple(axis, nd)
|
|
|
|
for ax in axis:
|
|
keepdim[ax] = 1
|
|
|
|
if len(axis) == 1:
|
|
kwargs['axis'] = axis[0]
|
|
else:
|
|
keep = set(range(nd)) - set(axis)
|
|
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
|
|
keepdim = tuple(keepdim)
|
|
else:
|
|
keepdim = (1,) * a.ndim
|
|
|
|
r = func(a, **kwargs)
|
|
return r, keepdim
|
|
|
|
|
|
def _median_dispatcher(
|
|
a, axis=None, out=None, overwrite_input=None, keepdims=None):
|
|
return (a, out)
|
|
|
|
|
|
@array_function_dispatch(_median_dispatcher)
|
|
def median(a, axis=None, out=None, overwrite_input=False, keepdims=False):
|
|
"""
|
|
Compute the median along the specified axis.
|
|
|
|
Returns the median of the array elements.
|
|
|
|
Parameters
|
|
----------
|
|
a : array_like
|
|
Input array or object that can be converted to an array.
|
|
axis : {int, sequence of int, None}, optional
|
|
Axis or axes along which the medians are computed. The default
|
|
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. If `overwrite_input` is ``True`` and `a` 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. If the input contains integers
|
|
or floats smaller than ``float64``, then the output data-type is
|
|
``np.float64``. Otherwise, the data-type of the output is the
|
|
same as that of the input. If `out` is specified, that array is
|
|
returned instead.
|
|
|
|
See Also
|
|
--------
|
|
mean, 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, and the average of the
|
|
two middle values of ``V_sorted`` when ``N`` is even.
|
|
|
|
Examples
|
|
--------
|
|
>>> a = np.array([[10, 7, 4], [3, 2, 1]])
|
|
>>> a
|
|
array([[10, 7, 4],
|
|
[ 3, 2, 1]])
|
|
>>> np.median(a)
|
|
3.5
|
|
>>> np.median(a, axis=0)
|
|
array([6.5, 4.5, 2.5])
|
|
>>> np.median(a, axis=1)
|
|
array([7., 2.])
|
|
>>> m = np.median(a, axis=0)
|
|
>>> out = np.zeros_like(m)
|
|
>>> np.median(a, axis=0, out=m)
|
|
array([6.5, 4.5, 2.5])
|
|
>>> m
|
|
array([6.5, 4.5, 2.5])
|
|
>>> b = a.copy()
|
|
>>> np.median(b, axis=1, overwrite_input=True)
|
|
array([7., 2.])
|
|
>>> assert not np.all(a==b)
|
|
>>> b = a.copy()
|
|
>>> np.median(b, axis=None, overwrite_input=True)
|
|
3.5
|
|
>>> 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)
|
|
|
|
# Set the partition indexes
|
|
if axis is None:
|
|
sz = a.size
|
|
else:
|
|
sz = a.shape[axis]
|
|
if sz % 2 == 0:
|
|
szh = sz // 2
|
|
kth = [szh - 1, szh]
|
|
else:
|
|
kth = [(sz - 1) // 2]
|
|
# Check if the array contains any nan's
|
|
if np.issubdtype(a.dtype, np.inexact):
|
|
kth.append(-1)
|
|
|
|
if overwrite_input:
|
|
if axis is None:
|
|
part = a.ravel()
|
|
part.partition(kth)
|
|
else:
|
|
a.partition(kth, axis=axis)
|
|
part = a
|
|
else:
|
|
part = partition(a, kth, axis=axis)
|
|
|
|
if part.shape == ():
|
|
# make 0-D arrays work
|
|
return part.item()
|
|
if axis is None:
|
|
axis = 0
|
|
|
|
indexer = [slice(None)] * part.ndim
|
|
index = part.shape[axis] // 2
|
|
if part.shape[axis] % 2 == 1:
|
|
# index with slice to allow mean (below) to work
|
|
indexer[axis] = slice(index, index+1)
|
|
else:
|
|
indexer[axis] = slice(index-1, index+1)
|
|
indexer = tuple(indexer)
|
|
|
|
# Check if the array contains any nan's
|
|
if np.issubdtype(a.dtype, np.inexact) and sz > 0:
|
|
# warn and return nans like mean would
|
|
rout = mean(part[indexer], axis=axis, out=out)
|
|
return np.lib.utils._median_nancheck(part, rout, axis, out)
|
|
else:
|
|
# if there are no nans
|
|
# Use mean in odd and even case to coerce data type
|
|
# and check, use out array.
|
|
return mean(part[indexer], axis=axis, out=out)
|
|
|
|
|
|
def _percentile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
|
|
interpolation=None, keepdims=None):
|
|
return (a, q, out)
|
|
|
|
|
|
@array_function_dispatch(_percentile_dispatcher)
|
|
def percentile(a, q, axis=None, out=None,
|
|
overwrite_input=False, interpolation='linear', keepdims=False):
|
|
"""
|
|
Compute the q-th percentile of the data along the specified axis.
|
|
|
|
Returns the q-th percentile(s) of the array elements.
|
|
|
|
Parameters
|
|
----------
|
|
a : array_like
|
|
Input array or object that can be converted to an array.
|
|
q : array_like of float
|
|
Percentile or sequence of percentiles to compute, which must be between
|
|
0 and 100 inclusive.
|
|
axis : {int, tuple of int, None}, optional
|
|
Axis or axes along which the percentiles are computed. The
|
|
default is to compute the percentile(s) along a flattened
|
|
version of the array.
|
|
|
|
.. versionchanged:: 1.9.0
|
|
A tuple of axes is supported
|
|
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 the input array `a` to be modified by intermediate
|
|
calculations, to save memory. In this case, the contents of the input
|
|
`a` after this function completes is undefined.
|
|
|
|
interpolation : {'linear', 'lower', 'higher', 'midpoint', 'nearest'}
|
|
This optional parameter specifies the interpolation method to
|
|
use when the desired percentile lies between two data points
|
|
``i < 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``.
|
|
|
|
.. 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 array `a`.
|
|
|
|
.. versionadded:: 1.9.0
|
|
|
|
Returns
|
|
-------
|
|
percentile : scalar or ndarray
|
|
If `q` is a single percentile and `axis=None`, then the result
|
|
is a scalar. If multiple percentiles are given, first axis of
|
|
the result corresponds to the percentiles. The other axes are
|
|
the axes that remain after the reduction of `a`. If the input
|
|
contains integers or floats smaller than ``float64``, the output
|
|
data-type is ``float64``. Otherwise, the output data-type is the
|
|
same as that of the input. If `out` is specified, that array is
|
|
returned instead.
|
|
|
|
See Also
|
|
--------
|
|
mean
|
|
median : equivalent to ``percentile(..., 50)``
|
|
nanpercentile
|
|
quantile : equivalent to percentile, except with q in the range [0, 1].
|
|
|
|
Notes
|
|
-----
|
|
Given a vector ``V`` of length ``N``, the q-th percentile of
|
|
``V`` is the value ``q/100`` of the way from the minimum to the
|
|
maximum 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 the location of ``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
|
|
array([[10, 7, 4],
|
|
[ 3, 2, 1]])
|
|
>>> np.percentile(a, 50)
|
|
3.5
|
|
>>> np.percentile(a, 50, axis=0)
|
|
array([6.5, 4.5, 2.5])
|
|
>>> np.percentile(a, 50, axis=1)
|
|
array([7., 2.])
|
|
>>> np.percentile(a, 50, axis=1, keepdims=True)
|
|
array([[7.],
|
|
[2.]])
|
|
|
|
>>> m = np.percentile(a, 50, axis=0)
|
|
>>> out = np.zeros_like(m)
|
|
>>> np.percentile(a, 50, axis=0, out=out)
|
|
array([6.5, 4.5, 2.5])
|
|
>>> m
|
|
array([6.5, 4.5, 2.5])
|
|
|
|
>>> b = a.copy()
|
|
>>> np.percentile(b, 50, axis=1, overwrite_input=True)
|
|
array([7., 2.])
|
|
>>> assert not np.all(a == b)
|
|
|
|
The different types of interpolation can be visualized graphically:
|
|
|
|
.. plot::
|
|
|
|
import matplotlib.pyplot as plt
|
|
|
|
a = np.arange(4)
|
|
p = np.linspace(0, 100, 6001)
|
|
ax = plt.gca()
|
|
lines = [
|
|
('linear', None),
|
|
('higher', '--'),
|
|
('lower', '--'),
|
|
('nearest', '-.'),
|
|
('midpoint', '-.'),
|
|
]
|
|
for interpolation, style in lines:
|
|
ax.plot(
|
|
p, np.percentile(a, p, interpolation=interpolation),
|
|
label=interpolation, linestyle=style)
|
|
ax.set(
|
|
title='Interpolation methods for list: ' + str(a),
|
|
xlabel='Percentile',
|
|
ylabel='List item returned',
|
|
yticks=a)
|
|
ax.legend()
|
|
plt.show()
|
|
|
|
"""
|
|
q = np.true_divide(q, 100)
|
|
q = asanyarray(q) # undo any decay that the ufunc performed (see gh-13105)
|
|
if not _quantile_is_valid(q):
|
|
raise ValueError("Percentiles must be in the range [0, 100]")
|
|
return _quantile_unchecked(
|
|
a, q, axis, out, overwrite_input, interpolation, keepdims)
|
|
|
|
|
|
def _quantile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
|
|
interpolation=None, keepdims=None):
|
|
return (a, q, out)
|
|
|
|
|
|
@array_function_dispatch(_quantile_dispatcher)
|
|
def quantile(a, q, axis=None, out=None,
|
|
overwrite_input=False, interpolation='linear', keepdims=False):
|
|
"""
|
|
Compute the q-th quantile of the data along the specified axis.
|
|
|
|
.. versionadded:: 1.15.0
|
|
|
|
Parameters
|
|
----------
|
|
a : array_like
|
|
Input array or object that can be converted to an array.
|
|
q : array_like of float
|
|
Quantile or sequence of quantiles to compute, which must be between
|
|
0 and 1 inclusive.
|
|
axis : {int, tuple of int, None}, optional
|
|
Axis or axes along which the quantiles are computed. The
|
|
default is to compute the quantile(s) along a flattened
|
|
version of the array.
|
|
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 the input array `a` to be modified by intermediate
|
|
calculations, to save memory. In this case, the contents of the input
|
|
`a` after this function completes is undefined.
|
|
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 < 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 array `a`.
|
|
|
|
Returns
|
|
-------
|
|
quantile : scalar or ndarray
|
|
If `q` is a single quantile and `axis=None`, then the result
|
|
is a scalar. If multiple quantiles are given, first axis of
|
|
the result corresponds to the quantiles. The other axes are
|
|
the axes that remain after the reduction of `a`. If the input
|
|
contains integers or floats smaller than ``float64``, the output
|
|
data-type is ``float64``. Otherwise, the output data-type is the
|
|
same as that of the input. If `out` is specified, that array is
|
|
returned instead.
|
|
|
|
See Also
|
|
--------
|
|
mean
|
|
percentile : equivalent to quantile, but with q in the range [0, 100].
|
|
median : equivalent to ``quantile(..., 0.5)``
|
|
nanquantile
|
|
|
|
Notes
|
|
-----
|
|
Given a vector ``V`` of length ``N``, the q-th quantile of
|
|
``V`` is the value ``q`` of the way from the minimum to the
|
|
maximum in a sorted copy of ``V``. The values and distances of
|
|
the two nearest neighbors as well as the `interpolation` parameter
|
|
will determine the quantile if the normalized ranking does not
|
|
match the location of ``q`` exactly. This function is the same as
|
|
the median if ``q=0.5``, the same as the minimum if ``q=0.0`` and the
|
|
same as the maximum if ``q=1.0``.
|
|
|
|
Examples
|
|
--------
|
|
>>> a = np.array([[10, 7, 4], [3, 2, 1]])
|
|
>>> a
|
|
array([[10, 7, 4],
|
|
[ 3, 2, 1]])
|
|
>>> np.quantile(a, 0.5)
|
|
3.5
|
|
>>> np.quantile(a, 0.5, axis=0)
|
|
array([6.5, 4.5, 2.5])
|
|
>>> np.quantile(a, 0.5, axis=1)
|
|
array([7., 2.])
|
|
>>> np.quantile(a, 0.5, axis=1, keepdims=True)
|
|
array([[7.],
|
|
[2.]])
|
|
>>> m = np.quantile(a, 0.5, axis=0)
|
|
>>> out = np.zeros_like(m)
|
|
>>> np.quantile(a, 0.5, axis=0, out=out)
|
|
array([6.5, 4.5, 2.5])
|
|
>>> m
|
|
array([6.5, 4.5, 2.5])
|
|
>>> b = a.copy()
|
|
>>> np.quantile(b, 0.5, axis=1, overwrite_input=True)
|
|
array([7., 2.])
|
|
>>> assert not np.all(a == b)
|
|
"""
|
|
q = np.asanyarray(q)
|
|
if not _quantile_is_valid(q):
|
|
raise ValueError("Quantiles must be in the range [0, 1]")
|
|
return _quantile_unchecked(
|
|
a, q, axis, out, overwrite_input, interpolation, keepdims)
|
|
|
|
|
|
def _quantile_unchecked(a, q, axis=None, out=None, overwrite_input=False,
|
|
interpolation='linear', keepdims=False):
|
|
"""Assumes that q is in [0, 1], and is an ndarray"""
|
|
r, k = _ureduce(a, func=_quantile_ureduce_func, q=q, axis=axis, out=out,
|
|
overwrite_input=overwrite_input,
|
|
interpolation=interpolation)
|
|
if keepdims:
|
|
return r.reshape(q.shape + k)
|
|
else:
|
|
return r
|
|
|
|
|
|
def _quantile_is_valid(q):
|
|
# avoid expensive reductions, relevant for arrays with < O(1000) elements
|
|
if q.ndim == 1 and q.size < 10:
|
|
for i in range(q.size):
|
|
if q[i] < 0.0 or q[i] > 1.0:
|
|
return False
|
|
else:
|
|
# faster than any()
|
|
if np.count_nonzero(q < 0.0) or np.count_nonzero(q > 1.0):
|
|
return False
|
|
return True
|
|
|
|
|
|
def _quantile_ureduce_func(a, q, axis=None, out=None, overwrite_input=False,
|
|
interpolation='linear', keepdims=False):
|
|
a = asarray(a)
|
|
if q.ndim == 0:
|
|
# Do not allow 0-d arrays because following code fails for scalar
|
|
zerod = True
|
|
q = q[None]
|
|
else:
|
|
zerod = False
|
|
|
|
# prepare a for partitioning
|
|
if overwrite_input:
|
|
if axis is None:
|
|
ap = a.ravel()
|
|
else:
|
|
ap = a
|
|
else:
|
|
if axis is None:
|
|
ap = a.flatten()
|
|
else:
|
|
ap = a.copy()
|
|
|
|
if axis is None:
|
|
axis = 0
|
|
|
|
Nx = ap.shape[axis]
|
|
indices = q * (Nx - 1)
|
|
|
|
# round fractional indices according to interpolation method
|
|
if interpolation == 'lower':
|
|
indices = floor(indices).astype(intp)
|
|
elif interpolation == 'higher':
|
|
indices = ceil(indices).astype(intp)
|
|
elif interpolation == 'midpoint':
|
|
indices = 0.5 * (floor(indices) + ceil(indices))
|
|
elif interpolation == 'nearest':
|
|
indices = around(indices).astype(intp)
|
|
elif interpolation == 'linear':
|
|
pass # keep index as fraction and interpolate
|
|
else:
|
|
raise ValueError(
|
|
"interpolation can only be 'linear', 'lower' 'higher', "
|
|
"'midpoint', or 'nearest'")
|
|
|
|
n = np.array(False, dtype=bool) # check for nan's flag
|
|
if indices.dtype == intp: # take the points along axis
|
|
# Check if the array contains any nan's
|
|
if np.issubdtype(a.dtype, np.inexact):
|
|
indices = concatenate((indices, [-1]))
|
|
|
|
ap.partition(indices, axis=axis)
|
|
# ensure axis with q-th is first
|
|
ap = np.moveaxis(ap, axis, 0)
|
|
axis = 0
|
|
|
|
# Check if the array contains any nan's
|
|
if np.issubdtype(a.dtype, np.inexact):
|
|
indices = indices[:-1]
|
|
n = np.isnan(ap[-1:, ...])
|
|
|
|
if zerod:
|
|
indices = indices[0]
|
|
r = take(ap, indices, axis=axis, out=out)
|
|
|
|
|
|
else: # weight the points above and below the indices
|
|
indices_below = floor(indices).astype(intp)
|
|
indices_above = indices_below + 1
|
|
indices_above[indices_above > Nx - 1] = Nx - 1
|
|
|
|
# Check if the array contains any nan's
|
|
if np.issubdtype(a.dtype, np.inexact):
|
|
indices_above = concatenate((indices_above, [-1]))
|
|
|
|
weights_above = indices - indices_below
|
|
weights_below = 1 - weights_above
|
|
|
|
weights_shape = [1, ] * ap.ndim
|
|
weights_shape[axis] = len(indices)
|
|
weights_below.shape = weights_shape
|
|
weights_above.shape = weights_shape
|
|
|
|
ap.partition(concatenate((indices_below, indices_above)), axis=axis)
|
|
|
|
# ensure axis with q-th is first
|
|
ap = np.moveaxis(ap, axis, 0)
|
|
weights_below = np.moveaxis(weights_below, axis, 0)
|
|
weights_above = np.moveaxis(weights_above, axis, 0)
|
|
axis = 0
|
|
|
|
# Check if the array contains any nan's
|
|
if np.issubdtype(a.dtype, np.inexact):
|
|
indices_above = indices_above[:-1]
|
|
n = np.isnan(ap[-1:, ...])
|
|
|
|
x1 = take(ap, indices_below, axis=axis) * weights_below
|
|
x2 = take(ap, indices_above, axis=axis) * weights_above
|
|
|
|
# ensure axis with q-th is first
|
|
x1 = np.moveaxis(x1, axis, 0)
|
|
x2 = np.moveaxis(x2, axis, 0)
|
|
|
|
if zerod:
|
|
x1 = x1.squeeze(0)
|
|
x2 = x2.squeeze(0)
|
|
|
|
if out is not None:
|
|
r = add(x1, x2, out=out)
|
|
else:
|
|
r = add(x1, x2)
|
|
|
|
if np.any(n):
|
|
if zerod:
|
|
if ap.ndim == 1:
|
|
if out is not None:
|
|
out[...] = a.dtype.type(np.nan)
|
|
r = out
|
|
else:
|
|
r = a.dtype.type(np.nan)
|
|
else:
|
|
r[..., n.squeeze(0)] = a.dtype.type(np.nan)
|
|
else:
|
|
if r.ndim == 1:
|
|
r[:] = a.dtype.type(np.nan)
|
|
else:
|
|
r[..., n.repeat(q.size, 0)] = a.dtype.type(np.nan)
|
|
|
|
return r
|
|
|
|
|
|
def _trapz_dispatcher(y, x=None, dx=None, axis=None):
|
|
return (y, x)
|
|
|
|
|
|
@array_function_dispatch(_trapz_dispatcher)
|
|
def trapz(y, x=None, dx=1.0, axis=-1):
|
|
"""
|
|
Integrate along the given axis using the composite trapezoidal rule.
|
|
|
|
Integrate `y` (`x`) along given axis.
|
|
|
|
Parameters
|
|
----------
|
|
y : array_like
|
|
Input array to integrate.
|
|
x : array_like, optional
|
|
The sample points corresponding to the `y` values. If `x` is None,
|
|
the sample points are assumed to be evenly spaced `dx` apart. The
|
|
default is None.
|
|
dx : scalar, optional
|
|
The spacing between sample points when `x` is None. The default is 1.
|
|
axis : int, optional
|
|
The axis along which to integrate.
|
|
|
|
Returns
|
|
-------
|
|
trapz : float
|
|
Definite integral as approximated by trapezoidal rule.
|
|
|
|
See Also
|
|
--------
|
|
sum, cumsum
|
|
|
|
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.
|
|
|
|
|
|
References
|
|
----------
|
|
.. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule
|
|
|
|
.. [2] Illustration image:
|
|
https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
|
|
|
|
Examples
|
|
--------
|
|
>>> np.trapz([1,2,3])
|
|
4.0
|
|
>>> np.trapz([1,2,3], x=[4,6,8])
|
|
8.0
|
|
>>> np.trapz([1,2,3], dx=2)
|
|
8.0
|
|
>>> a = np.arange(6).reshape(2, 3)
|
|
>>> a
|
|
array([[0, 1, 2],
|
|
[3, 4, 5]])
|
|
>>> np.trapz(a, axis=0)
|
|
array([1.5, 2.5, 3.5])
|
|
>>> np.trapz(a, axis=1)
|
|
array([2., 8.])
|
|
|
|
"""
|
|
y = asanyarray(y)
|
|
if x is None:
|
|
d = dx
|
|
else:
|
|
x = asanyarray(x)
|
|
if x.ndim == 1:
|
|
d = diff(x)
|
|
# reshape to correct shape
|
|
shape = [1]*y.ndim
|
|
shape[axis] = d.shape[0]
|
|
d = d.reshape(shape)
|
|
else:
|
|
d = diff(x, axis=axis)
|
|
nd = y.ndim
|
|
slice1 = [slice(None)]*nd
|
|
slice2 = [slice(None)]*nd
|
|
slice1[axis] = slice(1, None)
|
|
slice2[axis] = slice(None, -1)
|
|
try:
|
|
ret = (d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0).sum(axis)
|
|
except ValueError:
|
|
# Operations didn't work, cast to ndarray
|
|
d = np.asarray(d)
|
|
y = np.asarray(y)
|
|
ret = add.reduce(d * (y[tuple(slice1)]+y[tuple(slice2)])/2.0, axis)
|
|
return ret
|
|
|
|
|
|
def _meshgrid_dispatcher(*xi, **kwargs):
|
|
return xi
|
|
|
|
|
|
# Based on scitools meshgrid
|
|
@array_function_dispatch(_meshgrid_dispatcher)
|
|
def meshgrid(*xi, **kwargs):
|
|
"""
|
|
Return coordinate matrices from coordinate vectors.
|
|
|
|
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.
|
|
|
|
Parameters
|
|
----------
|
|
x1, x2,..., xn : array_like
|
|
1-D arrays representing the coordinates of a grid.
|
|
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.
|
|
|
|
.. 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.
|
|
|
|
.. versionadded:: 1.7.0
|
|
|
|
Returns
|
|
-------
|
|
X1, X2,..., XN : ndarray
|
|
For vectors `x1`, `x2`,..., 'xn' with lengths ``Ni=len(xi)`` ,
|
|
return ``(N1, N2, N3,...Nn)`` shaped arrays if indexing='ij'
|
|
or ``(N2, N1, N3,...Nn)`` shaped arrays if indexing='xy'
|
|
with the elements of `xi` repeated to fill the matrix along
|
|
the first dimension for `x1`, the second for `x2` and so on.
|
|
|
|
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::
|
|
|
|
xv, yv = np.meshgrid(x, y, sparse=False, indexing='ij')
|
|
for i in range(nx):
|
|
for j in range(ny):
|
|
# treat xv[i,j], yv[i,j]
|
|
|
|
xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy')
|
|
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
|
|
--------
|
|
index_tricks.mgrid : Construct a multi-dimensional "meshgrid"
|
|
using indexing notation.
|
|
index_tricks.ogrid : Construct an open multi-dimensional "meshgrid"
|
|
using indexing notation.
|
|
|
|
Examples
|
|
--------
|
|
>>> nx, ny = (3, 2)
|
|
>>> x = np.linspace(0, 1, nx)
|
|
>>> y = np.linspace(0, 1, ny)
|
|
>>> xv, yv = np.meshgrid(x, y)
|
|
>>> xv
|
|
array([[0. , 0.5, 1. ],
|
|
[0. , 0.5, 1. ]])
|
|
>>> yv
|
|
array([[0., 0., 0.],
|
|
[1., 1., 1.]])
|
|
>>> xv, yv = np.meshgrid(x, y, sparse=True) # make sparse output arrays
|
|
>>> xv
|
|
array([[0. , 0.5, 1. ]])
|
|
>>> yv
|
|
array([[0.],
|
|
[1.]])
|
|
|
|
`meshgrid` is very useful to evaluate functions on a grid.
|
|
|
|
>>> import matplotlib.pyplot as plt
|
|
>>> x = np.arange(-5, 5, 0.1)
|
|
>>> y = np.arange(-5, 5, 0.1)
|
|
>>> xx, yy = np.meshgrid(x, y, sparse=True)
|
|
>>> z = np.sin(xx**2 + yy**2) / (xx**2 + yy**2)
|
|
>>> h = plt.contourf(x,y,z)
|
|
>>> plt.show()
|
|
|
|
"""
|
|
ndim = len(xi)
|
|
|
|
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 indexing not in ['xy', 'ij']:
|
|
raise ValueError(
|
|
"Valid values for `indexing` are 'xy' and 'ij'.")
|
|
|
|
s0 = (1,) * ndim
|
|
output = [np.asanyarray(x).reshape(s0[:i] + (-1,) + s0[i + 1:])
|
|
for i, x in enumerate(xi)]
|
|
|
|
if indexing == 'xy' and ndim > 1:
|
|
# switch first and second axis
|
|
output[0].shape = (1, -1) + s0[2:]
|
|
output[1].shape = (-1, 1) + s0[2:]
|
|
|
|
if not sparse:
|
|
# Return the full N-D matrix (not only the 1-D vector)
|
|
output = np.broadcast_arrays(*output, subok=True)
|
|
|
|
if copy_:
|
|
output = [x.copy() for x in output]
|
|
|
|
return output
|
|
|
|
|
|
def _delete_dispatcher(arr, obj, axis=None):
|
|
return (arr, obj)
|
|
|
|
|
|
@array_function_dispatch(_delete_dispatcher)
|
|
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]`.
|
|
|
|
Parameters
|
|
----------
|
|
arr : array_like
|
|
Input array.
|
|
obj : slice, int or array of ints
|
|
Indicate indices of sub-arrays to remove along the specified axis.
|
|
axis : int, optional
|
|
The axis along which to delete the subarray defined by `obj`.
|
|
If `axis` is None, `obj` is applied to the flattened array.
|
|
|
|
Returns
|
|
-------
|
|
out : ndarray
|
|
A copy of `arr` with the elements specified by `obj` removed. Note
|
|
that `delete` does not occur in-place. If `axis` is None, `out` is
|
|
a flattened array.
|
|
|
|
See Also
|
|
--------
|
|
insert : Insert elements into an array.
|
|
append : Append elements at the end of an array.
|
|
|
|
Notes
|
|
-----
|
|
Often it is preferable to use a boolean mask. For example:
|
|
|
|
>>> arr = np.arange(12) + 1
|
|
>>> 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`.
|
|
|
|
Examples
|
|
--------
|
|
>>> arr = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]])
|
|
>>> arr
|
|
array([[ 1, 2, 3, 4],
|
|
[ 5, 6, 7, 8],
|
|
[ 9, 10, 11, 12]])
|
|
>>> np.delete(arr, 1, 0)
|
|
array([[ 1, 2, 3, 4],
|
|
[ 9, 10, 11, 12]])
|
|
|
|
>>> np.delete(arr, np.s_[::2], 1)
|
|
array([[ 2, 4],
|
|
[ 6, 8],
|
|
[10, 12]])
|
|
>>> np.delete(arr, [1,3,5], None)
|
|
array([ 1, 3, 5, 7, 8, 9, 10, 11, 12])
|
|
|
|
"""
|
|
wrap = None
|
|
if type(arr) is not ndarray:
|
|
try:
|
|
wrap = arr.__array_wrap__
|
|
except AttributeError:
|
|
pass
|
|
|
|
arr = asarray(arr)
|
|
ndim = arr.ndim
|
|
arrorder = 'F' if arr.flags.fnc else 'C'
|
|
if axis is None:
|
|
if ndim != 1:
|
|
arr = arr.ravel()
|
|
ndim = arr.ndim
|
|
axis = -1
|
|
|
|
if ndim == 0:
|
|
# 2013-09-24, 1.9
|
|
warnings.warn(
|
|
"in the future the special handling of scalars will be removed "
|
|
"from delete and raise an error", DeprecationWarning, stacklevel=3)
|
|
if wrap:
|
|
return wrap(arr)
|
|
else:
|
|
return arr.copy(order=arrorder)
|
|
|
|
axis = normalize_axis_index(axis, ndim)
|
|
|
|
slobj = [slice(None)]*ndim
|
|
N = arr.shape[axis]
|
|
newshape = list(arr.shape)
|
|
|
|
if isinstance(obj, slice):
|
|
start, stop, step = obj.indices(N)
|
|
xr = range(start, stop, step)
|
|
numtodel = len(xr)
|
|
|
|
if numtodel <= 0:
|
|
if wrap:
|
|
return wrap(arr.copy(order=arrorder))
|
|
else:
|
|
return arr.copy(order=arrorder)
|
|
|
|
# Invert if step is negative:
|
|
if step < 0:
|
|
step = -step
|
|
start = xr[-1]
|
|
stop = xr[0] + 1
|
|
|
|
newshape[axis] -= numtodel
|
|
new = empty(newshape, arr.dtype, arrorder)
|
|
# copy initial chunk
|
|
if start == 0:
|
|
pass
|
|
else:
|
|
slobj[axis] = slice(None, start)
|
|
new[tuple(slobj)] = arr[tuple(slobj)]
|
|
# copy end chunk
|
|
if stop == N:
|
|
pass
|
|
else:
|
|
slobj[axis] = slice(stop-numtodel, None)
|
|
slobj2 = [slice(None)]*ndim
|
|
slobj2[axis] = slice(stop, None)
|
|
new[tuple(slobj)] = arr[tuple(slobj2)]
|
|
# copy middle pieces
|
|
if step == 1:
|
|
pass
|
|
else: # use array indexing.
|
|
keep = ones(stop-start, dtype=bool)
|
|
keep[:stop-start:step] = False
|
|
slobj[axis] = slice(start, stop-numtodel)
|
|
slobj2 = [slice(None)]*ndim
|
|
slobj2[axis] = slice(start, stop)
|
|
arr = arr[tuple(slobj2)]
|
|
slobj2[axis] = keep
|
|
new[tuple(slobj)] = arr[tuple(slobj2)]
|
|
if wrap:
|
|
return wrap(new)
|
|
else:
|
|
return new
|
|
|
|
_obj = obj
|
|
obj = np.asarray(obj)
|
|
# After removing the special handling of booleans and out of
|
|
# bounds values, the conversion to the array can be removed.
|
|
if obj.dtype == bool:
|
|
warnings.warn("in the future insert will treat boolean arrays and "
|
|
"array-likes as boolean index instead of casting it "
|
|
"to integer", FutureWarning, stacklevel=3)
|
|
obj = obj.astype(intp)
|
|
if isinstance(_obj, (int, long, integer)):
|
|
# optimization for a single value
|
|
obj = obj.item()
|
|
if (obj < -N or obj >= N):
|
|
raise IndexError(
|
|
"index %i is out of bounds for axis %i with "
|
|
"size %i" % (obj, axis, N))
|
|
if (obj < 0):
|
|
obj += N
|
|
newshape[axis] -= 1
|
|
new = empty(newshape, arr.dtype, arrorder)
|
|
slobj[axis] = slice(None, obj)
|
|
new[tuple(slobj)] = arr[tuple(slobj)]
|
|
slobj[axis] = slice(obj, None)
|
|
slobj2 = [slice(None)]*ndim
|
|
slobj2[axis] = slice(obj+1, None)
|
|
new[tuple(slobj)] = arr[tuple(slobj2)]
|
|
else:
|
|
if obj.size == 0 and not isinstance(_obj, np.ndarray):
|
|
obj = obj.astype(intp)
|
|
if not np.can_cast(obj, intp, 'same_kind'):
|
|
# obj.size = 1 special case always failed and would just
|
|
# give superfluous warnings.
|
|
# 2013-09-24, 1.9
|
|
warnings.warn(
|
|
"using a non-integer array as obj in delete will result in an "
|
|
"error in the future", DeprecationWarning, stacklevel=3)
|
|
obj = obj.astype(intp)
|
|
keep = ones(N, dtype=bool)
|
|
|
|
# Test if there are out of bound indices, this is deprecated
|
|
inside_bounds = (obj < N) & (obj >= -N)
|
|
if not inside_bounds.all():
|
|
# 2013-09-24, 1.9
|
|
warnings.warn(
|
|
"in the future out of bounds indices will raise an error "
|
|
"instead of being ignored by `numpy.delete`.",
|
|
DeprecationWarning, stacklevel=3)
|
|
obj = obj[inside_bounds]
|
|
positive_indices = obj >= 0
|
|
if not positive_indices.all():
|
|
warnings.warn(
|
|
"in the future negative indices will not be ignored by "
|
|
"`numpy.delete`.", FutureWarning, stacklevel=3)
|
|
obj = obj[positive_indices]
|
|
|
|
keep[obj, ] = False
|
|
slobj[axis] = keep
|
|
new = arr[tuple(slobj)]
|
|
|
|
if wrap:
|
|
return wrap(new)
|
|
else:
|
|
return new
|
|
|
|
|
|
def _insert_dispatcher(arr, obj, values, axis=None):
|
|
return (arr, obj, values)
|
|
|
|
|
|
@array_function_dispatch(_insert_dispatcher)
|
|
def insert(arr, obj, values, axis=None):
|
|
"""
|
|
Insert values along the given axis before the given indices.
|
|
|
|
Parameters
|
|
----------
|
|
arr : array_like
|
|
Input array.
|
|
obj : int, slice or sequence of ints
|
|
Object that defines the index or indices before which `values` is
|
|
inserted.
|
|
|
|
.. 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).
|
|
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`.
|
|
`values` should be shaped so that ``arr[...,obj,...] = values``
|
|
is legal.
|
|
axis : int, optional
|
|
Axis along which to insert `values`. If `axis` is None then `arr`
|
|
is flattened first.
|
|
|
|
Returns
|
|
-------
|
|
out : ndarray
|
|
A copy of `arr` with `values` inserted. Note that `insert`
|
|
does not occur in-place: a new array is returned. If
|
|
`axis` is None, `out` is a flattened array.
|
|
|
|
See Also
|
|
--------
|
|
append : Append elements at the end of an array.
|
|
concatenate : Join a sequence of arrays along an existing axis.
|
|
delete : Delete elements from an array.
|
|
|
|
Notes
|
|
-----
|
|
Note that for higher dimensional inserts `obj=0` behaves very different
|
|
from `obj=[0]` just like `arr[:,0,:] = values` is different from
|
|
`arr[:,[0],:] = values`.
|
|
|
|
Examples
|
|
--------
|
|
>>> a = np.array([[1, 1], [2, 2], [3, 3]])
|
|
>>> a
|
|
array([[1, 1],
|
|
[2, 2],
|
|
[3, 3]])
|
|
>>> np.insert(a, 1, 5)
|
|
array([1, 5, 1, ..., 2, 3, 3])
|
|
>>> np.insert(a, 1, 5, axis=1)
|
|
array([[1, 5, 1],
|
|
[2, 5, 2],
|
|
[3, 5, 3]])
|
|
|
|
Difference between sequence and scalars:
|
|
|
|
>>> np.insert(a, [1], [[1],[2],[3]], axis=1)
|
|
array([[1, 1, 1],
|
|
[2, 2, 2],
|
|
[3, 3, 3]])
|
|
>>> np.array_equal(np.insert(a, 1, [1, 2, 3], axis=1),
|
|
... np.insert(a, [1], [[1],[2],[3]], axis=1))
|
|
True
|
|
|
|
>>> b = a.flatten()
|
|
>>> b
|
|
array([1, 1, 2, 2, 3, 3])
|
|
>>> np.insert(b, [2, 2], [5, 6])
|
|
array([1, 1, 5, ..., 2, 3, 3])
|
|
|
|
>>> np.insert(b, slice(2, 4), [5, 6])
|
|
array([1, 1, 5, ..., 2, 3, 3])
|
|
|
|
>>> np.insert(b, [2, 2], [7.13, False]) # type casting
|
|
array([1, 1, 7, ..., 2, 3, 3])
|
|
|
|
>>> x = np.arange(8).reshape(2, 4)
|
|
>>> idx = (1, 3)
|
|
>>> np.insert(x, idx, 999, axis=1)
|
|
array([[ 0, 999, 1, 2, 999, 3],
|
|
[ 4, 999, 5, 6, 999, 7]])
|
|
|
|
"""
|
|
wrap = None
|
|
if type(arr) is not ndarray:
|
|
try:
|
|
wrap = arr.__array_wrap__
|
|
except AttributeError:
|
|
pass
|
|
|
|
arr = asarray(arr)
|
|
ndim = arr.ndim
|
|
arrorder = 'F' if arr.flags.fnc else 'C'
|
|
if axis is None:
|
|
if ndim != 1:
|
|
arr = arr.ravel()
|
|
ndim = arr.ndim
|
|
axis = ndim - 1
|
|
elif ndim == 0:
|
|
# 2013-09-24, 1.9
|
|
warnings.warn(
|
|
"in the future the special handling of scalars will be removed "
|
|
"from insert and raise an error", DeprecationWarning, stacklevel=3)
|
|
arr = arr.copy(order=arrorder)
|
|
arr[...] = values
|
|
if wrap:
|
|
return wrap(arr)
|
|
else:
|
|
return arr
|
|
else:
|
|
axis = normalize_axis_index(axis, ndim)
|
|
slobj = [slice(None)]*ndim
|
|
N = arr.shape[axis]
|
|
newshape = list(arr.shape)
|
|
|
|
if isinstance(obj, slice):
|
|
# turn it into a range object
|
|
indices = arange(*obj.indices(N), **{'dtype': intp})
|
|
else:
|
|
# need to copy obj, because indices will be changed in-place
|
|
indices = np.array(obj)
|
|
if indices.dtype == bool:
|
|
# See also delete
|
|
warnings.warn(
|
|
"in the future insert will treat boolean arrays and "
|
|
"array-likes as a boolean index instead of casting it to "
|
|
"integer", FutureWarning, stacklevel=3)
|
|
indices = indices.astype(intp)
|
|
# Code after warning period:
|
|
#if obj.ndim != 1:
|
|
# raise ValueError('boolean array argument obj to insert '
|
|
# 'must be one dimensional')
|
|
#indices = np.flatnonzero(obj)
|
|
elif indices.ndim > 1:
|
|
raise ValueError(
|
|
"index array argument obj to insert must be one dimensional "
|
|
"or scalar")
|
|
if indices.size == 1:
|
|
index = indices.item()
|
|
if index < -N or index > N:
|
|
raise IndexError(
|
|
"index %i is out of bounds for axis %i with "
|
|
"size %i" % (obj, axis, N))
|
|
if (index < 0):
|
|
index += N
|
|
|
|
# There are some object array corner cases here, but we cannot avoid
|
|
# that:
|
|
values = array(values, copy=False, ndmin=arr.ndim, dtype=arr.dtype)
|
|
if indices.ndim == 0:
|
|
# broadcasting is very different here, since a[:,0,:] = ... behaves
|
|
# very different from a[:,[0],:] = ...! This changes values so that
|
|
# it works likes the second case. (here a[:,0:1,:])
|
|
values = np.moveaxis(values, 0, axis)
|
|
numnew = values.shape[axis]
|
|
newshape[axis] += numnew
|
|
new = empty(newshape, arr.dtype, arrorder)
|
|
slobj[axis] = slice(None, index)
|
|
new[tuple(slobj)] = arr[tuple(slobj)]
|
|
slobj[axis] = slice(index, index+numnew)
|
|
new[tuple(slobj)] = values
|
|
slobj[axis] = slice(index+numnew, None)
|
|
slobj2 = [slice(None)] * ndim
|
|
slobj2[axis] = slice(index, None)
|
|
new[tuple(slobj)] = arr[tuple(slobj2)]
|
|
if wrap:
|
|
return wrap(new)
|
|
return new
|
|
elif indices.size == 0 and not isinstance(obj, np.ndarray):
|
|
# Can safely cast the empty list to intp
|
|
indices = indices.astype(intp)
|
|
|
|
if not np.can_cast(indices, intp, 'same_kind'):
|
|
# 2013-09-24, 1.9
|
|
warnings.warn(
|
|
"using a non-integer array as obj in insert will result in an "
|
|
"error in the future", DeprecationWarning, stacklevel=3)
|
|
indices = indices.astype(intp)
|
|
|
|
indices[indices < 0] += N
|
|
|
|
numnew = len(indices)
|
|
order = indices.argsort(kind='mergesort') # stable sort
|
|
indices[order] += np.arange(numnew)
|
|
|
|
newshape[axis] += numnew
|
|
old_mask = ones(newshape[axis], dtype=bool)
|
|
old_mask[indices] = False
|
|
|
|
new = empty(newshape, arr.dtype, arrorder)
|
|
slobj2 = [slice(None)]*ndim
|
|
slobj[axis] = indices
|
|
slobj2[axis] = old_mask
|
|
new[tuple(slobj)] = values
|
|
new[tuple(slobj2)] = arr
|
|
|
|
if wrap:
|
|
return wrap(new)
|
|
return new
|
|
|
|
|
|
def _append_dispatcher(arr, values, axis=None):
|
|
return (arr, values)
|
|
|
|
|
|
@array_function_dispatch(_append_dispatcher)
|
|
def append(arr, values, axis=None):
|
|
"""
|
|
Append values to the end of an array.
|
|
|
|
Parameters
|
|
----------
|
|
arr : array_like
|
|
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.
|
|
axis : int, optional
|
|
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.
|
|
|
|
See Also
|
|
--------
|
|
insert : Insert elements into an array.
|
|
delete : Delete elements from an array.
|
|
|
|
Examples
|
|
--------
|
|
>>> np.append([1, 2, 3], [[4, 5, 6], [7, 8, 9]])
|
|
array([1, 2, 3, ..., 7, 8, 9])
|
|
|
|
When `axis` is specified, `values` must have the correct shape.
|
|
|
|
>>> np.append([[1, 2, 3], [4, 5, 6]], [[7, 8, 9]], axis=0)
|
|
array([[1, 2, 3],
|
|
[4, 5, 6],
|
|
[7, 8, 9]])
|
|
>>> np.append([[1, 2, 3], [4, 5, 6]], [7, 8, 9], axis=0)
|
|
Traceback (most recent call last):
|
|
...
|
|
ValueError: all the input arrays must have same number of dimensions
|
|
|
|
"""
|
|
arr = asanyarray(arr)
|
|
if axis is None:
|
|
if arr.ndim != 1:
|
|
arr = arr.ravel()
|
|
values = ravel(values)
|
|
axis = arr.ndim-1
|
|
return concatenate((arr, values), axis=axis)
|
|
|
|
|
|
def _digitize_dispatcher(x, bins, right=None):
|
|
return (x, bins)
|
|
|
|
|
|
@array_function_dispatch(_digitize_dispatcher)
|
|
def digitize(x, bins, right=False):
|
|
"""
|
|
Return the indices of the bins to which each value in input array belongs.
|
|
|
|
========= ============= ============================
|
|
`right` order of bins returned index `i` satisfies
|
|
========= ============= ============================
|
|
``False`` increasing ``bins[i-1] <= x < bins[i]``
|
|
``True`` increasing ``bins[i-1] < x <= bins[i]``
|
|
``False`` decreasing ``bins[i-1] > x >= bins[i]``
|
|
``True`` decreasing ``bins[i-1] >= x > bins[i]``
|
|
========= ============= ============================
|
|
|
|
If values in `x` are beyond the bounds of `bins`, 0 or ``len(bins)`` is
|
|
returned as appropriate.
|
|
|
|
Parameters
|
|
----------
|
|
x : array_like
|
|
Input array to be binned. Prior to NumPy 1.10.0, this array had to
|
|
be 1-dimensional, but can now have any shape.
|
|
bins : array_like
|
|
Array of bins. It has to be 1-dimensional and monotonic.
|
|
right : bool, optional
|
|
Indicating whether the intervals include the right or the left bin
|
|
edge. Default behavior is (right==False) indicating that the interval
|
|
does not include the right edge. The left bin end is open in this
|
|
case, i.e., bins[i-1] <= x < bins[i] is the default behavior for
|
|
monotonically increasing bins.
|
|
|
|
Returns
|
|
-------
|
|
indices : ndarray of ints
|
|
Output array of indices, of same shape as `x`.
|
|
|
|
Raises
|
|
------
|
|
ValueError
|
|
If `bins` is not monotonic.
|
|
TypeError
|
|
If the type of the input is complex.
|
|
|
|
See Also
|
|
--------
|
|
bincount, histogram, unique, searchsorted
|
|
|
|
Notes
|
|
-----
|
|
If values in `x` are such that they fall outside the bin range,
|
|
attempting to index `bins` with the indices that `digitize` returns
|
|
will result in an IndexError.
|
|
|
|
.. versionadded:: 1.10.0
|
|
|
|
`np.digitize` is implemented in terms of `np.searchsorted`. This means
|
|
that a binary search is used to bin the values, which scales much better
|
|
for larger number of bins than the previous linear search. It also removes
|
|
the requirement for the input array to be 1-dimensional.
|
|
|
|
For monotonically _increasing_ `bins`, the following are equivalent::
|
|
|
|
np.digitize(x, bins, right=True)
|
|
np.searchsorted(bins, x, side='left')
|
|
|
|
Note that as the order of the arguments are reversed, the side must be too.
|
|
The `searchsorted` call is marginally faster, as it does not do any
|
|
monotonicity checks. Perhaps more importantly, it supports all dtypes.
|
|
|
|
Examples
|
|
--------
|
|
>>> x = np.array([0.2, 6.4, 3.0, 1.6])
|
|
>>> bins = np.array([0.0, 1.0, 2.5, 4.0, 10.0])
|
|
>>> inds = np.digitize(x, bins)
|
|
>>> inds
|
|
array([1, 4, 3, 2])
|
|
>>> for n in range(x.size):
|
|
... print(bins[inds[n]-1], "<=", x[n], "<", bins[inds[n]])
|
|
...
|
|
0.0 <= 0.2 < 1.0
|
|
4.0 <= 6.4 < 10.0
|
|
2.5 <= 3.0 < 4.0
|
|
1.0 <= 1.6 < 2.5
|
|
|
|
>>> x = np.array([1.2, 10.0, 12.4, 15.5, 20.])
|
|
>>> bins = np.array([0, 5, 10, 15, 20])
|
|
>>> np.digitize(x,bins,right=True)
|
|
array([1, 2, 3, 4, 4])
|
|
>>> np.digitize(x,bins,right=False)
|
|
array([1, 3, 3, 4, 5])
|
|
"""
|
|
x = _nx.asarray(x)
|
|
bins = _nx.asarray(bins)
|
|
|
|
# here for compatibility, searchsorted below is happy to take this
|
|
if np.issubdtype(x.dtype, _nx.complexfloating):
|
|
raise TypeError("x may not be complex")
|
|
|
|
mono = _monotonicity(bins)
|
|
if mono == 0:
|
|
raise ValueError("bins must be monotonically increasing or decreasing")
|
|
|
|
# this is backwards because the arguments below are swapped
|
|
side = 'left' if right else 'right'
|
|
if mono == -1:
|
|
# reverse the bins, and invert the results
|
|
return len(bins) - _nx.searchsorted(bins[::-1], x, side=side)
|
|
else:
|
|
return _nx.searchsorted(bins, x, side=side)
|