Chris@87: from __future__ import division, absolute_import, print_function Chris@87: Chris@87: import warnings Chris@87: import sys Chris@87: import collections Chris@87: import operator Chris@87: Chris@87: import numpy as np Chris@87: import numpy.core.numeric as _nx Chris@87: from numpy.core import linspace, atleast_1d, atleast_2d Chris@87: from numpy.core.numeric import ( Chris@87: ones, zeros, arange, concatenate, array, asarray, asanyarray, empty, Chris@87: empty_like, ndarray, around, floor, ceil, take, dot, where, intp, Chris@87: integer, isscalar Chris@87: ) Chris@87: from numpy.core.umath import ( Chris@87: pi, multiply, add, arctan2, frompyfunc, cos, less_equal, sqrt, sin, Chris@87: mod, exp, log10 Chris@87: ) Chris@87: from numpy.core.fromnumeric import ( Chris@87: ravel, nonzero, sort, partition, mean Chris@87: ) Chris@87: from numpy.core.numerictypes import typecodes, number Chris@87: from numpy.lib.twodim_base import diag Chris@87: from .utils import deprecate Chris@87: from ._compiled_base import _insert, add_docstring Chris@87: from ._compiled_base import digitize, bincount, interp as compiled_interp Chris@87: from ._compiled_base import add_newdoc_ufunc Chris@87: from numpy.compat import long Chris@87: Chris@87: # Force range to be a generator, for np.delete's usage. Chris@87: if sys.version_info[0] < 3: Chris@87: range = xrange Chris@87: Chris@87: Chris@87: __all__ = [ Chris@87: 'select', 'piecewise', 'trim_zeros', 'copy', 'iterable', 'percentile', Chris@87: 'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'disp', Chris@87: 'extract', 'place', 'vectorize', 'asarray_chkfinite', 'average', Chris@87: 'histogram', 'histogramdd', 'bincount', 'digitize', 'cov', 'corrcoef', Chris@87: 'msort', 'median', 'sinc', 'hamming', 'hanning', 'bartlett', Chris@87: 'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc', 'add_docstring', Chris@87: 'meshgrid', 'delete', 'insert', 'append', 'interp', 'add_newdoc_ufunc' Chris@87: ] Chris@87: Chris@87: Chris@87: def iterable(y): Chris@87: """ Chris@87: Check whether or not an object can be iterated over. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: y : object Chris@87: Input object. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: b : {0, 1} Chris@87: Return 1 if the object has an iterator method or is a sequence, Chris@87: and 0 otherwise. Chris@87: Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> np.iterable([1, 2, 3]) Chris@87: 1 Chris@87: >>> np.iterable(2) Chris@87: 0 Chris@87: Chris@87: """ Chris@87: try: Chris@87: iter(y) Chris@87: except: Chris@87: return 0 Chris@87: return 1 Chris@87: Chris@87: Chris@87: def histogram(a, bins=10, range=None, normed=False, weights=None, Chris@87: density=None): Chris@87: """ Chris@87: Compute the histogram of a set of data. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: a : array_like Chris@87: Input data. The histogram is computed over the flattened array. Chris@87: bins : int or sequence of scalars, optional Chris@87: If `bins` is an int, it defines the number of equal-width Chris@87: bins in the given range (10, by default). If `bins` is a sequence, Chris@87: it defines the bin edges, including the rightmost edge, allowing Chris@87: for non-uniform bin widths. Chris@87: range : (float, float), optional Chris@87: The lower and upper range of the bins. If not provided, range Chris@87: is simply ``(a.min(), a.max())``. Values outside the range are Chris@87: ignored. Chris@87: normed : bool, optional Chris@87: This keyword is deprecated in Numpy 1.6 due to confusing/buggy Chris@87: behavior. It will be removed in Numpy 2.0. Use the density keyword Chris@87: instead. Chris@87: If False, the result will contain the number of samples Chris@87: in each bin. If True, the result is the value of the Chris@87: probability *density* function at the bin, normalized such that Chris@87: the *integral* over the range is 1. Note that this latter behavior is Chris@87: known to be buggy with unequal bin widths; use `density` instead. Chris@87: weights : array_like, optional Chris@87: An array of weights, of the same shape as `a`. Each value in `a` Chris@87: only contributes its associated weight towards the bin count Chris@87: (instead of 1). If `normed` is True, the weights are normalized, Chris@87: so that the integral of the density over the range remains 1 Chris@87: density : bool, optional Chris@87: If False, the result will contain the number of samples Chris@87: in each bin. If True, the result is the value of the Chris@87: probability *density* function at the bin, normalized such that Chris@87: the *integral* over the range is 1. Note that the sum of the Chris@87: histogram values will not be equal to 1 unless bins of unity Chris@87: width are chosen; it is not a probability *mass* function. Chris@87: Overrides the `normed` keyword if given. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: hist : array Chris@87: The values of the histogram. See `normed` and `weights` for a Chris@87: description of the possible semantics. Chris@87: bin_edges : array of dtype float Chris@87: Return the bin edges ``(length(hist)+1)``. Chris@87: Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: histogramdd, bincount, searchsorted, digitize Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: All but the last (righthand-most) bin is half-open. In other words, if Chris@87: `bins` is:: Chris@87: Chris@87: [1, 2, 3, 4] Chris@87: Chris@87: then the first bin is ``[1, 2)`` (including 1, but excluding 2) and the Chris@87: second ``[2, 3)``. The last bin, however, is ``[3, 4]``, which *includes* Chris@87: 4. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> np.histogram([1, 2, 1], bins=[0, 1, 2, 3]) Chris@87: (array([0, 2, 1]), array([0, 1, 2, 3])) Chris@87: >>> np.histogram(np.arange(4), bins=np.arange(5), density=True) Chris@87: (array([ 0.25, 0.25, 0.25, 0.25]), array([0, 1, 2, 3, 4])) Chris@87: >>> np.histogram([[1, 2, 1], [1, 0, 1]], bins=[0,1,2,3]) Chris@87: (array([1, 4, 1]), array([0, 1, 2, 3])) Chris@87: Chris@87: >>> a = np.arange(5) Chris@87: >>> hist, bin_edges = np.histogram(a, density=True) Chris@87: >>> hist Chris@87: array([ 0.5, 0. , 0.5, 0. , 0. , 0.5, 0. , 0.5, 0. , 0.5]) Chris@87: >>> hist.sum() Chris@87: 2.4999999999999996 Chris@87: >>> np.sum(hist*np.diff(bin_edges)) Chris@87: 1.0 Chris@87: Chris@87: """ Chris@87: Chris@87: a = asarray(a) Chris@87: if weights is not None: Chris@87: weights = asarray(weights) Chris@87: if np.any(weights.shape != a.shape): Chris@87: raise ValueError( Chris@87: 'weights should have the same shape as a.') Chris@87: weights = weights.ravel() Chris@87: a = a.ravel() Chris@87: Chris@87: if (range is not None): Chris@87: mn, mx = range Chris@87: if (mn > mx): Chris@87: raise AttributeError( Chris@87: 'max must be larger than min in range parameter.') Chris@87: Chris@87: if not iterable(bins): Chris@87: if np.isscalar(bins) and bins < 1: Chris@87: raise ValueError( Chris@87: '`bins` should be a positive integer.') Chris@87: if range is None: Chris@87: if a.size == 0: Chris@87: # handle empty arrays. Can't determine range, so use 0-1. Chris@87: range = (0, 1) Chris@87: else: Chris@87: range = (a.min(), a.max()) Chris@87: mn, mx = [mi + 0.0 for mi in range] Chris@87: if mn == mx: Chris@87: mn -= 0.5 Chris@87: mx += 0.5 Chris@87: bins = linspace(mn, mx, bins + 1, endpoint=True) Chris@87: else: Chris@87: bins = asarray(bins) Chris@87: if (np.diff(bins) < 0).any(): Chris@87: raise AttributeError( Chris@87: 'bins must increase monotonically.') Chris@87: Chris@87: # Histogram is an integer or a float array depending on the weights. Chris@87: if weights is None: Chris@87: ntype = int Chris@87: else: Chris@87: ntype = weights.dtype Chris@87: n = np.zeros(bins.shape, ntype) Chris@87: Chris@87: block = 65536 Chris@87: if weights is None: Chris@87: for i in arange(0, len(a), block): Chris@87: sa = sort(a[i:i+block]) Chris@87: n += np.r_[sa.searchsorted(bins[:-1], 'left'), Chris@87: sa.searchsorted(bins[-1], 'right')] Chris@87: else: Chris@87: zero = array(0, dtype=ntype) Chris@87: for i in arange(0, len(a), block): Chris@87: tmp_a = a[i:i+block] Chris@87: tmp_w = weights[i:i+block] Chris@87: sorting_index = np.argsort(tmp_a) Chris@87: sa = tmp_a[sorting_index] Chris@87: sw = tmp_w[sorting_index] Chris@87: cw = np.concatenate(([zero, ], sw.cumsum())) Chris@87: bin_index = np.r_[sa.searchsorted(bins[:-1], 'left'), Chris@87: sa.searchsorted(bins[-1], 'right')] Chris@87: n += cw[bin_index] Chris@87: Chris@87: n = np.diff(n) Chris@87: Chris@87: if density is not None: Chris@87: if density: Chris@87: db = array(np.diff(bins), float) Chris@87: return n/db/n.sum(), bins Chris@87: else: Chris@87: return n, bins Chris@87: else: Chris@87: # deprecated, buggy behavior. Remove for Numpy 2.0 Chris@87: if normed: Chris@87: db = array(np.diff(bins), float) Chris@87: return n/(n*db).sum(), bins Chris@87: else: Chris@87: return n, bins Chris@87: Chris@87: Chris@87: def histogramdd(sample, bins=10, range=None, normed=False, weights=None): Chris@87: """ Chris@87: Compute the multidimensional histogram of some data. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: sample : array_like Chris@87: The data to be histogrammed. It must be an (N,D) array or data Chris@87: that can be converted to such. The rows of the resulting array Chris@87: are the coordinates of points in a D dimensional polytope. Chris@87: bins : sequence or int, optional Chris@87: The bin specification: Chris@87: Chris@87: * A sequence of arrays describing the bin edges along each dimension. Chris@87: * The number of bins for each dimension (nx, ny, ... =bins) Chris@87: * The number of bins for all dimensions (nx=ny=...=bins). Chris@87: Chris@87: range : sequence, optional Chris@87: A sequence of lower and upper bin edges to be used if the edges are Chris@87: not given explicitly in `bins`. Defaults to the minimum and maximum Chris@87: values along each dimension. Chris@87: normed : bool, optional Chris@87: If False, returns the number of samples in each bin. If True, Chris@87: returns the bin density ``bin_count / sample_count / bin_volume``. Chris@87: weights : array_like (N,), optional Chris@87: An array of values `w_i` weighing each sample `(x_i, y_i, z_i, ...)`. Chris@87: Weights are normalized to 1 if normed is True. If normed is False, Chris@87: the values of the returned histogram are equal to the sum of the Chris@87: weights belonging to the samples falling into each bin. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: H : ndarray Chris@87: The multidimensional histogram of sample x. See normed and weights Chris@87: for the different possible semantics. Chris@87: edges : list Chris@87: A list of D arrays describing the bin edges for each dimension. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: histogram: 1-D histogram Chris@87: histogram2d: 2-D histogram Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> r = np.random.randn(100,3) Chris@87: >>> H, edges = np.histogramdd(r, bins = (5, 8, 4)) Chris@87: >>> H.shape, edges[0].size, edges[1].size, edges[2].size Chris@87: ((5, 8, 4), 6, 9, 5) Chris@87: Chris@87: """ Chris@87: Chris@87: try: Chris@87: # Sample is an ND-array. Chris@87: N, D = sample.shape Chris@87: except (AttributeError, ValueError): Chris@87: # Sample is a sequence of 1D arrays. Chris@87: sample = atleast_2d(sample).T Chris@87: N, D = sample.shape Chris@87: Chris@87: nbin = empty(D, int) Chris@87: edges = D*[None] Chris@87: dedges = D*[None] Chris@87: if weights is not None: Chris@87: weights = asarray(weights) Chris@87: Chris@87: try: Chris@87: M = len(bins) Chris@87: if M != D: Chris@87: raise AttributeError( Chris@87: 'The dimension of bins must be equal to the dimension of the ' Chris@87: ' sample x.') Chris@87: except TypeError: Chris@87: # bins is an integer Chris@87: bins = D*[bins] Chris@87: Chris@87: # Select range for each dimension Chris@87: # Used only if number of bins is given. Chris@87: if range is None: Chris@87: # Handle empty input. Range can't be determined in that case, use 0-1. Chris@87: if N == 0: Chris@87: smin = zeros(D) Chris@87: smax = ones(D) Chris@87: else: Chris@87: smin = atleast_1d(array(sample.min(0), float)) Chris@87: smax = atleast_1d(array(sample.max(0), float)) Chris@87: else: Chris@87: smin = zeros(D) Chris@87: smax = zeros(D) Chris@87: for i in arange(D): Chris@87: smin[i], smax[i] = range[i] Chris@87: Chris@87: # Make sure the bins have a finite width. Chris@87: for i in arange(len(smin)): Chris@87: if smin[i] == smax[i]: Chris@87: smin[i] = smin[i] - .5 Chris@87: smax[i] = smax[i] + .5 Chris@87: Chris@87: # avoid rounding issues for comparisons when dealing with inexact types Chris@87: if np.issubdtype(sample.dtype, np.inexact): Chris@87: edge_dt = sample.dtype Chris@87: else: Chris@87: edge_dt = float Chris@87: # Create edge arrays Chris@87: for i in arange(D): Chris@87: if isscalar(bins[i]): Chris@87: if bins[i] < 1: Chris@87: raise ValueError( Chris@87: "Element at index %s in `bins` should be a positive " Chris@87: "integer." % i) Chris@87: nbin[i] = bins[i] + 2 # +2 for outlier bins Chris@87: edges[i] = linspace(smin[i], smax[i], nbin[i]-1, dtype=edge_dt) Chris@87: else: Chris@87: edges[i] = asarray(bins[i], edge_dt) Chris@87: nbin[i] = len(edges[i]) + 1 # +1 for outlier bins Chris@87: dedges[i] = diff(edges[i]) Chris@87: if np.any(np.asarray(dedges[i]) <= 0): Chris@87: raise ValueError( Chris@87: "Found bin edge of size <= 0. Did you specify `bins` with" Chris@87: "non-monotonic sequence?") Chris@87: Chris@87: nbin = asarray(nbin) Chris@87: Chris@87: # Handle empty input. Chris@87: if N == 0: Chris@87: return np.zeros(nbin-2), edges Chris@87: Chris@87: # Compute the bin number each sample falls into. Chris@87: Ncount = {} Chris@87: for i in arange(D): Chris@87: Ncount[i] = digitize(sample[:, i], edges[i]) Chris@87: Chris@87: # Using digitize, values that fall on an edge are put in the right bin. Chris@87: # For the rightmost bin, we want values equal to the right edge to be Chris@87: # counted in the last bin, and not as an outlier. Chris@87: for i in arange(D): Chris@87: # Rounding precision Chris@87: mindiff = dedges[i].min() Chris@87: if not np.isinf(mindiff): Chris@87: decimal = int(-log10(mindiff)) + 6 Chris@87: # Find which points are on the rightmost edge. Chris@87: not_smaller_than_edge = (sample[:, i] >= edges[i][-1]) Chris@87: on_edge = (around(sample[:, i], decimal) == Chris@87: around(edges[i][-1], decimal)) Chris@87: # Shift these points one bin to the left. Chris@87: Ncount[i][where(on_edge & not_smaller_than_edge)[0]] -= 1 Chris@87: Chris@87: # Flattened histogram matrix (1D) Chris@87: # Reshape is used so that overlarge arrays Chris@87: # will raise an error. Chris@87: hist = zeros(nbin, float).reshape(-1) Chris@87: Chris@87: # Compute the sample indices in the flattened histogram matrix. Chris@87: ni = nbin.argsort() Chris@87: xy = zeros(N, int) Chris@87: for i in arange(0, D-1): Chris@87: xy += Ncount[ni[i]] * nbin[ni[i+1:]].prod() Chris@87: xy += Ncount[ni[-1]] Chris@87: Chris@87: # Compute the number of repetitions in xy and assign it to the Chris@87: # flattened histmat. Chris@87: if len(xy) == 0: Chris@87: return zeros(nbin-2, int), edges Chris@87: Chris@87: flatcount = bincount(xy, weights) Chris@87: a = arange(len(flatcount)) Chris@87: hist[a] = flatcount Chris@87: Chris@87: # Shape into a proper matrix Chris@87: hist = hist.reshape(sort(nbin)) Chris@87: for i in arange(nbin.size): Chris@87: j = ni.argsort()[i] Chris@87: hist = hist.swapaxes(i, j) Chris@87: ni[i], ni[j] = ni[j], ni[i] Chris@87: Chris@87: # Remove outliers (indices 0 and -1 for each dimension). Chris@87: core = D*[slice(1, -1)] Chris@87: hist = hist[core] Chris@87: Chris@87: # Normalize if normed is True Chris@87: if normed: Chris@87: s = hist.sum() Chris@87: for i in arange(D): Chris@87: shape = ones(D, int) Chris@87: shape[i] = nbin[i] - 2 Chris@87: hist = hist / dedges[i].reshape(shape) Chris@87: hist /= s Chris@87: Chris@87: if (hist.shape != nbin - 2).any(): Chris@87: raise RuntimeError( Chris@87: "Internal Shape Error") Chris@87: return hist, edges Chris@87: Chris@87: Chris@87: def average(a, axis=None, weights=None, returned=False): Chris@87: """ Chris@87: Compute the weighted average along the specified axis. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: a : array_like Chris@87: Array containing data to be averaged. If `a` is not an array, a Chris@87: conversion is attempted. Chris@87: axis : int, optional Chris@87: Axis along which to average `a`. If `None`, averaging is done over Chris@87: the flattened array. Chris@87: weights : array_like, optional Chris@87: An array of weights associated with the values in `a`. Each value in Chris@87: `a` contributes to the average according to its associated weight. Chris@87: The weights array can either be 1-D (in which case its length must be Chris@87: the size of `a` along the given axis) or of the same shape as `a`. Chris@87: If `weights=None`, then all data in `a` are assumed to have a Chris@87: weight equal to one. Chris@87: returned : bool, optional Chris@87: Default is `False`. If `True`, the tuple (`average`, `sum_of_weights`) Chris@87: is returned, otherwise only the average is returned. Chris@87: If `weights=None`, `sum_of_weights` is equivalent to the number of Chris@87: elements over which the average is taken. Chris@87: Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: average, [sum_of_weights] : {array_type, double} Chris@87: Return the average along the specified axis. When returned is `True`, Chris@87: return a tuple with the average as the first element and the sum Chris@87: of the weights as the second element. The return type is `Float` Chris@87: if `a` is of integer type, otherwise it is of the same type as `a`. Chris@87: `sum_of_weights` is of the same type as `average`. Chris@87: Chris@87: Raises Chris@87: ------ Chris@87: ZeroDivisionError Chris@87: When all weights along axis are zero. See `numpy.ma.average` for a Chris@87: version robust to this type of error. Chris@87: TypeError Chris@87: When the length of 1D `weights` is not the same as the shape of `a` Chris@87: along axis. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: mean Chris@87: Chris@87: ma.average : average for masked arrays -- useful if your data contains Chris@87: "missing" values Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> data = range(1,5) Chris@87: >>> data Chris@87: [1, 2, 3, 4] Chris@87: >>> np.average(data) Chris@87: 2.5 Chris@87: >>> np.average(range(1,11), weights=range(10,0,-1)) Chris@87: 4.0 Chris@87: Chris@87: >>> data = np.arange(6).reshape((3,2)) Chris@87: >>> data Chris@87: array([[0, 1], Chris@87: [2, 3], Chris@87: [4, 5]]) Chris@87: >>> np.average(data, axis=1, weights=[1./4, 3./4]) Chris@87: array([ 0.75, 2.75, 4.75]) Chris@87: >>> np.average(data, weights=[1./4, 3./4]) Chris@87: Traceback (most recent call last): Chris@87: ... Chris@87: TypeError: Axis must be specified when shapes of a and weights differ. Chris@87: Chris@87: """ Chris@87: if not isinstance(a, np.matrix): Chris@87: a = np.asarray(a) Chris@87: Chris@87: if weights is None: Chris@87: avg = a.mean(axis) Chris@87: scl = avg.dtype.type(a.size/avg.size) Chris@87: else: Chris@87: a = a + 0.0 Chris@87: wgt = np.array(weights, dtype=a.dtype, copy=0) Chris@87: Chris@87: # Sanity checks Chris@87: if a.shape != wgt.shape: Chris@87: if axis is None: Chris@87: raise TypeError( Chris@87: "Axis must be specified when shapes of a and weights " Chris@87: "differ.") Chris@87: if wgt.ndim != 1: Chris@87: raise TypeError( Chris@87: "1D weights expected when shapes of a and weights differ.") Chris@87: if wgt.shape[0] != a.shape[axis]: Chris@87: raise ValueError( Chris@87: "Length of weights not compatible with specified axis.") Chris@87: Chris@87: # setup wgt to broadcast along axis Chris@87: wgt = np.array(wgt, copy=0, ndmin=a.ndim).swapaxes(-1, axis) Chris@87: Chris@87: scl = wgt.sum(axis=axis) Chris@87: if (scl == 0.0).any(): Chris@87: raise ZeroDivisionError( Chris@87: "Weights sum to zero, can't be normalized") Chris@87: Chris@87: avg = np.multiply(a, wgt).sum(axis)/scl Chris@87: Chris@87: if returned: Chris@87: scl = np.multiply(avg, 0) + scl Chris@87: return avg, scl Chris@87: else: Chris@87: return avg Chris@87: Chris@87: Chris@87: def asarray_chkfinite(a, dtype=None, order=None): Chris@87: """ Chris@87: Convert the input to an array, checking for NaNs or Infs. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: a : array_like Chris@87: Input data, in any form that can be converted to an array. This Chris@87: includes lists, lists of tuples, tuples, tuples of tuples, tuples Chris@87: of lists and ndarrays. Success requires no NaNs or Infs. Chris@87: dtype : data-type, optional Chris@87: By default, the data-type is inferred from the input data. Chris@87: order : {'C', 'F'}, optional Chris@87: Whether to use row-major ('C') or column-major ('FORTRAN') memory Chris@87: representation. Defaults to 'C'. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: Array interpretation of `a`. No copy is performed if the input Chris@87: is already an ndarray. If `a` is a subclass of ndarray, a base Chris@87: class ndarray is returned. Chris@87: Chris@87: Raises Chris@87: ------ Chris@87: ValueError Chris@87: Raises ValueError if `a` contains NaN (Not a Number) or Inf (Infinity). Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: asarray : Create and array. Chris@87: asanyarray : Similar function which passes through subclasses. Chris@87: ascontiguousarray : Convert input to a contiguous array. Chris@87: asfarray : Convert input to a floating point ndarray. Chris@87: asfortranarray : Convert input to an ndarray with column-major Chris@87: memory order. Chris@87: fromiter : Create an array from an iterator. Chris@87: fromfunction : Construct an array by executing a function on grid Chris@87: positions. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: Convert a list into an array. If all elements are finite Chris@87: ``asarray_chkfinite`` is identical to ``asarray``. Chris@87: Chris@87: >>> a = [1, 2] Chris@87: >>> np.asarray_chkfinite(a, dtype=float) Chris@87: array([1., 2.]) Chris@87: Chris@87: Raises ValueError if array_like contains Nans or Infs. Chris@87: Chris@87: >>> a = [1, 2, np.inf] Chris@87: >>> try: Chris@87: ... np.asarray_chkfinite(a) Chris@87: ... except ValueError: Chris@87: ... print 'ValueError' Chris@87: ... Chris@87: ValueError Chris@87: Chris@87: """ Chris@87: a = asarray(a, dtype=dtype, order=order) Chris@87: if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all(): Chris@87: raise ValueError( Chris@87: "array must not contain infs or NaNs") Chris@87: return a Chris@87: Chris@87: Chris@87: def piecewise(x, condlist, funclist, *args, **kw): Chris@87: """ Chris@87: Evaluate a piecewise-defined function. Chris@87: Chris@87: Given a set of conditions and corresponding functions, evaluate each Chris@87: function on the input data wherever its condition is true. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x : ndarray Chris@87: The input domain. Chris@87: condlist : list of bool arrays Chris@87: Each boolean array corresponds to a function in `funclist`. Wherever Chris@87: `condlist[i]` is True, `funclist[i](x)` is used as the output value. Chris@87: Chris@87: Each boolean array in `condlist` selects a piece of `x`, Chris@87: and should therefore be of the same shape as `x`. Chris@87: Chris@87: The length of `condlist` must correspond to that of `funclist`. Chris@87: If one extra function is given, i.e. if Chris@87: ``len(funclist) - len(condlist) == 1``, then that extra function Chris@87: is the default value, used wherever all conditions are false. Chris@87: funclist : list of callables, f(x,*args,**kw), or scalars Chris@87: Each function is evaluated over `x` wherever its corresponding Chris@87: condition is True. It should take an array as input and give an array Chris@87: or a scalar value as output. If, instead of a callable, Chris@87: a scalar is provided then a constant function (``lambda x: scalar``) is Chris@87: assumed. Chris@87: args : tuple, optional Chris@87: Any further arguments given to `piecewise` are passed to the functions Chris@87: upon execution, i.e., if called ``piecewise(..., ..., 1, 'a')``, then Chris@87: each function is called as ``f(x, 1, 'a')``. Chris@87: kw : dict, optional Chris@87: Keyword arguments used in calling `piecewise` are passed to the Chris@87: functions upon execution, i.e., if called Chris@87: ``piecewise(..., ..., lambda=1)``, then each function is called as Chris@87: ``f(x, lambda=1)``. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: The output is the same shape and type as x and is found by Chris@87: calling the functions in `funclist` on the appropriate portions of `x`, Chris@87: as defined by the boolean arrays in `condlist`. Portions not covered Chris@87: by any condition have a default value of 0. Chris@87: Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: choose, select, where Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: This is similar to choose or select, except that functions are Chris@87: evaluated on elements of `x` that satisfy the corresponding condition from Chris@87: `condlist`. Chris@87: Chris@87: The result is:: Chris@87: Chris@87: |-- Chris@87: |funclist[0](x[condlist[0]]) Chris@87: out = |funclist[1](x[condlist[1]]) Chris@87: |... Chris@87: |funclist[n2](x[condlist[n2]]) Chris@87: |-- Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: Define the sigma function, which is -1 for ``x < 0`` and +1 for ``x >= 0``. Chris@87: Chris@87: >>> x = np.linspace(-2.5, 2.5, 6) Chris@87: >>> np.piecewise(x, [x < 0, x >= 0], [-1, 1]) Chris@87: array([-1., -1., -1., 1., 1., 1.]) Chris@87: Chris@87: Define the absolute value, which is ``-x`` for ``x <0`` and ``x`` for Chris@87: ``x >= 0``. Chris@87: Chris@87: >>> np.piecewise(x, [x < 0, x >= 0], [lambda x: -x, lambda x: x]) Chris@87: array([ 2.5, 1.5, 0.5, 0.5, 1.5, 2.5]) Chris@87: Chris@87: """ Chris@87: x = asanyarray(x) Chris@87: n2 = len(funclist) Chris@87: if (isscalar(condlist) or not (isinstance(condlist[0], list) or Chris@87: isinstance(condlist[0], ndarray))): Chris@87: condlist = [condlist] Chris@87: condlist = array(condlist, dtype=bool) Chris@87: n = len(condlist) Chris@87: # This is a hack to work around problems with NumPy's Chris@87: # handling of 0-d arrays and boolean indexing with Chris@87: # numpy.bool_ scalars Chris@87: zerod = False Chris@87: if x.ndim == 0: Chris@87: x = x[None] Chris@87: zerod = True Chris@87: if condlist.shape[-1] != 1: Chris@87: condlist = condlist.T Chris@87: if n == n2 - 1: # compute the "otherwise" condition. Chris@87: totlist = np.logical_or.reduce(condlist, axis=0) Chris@87: condlist = np.vstack([condlist, ~totlist]) Chris@87: n += 1 Chris@87: if (n != n2): Chris@87: raise ValueError( Chris@87: "function list and condition list must be the same") Chris@87: Chris@87: y = zeros(x.shape, x.dtype) Chris@87: for k in range(n): Chris@87: item = funclist[k] Chris@87: if not isinstance(item, collections.Callable): Chris@87: y[condlist[k]] = item Chris@87: else: Chris@87: vals = x[condlist[k]] Chris@87: if vals.size > 0: Chris@87: y[condlist[k]] = item(vals, *args, **kw) Chris@87: if zerod: Chris@87: y = y.squeeze() Chris@87: return y Chris@87: Chris@87: Chris@87: def select(condlist, choicelist, default=0): Chris@87: """ Chris@87: Return an array drawn from elements in choicelist, depending on conditions. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: condlist : list of bool ndarrays Chris@87: The list of conditions which determine from which array in `choicelist` Chris@87: the output elements are taken. When multiple conditions are satisfied, Chris@87: the first one encountered in `condlist` is used. Chris@87: choicelist : list of ndarrays Chris@87: The list of arrays from which the output elements are taken. It has Chris@87: to be of the same length as `condlist`. Chris@87: default : scalar, optional Chris@87: The element inserted in `output` when all conditions evaluate to False. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: output : ndarray Chris@87: The output at position m is the m-th element of the array in Chris@87: `choicelist` where the m-th element of the corresponding array in Chris@87: `condlist` is True. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: where : Return elements from one of two arrays depending on condition. Chris@87: take, choose, compress, diag, diagonal Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> x = np.arange(10) Chris@87: >>> condlist = [x<3, x>5] Chris@87: >>> choicelist = [x, x**2] Chris@87: >>> np.select(condlist, choicelist) Chris@87: array([ 0, 1, 2, 0, 0, 0, 36, 49, 64, 81]) Chris@87: Chris@87: """ Chris@87: # Check the size of condlist and choicelist are the same, or abort. Chris@87: if len(condlist) != len(choicelist): Chris@87: raise ValueError( Chris@87: 'list of cases must be same length as list of conditions') Chris@87: Chris@87: # Now that the dtype is known, handle the deprecated select([], []) case Chris@87: if len(condlist) == 0: Chris@87: warnings.warn("select with an empty condition list is not possible" Chris@87: "and will be deprecated", Chris@87: DeprecationWarning) Chris@87: return np.asarray(default)[()] Chris@87: Chris@87: choicelist = [np.asarray(choice) for choice in choicelist] Chris@87: choicelist.append(np.asarray(default)) Chris@87: Chris@87: # need to get the result type before broadcasting for correct scalar Chris@87: # behaviour Chris@87: dtype = np.result_type(*choicelist) Chris@87: Chris@87: # Convert conditions to arrays and broadcast conditions and choices Chris@87: # as the shape is needed for the result. Doing it seperatly optimizes Chris@87: # for example when all choices are scalars. Chris@87: condlist = np.broadcast_arrays(*condlist) Chris@87: choicelist = np.broadcast_arrays(*choicelist) Chris@87: Chris@87: # If cond array is not an ndarray in boolean format or scalar bool, abort. Chris@87: deprecated_ints = False Chris@87: for i in range(len(condlist)): Chris@87: cond = condlist[i] Chris@87: if cond.dtype.type is not np.bool_: Chris@87: if np.issubdtype(cond.dtype, np.integer): Chris@87: # A previous implementation accepted int ndarrays accidentally. Chris@87: # Supported here deliberately, but deprecated. Chris@87: condlist[i] = condlist[i].astype(bool) Chris@87: deprecated_ints = True Chris@87: else: Chris@87: raise ValueError( Chris@87: 'invalid entry in choicelist: should be boolean ndarray') Chris@87: Chris@87: if deprecated_ints: Chris@87: msg = "select condlists containing integer ndarrays is deprecated " \ Chris@87: "and will be removed in the future. Use `.astype(bool)` to " \ Chris@87: "convert to bools." Chris@87: warnings.warn(msg, DeprecationWarning) Chris@87: Chris@87: if choicelist[0].ndim == 0: Chris@87: # This may be common, so avoid the call. Chris@87: result_shape = condlist[0].shape Chris@87: else: Chris@87: result_shape = np.broadcast_arrays(condlist[0], choicelist[0])[0].shape Chris@87: Chris@87: result = np.full(result_shape, choicelist[-1], dtype) Chris@87: Chris@87: # Use np.copyto to burn each choicelist array onto result, using the Chris@87: # corresponding condlist as a boolean mask. This is done in reverse Chris@87: # order since the first choice should take precedence. Chris@87: choicelist = choicelist[-2::-1] Chris@87: condlist = condlist[::-1] Chris@87: for choice, cond in zip(choicelist, condlist): Chris@87: np.copyto(result, choice, where=cond) Chris@87: Chris@87: return result Chris@87: Chris@87: Chris@87: def copy(a, order='K'): Chris@87: """ Chris@87: Return an array copy of the given object. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: a : array_like Chris@87: Input data. Chris@87: order : {'C', 'F', 'A', 'K'}, optional Chris@87: Controls the memory layout of the copy. 'C' means C-order, Chris@87: 'F' means F-order, 'A' means 'F' if `a` is Fortran contiguous, Chris@87: 'C' otherwise. 'K' means match the layout of `a` as closely Chris@87: as possible. (Note that this function and :meth:ndarray.copy are very Chris@87: similar, but have different default values for their order= Chris@87: arguments.) Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: arr : ndarray Chris@87: Array interpretation of `a`. Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: This is equivalent to Chris@87: Chris@87: >>> np.array(a, copy=True) #doctest: +SKIP Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: Create an array x, with a reference y and a copy z: Chris@87: Chris@87: >>> x = np.array([1, 2, 3]) Chris@87: >>> y = x Chris@87: >>> z = np.copy(x) Chris@87: Chris@87: Note that, when we modify x, y changes, but not z: Chris@87: Chris@87: >>> x[0] = 10 Chris@87: >>> x[0] == y[0] Chris@87: True Chris@87: >>> x[0] == z[0] Chris@87: False Chris@87: Chris@87: """ Chris@87: return array(a, order=order, copy=True) Chris@87: Chris@87: # Basic operations Chris@87: Chris@87: Chris@87: def gradient(f, *varargs, **kwargs): Chris@87: """ Chris@87: Return the gradient of an N-dimensional array. Chris@87: Chris@87: The gradient is computed using second order accurate central differences Chris@87: in the interior and either first differences or second order accurate Chris@87: one-sides (forward or backwards) differences at the boundaries. The Chris@87: returned gradient hence has the same shape as the input array. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: f : array_like Chris@87: An N-dimensional array containing samples of a scalar function. Chris@87: varargs : list of scalar, optional Chris@87: N scalars specifying the sample distances for each dimension, Chris@87: i.e. `dx`, `dy`, `dz`, ... Default distance: 1. Chris@87: edge_order : {1, 2}, optional Chris@87: Gradient is calculated using N\ :sup:`th` order accurate differences Chris@87: at the boundaries. Default: 1. Chris@87: Chris@87: .. versionadded:: 1.9.1 Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: gradient : ndarray Chris@87: N arrays of the same shape as `f` giving the derivative of `f` with Chris@87: respect to each dimension. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> x = np.array([1, 2, 4, 7, 11, 16], dtype=np.float) Chris@87: >>> np.gradient(x) Chris@87: array([ 1. , 1.5, 2.5, 3.5, 4.5, 5. ]) Chris@87: >>> np.gradient(x, 2) Chris@87: array([ 0.5 , 0.75, 1.25, 1.75, 2.25, 2.5 ]) Chris@87: Chris@87: >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=np.float)) Chris@87: [array([[ 2., 2., -1.], Chris@87: [ 2., 2., -1.]]), array([[ 1. , 2.5, 4. ], Chris@87: [ 1. , 1. , 1. ]])] Chris@87: Chris@87: >>> x = np.array([0, 1, 2, 3, 4]) Chris@87: >>> dx = np.gradient(x) Chris@87: >>> y = x**2 Chris@87: >>> np.gradient(y, dx, edge_order=2) Chris@87: array([-0., 2., 4., 6., 8.]) Chris@87: """ Chris@87: f = np.asanyarray(f) Chris@87: N = len(f.shape) # number of dimensions Chris@87: n = len(varargs) Chris@87: if n == 0: Chris@87: dx = [1.0]*N Chris@87: elif n == 1: Chris@87: dx = [varargs[0]]*N Chris@87: elif n == N: Chris@87: dx = list(varargs) Chris@87: else: Chris@87: raise SyntaxError( Chris@87: "invalid number of arguments") Chris@87: Chris@87: edge_order = kwargs.pop('edge_order', 1) Chris@87: if kwargs: Chris@87: raise TypeError('"{}" are not valid keyword arguments.'.format( Chris@87: '", "'.join(kwargs.keys()))) Chris@87: if edge_order > 2: Chris@87: raise ValueError("'edge_order' greater than 2 not supported") Chris@87: Chris@87: # use central differences on interior and one-sided differences on the Chris@87: # endpoints. This preserves second order-accuracy over the full domain. Chris@87: Chris@87: outvals = [] Chris@87: Chris@87: # create slice objects --- initially all are [:, :, ..., :] Chris@87: slice1 = [slice(None)]*N Chris@87: slice2 = [slice(None)]*N Chris@87: slice3 = [slice(None)]*N Chris@87: slice4 = [slice(None)]*N Chris@87: Chris@87: otype = f.dtype.char Chris@87: if otype not in ['f', 'd', 'F', 'D', 'm', 'M']: Chris@87: otype = 'd' Chris@87: Chris@87: # Difference of datetime64 elements results in timedelta64 Chris@87: if otype == 'M': Chris@87: # Need to use the full dtype name because it contains unit information Chris@87: otype = f.dtype.name.replace('datetime', 'timedelta') Chris@87: elif otype == 'm': Chris@87: # Needs to keep the specific units, can't be a general unit Chris@87: otype = f.dtype Chris@87: Chris@87: # Convert datetime64 data into ints. Make dummy variable `y` Chris@87: # that is a view of ints if the data is datetime64, otherwise Chris@87: # just set y equal to the the array `f`. Chris@87: if f.dtype.char in ["M", "m"]: Chris@87: y = f.view('int64') Chris@87: else: Chris@87: y = f Chris@87: Chris@87: for axis in range(N): Chris@87: Chris@87: if y.shape[axis] < 2: Chris@87: raise ValueError( Chris@87: "Shape of array too small to calculate a numerical gradient, " Chris@87: "at least two elements are required.") Chris@87: Chris@87: # Numerical differentiation: 1st order edges, 2nd order interior Chris@87: if y.shape[axis] == 2 or edge_order == 1: Chris@87: # Use first order differences for time data Chris@87: out = np.empty_like(y, dtype=otype) Chris@87: Chris@87: slice1[axis] = slice(1, -1) Chris@87: slice2[axis] = slice(2, None) Chris@87: slice3[axis] = slice(None, -2) Chris@87: # 1D equivalent -- out[1:-1] = (y[2:] - y[:-2])/2.0 Chris@87: out[slice1] = (y[slice2] - y[slice3])/2.0 Chris@87: Chris@87: slice1[axis] = 0 Chris@87: slice2[axis] = 1 Chris@87: slice3[axis] = 0 Chris@87: # 1D equivalent -- out[0] = (y[1] - y[0]) Chris@87: out[slice1] = (y[slice2] - y[slice3]) Chris@87: Chris@87: slice1[axis] = -1 Chris@87: slice2[axis] = -1 Chris@87: slice3[axis] = -2 Chris@87: # 1D equivalent -- out[-1] = (y[-1] - y[-2]) Chris@87: out[slice1] = (y[slice2] - y[slice3]) Chris@87: Chris@87: # Numerical differentiation: 2st order edges, 2nd order interior Chris@87: else: Chris@87: # Use second order differences where possible Chris@87: out = np.empty_like(y, dtype=otype) Chris@87: Chris@87: slice1[axis] = slice(1, -1) Chris@87: slice2[axis] = slice(2, None) Chris@87: slice3[axis] = slice(None, -2) Chris@87: # 1D equivalent -- out[1:-1] = (y[2:] - y[:-2])/2.0 Chris@87: out[slice1] = (y[slice2] - y[slice3])/2.0 Chris@87: Chris@87: slice1[axis] = 0 Chris@87: slice2[axis] = 0 Chris@87: slice3[axis] = 1 Chris@87: slice4[axis] = 2 Chris@87: # 1D equivalent -- out[0] = -(3*y[0] - 4*y[1] + y[2]) / 2.0 Chris@87: out[slice1] = -(3.0*y[slice2] - 4.0*y[slice3] + y[slice4])/2.0 Chris@87: Chris@87: slice1[axis] = -1 Chris@87: slice2[axis] = -1 Chris@87: slice3[axis] = -2 Chris@87: slice4[axis] = -3 Chris@87: # 1D equivalent -- out[-1] = (3*y[-1] - 4*y[-2] + y[-3]) Chris@87: out[slice1] = (3.0*y[slice2] - 4.0*y[slice3] + y[slice4])/2.0 Chris@87: Chris@87: # divide by step size Chris@87: out /= dx[axis] Chris@87: outvals.append(out) Chris@87: Chris@87: # reset the slice object in this dimension to ":" Chris@87: slice1[axis] = slice(None) Chris@87: slice2[axis] = slice(None) Chris@87: slice3[axis] = slice(None) Chris@87: slice4[axis] = slice(None) Chris@87: Chris@87: if N == 1: Chris@87: return outvals[0] Chris@87: else: Chris@87: return outvals Chris@87: Chris@87: Chris@87: def diff(a, n=1, axis=-1): Chris@87: """ Chris@87: Calculate the n-th order discrete difference along given axis. Chris@87: Chris@87: The first order difference is given by ``out[n] = a[n+1] - a[n]`` along Chris@87: the given axis, higher order differences are calculated by using `diff` Chris@87: recursively. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: a : array_like Chris@87: Input array Chris@87: n : int, optional Chris@87: The number of times values are differenced. Chris@87: axis : int, optional Chris@87: The axis along which the difference is taken, default is the last axis. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: diff : ndarray Chris@87: The `n` order differences. The shape of the output is the same as `a` Chris@87: except along `axis` where the dimension is smaller by `n`. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: gradient, ediff1d, cumsum Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> x = np.array([1, 2, 4, 7, 0]) Chris@87: >>> np.diff(x) Chris@87: array([ 1, 2, 3, -7]) Chris@87: >>> np.diff(x, n=2) Chris@87: array([ 1, 1, -10]) Chris@87: Chris@87: >>> x = np.array([[1, 3, 6, 10], [0, 5, 6, 8]]) Chris@87: >>> np.diff(x) Chris@87: array([[2, 3, 4], Chris@87: [5, 1, 2]]) Chris@87: >>> np.diff(x, axis=0) Chris@87: array([[-1, 2, 0, -2]]) Chris@87: Chris@87: """ Chris@87: if n == 0: Chris@87: return a Chris@87: if n < 0: Chris@87: raise ValueError( Chris@87: "order must be non-negative but got " + repr(n)) Chris@87: a = asanyarray(a) Chris@87: nd = len(a.shape) Chris@87: slice1 = [slice(None)]*nd Chris@87: slice2 = [slice(None)]*nd Chris@87: slice1[axis] = slice(1, None) Chris@87: slice2[axis] = slice(None, -1) Chris@87: slice1 = tuple(slice1) Chris@87: slice2 = tuple(slice2) Chris@87: if n > 1: Chris@87: return diff(a[slice1]-a[slice2], n-1, axis=axis) Chris@87: else: Chris@87: return a[slice1]-a[slice2] Chris@87: Chris@87: Chris@87: def interp(x, xp, fp, left=None, right=None): Chris@87: """ Chris@87: One-dimensional linear interpolation. Chris@87: Chris@87: Returns the one-dimensional piecewise linear interpolant to a function Chris@87: with given values at discrete data-points. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x : array_like Chris@87: The x-coordinates of the interpolated values. Chris@87: Chris@87: xp : 1-D sequence of floats Chris@87: The x-coordinates of the data points, must be increasing. Chris@87: Chris@87: fp : 1-D sequence of floats Chris@87: The y-coordinates of the data points, same length as `xp`. Chris@87: Chris@87: left : float, optional Chris@87: Value to return for `x < xp[0]`, default is `fp[0]`. Chris@87: Chris@87: right : float, optional Chris@87: Value to return for `x > xp[-1]`, default is `fp[-1]`. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: y : {float, ndarray} Chris@87: The interpolated values, same shape as `x`. Chris@87: Chris@87: Raises Chris@87: ------ Chris@87: ValueError Chris@87: If `xp` and `fp` have different length Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: Does not check that the x-coordinate sequence `xp` is increasing. Chris@87: If `xp` is not increasing, the results are nonsense. Chris@87: A simple check for increasing is:: Chris@87: Chris@87: np.all(np.diff(xp) > 0) Chris@87: Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> xp = [1, 2, 3] Chris@87: >>> fp = [3, 2, 0] Chris@87: >>> np.interp(2.5, xp, fp) Chris@87: 1.0 Chris@87: >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp) Chris@87: array([ 3. , 3. , 2.5 , 0.56, 0. ]) Chris@87: >>> UNDEF = -99.0 Chris@87: >>> np.interp(3.14, xp, fp, right=UNDEF) Chris@87: -99.0 Chris@87: Chris@87: Plot an interpolant to the sine function: Chris@87: Chris@87: >>> x = np.linspace(0, 2*np.pi, 10) Chris@87: >>> y = np.sin(x) Chris@87: >>> xvals = np.linspace(0, 2*np.pi, 50) Chris@87: >>> yinterp = np.interp(xvals, x, y) Chris@87: >>> import matplotlib.pyplot as plt Chris@87: >>> plt.plot(x, y, 'o') Chris@87: [] Chris@87: >>> plt.plot(xvals, yinterp, '-x') Chris@87: [] Chris@87: >>> plt.show() Chris@87: Chris@87: """ Chris@87: if isinstance(x, (float, int, number)): Chris@87: return compiled_interp([x], xp, fp, left, right).item() Chris@87: elif isinstance(x, np.ndarray) and x.ndim == 0: Chris@87: return compiled_interp([x], xp, fp, left, right).item() Chris@87: else: Chris@87: return compiled_interp(x, xp, fp, left, right) Chris@87: Chris@87: Chris@87: def angle(z, deg=0): Chris@87: """ Chris@87: Return the angle of the complex argument. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: z : array_like Chris@87: A complex number or sequence of complex numbers. Chris@87: deg : bool, optional Chris@87: Return angle in degrees if True, radians if False (default). Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: angle : {ndarray, scalar} Chris@87: The counterclockwise angle from the positive real axis on Chris@87: the complex plane, with dtype as numpy.float64. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: arctan2 Chris@87: absolute Chris@87: Chris@87: Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> np.angle([1.0, 1.0j, 1+1j]) # in radians Chris@87: array([ 0. , 1.57079633, 0.78539816]) Chris@87: >>> np.angle(1+1j, deg=True) # in degrees Chris@87: 45.0 Chris@87: Chris@87: """ Chris@87: if deg: Chris@87: fact = 180/pi Chris@87: else: Chris@87: fact = 1.0 Chris@87: z = asarray(z) Chris@87: if (issubclass(z.dtype.type, _nx.complexfloating)): Chris@87: zimag = z.imag Chris@87: zreal = z.real Chris@87: else: Chris@87: zimag = 0 Chris@87: zreal = z Chris@87: return arctan2(zimag, zreal) * fact Chris@87: Chris@87: Chris@87: def unwrap(p, discont=pi, axis=-1): Chris@87: """ Chris@87: Unwrap by changing deltas between values to 2*pi complement. Chris@87: Chris@87: Unwrap radian phase `p` by changing absolute jumps greater than Chris@87: `discont` to their 2*pi complement along the given axis. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: p : array_like Chris@87: Input array. Chris@87: discont : float, optional Chris@87: Maximum discontinuity between values, default is ``pi``. Chris@87: axis : int, optional Chris@87: Axis along which unwrap will operate, default is the last axis. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: Output array. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: rad2deg, deg2rad Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: If the discontinuity in `p` is smaller than ``pi``, but larger than Chris@87: `discont`, no unwrapping is done because taking the 2*pi complement Chris@87: would only make the discontinuity larger. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> phase = np.linspace(0, np.pi, num=5) Chris@87: >>> phase[3:] += np.pi Chris@87: >>> phase Chris@87: array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531]) Chris@87: >>> np.unwrap(phase) Chris@87: array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ]) Chris@87: Chris@87: """ Chris@87: p = asarray(p) Chris@87: nd = len(p.shape) Chris@87: dd = diff(p, axis=axis) Chris@87: slice1 = [slice(None, None)]*nd # full slices Chris@87: slice1[axis] = slice(1, None) Chris@87: ddmod = mod(dd + pi, 2*pi) - pi Chris@87: _nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0)) Chris@87: ph_correct = ddmod - dd Chris@87: _nx.copyto(ph_correct, 0, where=abs(dd) < discont) Chris@87: up = array(p, copy=True, dtype='d') Chris@87: up[slice1] = p[slice1] + ph_correct.cumsum(axis) Chris@87: return up Chris@87: Chris@87: Chris@87: def sort_complex(a): Chris@87: """ Chris@87: Sort a complex array using the real part first, then the imaginary part. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: a : array_like Chris@87: Input array Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : complex ndarray Chris@87: Always returns a sorted complex array. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> np.sort_complex([5, 3, 6, 2, 1]) Chris@87: array([ 1.+0.j, 2.+0.j, 3.+0.j, 5.+0.j, 6.+0.j]) Chris@87: Chris@87: >>> np.sort_complex([1 + 2j, 2 - 1j, 3 - 2j, 3 - 3j, 3 + 5j]) Chris@87: array([ 1.+2.j, 2.-1.j, 3.-3.j, 3.-2.j, 3.+5.j]) Chris@87: Chris@87: """ Chris@87: b = array(a, copy=True) Chris@87: b.sort() Chris@87: if not issubclass(b.dtype.type, _nx.complexfloating): Chris@87: if b.dtype.char in 'bhBH': Chris@87: return b.astype('F') Chris@87: elif b.dtype.char == 'g': Chris@87: return b.astype('G') Chris@87: else: Chris@87: return b.astype('D') Chris@87: else: Chris@87: return b Chris@87: Chris@87: Chris@87: def trim_zeros(filt, trim='fb'): Chris@87: """ Chris@87: Trim the leading and/or trailing zeros from a 1-D array or sequence. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: filt : 1-D array or sequence Chris@87: Input array. Chris@87: trim : str, optional Chris@87: A string with 'f' representing trim from front and 'b' to trim from Chris@87: back. Default is 'fb', trim zeros from both front and back of the Chris@87: array. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: trimmed : 1-D array or sequence Chris@87: The result of trimming the input. The input data type is preserved. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> a = np.array((0, 0, 0, 1, 2, 3, 0, 2, 1, 0)) Chris@87: >>> np.trim_zeros(a) Chris@87: array([1, 2, 3, 0, 2, 1]) Chris@87: Chris@87: >>> np.trim_zeros(a, 'b') Chris@87: array([0, 0, 0, 1, 2, 3, 0, 2, 1]) Chris@87: Chris@87: The input data type is preserved, list/tuple in means list/tuple out. Chris@87: Chris@87: >>> np.trim_zeros([0, 1, 2, 0]) Chris@87: [1, 2] Chris@87: Chris@87: """ Chris@87: first = 0 Chris@87: trim = trim.upper() Chris@87: if 'F' in trim: Chris@87: for i in filt: Chris@87: if i != 0.: Chris@87: break Chris@87: else: Chris@87: first = first + 1 Chris@87: last = len(filt) Chris@87: if 'B' in trim: Chris@87: for i in filt[::-1]: Chris@87: if i != 0.: Chris@87: break Chris@87: else: Chris@87: last = last - 1 Chris@87: return filt[first:last] Chris@87: Chris@87: Chris@87: @deprecate Chris@87: def unique(x): Chris@87: """ Chris@87: This function is deprecated. Use numpy.lib.arraysetops.unique() Chris@87: instead. Chris@87: """ Chris@87: try: Chris@87: tmp = x.flatten() Chris@87: if tmp.size == 0: Chris@87: return tmp Chris@87: tmp.sort() Chris@87: idx = concatenate(([True], tmp[1:] != tmp[:-1])) Chris@87: return tmp[idx] Chris@87: except AttributeError: Chris@87: items = sorted(set(x)) Chris@87: return asarray(items) Chris@87: Chris@87: Chris@87: def extract(condition, arr): Chris@87: """ Chris@87: Return the elements of an array that satisfy some condition. Chris@87: Chris@87: This is equivalent to ``np.compress(ravel(condition), ravel(arr))``. If Chris@87: `condition` is boolean ``np.extract`` is equivalent to ``arr[condition]``. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: condition : array_like Chris@87: An array whose nonzero or True entries indicate the elements of `arr` Chris@87: to extract. Chris@87: arr : array_like Chris@87: Input array of the same size as `condition`. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: extract : ndarray Chris@87: Rank 1 array of values from `arr` where `condition` is True. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: take, put, copyto, compress Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> arr = np.arange(12).reshape((3, 4)) Chris@87: >>> arr Chris@87: array([[ 0, 1, 2, 3], Chris@87: [ 4, 5, 6, 7], Chris@87: [ 8, 9, 10, 11]]) Chris@87: >>> condition = np.mod(arr, 3)==0 Chris@87: >>> condition Chris@87: array([[ True, False, False, True], Chris@87: [False, False, True, False], Chris@87: [False, True, False, False]], dtype=bool) Chris@87: >>> np.extract(condition, arr) Chris@87: array([0, 3, 6, 9]) Chris@87: Chris@87: Chris@87: If `condition` is boolean: Chris@87: Chris@87: >>> arr[condition] Chris@87: array([0, 3, 6, 9]) Chris@87: Chris@87: """ Chris@87: return _nx.take(ravel(arr), nonzero(ravel(condition))[0]) Chris@87: Chris@87: Chris@87: def place(arr, mask, vals): Chris@87: """ Chris@87: Change elements of an array based on conditional and input values. Chris@87: Chris@87: Similar to ``np.copyto(arr, vals, where=mask)``, the difference is that Chris@87: `place` uses the first N elements of `vals`, where N is the number of Chris@87: True values in `mask`, while `copyto` uses the elements where `mask` Chris@87: is True. Chris@87: Chris@87: Note that `extract` does the exact opposite of `place`. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: arr : array_like Chris@87: Array to put data into. Chris@87: mask : array_like Chris@87: Boolean mask array. Must have the same size as `a`. Chris@87: vals : 1-D sequence Chris@87: Values to put into `a`. Only the first N elements are used, where Chris@87: N is the number of True values in `mask`. If `vals` is smaller Chris@87: than N it will be repeated. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: copyto, put, take, extract Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> arr = np.arange(6).reshape(2, 3) Chris@87: >>> np.place(arr, arr>2, [44, 55]) Chris@87: >>> arr Chris@87: array([[ 0, 1, 2], Chris@87: [44, 55, 44]]) Chris@87: Chris@87: """ Chris@87: return _insert(arr, mask, vals) Chris@87: Chris@87: Chris@87: def disp(mesg, device=None, linefeed=True): Chris@87: """ Chris@87: Display a message on a device. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: mesg : str Chris@87: Message to display. Chris@87: device : object Chris@87: Device to write message. If None, defaults to ``sys.stdout`` which is Chris@87: very similar to ``print``. `device` needs to have ``write()`` and Chris@87: ``flush()`` methods. Chris@87: linefeed : bool, optional Chris@87: Option whether to print a line feed or not. Defaults to True. Chris@87: Chris@87: Raises Chris@87: ------ Chris@87: AttributeError Chris@87: If `device` does not have a ``write()`` or ``flush()`` method. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: Besides ``sys.stdout``, a file-like object can also be used as it has Chris@87: both required methods: Chris@87: Chris@87: >>> from StringIO import StringIO Chris@87: >>> buf = StringIO() Chris@87: >>> np.disp('"Display" in a file', device=buf) Chris@87: >>> buf.getvalue() Chris@87: '"Display" in a file\\n' Chris@87: Chris@87: """ Chris@87: if device is None: Chris@87: device = sys.stdout Chris@87: if linefeed: Chris@87: device.write('%s\n' % mesg) Chris@87: else: Chris@87: device.write('%s' % mesg) Chris@87: device.flush() Chris@87: return Chris@87: Chris@87: Chris@87: class vectorize(object): Chris@87: """ Chris@87: vectorize(pyfunc, otypes='', doc=None, excluded=None, cache=False) Chris@87: Chris@87: Generalized function class. Chris@87: Chris@87: Define a vectorized function which takes a nested sequence Chris@87: of objects or numpy arrays as inputs and returns a Chris@87: numpy array as output. The vectorized function evaluates `pyfunc` over Chris@87: successive tuples of the input arrays like the python map function, Chris@87: except it uses the broadcasting rules of numpy. Chris@87: Chris@87: The data type of the output of `vectorized` is determined by calling Chris@87: the function with the first element of the input. This can be avoided Chris@87: by specifying the `otypes` argument. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: pyfunc : callable Chris@87: A python function or method. Chris@87: otypes : str or list of dtypes, optional Chris@87: The output data type. It must be specified as either a string of Chris@87: typecode characters or a list of data type specifiers. There should Chris@87: be one data type specifier for each output. Chris@87: doc : str, optional Chris@87: The docstring for the function. If `None`, the docstring will be the Chris@87: ``pyfunc.__doc__``. Chris@87: excluded : set, optional Chris@87: Set of strings or integers representing the positional or keyword Chris@87: arguments for which the function will not be vectorized. These will be Chris@87: passed directly to `pyfunc` unmodified. Chris@87: Chris@87: .. versionadded:: 1.7.0 Chris@87: Chris@87: cache : bool, optional Chris@87: If `True`, then cache the first function call that determines the number Chris@87: of outputs if `otypes` is not provided. Chris@87: Chris@87: .. versionadded:: 1.7.0 Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: vectorized : callable Chris@87: Vectorized function. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> def myfunc(a, b): Chris@87: ... "Return a-b if a>b, otherwise return a+b" Chris@87: ... if a > b: Chris@87: ... return a - b Chris@87: ... else: Chris@87: ... return a + b Chris@87: Chris@87: >>> vfunc = np.vectorize(myfunc) Chris@87: >>> vfunc([1, 2, 3, 4], 2) Chris@87: array([3, 4, 1, 2]) Chris@87: Chris@87: The docstring is taken from the input function to `vectorize` unless it Chris@87: is specified Chris@87: Chris@87: >>> vfunc.__doc__ Chris@87: 'Return a-b if a>b, otherwise return a+b' Chris@87: >>> vfunc = np.vectorize(myfunc, doc='Vectorized `myfunc`') Chris@87: >>> vfunc.__doc__ Chris@87: 'Vectorized `myfunc`' Chris@87: Chris@87: The output type is determined by evaluating the first element of the input, Chris@87: unless it is specified Chris@87: Chris@87: >>> out = vfunc([1, 2, 3, 4], 2) Chris@87: >>> type(out[0]) Chris@87: Chris@87: >>> vfunc = np.vectorize(myfunc, otypes=[np.float]) Chris@87: >>> out = vfunc([1, 2, 3, 4], 2) Chris@87: >>> type(out[0]) Chris@87: Chris@87: Chris@87: The `excluded` argument can be used to prevent vectorizing over certain Chris@87: arguments. This can be useful for array-like arguments of a fixed length Chris@87: such as the coefficients for a polynomial as in `polyval`: Chris@87: Chris@87: >>> def mypolyval(p, x): Chris@87: ... _p = list(p) Chris@87: ... res = _p.pop(0) Chris@87: ... while _p: Chris@87: ... res = res*x + _p.pop(0) Chris@87: ... return res Chris@87: >>> vpolyval = np.vectorize(mypolyval, excluded=['p']) Chris@87: >>> vpolyval(p=[1, 2, 3], x=[0, 1]) Chris@87: array([3, 6]) Chris@87: Chris@87: Positional arguments may also be excluded by specifying their position: Chris@87: Chris@87: >>> vpolyval.excluded.add(0) Chris@87: >>> vpolyval([1, 2, 3], x=[0, 1]) Chris@87: array([3, 6]) Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: The `vectorize` function is provided primarily for convenience, not for Chris@87: performance. The implementation is essentially a for loop. Chris@87: Chris@87: If `otypes` is not specified, then a call to the function with the Chris@87: first argument will be used to determine the number of outputs. The Chris@87: results of this call will be cached if `cache` is `True` to prevent Chris@87: calling the function twice. However, to implement the cache, the Chris@87: original function must be wrapped which will slow down subsequent Chris@87: calls, so only do this if your function is expensive. Chris@87: Chris@87: The new keyword argument interface and `excluded` argument support Chris@87: further degrades performance. Chris@87: Chris@87: """ Chris@87: Chris@87: def __init__(self, pyfunc, otypes='', doc=None, excluded=None, Chris@87: cache=False): Chris@87: self.pyfunc = pyfunc Chris@87: self.cache = cache Chris@87: self._ufunc = None # Caching to improve default performance Chris@87: Chris@87: if doc is None: Chris@87: self.__doc__ = pyfunc.__doc__ Chris@87: else: Chris@87: self.__doc__ = doc Chris@87: Chris@87: if isinstance(otypes, str): Chris@87: self.otypes = otypes Chris@87: for char in self.otypes: Chris@87: if char not in typecodes['All']: Chris@87: raise ValueError( Chris@87: "Invalid otype specified: %s" % (char,)) Chris@87: elif iterable(otypes): Chris@87: self.otypes = ''.join([_nx.dtype(x).char for x in otypes]) Chris@87: else: Chris@87: raise ValueError( Chris@87: "Invalid otype specification") Chris@87: Chris@87: # Excluded variable support Chris@87: if excluded is None: Chris@87: excluded = set() Chris@87: self.excluded = set(excluded) Chris@87: Chris@87: def __call__(self, *args, **kwargs): Chris@87: """ Chris@87: Return arrays with the results of `pyfunc` broadcast (vectorized) over Chris@87: `args` and `kwargs` not in `excluded`. Chris@87: """ Chris@87: excluded = self.excluded Chris@87: if not kwargs and not excluded: Chris@87: func = self.pyfunc Chris@87: vargs = args Chris@87: else: Chris@87: # The wrapper accepts only positional arguments: we use `names` and Chris@87: # `inds` to mutate `the_args` and `kwargs` to pass to the original Chris@87: # function. Chris@87: nargs = len(args) Chris@87: Chris@87: names = [_n for _n in kwargs if _n not in excluded] Chris@87: inds = [_i for _i in range(nargs) if _i not in excluded] Chris@87: the_args = list(args) Chris@87: Chris@87: def func(*vargs): Chris@87: for _n, _i in enumerate(inds): Chris@87: the_args[_i] = vargs[_n] Chris@87: kwargs.update(zip(names, vargs[len(inds):])) Chris@87: return self.pyfunc(*the_args, **kwargs) Chris@87: Chris@87: vargs = [args[_i] for _i in inds] Chris@87: vargs.extend([kwargs[_n] for _n in names]) Chris@87: Chris@87: return self._vectorize_call(func=func, args=vargs) Chris@87: Chris@87: def _get_ufunc_and_otypes(self, func, args): Chris@87: """Return (ufunc, otypes).""" Chris@87: # frompyfunc will fail if args is empty Chris@87: if not args: Chris@87: raise ValueError('args can not be empty') Chris@87: Chris@87: if self.otypes: Chris@87: otypes = self.otypes Chris@87: nout = len(otypes) Chris@87: Chris@87: # Note logic here: We only *use* self._ufunc if func is self.pyfunc Chris@87: # even though we set self._ufunc regardless. Chris@87: if func is self.pyfunc and self._ufunc is not None: Chris@87: ufunc = self._ufunc Chris@87: else: Chris@87: ufunc = self._ufunc = frompyfunc(func, len(args), nout) Chris@87: else: Chris@87: # Get number of outputs and output types by calling the function on Chris@87: # the first entries of args. We also cache the result to prevent Chris@87: # the subsequent call when the ufunc is evaluated. Chris@87: # Assumes that ufunc first evaluates the 0th elements in the input Chris@87: # arrays (the input values are not checked to ensure this) Chris@87: inputs = [asarray(_a).flat[0] for _a in args] Chris@87: outputs = func(*inputs) Chris@87: Chris@87: # Performance note: profiling indicates that -- for simple Chris@87: # functions at least -- this wrapping can almost double the Chris@87: # execution time. Chris@87: # Hence we make it optional. Chris@87: if self.cache: Chris@87: _cache = [outputs] Chris@87: Chris@87: def _func(*vargs): Chris@87: if _cache: Chris@87: return _cache.pop() Chris@87: else: Chris@87: return func(*vargs) Chris@87: else: Chris@87: _func = func Chris@87: Chris@87: if isinstance(outputs, tuple): Chris@87: nout = len(outputs) Chris@87: else: Chris@87: nout = 1 Chris@87: outputs = (outputs,) Chris@87: Chris@87: otypes = ''.join([asarray(outputs[_k]).dtype.char Chris@87: for _k in range(nout)]) Chris@87: Chris@87: # Performance note: profiling indicates that creating the ufunc is Chris@87: # not a significant cost compared with wrapping so it seems not Chris@87: # worth trying to cache this. Chris@87: ufunc = frompyfunc(_func, len(args), nout) Chris@87: Chris@87: return ufunc, otypes Chris@87: Chris@87: def _vectorize_call(self, func, args): Chris@87: """Vectorized call to `func` over positional `args`.""" Chris@87: if not args: Chris@87: _res = func() Chris@87: else: Chris@87: ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args) Chris@87: Chris@87: # Convert args to object arrays first Chris@87: inputs = [array(_a, copy=False, subok=True, dtype=object) Chris@87: for _a in args] Chris@87: Chris@87: outputs = ufunc(*inputs) Chris@87: Chris@87: if ufunc.nout == 1: Chris@87: _res = array(outputs, Chris@87: copy=False, subok=True, dtype=otypes[0]) Chris@87: else: Chris@87: _res = tuple([array(_x, copy=False, subok=True, dtype=_t) Chris@87: for _x, _t in zip(outputs, otypes)]) Chris@87: return _res Chris@87: Chris@87: Chris@87: def cov(m, y=None, rowvar=1, bias=0, ddof=None): Chris@87: """ Chris@87: Estimate a covariance matrix, given data. Chris@87: Chris@87: Covariance indicates the level to which two variables vary together. Chris@87: If we examine N-dimensional samples, :math:`X = [x_1, x_2, ... x_N]^T`, Chris@87: then the covariance matrix element :math:`C_{ij}` is the covariance of Chris@87: :math:`x_i` and :math:`x_j`. The element :math:`C_{ii}` is the variance Chris@87: of :math:`x_i`. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: m : array_like Chris@87: A 1-D or 2-D array containing multiple variables and observations. Chris@87: Each row of `m` represents a variable, and each column a single Chris@87: observation of all those variables. Also see `rowvar` below. Chris@87: y : array_like, optional Chris@87: An additional set of variables and observations. `y` has the same Chris@87: form as that of `m`. Chris@87: rowvar : int, optional Chris@87: If `rowvar` is non-zero (default), then each row represents a Chris@87: variable, with observations in the columns. Otherwise, the relationship Chris@87: is transposed: each column represents a variable, while the rows Chris@87: contain observations. Chris@87: bias : int, optional Chris@87: Default normalization is by ``(N - 1)``, where ``N`` is the number of Chris@87: observations given (unbiased estimate). If `bias` is 1, then Chris@87: normalization is by ``N``. These values can be overridden by using Chris@87: the keyword ``ddof`` in numpy versions >= 1.5. Chris@87: ddof : int, optional Chris@87: .. versionadded:: 1.5 Chris@87: If not ``None`` normalization is by ``(N - ddof)``, where ``N`` is Chris@87: the number of observations; this overrides the value implied by Chris@87: ``bias``. The default value is ``None``. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: The covariance matrix of the variables. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: corrcoef : Normalized covariance matrix Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: Consider two variables, :math:`x_0` and :math:`x_1`, which Chris@87: correlate perfectly, but in opposite directions: Chris@87: Chris@87: >>> x = np.array([[0, 2], [1, 1], [2, 0]]).T Chris@87: >>> x Chris@87: array([[0, 1, 2], Chris@87: [2, 1, 0]]) Chris@87: Chris@87: Note how :math:`x_0` increases while :math:`x_1` decreases. The covariance Chris@87: matrix shows this clearly: Chris@87: Chris@87: >>> np.cov(x) Chris@87: array([[ 1., -1.], Chris@87: [-1., 1.]]) Chris@87: Chris@87: Note that element :math:`C_{0,1}`, which shows the correlation between Chris@87: :math:`x_0` and :math:`x_1`, is negative. Chris@87: Chris@87: Further, note how `x` and `y` are combined: Chris@87: Chris@87: >>> x = [-2.1, -1, 4.3] Chris@87: >>> y = [3, 1.1, 0.12] Chris@87: >>> X = np.vstack((x,y)) Chris@87: >>> print np.cov(X) Chris@87: [[ 11.71 -4.286 ] Chris@87: [ -4.286 2.14413333]] Chris@87: >>> print np.cov(x, y) Chris@87: [[ 11.71 -4.286 ] Chris@87: [ -4.286 2.14413333]] Chris@87: >>> print np.cov(x) Chris@87: 11.71 Chris@87: Chris@87: """ Chris@87: # Check inputs Chris@87: if ddof is not None and ddof != int(ddof): Chris@87: raise ValueError( Chris@87: "ddof must be integer") Chris@87: Chris@87: # Handles complex arrays too Chris@87: m = np.asarray(m) Chris@87: if y is None: Chris@87: dtype = np.result_type(m, np.float64) Chris@87: else: Chris@87: y = np.asarray(y) Chris@87: dtype = np.result_type(m, y, np.float64) Chris@87: X = array(m, ndmin=2, dtype=dtype) Chris@87: Chris@87: if X.shape[0] == 1: Chris@87: rowvar = 1 Chris@87: if rowvar: Chris@87: N = X.shape[1] Chris@87: axis = 0 Chris@87: else: Chris@87: N = X.shape[0] Chris@87: axis = 1 Chris@87: Chris@87: # check ddof Chris@87: if ddof is None: Chris@87: if bias == 0: Chris@87: ddof = 1 Chris@87: else: Chris@87: ddof = 0 Chris@87: fact = float(N - ddof) Chris@87: if fact <= 0: Chris@87: warnings.warn("Degrees of freedom <= 0 for slice", RuntimeWarning) Chris@87: fact = 0.0 Chris@87: Chris@87: if y is not None: Chris@87: y = array(y, copy=False, ndmin=2, dtype=dtype) Chris@87: X = concatenate((X, y), axis) Chris@87: Chris@87: X -= X.mean(axis=1-axis, keepdims=True) Chris@87: if not rowvar: Chris@87: return (dot(X.T, X.conj()) / fact).squeeze() Chris@87: else: Chris@87: return (dot(X, X.T.conj()) / fact).squeeze() Chris@87: Chris@87: Chris@87: def corrcoef(x, y=None, rowvar=1, bias=0, ddof=None): Chris@87: """ Chris@87: Return correlation coefficients. Chris@87: Chris@87: Please refer to the documentation for `cov` for more detail. The Chris@87: relationship between the correlation coefficient matrix, `P`, and the Chris@87: covariance matrix, `C`, is Chris@87: Chris@87: .. math:: P_{ij} = \\frac{ C_{ij} } { \\sqrt{ C_{ii} * C_{jj} } } Chris@87: Chris@87: The values of `P` are between -1 and 1, inclusive. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x : array_like Chris@87: A 1-D or 2-D array containing multiple variables and observations. Chris@87: Each row of `m` represents a variable, and each column a single Chris@87: observation of all those variables. Also see `rowvar` below. Chris@87: y : array_like, optional Chris@87: An additional set of variables and observations. `y` has the same Chris@87: shape as `m`. Chris@87: rowvar : int, optional Chris@87: If `rowvar` is non-zero (default), then each row represents a Chris@87: variable, with observations in the columns. Otherwise, the relationship Chris@87: is transposed: each column represents a variable, while the rows Chris@87: contain observations. Chris@87: bias : int, optional Chris@87: Default normalization is by ``(N - 1)``, where ``N`` is the number of Chris@87: observations (unbiased estimate). If `bias` is 1, then Chris@87: normalization is by ``N``. These values can be overridden by using Chris@87: the keyword ``ddof`` in numpy versions >= 1.5. Chris@87: ddof : {None, int}, optional Chris@87: .. versionadded:: 1.5 Chris@87: If not ``None`` normalization is by ``(N - ddof)``, where ``N`` is Chris@87: the number of observations; this overrides the value implied by Chris@87: ``bias``. The default value is ``None``. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: The correlation coefficient matrix of the variables. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: cov : Covariance matrix Chris@87: Chris@87: """ Chris@87: c = cov(x, y, rowvar, bias, ddof) Chris@87: try: Chris@87: d = diag(c) Chris@87: except ValueError: # scalar covariance Chris@87: # nan if incorrect value (nan, inf, 0), 1 otherwise Chris@87: return c / c Chris@87: return c / sqrt(multiply.outer(d, d)) Chris@87: Chris@87: Chris@87: def blackman(M): Chris@87: """ Chris@87: Return the Blackman window. Chris@87: Chris@87: The Blackman window is a taper formed by using the first three Chris@87: terms of a summation of cosines. It was designed to have close to the Chris@87: minimal leakage possible. It is close to optimal, only slightly worse Chris@87: than a Kaiser window. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: M : int Chris@87: Number of points in the output window. If zero or less, an empty Chris@87: array is returned. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: The window, with the maximum value normalized to one (the value one Chris@87: appears only if the number of samples is odd). Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: bartlett, hamming, hanning, kaiser Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: The Blackman window is defined as Chris@87: Chris@87: .. math:: w(n) = 0.42 - 0.5 \\cos(2\\pi n/M) + 0.08 \\cos(4\\pi n/M) Chris@87: Chris@87: Most references to the Blackman window come from the signal processing Chris@87: literature, where it is used as one of many windowing functions for Chris@87: smoothing values. It is also known as an apodization (which means Chris@87: "removing the foot", i.e. smoothing discontinuities at the beginning Chris@87: and end of the sampled signal) or tapering function. It is known as a Chris@87: "near optimal" tapering function, almost as good (by some measures) Chris@87: as the kaiser window. Chris@87: Chris@87: References Chris@87: ---------- Chris@87: Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra, Chris@87: Dover Publications, New York. Chris@87: Chris@87: Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing. Chris@87: Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> np.blackman(12) Chris@87: array([ -1.38777878e-17, 3.26064346e-02, 1.59903635e-01, Chris@87: 4.14397981e-01, 7.36045180e-01, 9.67046769e-01, Chris@87: 9.67046769e-01, 7.36045180e-01, 4.14397981e-01, Chris@87: 1.59903635e-01, 3.26064346e-02, -1.38777878e-17]) Chris@87: Chris@87: Chris@87: Plot the window and the frequency response: Chris@87: Chris@87: >>> from numpy.fft import fft, fftshift Chris@87: >>> window = np.blackman(51) Chris@87: >>> plt.plot(window) Chris@87: [] Chris@87: >>> plt.title("Blackman window") Chris@87: Chris@87: >>> plt.ylabel("Amplitude") Chris@87: Chris@87: >>> plt.xlabel("Sample") Chris@87: Chris@87: >>> plt.show() Chris@87: Chris@87: >>> plt.figure() Chris@87: Chris@87: >>> A = fft(window, 2048) / 25.5 Chris@87: >>> mag = np.abs(fftshift(A)) Chris@87: >>> freq = np.linspace(-0.5, 0.5, len(A)) Chris@87: >>> response = 20 * np.log10(mag) Chris@87: >>> response = np.clip(response, -100, 100) Chris@87: >>> plt.plot(freq, response) Chris@87: [] Chris@87: >>> plt.title("Frequency response of Blackman window") Chris@87: Chris@87: >>> plt.ylabel("Magnitude [dB]") Chris@87: Chris@87: >>> plt.xlabel("Normalized frequency [cycles per sample]") Chris@87: Chris@87: >>> plt.axis('tight') Chris@87: (-0.5, 0.5, -100.0, ...) Chris@87: >>> plt.show() Chris@87: Chris@87: """ Chris@87: if M < 1: Chris@87: return array([]) Chris@87: if M == 1: Chris@87: return ones(1, float) Chris@87: n = arange(0, M) Chris@87: return 0.42 - 0.5*cos(2.0*pi*n/(M-1)) + 0.08*cos(4.0*pi*n/(M-1)) Chris@87: Chris@87: Chris@87: def bartlett(M): Chris@87: """ Chris@87: Return the Bartlett window. Chris@87: Chris@87: The Bartlett window is very similar to a triangular window, except Chris@87: that the end points are at zero. It is often used in signal Chris@87: processing for tapering a signal, without generating too much Chris@87: ripple in the frequency domain. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: M : int Chris@87: Number of points in the output window. If zero or less, an Chris@87: empty array is returned. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : array Chris@87: The triangular window, with the maximum value normalized to one Chris@87: (the value one appears only if the number of samples is odd), with Chris@87: the first and last samples equal to zero. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: blackman, hamming, hanning, kaiser Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: The Bartlett window is defined as Chris@87: Chris@87: .. math:: w(n) = \\frac{2}{M-1} \\left( Chris@87: \\frac{M-1}{2} - \\left|n - \\frac{M-1}{2}\\right| Chris@87: \\right) Chris@87: Chris@87: Most references to the Bartlett window come from the signal Chris@87: processing literature, where it is used as one of many windowing Chris@87: functions for smoothing values. Note that convolution with this Chris@87: window produces linear interpolation. It is also known as an Chris@87: apodization (which means"removing the foot", i.e. smoothing Chris@87: discontinuities at the beginning and end of the sampled signal) or Chris@87: tapering function. The fourier transform of the Bartlett is the product Chris@87: of two sinc functions. Chris@87: Note the excellent discussion in Kanasewich. Chris@87: Chris@87: References Chris@87: ---------- Chris@87: .. [1] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra", Chris@87: Biometrika 37, 1-16, 1950. Chris@87: .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", Chris@87: The University of Alberta Press, 1975, pp. 109-110. Chris@87: .. [3] A.V. Oppenheim and R.W. Schafer, "Discrete-Time Signal Chris@87: Processing", Prentice-Hall, 1999, pp. 468-471. Chris@87: .. [4] Wikipedia, "Window function", Chris@87: http://en.wikipedia.org/wiki/Window_function Chris@87: .. [5] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, Chris@87: "Numerical Recipes", Cambridge University Press, 1986, page 429. Chris@87: Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> np.bartlett(12) Chris@87: array([ 0. , 0.18181818, 0.36363636, 0.54545455, 0.72727273, Chris@87: 0.90909091, 0.90909091, 0.72727273, 0.54545455, 0.36363636, Chris@87: 0.18181818, 0. ]) Chris@87: Chris@87: Plot the window and its frequency response (requires SciPy and matplotlib): Chris@87: Chris@87: >>> from numpy.fft import fft, fftshift Chris@87: >>> window = np.bartlett(51) Chris@87: >>> plt.plot(window) Chris@87: [] Chris@87: >>> plt.title("Bartlett window") Chris@87: Chris@87: >>> plt.ylabel("Amplitude") Chris@87: Chris@87: >>> plt.xlabel("Sample") Chris@87: Chris@87: >>> plt.show() Chris@87: Chris@87: >>> plt.figure() Chris@87: Chris@87: >>> A = fft(window, 2048) / 25.5 Chris@87: >>> mag = np.abs(fftshift(A)) Chris@87: >>> freq = np.linspace(-0.5, 0.5, len(A)) Chris@87: >>> response = 20 * np.log10(mag) Chris@87: >>> response = np.clip(response, -100, 100) Chris@87: >>> plt.plot(freq, response) Chris@87: [] Chris@87: >>> plt.title("Frequency response of Bartlett window") Chris@87: Chris@87: >>> plt.ylabel("Magnitude [dB]") Chris@87: Chris@87: >>> plt.xlabel("Normalized frequency [cycles per sample]") Chris@87: Chris@87: >>> plt.axis('tight') Chris@87: (-0.5, 0.5, -100.0, ...) Chris@87: >>> plt.show() Chris@87: Chris@87: """ Chris@87: if M < 1: Chris@87: return array([]) Chris@87: if M == 1: Chris@87: return ones(1, float) Chris@87: n = arange(0, M) Chris@87: return where(less_equal(n, (M-1)/2.0), 2.0*n/(M-1), 2.0 - 2.0*n/(M-1)) Chris@87: Chris@87: Chris@87: def hanning(M): Chris@87: """ Chris@87: Return the Hanning window. Chris@87: Chris@87: The Hanning window is a taper formed by using a weighted cosine. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: M : int Chris@87: Number of points in the output window. If zero or less, an Chris@87: empty array is returned. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray, shape(M,) Chris@87: The window, with the maximum value normalized to one (the value Chris@87: one appears only if `M` is odd). Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: bartlett, blackman, hamming, kaiser Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: The Hanning window is defined as Chris@87: Chris@87: .. math:: w(n) = 0.5 - 0.5cos\\left(\\frac{2\\pi{n}}{M-1}\\right) Chris@87: \\qquad 0 \\leq n \\leq M-1 Chris@87: Chris@87: The Hanning was named for Julius van Hann, an Austrian meteorologist. Chris@87: It is also known as the Cosine Bell. Some authors prefer that it be Chris@87: called a Hann window, to help avoid confusion with the very similar Chris@87: Hamming window. Chris@87: Chris@87: Most references to the Hanning window come from the signal processing Chris@87: literature, where it is used as one of many windowing functions for Chris@87: smoothing values. It is also known as an apodization (which means Chris@87: "removing the foot", i.e. smoothing discontinuities at the beginning Chris@87: and end of the sampled signal) or tapering function. Chris@87: Chris@87: References Chris@87: ---------- Chris@87: .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power Chris@87: spectra, Dover Publications, New York. Chris@87: .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", Chris@87: The University of Alberta Press, 1975, pp. 106-108. Chris@87: .. [3] Wikipedia, "Window function", Chris@87: http://en.wikipedia.org/wiki/Window_function Chris@87: .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, Chris@87: "Numerical Recipes", Cambridge University Press, 1986, page 425. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> np.hanning(12) Chris@87: array([ 0. , 0.07937323, 0.29229249, 0.57115742, 0.82743037, Chris@87: 0.97974649, 0.97974649, 0.82743037, 0.57115742, 0.29229249, Chris@87: 0.07937323, 0. ]) Chris@87: Chris@87: Plot the window and its frequency response: Chris@87: Chris@87: >>> from numpy.fft import fft, fftshift Chris@87: >>> window = np.hanning(51) Chris@87: >>> plt.plot(window) Chris@87: [] Chris@87: >>> plt.title("Hann window") Chris@87: Chris@87: >>> plt.ylabel("Amplitude") Chris@87: Chris@87: >>> plt.xlabel("Sample") Chris@87: Chris@87: >>> plt.show() Chris@87: Chris@87: >>> plt.figure() Chris@87: Chris@87: >>> A = fft(window, 2048) / 25.5 Chris@87: >>> mag = np.abs(fftshift(A)) Chris@87: >>> freq = np.linspace(-0.5, 0.5, len(A)) Chris@87: >>> response = 20 * np.log10(mag) Chris@87: >>> response = np.clip(response, -100, 100) Chris@87: >>> plt.plot(freq, response) Chris@87: [] Chris@87: >>> plt.title("Frequency response of the Hann window") Chris@87: Chris@87: >>> plt.ylabel("Magnitude [dB]") Chris@87: Chris@87: >>> plt.xlabel("Normalized frequency [cycles per sample]") Chris@87: Chris@87: >>> plt.axis('tight') Chris@87: (-0.5, 0.5, -100.0, ...) Chris@87: >>> plt.show() Chris@87: Chris@87: """ Chris@87: if M < 1: Chris@87: return array([]) Chris@87: if M == 1: Chris@87: return ones(1, float) Chris@87: n = arange(0, M) Chris@87: return 0.5 - 0.5*cos(2.0*pi*n/(M-1)) Chris@87: Chris@87: Chris@87: def hamming(M): Chris@87: """ Chris@87: Return the Hamming window. Chris@87: Chris@87: The Hamming window is a taper formed by using a weighted cosine. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: M : int Chris@87: Number of points in the output window. If zero or less, an Chris@87: empty array is returned. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: The window, with the maximum value normalized to one (the value Chris@87: one appears only if the number of samples is odd). Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: bartlett, blackman, hanning, kaiser Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: The Hamming window is defined as Chris@87: Chris@87: .. math:: w(n) = 0.54 - 0.46cos\\left(\\frac{2\\pi{n}}{M-1}\\right) Chris@87: \\qquad 0 \\leq n \\leq M-1 Chris@87: Chris@87: The Hamming was named for R. W. Hamming, an associate of J. W. Tukey Chris@87: and is described in Blackman and Tukey. It was recommended for Chris@87: smoothing the truncated autocovariance function in the time domain. Chris@87: Most references to the Hamming window come from the signal processing Chris@87: literature, where it is used as one of many windowing functions for Chris@87: smoothing values. It is also known as an apodization (which means Chris@87: "removing the foot", i.e. smoothing discontinuities at the beginning Chris@87: and end of the sampled signal) or tapering function. Chris@87: Chris@87: References Chris@87: ---------- Chris@87: .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power Chris@87: spectra, Dover Publications, New York. Chris@87: .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The Chris@87: University of Alberta Press, 1975, pp. 109-110. Chris@87: .. [3] Wikipedia, "Window function", Chris@87: http://en.wikipedia.org/wiki/Window_function Chris@87: .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, Chris@87: "Numerical Recipes", Cambridge University Press, 1986, page 425. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> np.hamming(12) Chris@87: array([ 0.08 , 0.15302337, 0.34890909, 0.60546483, 0.84123594, Chris@87: 0.98136677, 0.98136677, 0.84123594, 0.60546483, 0.34890909, Chris@87: 0.15302337, 0.08 ]) Chris@87: Chris@87: Plot the window and the frequency response: Chris@87: Chris@87: >>> from numpy.fft import fft, fftshift Chris@87: >>> window = np.hamming(51) Chris@87: >>> plt.plot(window) Chris@87: [] Chris@87: >>> plt.title("Hamming window") Chris@87: Chris@87: >>> plt.ylabel("Amplitude") Chris@87: Chris@87: >>> plt.xlabel("Sample") Chris@87: Chris@87: >>> plt.show() Chris@87: Chris@87: >>> plt.figure() Chris@87: Chris@87: >>> A = fft(window, 2048) / 25.5 Chris@87: >>> mag = np.abs(fftshift(A)) Chris@87: >>> freq = np.linspace(-0.5, 0.5, len(A)) Chris@87: >>> response = 20 * np.log10(mag) Chris@87: >>> response = np.clip(response, -100, 100) Chris@87: >>> plt.plot(freq, response) Chris@87: [] Chris@87: >>> plt.title("Frequency response of Hamming window") Chris@87: Chris@87: >>> plt.ylabel("Magnitude [dB]") Chris@87: Chris@87: >>> plt.xlabel("Normalized frequency [cycles per sample]") Chris@87: Chris@87: >>> plt.axis('tight') Chris@87: (-0.5, 0.5, -100.0, ...) Chris@87: >>> plt.show() Chris@87: Chris@87: """ Chris@87: if M < 1: Chris@87: return array([]) Chris@87: if M == 1: Chris@87: return ones(1, float) Chris@87: n = arange(0, M) Chris@87: return 0.54 - 0.46*cos(2.0*pi*n/(M-1)) Chris@87: Chris@87: ## Code from cephes for i0 Chris@87: Chris@87: _i0A = [ Chris@87: -4.41534164647933937950E-18, Chris@87: 3.33079451882223809783E-17, Chris@87: -2.43127984654795469359E-16, Chris@87: 1.71539128555513303061E-15, Chris@87: -1.16853328779934516808E-14, Chris@87: 7.67618549860493561688E-14, Chris@87: -4.85644678311192946090E-13, Chris@87: 2.95505266312963983461E-12, Chris@87: -1.72682629144155570723E-11, Chris@87: 9.67580903537323691224E-11, Chris@87: -5.18979560163526290666E-10, Chris@87: 2.65982372468238665035E-9, Chris@87: -1.30002500998624804212E-8, Chris@87: 6.04699502254191894932E-8, Chris@87: -2.67079385394061173391E-7, Chris@87: 1.11738753912010371815E-6, Chris@87: -4.41673835845875056359E-6, Chris@87: 1.64484480707288970893E-5, Chris@87: -5.75419501008210370398E-5, Chris@87: 1.88502885095841655729E-4, Chris@87: -5.76375574538582365885E-4, Chris@87: 1.63947561694133579842E-3, Chris@87: -4.32430999505057594430E-3, Chris@87: 1.05464603945949983183E-2, Chris@87: -2.37374148058994688156E-2, Chris@87: 4.93052842396707084878E-2, Chris@87: -9.49010970480476444210E-2, Chris@87: 1.71620901522208775349E-1, Chris@87: -3.04682672343198398683E-1, Chris@87: 6.76795274409476084995E-1 Chris@87: ] Chris@87: Chris@87: _i0B = [ Chris@87: -7.23318048787475395456E-18, Chris@87: -4.83050448594418207126E-18, Chris@87: 4.46562142029675999901E-17, Chris@87: 3.46122286769746109310E-17, Chris@87: -2.82762398051658348494E-16, Chris@87: -3.42548561967721913462E-16, Chris@87: 1.77256013305652638360E-15, Chris@87: 3.81168066935262242075E-15, Chris@87: -9.55484669882830764870E-15, Chris@87: -4.15056934728722208663E-14, Chris@87: 1.54008621752140982691E-14, Chris@87: 3.85277838274214270114E-13, Chris@87: 7.18012445138366623367E-13, Chris@87: -1.79417853150680611778E-12, Chris@87: -1.32158118404477131188E-11, Chris@87: -3.14991652796324136454E-11, Chris@87: 1.18891471078464383424E-11, Chris@87: 4.94060238822496958910E-10, Chris@87: 3.39623202570838634515E-9, Chris@87: 2.26666899049817806459E-8, Chris@87: 2.04891858946906374183E-7, Chris@87: 2.89137052083475648297E-6, Chris@87: 6.88975834691682398426E-5, Chris@87: 3.36911647825569408990E-3, Chris@87: 8.04490411014108831608E-1 Chris@87: ] Chris@87: Chris@87: Chris@87: def _chbevl(x, vals): Chris@87: b0 = vals[0] Chris@87: b1 = 0.0 Chris@87: Chris@87: for i in range(1, len(vals)): Chris@87: b2 = b1 Chris@87: b1 = b0 Chris@87: b0 = x*b1 - b2 + vals[i] Chris@87: Chris@87: return 0.5*(b0 - b2) Chris@87: Chris@87: Chris@87: def _i0_1(x): Chris@87: return exp(x) * _chbevl(x/2.0-2, _i0A) Chris@87: Chris@87: Chris@87: def _i0_2(x): Chris@87: return exp(x) * _chbevl(32.0/x - 2.0, _i0B) / sqrt(x) Chris@87: Chris@87: Chris@87: def i0(x): Chris@87: """ Chris@87: Modified Bessel function of the first kind, order 0. Chris@87: Chris@87: Usually denoted :math:`I_0`. This function does broadcast, but will *not* Chris@87: "up-cast" int dtype arguments unless accompanied by at least one float or Chris@87: complex dtype argument (see Raises below). Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x : array_like, dtype float or complex Chris@87: Argument of the Bessel function. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray, shape = x.shape, dtype = x.dtype Chris@87: The modified Bessel function evaluated at each of the elements of `x`. Chris@87: Chris@87: Raises Chris@87: ------ Chris@87: TypeError: array cannot be safely cast to required type Chris@87: If argument consists exclusively of int dtypes. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: scipy.special.iv, scipy.special.ive Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: We use the algorithm published by Clenshaw [1]_ and referenced by Chris@87: Abramowitz and Stegun [2]_, for which the function domain is Chris@87: partitioned into the two intervals [0,8] and (8,inf), and Chebyshev Chris@87: polynomial expansions are employed in each interval. Relative error on Chris@87: the domain [0,30] using IEEE arithmetic is documented [3]_ as having a Chris@87: peak of 5.8e-16 with an rms of 1.4e-16 (n = 30000). Chris@87: Chris@87: References Chris@87: ---------- Chris@87: .. [1] C. W. Clenshaw, "Chebyshev series for mathematical functions", in Chris@87: *National Physical Laboratory Mathematical Tables*, vol. 5, London: Chris@87: Her Majesty's Stationery Office, 1962. Chris@87: .. [2] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical Chris@87: Functions*, 10th printing, New York: Dover, 1964, pp. 379. Chris@87: http://www.math.sfu.ca/~cbm/aands/page_379.htm Chris@87: .. [3] http://kobesearch.cpan.org/htdocs/Math-Cephes/Math/Cephes.html Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> np.i0([0.]) Chris@87: array(1.0) Chris@87: >>> np.i0([0., 1. + 2j]) Chris@87: array([ 1.00000000+0.j , 0.18785373+0.64616944j]) Chris@87: Chris@87: """ Chris@87: x = atleast_1d(x).copy() Chris@87: y = empty_like(x) Chris@87: ind = (x < 0) Chris@87: x[ind] = -x[ind] Chris@87: ind = (x <= 8.0) Chris@87: y[ind] = _i0_1(x[ind]) Chris@87: ind2 = ~ind Chris@87: y[ind2] = _i0_2(x[ind2]) Chris@87: return y.squeeze() Chris@87: Chris@87: ## End of cephes code for i0 Chris@87: Chris@87: Chris@87: def kaiser(M, beta): Chris@87: """ Chris@87: Return the Kaiser window. Chris@87: Chris@87: The Kaiser window is a taper formed by using a Bessel function. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: M : int Chris@87: Number of points in the output window. If zero or less, an Chris@87: empty array is returned. Chris@87: beta : float Chris@87: Shape parameter for window. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : array Chris@87: The window, with the maximum value normalized to one (the value Chris@87: one appears only if the number of samples is odd). Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: bartlett, blackman, hamming, hanning Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: The Kaiser window is defined as Chris@87: Chris@87: .. math:: w(n) = I_0\\left( \\beta \\sqrt{1-\\frac{4n^2}{(M-1)^2}} Chris@87: \\right)/I_0(\\beta) Chris@87: Chris@87: with Chris@87: Chris@87: .. math:: \\quad -\\frac{M-1}{2} \\leq n \\leq \\frac{M-1}{2}, Chris@87: Chris@87: where :math:`I_0` is the modified zeroth-order Bessel function. Chris@87: Chris@87: The Kaiser was named for Jim Kaiser, who discovered a simple Chris@87: approximation to the DPSS window based on Bessel functions. The Kaiser Chris@87: window is a very good approximation to the Digital Prolate Spheroidal Chris@87: Sequence, or Slepian window, which is the transform which maximizes the Chris@87: energy in the main lobe of the window relative to total energy. Chris@87: Chris@87: The Kaiser can approximate many other windows by varying the beta Chris@87: parameter. Chris@87: Chris@87: ==== ======================= Chris@87: beta Window shape Chris@87: ==== ======================= Chris@87: 0 Rectangular Chris@87: 5 Similar to a Hamming Chris@87: 6 Similar to a Hanning Chris@87: 8.6 Similar to a Blackman Chris@87: ==== ======================= Chris@87: Chris@87: A beta value of 14 is probably a good starting point. Note that as beta Chris@87: gets large, the window narrows, and so the number of samples needs to be Chris@87: large enough to sample the increasingly narrow spike, otherwise NaNs will Chris@87: get returned. Chris@87: Chris@87: Most references to the Kaiser window come from the signal processing Chris@87: literature, where it is used as one of many windowing functions for Chris@87: smoothing values. It is also known as an apodization (which means Chris@87: "removing the foot", i.e. smoothing discontinuities at the beginning Chris@87: and end of the sampled signal) or tapering function. Chris@87: Chris@87: References Chris@87: ---------- Chris@87: .. [1] J. F. Kaiser, "Digital Filters" - Ch 7 in "Systems analysis by Chris@87: digital computer", Editors: F.F. Kuo and J.F. Kaiser, p 218-285. Chris@87: John Wiley and Sons, New York, (1966). Chris@87: .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The Chris@87: University of Alberta Press, 1975, pp. 177-178. Chris@87: .. [3] Wikipedia, "Window function", Chris@87: http://en.wikipedia.org/wiki/Window_function Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> np.kaiser(12, 14) Chris@87: array([ 7.72686684e-06, 3.46009194e-03, 4.65200189e-02, Chris@87: 2.29737120e-01, 5.99885316e-01, 9.45674898e-01, Chris@87: 9.45674898e-01, 5.99885316e-01, 2.29737120e-01, Chris@87: 4.65200189e-02, 3.46009194e-03, 7.72686684e-06]) Chris@87: Chris@87: Chris@87: Plot the window and the frequency response: Chris@87: Chris@87: >>> from numpy.fft import fft, fftshift Chris@87: >>> window = np.kaiser(51, 14) Chris@87: >>> plt.plot(window) Chris@87: [] Chris@87: >>> plt.title("Kaiser window") Chris@87: Chris@87: >>> plt.ylabel("Amplitude") Chris@87: Chris@87: >>> plt.xlabel("Sample") Chris@87: Chris@87: >>> plt.show() Chris@87: Chris@87: >>> plt.figure() Chris@87: Chris@87: >>> A = fft(window, 2048) / 25.5 Chris@87: >>> mag = np.abs(fftshift(A)) Chris@87: >>> freq = np.linspace(-0.5, 0.5, len(A)) Chris@87: >>> response = 20 * np.log10(mag) Chris@87: >>> response = np.clip(response, -100, 100) Chris@87: >>> plt.plot(freq, response) Chris@87: [] Chris@87: >>> plt.title("Frequency response of Kaiser window") Chris@87: Chris@87: >>> plt.ylabel("Magnitude [dB]") Chris@87: Chris@87: >>> plt.xlabel("Normalized frequency [cycles per sample]") Chris@87: Chris@87: >>> plt.axis('tight') Chris@87: (-0.5, 0.5, -100.0, ...) Chris@87: >>> plt.show() Chris@87: Chris@87: """ Chris@87: from numpy.dual import i0 Chris@87: if M == 1: Chris@87: return np.array([1.]) Chris@87: n = arange(0, M) Chris@87: alpha = (M-1)/2.0 Chris@87: return i0(beta * sqrt(1-((n-alpha)/alpha)**2.0))/i0(float(beta)) Chris@87: Chris@87: Chris@87: def sinc(x): Chris@87: """ Chris@87: Return the sinc function. Chris@87: Chris@87: The sinc function is :math:`\\sin(\\pi x)/(\\pi x)`. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x : ndarray Chris@87: Array (possibly multi-dimensional) of values for which to to Chris@87: calculate ``sinc(x)``. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: ``sinc(x)``, which has the same shape as the input. Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: ``sinc(0)`` is the limit value 1. Chris@87: Chris@87: The name sinc is short for "sine cardinal" or "sinus cardinalis". Chris@87: Chris@87: The sinc function is used in various signal processing applications, Chris@87: including in anti-aliasing, in the construction of a Lanczos resampling Chris@87: filter, and in interpolation. Chris@87: Chris@87: For bandlimited interpolation of discrete-time signals, the ideal Chris@87: interpolation kernel is proportional to the sinc function. Chris@87: Chris@87: References Chris@87: ---------- Chris@87: .. [1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram Web Chris@87: Resource. http://mathworld.wolfram.com/SincFunction.html Chris@87: .. [2] Wikipedia, "Sinc function", Chris@87: http://en.wikipedia.org/wiki/Sinc_function Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> x = np.linspace(-4, 4, 41) Chris@87: >>> np.sinc(x) Chris@87: array([ -3.89804309e-17, -4.92362781e-02, -8.40918587e-02, Chris@87: -8.90384387e-02, -5.84680802e-02, 3.89804309e-17, Chris@87: 6.68206631e-02, 1.16434881e-01, 1.26137788e-01, Chris@87: 8.50444803e-02, -3.89804309e-17, -1.03943254e-01, Chris@87: -1.89206682e-01, -2.16236208e-01, -1.55914881e-01, Chris@87: 3.89804309e-17, 2.33872321e-01, 5.04551152e-01, Chris@87: 7.56826729e-01, 9.35489284e-01, 1.00000000e+00, Chris@87: 9.35489284e-01, 7.56826729e-01, 5.04551152e-01, Chris@87: 2.33872321e-01, 3.89804309e-17, -1.55914881e-01, Chris@87: -2.16236208e-01, -1.89206682e-01, -1.03943254e-01, Chris@87: -3.89804309e-17, 8.50444803e-02, 1.26137788e-01, Chris@87: 1.16434881e-01, 6.68206631e-02, 3.89804309e-17, Chris@87: -5.84680802e-02, -8.90384387e-02, -8.40918587e-02, Chris@87: -4.92362781e-02, -3.89804309e-17]) Chris@87: Chris@87: >>> plt.plot(x, np.sinc(x)) Chris@87: [] Chris@87: >>> plt.title("Sinc Function") Chris@87: Chris@87: >>> plt.ylabel("Amplitude") Chris@87: Chris@87: >>> plt.xlabel("X") Chris@87: Chris@87: >>> plt.show() Chris@87: Chris@87: It works in 2-D as well: Chris@87: Chris@87: >>> x = np.linspace(-4, 4, 401) Chris@87: >>> xx = np.outer(x, x) Chris@87: >>> plt.imshow(np.sinc(xx)) Chris@87: Chris@87: Chris@87: """ Chris@87: x = np.asanyarray(x) Chris@87: y = pi * where(x == 0, 1.0e-20, x) Chris@87: return sin(y)/y Chris@87: Chris@87: Chris@87: def msort(a): Chris@87: """ Chris@87: Return a copy of an array sorted along the first axis. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: a : array_like Chris@87: Array to be sorted. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: sorted_array : ndarray Chris@87: Array of the same type and shape as `a`. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: sort Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: ``np.msort(a)`` is equivalent to ``np.sort(a, axis=0)``. Chris@87: Chris@87: """ Chris@87: b = array(a, subok=True, copy=True) Chris@87: b.sort(0) Chris@87: return b Chris@87: Chris@87: Chris@87: def _ureduce(a, func, **kwargs): Chris@87: """ Chris@87: Internal Function. Chris@87: Call `func` with `a` as first argument swapping the axes to use extended Chris@87: axis on functions that don't support it natively. Chris@87: Chris@87: Returns result and a.shape with axis dims set to 1. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: a : array_like Chris@87: Input array or object that can be converted to an array. Chris@87: func : callable Chris@87: Reduction function Kapable of receiving an axis argument. Chris@87: It is is called with `a` as first argument followed by `kwargs`. Chris@87: kwargs : keyword arguments Chris@87: additional keyword arguments to pass to `func`. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: result : tuple Chris@87: Result of func(a, **kwargs) and a.shape with axis dims set to 1 Chris@87: which can be used to reshape the result to the same shape a ufunc with Chris@87: keepdims=True would produce. Chris@87: Chris@87: """ Chris@87: a = np.asanyarray(a) Chris@87: axis = kwargs.get('axis', None) Chris@87: if axis is not None: Chris@87: keepdim = list(a.shape) Chris@87: nd = a.ndim Chris@87: try: Chris@87: axis = operator.index(axis) Chris@87: if axis >= nd or axis < -nd: Chris@87: raise IndexError("axis %d out of bounds (%d)" % (axis, a.ndim)) Chris@87: keepdim[axis] = 1 Chris@87: except TypeError: Chris@87: sax = set() Chris@87: for x in axis: Chris@87: if x >= nd or x < -nd: Chris@87: raise IndexError("axis %d out of bounds (%d)" % (x, nd)) Chris@87: if x in sax: Chris@87: raise ValueError("duplicate value in axis") Chris@87: sax.add(x % nd) Chris@87: keepdim[x] = 1 Chris@87: keep = sax.symmetric_difference(frozenset(range(nd))) Chris@87: nkeep = len(keep) Chris@87: # swap axis that should not be reduced to front Chris@87: for i, s in enumerate(sorted(keep)): Chris@87: a = a.swapaxes(i, s) Chris@87: # merge reduced axis Chris@87: a = a.reshape(a.shape[:nkeep] + (-1,)) Chris@87: kwargs['axis'] = -1 Chris@87: else: Chris@87: keepdim = [1] * a.ndim Chris@87: Chris@87: r = func(a, **kwargs) Chris@87: return r, keepdim Chris@87: Chris@87: Chris@87: def median(a, axis=None, out=None, overwrite_input=False, keepdims=False): Chris@87: """ Chris@87: Compute the median along the specified axis. Chris@87: Chris@87: Returns the median of the array elements. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: a : array_like Chris@87: Input array or object that can be converted to an array. Chris@87: axis : int or sequence of int, optional Chris@87: Axis along which the medians are computed. The default (axis=None) Chris@87: is to compute the median along a flattened version of the array. Chris@87: A sequence of axes is supported since version 1.9.0. Chris@87: out : ndarray, optional Chris@87: Alternative output array in which to place the result. It must have Chris@87: the same shape and buffer length as the expected output, but the Chris@87: type (of the output) will be cast if necessary. Chris@87: overwrite_input : bool, optional Chris@87: If True, then allow use of memory of input array (a) for Chris@87: calculations. The input array will be modified by the call to Chris@87: median. This will save memory when you do not need to preserve the Chris@87: contents of the input array. Treat the input as undefined, but it Chris@87: will probably be fully or partially sorted. Default is False. Note Chris@87: that, if `overwrite_input` is True and the input is not already an Chris@87: ndarray, an error will be raised. Chris@87: keepdims : bool, optional Chris@87: If this is set to True, the axes which are reduced are left Chris@87: in the result as dimensions with size one. With this option, Chris@87: the result will broadcast correctly against the original `arr`. Chris@87: Chris@87: .. versionadded:: 1.9.0 Chris@87: Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: median : ndarray Chris@87: A new array holding the result (unless `out` is specified, in which Chris@87: case that array is returned instead). If the input contains Chris@87: integers, or floats of smaller precision than 64, then the output Chris@87: data-type is float64. Otherwise, the output data-type is the same Chris@87: as that of the input. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: mean, percentile Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: Given a vector V of length N, the median of V is the middle value of Chris@87: a sorted copy of V, ``V_sorted`` - i.e., ``V_sorted[(N-1)/2]``, when N is Chris@87: odd. When N is even, it is the average of the two middle values of Chris@87: ``V_sorted``. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> a = np.array([[10, 7, 4], [3, 2, 1]]) Chris@87: >>> a Chris@87: array([[10, 7, 4], Chris@87: [ 3, 2, 1]]) Chris@87: >>> np.median(a) Chris@87: 3.5 Chris@87: >>> np.median(a, axis=0) Chris@87: array([ 6.5, 4.5, 2.5]) Chris@87: >>> np.median(a, axis=1) Chris@87: array([ 7., 2.]) Chris@87: >>> m = np.median(a, axis=0) Chris@87: >>> out = np.zeros_like(m) Chris@87: >>> np.median(a, axis=0, out=m) Chris@87: array([ 6.5, 4.5, 2.5]) Chris@87: >>> m Chris@87: array([ 6.5, 4.5, 2.5]) Chris@87: >>> b = a.copy() Chris@87: >>> np.median(b, axis=1, overwrite_input=True) Chris@87: array([ 7., 2.]) Chris@87: >>> assert not np.all(a==b) Chris@87: >>> b = a.copy() Chris@87: >>> np.median(b, axis=None, overwrite_input=True) Chris@87: 3.5 Chris@87: >>> assert not np.all(a==b) Chris@87: Chris@87: """ Chris@87: r, k = _ureduce(a, func=_median, axis=axis, out=out, Chris@87: overwrite_input=overwrite_input) Chris@87: if keepdims: Chris@87: return r.reshape(k) Chris@87: else: Chris@87: return r Chris@87: Chris@87: def _median(a, axis=None, out=None, overwrite_input=False): Chris@87: # can't be reasonably be implemented in terms of percentile as we have to Chris@87: # call mean to not break astropy Chris@87: a = np.asanyarray(a) Chris@87: if axis is not None and axis >= a.ndim: Chris@87: raise IndexError( Chris@87: "axis %d out of bounds (%d)" % (axis, a.ndim)) Chris@87: Chris@87: if overwrite_input: Chris@87: if axis is None: Chris@87: part = a.ravel() Chris@87: sz = part.size Chris@87: if sz % 2 == 0: Chris@87: szh = sz // 2 Chris@87: part.partition((szh - 1, szh)) Chris@87: else: Chris@87: part.partition((sz - 1) // 2) Chris@87: else: Chris@87: sz = a.shape[axis] Chris@87: if sz % 2 == 0: Chris@87: szh = sz // 2 Chris@87: a.partition((szh - 1, szh), axis=axis) Chris@87: else: Chris@87: a.partition((sz - 1) // 2, axis=axis) Chris@87: part = a Chris@87: else: Chris@87: if axis is None: Chris@87: sz = a.size Chris@87: else: Chris@87: sz = a.shape[axis] Chris@87: if sz % 2 == 0: Chris@87: part = partition(a, ((sz // 2) - 1, sz // 2), axis=axis) Chris@87: else: Chris@87: part = partition(a, (sz - 1) // 2, axis=axis) Chris@87: if part.shape == (): Chris@87: # make 0-D arrays work Chris@87: return part.item() Chris@87: if axis is None: Chris@87: axis = 0 Chris@87: indexer = [slice(None)] * part.ndim Chris@87: index = part.shape[axis] // 2 Chris@87: if part.shape[axis] % 2 == 1: Chris@87: # index with slice to allow mean (below) to work Chris@87: indexer[axis] = slice(index, index+1) Chris@87: else: Chris@87: indexer[axis] = slice(index-1, index+1) Chris@87: # Use mean in odd and even case to coerce data type Chris@87: # and check, use out array. Chris@87: return mean(part[indexer], axis=axis, out=out) Chris@87: Chris@87: Chris@87: def percentile(a, q, axis=None, out=None, Chris@87: overwrite_input=False, interpolation='linear', keepdims=False): Chris@87: """ Chris@87: Compute the qth percentile of the data along the specified axis. Chris@87: Chris@87: Returns the qth percentile of the array elements. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: a : array_like Chris@87: Input array or object that can be converted to an array. Chris@87: q : float in range of [0,100] (or sequence of floats) Chris@87: Percentile to compute which must be between 0 and 100 inclusive. Chris@87: axis : int or sequence of int, optional Chris@87: Axis along which the percentiles are computed. The default (None) Chris@87: is to compute the percentiles along a flattened version of the array. Chris@87: A sequence of axes is supported since version 1.9.0. Chris@87: out : ndarray, optional Chris@87: Alternative output array in which to place the result. It must Chris@87: have the same shape and buffer length as the expected output, Chris@87: but the type (of the output) will be cast if necessary. Chris@87: overwrite_input : bool, optional Chris@87: If True, then allow use of memory of input array `a` for Chris@87: calculations. The input array will be modified by the call to Chris@87: percentile. This will save memory when you do not need to preserve Chris@87: the contents of the input array. In this case you should not make Chris@87: any assumptions about the content of the passed in array `a` after Chris@87: this function completes -- treat it as undefined. Default is False. Chris@87: Note that, if the `a` input is not already an array this parameter Chris@87: will have no effect, `a` will be converted to an array internally Chris@87: regardless of the value of this parameter. Chris@87: interpolation : {'linear', 'lower', 'higher', 'midpoint', 'nearest'} Chris@87: This optional parameter specifies the interpolation method to use, Chris@87: when the desired quantile lies between two data points `i` and `j`: Chris@87: * linear: `i + (j - i) * fraction`, where `fraction` is the Chris@87: fractional part of the index surrounded by `i` and `j`. Chris@87: * lower: `i`. Chris@87: * higher: `j`. Chris@87: * nearest: `i` or `j` whichever is nearest. Chris@87: * midpoint: (`i` + `j`) / 2. Chris@87: Chris@87: .. versionadded:: 1.9.0 Chris@87: keepdims : bool, optional Chris@87: If this is set to True, the axes which are reduced are left Chris@87: in the result as dimensions with size one. With this option, Chris@87: the result will broadcast correctly against the original `arr`. Chris@87: Chris@87: .. versionadded:: 1.9.0 Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: percentile : scalar or ndarray Chris@87: If a single percentile `q` is given and axis=None a scalar is Chris@87: returned. If multiple percentiles `q` are given an array holding Chris@87: the result is returned. The results are listed in the first axis. Chris@87: (If `out` is specified, in which case that array is returned Chris@87: instead). If the input contains integers, or floats of smaller Chris@87: precision than 64, then the output data-type is float64. Otherwise, Chris@87: the output data-type is the same as that of the input. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: mean, median Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: Given a vector V of length N, the q-th percentile of V is the q-th ranked Chris@87: value in a sorted copy of V. The values and distances of the two Chris@87: nearest neighbors as well as the `interpolation` parameter will Chris@87: determine the percentile if the normalized ranking does not match q Chris@87: exactly. This function is the same as the median if ``q=50``, the same Chris@87: as the minimum if ``q=0`` and the same as the maximum if ``q=100``. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> a = np.array([[10, 7, 4], [3, 2, 1]]) Chris@87: >>> a Chris@87: array([[10, 7, 4], Chris@87: [ 3, 2, 1]]) Chris@87: >>> np.percentile(a, 50) Chris@87: array([ 3.5]) Chris@87: >>> np.percentile(a, 50, axis=0) Chris@87: array([[ 6.5, 4.5, 2.5]]) Chris@87: >>> np.percentile(a, 50, axis=1) Chris@87: array([[ 7.], Chris@87: [ 2.]]) Chris@87: Chris@87: >>> m = np.percentile(a, 50, axis=0) Chris@87: >>> out = np.zeros_like(m) Chris@87: >>> np.percentile(a, 50, axis=0, out=m) Chris@87: array([[ 6.5, 4.5, 2.5]]) Chris@87: >>> m Chris@87: array([[ 6.5, 4.5, 2.5]]) Chris@87: Chris@87: >>> b = a.copy() Chris@87: >>> np.percentile(b, 50, axis=1, overwrite_input=True) Chris@87: array([[ 7.], Chris@87: [ 2.]]) Chris@87: >>> assert not np.all(a==b) Chris@87: >>> b = a.copy() Chris@87: >>> np.percentile(b, 50, axis=None, overwrite_input=True) Chris@87: array([ 3.5]) Chris@87: Chris@87: """ Chris@87: q = array(q, dtype=np.float64, copy=True) Chris@87: r, k = _ureduce(a, func=_percentile, q=q, axis=axis, out=out, Chris@87: overwrite_input=overwrite_input, Chris@87: interpolation=interpolation) Chris@87: if keepdims: Chris@87: if q.ndim == 0: Chris@87: return r.reshape(k) Chris@87: else: Chris@87: return r.reshape([len(q)] + k) Chris@87: else: Chris@87: return r Chris@87: Chris@87: Chris@87: def _percentile(a, q, axis=None, out=None, Chris@87: overwrite_input=False, interpolation='linear', keepdims=False): Chris@87: a = asarray(a) Chris@87: if q.ndim == 0: Chris@87: # Do not allow 0-d arrays because following code fails for scalar Chris@87: zerod = True Chris@87: q = q[None] Chris@87: else: Chris@87: zerod = False Chris@87: Chris@87: # avoid expensive reductions, relevant for arrays with < O(1000) elements Chris@87: if q.size < 10: Chris@87: for i in range(q.size): Chris@87: if q[i] < 0. or q[i] > 100.: Chris@87: raise ValueError("Percentiles must be in the range [0,100]") Chris@87: q[i] /= 100. Chris@87: else: Chris@87: # faster than any() Chris@87: if np.count_nonzero(q < 0.) or np.count_nonzero(q > 100.): Chris@87: raise ValueError("Percentiles must be in the range [0,100]") Chris@87: q /= 100. Chris@87: Chris@87: # prepare a for partioning Chris@87: if overwrite_input: Chris@87: if axis is None: Chris@87: ap = a.ravel() Chris@87: else: Chris@87: ap = a Chris@87: else: Chris@87: if axis is None: Chris@87: ap = a.flatten() Chris@87: else: Chris@87: ap = a.copy() Chris@87: Chris@87: if axis is None: Chris@87: axis = 0 Chris@87: Chris@87: Nx = ap.shape[axis] Chris@87: indices = q * (Nx - 1) Chris@87: Chris@87: # round fractional indices according to interpolation method Chris@87: if interpolation == 'lower': Chris@87: indices = floor(indices).astype(intp) Chris@87: elif interpolation == 'higher': Chris@87: indices = ceil(indices).astype(intp) Chris@87: elif interpolation == 'midpoint': Chris@87: indices = floor(indices) + 0.5 Chris@87: elif interpolation == 'nearest': Chris@87: indices = around(indices).astype(intp) Chris@87: elif interpolation == 'linear': Chris@87: pass # keep index as fraction and interpolate Chris@87: else: Chris@87: raise ValueError( Chris@87: "interpolation can only be 'linear', 'lower' 'higher', " Chris@87: "'midpoint', or 'nearest'") Chris@87: Chris@87: if indices.dtype == intp: # take the points along axis Chris@87: ap.partition(indices, axis=axis) Chris@87: # ensure axis with qth is first Chris@87: ap = np.rollaxis(ap, axis, 0) Chris@87: axis = 0 Chris@87: Chris@87: if zerod: Chris@87: indices = indices[0] Chris@87: r = take(ap, indices, axis=axis, out=out) Chris@87: else: # weight the points above and below the indices Chris@87: indices_below = floor(indices).astype(intp) Chris@87: indices_above = indices_below + 1 Chris@87: indices_above[indices_above > Nx - 1] = Nx - 1 Chris@87: Chris@87: weights_above = indices - indices_below Chris@87: weights_below = 1.0 - weights_above Chris@87: Chris@87: weights_shape = [1, ] * ap.ndim Chris@87: weights_shape[axis] = len(indices) Chris@87: weights_below.shape = weights_shape Chris@87: weights_above.shape = weights_shape Chris@87: Chris@87: ap.partition(concatenate((indices_below, indices_above)), axis=axis) Chris@87: x1 = take(ap, indices_below, axis=axis) * weights_below Chris@87: x2 = take(ap, indices_above, axis=axis) * weights_above Chris@87: Chris@87: # ensure axis with qth is first Chris@87: x1 = np.rollaxis(x1, axis, 0) Chris@87: x2 = np.rollaxis(x2, axis, 0) Chris@87: Chris@87: if zerod: Chris@87: x1 = x1.squeeze(0) Chris@87: x2 = x2.squeeze(0) Chris@87: Chris@87: if out is not None: Chris@87: r = add(x1, x2, out=out) Chris@87: else: Chris@87: r = add(x1, x2) Chris@87: Chris@87: return r Chris@87: Chris@87: Chris@87: def trapz(y, x=None, dx=1.0, axis=-1): Chris@87: """ Chris@87: Integrate along the given axis using the composite trapezoidal rule. Chris@87: Chris@87: Integrate `y` (`x`) along given axis. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: y : array_like Chris@87: Input array to integrate. Chris@87: x : array_like, optional Chris@87: If `x` is None, then spacing between all `y` elements is `dx`. Chris@87: dx : scalar, optional Chris@87: If `x` is None, spacing given by `dx` is assumed. Default is 1. Chris@87: axis : int, optional Chris@87: Specify the axis. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: trapz : float Chris@87: Definite integral as approximated by trapezoidal rule. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: sum, cumsum Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: Image [2]_ illustrates trapezoidal rule -- y-axis locations of points Chris@87: will be taken from `y` array, by default x-axis distances between Chris@87: points will be 1.0, alternatively they can be provided with `x` array Chris@87: or with `dx` scalar. Return value will be equal to combined area under Chris@87: the red lines. Chris@87: Chris@87: Chris@87: References Chris@87: ---------- Chris@87: .. [1] Wikipedia page: http://en.wikipedia.org/wiki/Trapezoidal_rule Chris@87: Chris@87: .. [2] Illustration image: Chris@87: http://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> np.trapz([1,2,3]) Chris@87: 4.0 Chris@87: >>> np.trapz([1,2,3], x=[4,6,8]) Chris@87: 8.0 Chris@87: >>> np.trapz([1,2,3], dx=2) Chris@87: 8.0 Chris@87: >>> a = np.arange(6).reshape(2, 3) Chris@87: >>> a Chris@87: array([[0, 1, 2], Chris@87: [3, 4, 5]]) Chris@87: >>> np.trapz(a, axis=0) Chris@87: array([ 1.5, 2.5, 3.5]) Chris@87: >>> np.trapz(a, axis=1) Chris@87: array([ 2., 8.]) Chris@87: Chris@87: """ Chris@87: y = asanyarray(y) Chris@87: if x is None: Chris@87: d = dx Chris@87: else: Chris@87: x = asanyarray(x) Chris@87: if x.ndim == 1: Chris@87: d = diff(x) Chris@87: # reshape to correct shape Chris@87: shape = [1]*y.ndim Chris@87: shape[axis] = d.shape[0] Chris@87: d = d.reshape(shape) Chris@87: else: Chris@87: d = diff(x, axis=axis) Chris@87: nd = len(y.shape) Chris@87: slice1 = [slice(None)]*nd Chris@87: slice2 = [slice(None)]*nd Chris@87: slice1[axis] = slice(1, None) Chris@87: slice2[axis] = slice(None, -1) Chris@87: try: Chris@87: ret = (d * (y[slice1] + y[slice2]) / 2.0).sum(axis) Chris@87: except ValueError: Chris@87: # Operations didn't work, cast to ndarray Chris@87: d = np.asarray(d) Chris@87: y = np.asarray(y) Chris@87: ret = add.reduce(d * (y[slice1]+y[slice2])/2.0, axis) Chris@87: return ret Chris@87: Chris@87: Chris@87: #always succeed Chris@87: def add_newdoc(place, obj, doc): Chris@87: """Adds documentation to obj which is in module place. Chris@87: Chris@87: If doc is a string add it to obj as a docstring Chris@87: Chris@87: If doc is a tuple, then the first element is interpreted as Chris@87: an attribute of obj and the second as the docstring Chris@87: (method, docstring) Chris@87: Chris@87: If doc is a list, then each element of the list should be a Chris@87: sequence of length two --> [(method1, docstring1), Chris@87: (method2, docstring2), ...] Chris@87: Chris@87: This routine never raises an error. Chris@87: Chris@87: This routine cannot modify read-only docstrings, as appear Chris@87: in new-style classes or built-in functions. Because this Chris@87: routine never raises an error the caller must check manually Chris@87: that the docstrings were changed. Chris@87: """ Chris@87: try: Chris@87: new = getattr(__import__(place, globals(), {}, [obj]), obj) Chris@87: if isinstance(doc, str): Chris@87: add_docstring(new, doc.strip()) Chris@87: elif isinstance(doc, tuple): Chris@87: add_docstring(getattr(new, doc[0]), doc[1].strip()) Chris@87: elif isinstance(doc, list): Chris@87: for val in doc: Chris@87: add_docstring(getattr(new, val[0]), val[1].strip()) Chris@87: except: Chris@87: pass Chris@87: Chris@87: Chris@87: # Based on scitools meshgrid Chris@87: def meshgrid(*xi, **kwargs): Chris@87: """ Chris@87: Return coordinate matrices from coordinate vectors. Chris@87: Chris@87: Make N-D coordinate arrays for vectorized evaluations of Chris@87: N-D scalar/vector fields over N-D grids, given Chris@87: one-dimensional coordinate arrays x1, x2,..., xn. Chris@87: Chris@87: .. versionchanged:: 1.9 Chris@87: 1-D and 0-D cases are allowed. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: x1, x2,..., xn : array_like Chris@87: 1-D arrays representing the coordinates of a grid. Chris@87: indexing : {'xy', 'ij'}, optional Chris@87: Cartesian ('xy', default) or matrix ('ij') indexing of output. Chris@87: See Notes for more details. Chris@87: Chris@87: .. versionadded:: 1.7.0 Chris@87: sparse : bool, optional Chris@87: If True a sparse grid is returned in order to conserve memory. Chris@87: Default is False. Chris@87: Chris@87: .. versionadded:: 1.7.0 Chris@87: copy : bool, optional Chris@87: If False, a view into the original arrays are returned in order to Chris@87: conserve memory. Default is True. Please note that Chris@87: ``sparse=False, copy=False`` will likely return non-contiguous Chris@87: arrays. Furthermore, more than one element of a broadcast array Chris@87: may refer to a single memory location. If you need to write to the Chris@87: arrays, make copies first. Chris@87: Chris@87: .. versionadded:: 1.7.0 Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: X1, X2,..., XN : ndarray Chris@87: For vectors `x1`, `x2`,..., 'xn' with lengths ``Ni=len(xi)`` , Chris@87: return ``(N1, N2, N3,...Nn)`` shaped arrays if indexing='ij' Chris@87: or ``(N2, N1, N3,...Nn)`` shaped arrays if indexing='xy' Chris@87: with the elements of `xi` repeated to fill the matrix along Chris@87: the first dimension for `x1`, the second for `x2` and so on. Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: This function supports both indexing conventions through the indexing Chris@87: keyword argument. Giving the string 'ij' returns a meshgrid with Chris@87: matrix indexing, while 'xy' returns a meshgrid with Cartesian indexing. Chris@87: In the 2-D case with inputs of length M and N, the outputs are of shape Chris@87: (N, M) for 'xy' indexing and (M, N) for 'ij' indexing. In the 3-D case Chris@87: with inputs of length M, N and P, outputs are of shape (N, M, P) for Chris@87: 'xy' indexing and (M, N, P) for 'ij' indexing. The difference is Chris@87: illustrated by the following code snippet:: Chris@87: Chris@87: xv, yv = meshgrid(x, y, sparse=False, indexing='ij') Chris@87: for i in range(nx): Chris@87: for j in range(ny): Chris@87: # treat xv[i,j], yv[i,j] Chris@87: Chris@87: xv, yv = meshgrid(x, y, sparse=False, indexing='xy') Chris@87: for i in range(nx): Chris@87: for j in range(ny): Chris@87: # treat xv[j,i], yv[j,i] Chris@87: Chris@87: In the 1-D and 0-D case, the indexing and sparse keywords have no effect. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: index_tricks.mgrid : Construct a multi-dimensional "meshgrid" Chris@87: using indexing notation. Chris@87: index_tricks.ogrid : Construct an open multi-dimensional "meshgrid" Chris@87: using indexing notation. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> nx, ny = (3, 2) Chris@87: >>> x = np.linspace(0, 1, nx) Chris@87: >>> y = np.linspace(0, 1, ny) Chris@87: >>> xv, yv = meshgrid(x, y) Chris@87: >>> xv Chris@87: array([[ 0. , 0.5, 1. ], Chris@87: [ 0. , 0.5, 1. ]]) Chris@87: >>> yv Chris@87: array([[ 0., 0., 0.], Chris@87: [ 1., 1., 1.]]) Chris@87: >>> xv, yv = meshgrid(x, y, sparse=True) # make sparse output arrays Chris@87: >>> xv Chris@87: array([[ 0. , 0.5, 1. ]]) Chris@87: >>> yv Chris@87: array([[ 0.], Chris@87: [ 1.]]) Chris@87: Chris@87: `meshgrid` is very useful to evaluate functions on a grid. Chris@87: Chris@87: >>> x = np.arange(-5, 5, 0.1) Chris@87: >>> y = np.arange(-5, 5, 0.1) Chris@87: >>> xx, yy = meshgrid(x, y, sparse=True) Chris@87: >>> z = np.sin(xx**2 + yy**2) / (xx**2 + yy**2) Chris@87: >>> h = plt.contourf(x,y,z) Chris@87: Chris@87: """ Chris@87: ndim = len(xi) Chris@87: Chris@87: copy_ = kwargs.pop('copy', True) Chris@87: sparse = kwargs.pop('sparse', False) Chris@87: indexing = kwargs.pop('indexing', 'xy') Chris@87: Chris@87: if kwargs: Chris@87: raise TypeError("meshgrid() got an unexpected keyword argument '%s'" Chris@87: % (list(kwargs)[0],)) Chris@87: Chris@87: if indexing not in ['xy', 'ij']: Chris@87: raise ValueError( Chris@87: "Valid values for `indexing` are 'xy' and 'ij'.") Chris@87: Chris@87: s0 = (1,) * ndim Chris@87: output = [np.asanyarray(x).reshape(s0[:i] + (-1,) + s0[i + 1::]) Chris@87: for i, x in enumerate(xi)] Chris@87: Chris@87: shape = [x.size for x in output] Chris@87: Chris@87: if indexing == 'xy' and ndim > 1: Chris@87: # switch first and second axis Chris@87: output[0].shape = (1, -1) + (1,)*(ndim - 2) Chris@87: output[1].shape = (-1, 1) + (1,)*(ndim - 2) Chris@87: shape[0], shape[1] = shape[1], shape[0] Chris@87: Chris@87: if sparse: Chris@87: if copy_: Chris@87: return [x.copy() for x in output] Chris@87: else: Chris@87: return output Chris@87: else: Chris@87: # Return the full N-D matrix (not only the 1-D vector) Chris@87: if copy_: Chris@87: mult_fact = np.ones(shape, dtype=int) Chris@87: return [x * mult_fact for x in output] Chris@87: else: Chris@87: return np.broadcast_arrays(*output) Chris@87: Chris@87: Chris@87: def delete(arr, obj, axis=None): Chris@87: """ Chris@87: Return a new array with sub-arrays along an axis deleted. For a one Chris@87: dimensional array, this returns those entries not returned by Chris@87: `arr[obj]`. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: arr : array_like Chris@87: Input array. Chris@87: obj : slice, int or array of ints Chris@87: Indicate which sub-arrays to remove. Chris@87: axis : int, optional Chris@87: The axis along which to delete the subarray defined by `obj`. Chris@87: If `axis` is None, `obj` is applied to the flattened array. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: A copy of `arr` with the elements specified by `obj` removed. Note Chris@87: that `delete` does not occur in-place. If `axis` is None, `out` is Chris@87: a flattened array. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: insert : Insert elements into an array. Chris@87: append : Append elements at the end of an array. Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: Often it is preferable to use a boolean mask. For example: Chris@87: Chris@87: >>> mask = np.ones(len(arr), dtype=bool) Chris@87: >>> mask[[0,2,4]] = False Chris@87: >>> result = arr[mask,...] Chris@87: Chris@87: Is equivalent to `np.delete(arr, [0,2,4], axis=0)`, but allows further Chris@87: use of `mask`. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> arr = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]]) Chris@87: >>> arr Chris@87: array([[ 1, 2, 3, 4], Chris@87: [ 5, 6, 7, 8], Chris@87: [ 9, 10, 11, 12]]) Chris@87: >>> np.delete(arr, 1, 0) Chris@87: array([[ 1, 2, 3, 4], Chris@87: [ 9, 10, 11, 12]]) Chris@87: Chris@87: >>> np.delete(arr, np.s_[::2], 1) Chris@87: array([[ 2, 4], Chris@87: [ 6, 8], Chris@87: [10, 12]]) Chris@87: >>> np.delete(arr, [1,3,5], None) Chris@87: array([ 1, 3, 5, 7, 8, 9, 10, 11, 12]) Chris@87: Chris@87: """ Chris@87: wrap = None Chris@87: if type(arr) is not ndarray: Chris@87: try: Chris@87: wrap = arr.__array_wrap__ Chris@87: except AttributeError: Chris@87: pass Chris@87: Chris@87: arr = asarray(arr) Chris@87: ndim = arr.ndim Chris@87: if axis is None: Chris@87: if ndim != 1: Chris@87: arr = arr.ravel() Chris@87: ndim = arr.ndim Chris@87: axis = ndim - 1 Chris@87: if ndim == 0: Chris@87: warnings.warn( Chris@87: "in the future the special handling of scalars will be removed " Chris@87: "from delete and raise an error", DeprecationWarning) Chris@87: if wrap: Chris@87: return wrap(arr) Chris@87: else: Chris@87: return arr.copy() Chris@87: Chris@87: slobj = [slice(None)]*ndim Chris@87: N = arr.shape[axis] Chris@87: newshape = list(arr.shape) Chris@87: Chris@87: if isinstance(obj, slice): Chris@87: start, stop, step = obj.indices(N) Chris@87: xr = range(start, stop, step) Chris@87: numtodel = len(xr) Chris@87: Chris@87: if numtodel <= 0: Chris@87: if wrap: Chris@87: return wrap(arr.copy()) Chris@87: else: Chris@87: return arr.copy() Chris@87: Chris@87: # Invert if step is negative: Chris@87: if step < 0: Chris@87: step = -step Chris@87: start = xr[-1] Chris@87: stop = xr[0] + 1 Chris@87: Chris@87: newshape[axis] -= numtodel Chris@87: new = empty(newshape, arr.dtype, arr.flags.fnc) Chris@87: # copy initial chunk Chris@87: if start == 0: Chris@87: pass Chris@87: else: Chris@87: slobj[axis] = slice(None, start) Chris@87: new[slobj] = arr[slobj] Chris@87: # copy end chunck Chris@87: if stop == N: Chris@87: pass Chris@87: else: Chris@87: slobj[axis] = slice(stop-numtodel, None) Chris@87: slobj2 = [slice(None)]*ndim Chris@87: slobj2[axis] = slice(stop, None) Chris@87: new[slobj] = arr[slobj2] Chris@87: # copy middle pieces Chris@87: if step == 1: Chris@87: pass Chris@87: else: # use array indexing. Chris@87: keep = ones(stop-start, dtype=bool) Chris@87: keep[:stop-start:step] = False Chris@87: slobj[axis] = slice(start, stop-numtodel) Chris@87: slobj2 = [slice(None)]*ndim Chris@87: slobj2[axis] = slice(start, stop) Chris@87: arr = arr[slobj2] Chris@87: slobj2[axis] = keep Chris@87: new[slobj] = arr[slobj2] Chris@87: if wrap: Chris@87: return wrap(new) Chris@87: else: Chris@87: return new Chris@87: Chris@87: _obj = obj Chris@87: obj = np.asarray(obj) Chris@87: # After removing the special handling of booleans and out of Chris@87: # bounds values, the conversion to the array can be removed. Chris@87: if obj.dtype == bool: Chris@87: warnings.warn( Chris@87: "in the future insert will treat boolean arrays and array-likes " Chris@87: "as boolean index instead of casting it to integer", FutureWarning) Chris@87: obj = obj.astype(intp) Chris@87: if isinstance(_obj, (int, long, integer)): Chris@87: # optimization for a single value Chris@87: obj = obj.item() Chris@87: if (obj < -N or obj >= N): Chris@87: raise IndexError( Chris@87: "index %i is out of bounds for axis %i with " Chris@87: "size %i" % (obj, axis, N)) Chris@87: if (obj < 0): Chris@87: obj += N Chris@87: newshape[axis] -= 1 Chris@87: new = empty(newshape, arr.dtype, arr.flags.fnc) Chris@87: slobj[axis] = slice(None, obj) Chris@87: new[slobj] = arr[slobj] Chris@87: slobj[axis] = slice(obj, None) Chris@87: slobj2 = [slice(None)]*ndim Chris@87: slobj2[axis] = slice(obj+1, None) Chris@87: new[slobj] = arr[slobj2] Chris@87: else: Chris@87: if obj.size == 0 and not isinstance(_obj, np.ndarray): Chris@87: obj = obj.astype(intp) Chris@87: if not np.can_cast(obj, intp, 'same_kind'): Chris@87: # obj.size = 1 special case always failed and would just Chris@87: # give superfluous warnings. Chris@87: warnings.warn( Chris@87: "using a non-integer array as obj in delete will result in an " Chris@87: "error in the future", DeprecationWarning) Chris@87: obj = obj.astype(intp) Chris@87: keep = ones(N, dtype=bool) Chris@87: Chris@87: # Test if there are out of bound indices, this is deprecated Chris@87: inside_bounds = (obj < N) & (obj >= -N) Chris@87: if not inside_bounds.all(): Chris@87: warnings.warn( Chris@87: "in the future out of bounds indices will raise an error " Chris@87: "instead of being ignored by `numpy.delete`.", Chris@87: DeprecationWarning) Chris@87: obj = obj[inside_bounds] Chris@87: positive_indices = obj >= 0 Chris@87: if not positive_indices.all(): Chris@87: warnings.warn( Chris@87: "in the future negative indices will not be ignored by " Chris@87: "`numpy.delete`.", FutureWarning) Chris@87: obj = obj[positive_indices] Chris@87: Chris@87: keep[obj, ] = False Chris@87: slobj[axis] = keep Chris@87: new = arr[slobj] Chris@87: Chris@87: if wrap: Chris@87: return wrap(new) Chris@87: else: Chris@87: return new Chris@87: Chris@87: Chris@87: def insert(arr, obj, values, axis=None): Chris@87: """ Chris@87: Insert values along the given axis before the given indices. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: arr : array_like Chris@87: Input array. Chris@87: obj : int, slice or sequence of ints Chris@87: Object that defines the index or indices before which `values` is Chris@87: inserted. Chris@87: Chris@87: .. versionadded:: 1.8.0 Chris@87: Chris@87: Support for multiple insertions when `obj` is a single scalar or a Chris@87: sequence with one element (similar to calling insert multiple Chris@87: times). Chris@87: values : array_like Chris@87: Values to insert into `arr`. If the type of `values` is different Chris@87: from that of `arr`, `values` is converted to the type of `arr`. Chris@87: `values` should be shaped so that ``arr[...,obj,...] = values`` Chris@87: is legal. Chris@87: axis : int, optional Chris@87: Axis along which to insert `values`. If `axis` is None then `arr` Chris@87: is flattened first. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: out : ndarray Chris@87: A copy of `arr` with `values` inserted. Note that `insert` Chris@87: does not occur in-place: a new array is returned. If Chris@87: `axis` is None, `out` is a flattened array. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: append : Append elements at the end of an array. Chris@87: concatenate : Join a sequence of arrays together. Chris@87: delete : Delete elements from an array. Chris@87: Chris@87: Notes Chris@87: ----- Chris@87: Note that for higher dimensional inserts `obj=0` behaves very different Chris@87: from `obj=[0]` just like `arr[:,0,:] = values` is different from Chris@87: `arr[:,[0],:] = values`. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> a = np.array([[1, 1], [2, 2], [3, 3]]) Chris@87: >>> a Chris@87: array([[1, 1], Chris@87: [2, 2], Chris@87: [3, 3]]) Chris@87: >>> np.insert(a, 1, 5) Chris@87: array([1, 5, 1, 2, 2, 3, 3]) Chris@87: >>> np.insert(a, 1, 5, axis=1) Chris@87: array([[1, 5, 1], Chris@87: [2, 5, 2], Chris@87: [3, 5, 3]]) Chris@87: Chris@87: Difference between sequence and scalars: Chris@87: >>> np.insert(a, [1], [[1],[2],[3]], axis=1) Chris@87: array([[1, 1, 1], Chris@87: [2, 2, 2], Chris@87: [3, 3, 3]]) Chris@87: >>> np.array_equal(np.insert(a, 1, [1, 2, 3], axis=1), Chris@87: ... np.insert(a, [1], [[1],[2],[3]], axis=1)) Chris@87: True Chris@87: Chris@87: >>> b = a.flatten() Chris@87: >>> b Chris@87: array([1, 1, 2, 2, 3, 3]) Chris@87: >>> np.insert(b, [2, 2], [5, 6]) Chris@87: array([1, 1, 5, 6, 2, 2, 3, 3]) Chris@87: Chris@87: >>> np.insert(b, slice(2, 4), [5, 6]) Chris@87: array([1, 1, 5, 2, 6, 2, 3, 3]) Chris@87: Chris@87: >>> np.insert(b, [2, 2], [7.13, False]) # type casting Chris@87: array([1, 1, 7, 0, 2, 2, 3, 3]) Chris@87: Chris@87: >>> x = np.arange(8).reshape(2, 4) Chris@87: >>> idx = (1, 3) Chris@87: >>> np.insert(x, idx, 999, axis=1) Chris@87: array([[ 0, 999, 1, 2, 999, 3], Chris@87: [ 4, 999, 5, 6, 999, 7]]) Chris@87: Chris@87: """ Chris@87: wrap = None Chris@87: if type(arr) is not ndarray: Chris@87: try: Chris@87: wrap = arr.__array_wrap__ Chris@87: except AttributeError: Chris@87: pass Chris@87: Chris@87: arr = asarray(arr) Chris@87: ndim = arr.ndim Chris@87: if axis is None: Chris@87: if ndim != 1: Chris@87: arr = arr.ravel() Chris@87: ndim = arr.ndim Chris@87: axis = ndim - 1 Chris@87: else: Chris@87: if ndim > 0 and (axis < -ndim or axis >= ndim): Chris@87: raise IndexError( Chris@87: "axis %i is out of bounds for an array of " Chris@87: "dimension %i" % (axis, ndim)) Chris@87: if (axis < 0): Chris@87: axis += ndim Chris@87: if (ndim == 0): Chris@87: warnings.warn( Chris@87: "in the future the special handling of scalars will be removed " Chris@87: "from insert and raise an error", DeprecationWarning) Chris@87: arr = arr.copy() Chris@87: arr[...] = values Chris@87: if wrap: Chris@87: return wrap(arr) Chris@87: else: Chris@87: return arr Chris@87: slobj = [slice(None)]*ndim Chris@87: N = arr.shape[axis] Chris@87: newshape = list(arr.shape) Chris@87: Chris@87: if isinstance(obj, slice): Chris@87: # turn it into a range object Chris@87: indices = arange(*obj.indices(N), **{'dtype': intp}) Chris@87: else: Chris@87: # need to copy obj, because indices will be changed in-place Chris@87: indices = np.array(obj) Chris@87: if indices.dtype == bool: Chris@87: # See also delete Chris@87: warnings.warn( Chris@87: "in the future insert will treat boolean arrays and " Chris@87: "array-likes as a boolean index instead of casting it to " Chris@87: "integer", FutureWarning) Chris@87: indices = indices.astype(intp) Chris@87: # Code after warning period: Chris@87: #if obj.ndim != 1: Chris@87: # raise ValueError('boolean array argument obj to insert ' Chris@87: # 'must be one dimensional') Chris@87: #indices = np.flatnonzero(obj) Chris@87: elif indices.ndim > 1: Chris@87: raise ValueError( Chris@87: "index array argument obj to insert must be one dimensional " Chris@87: "or scalar") Chris@87: if indices.size == 1: Chris@87: index = indices.item() Chris@87: if index < -N or index > N: Chris@87: raise IndexError( Chris@87: "index %i is out of bounds for axis %i with " Chris@87: "size %i" % (obj, axis, N)) Chris@87: if (index < 0): Chris@87: index += N Chris@87: Chris@87: # There are some object array corner cases here, but we cannot avoid Chris@87: # that: Chris@87: values = array(values, copy=False, ndmin=arr.ndim, dtype=arr.dtype) Chris@87: if indices.ndim == 0: Chris@87: # broadcasting is very different here, since a[:,0,:] = ... behaves Chris@87: # very different from a[:,[0],:] = ...! This changes values so that Chris@87: # it works likes the second case. (here a[:,0:1,:]) Chris@87: values = np.rollaxis(values, 0, (axis % values.ndim) + 1) Chris@87: numnew = values.shape[axis] Chris@87: newshape[axis] += numnew Chris@87: new = empty(newshape, arr.dtype, arr.flags.fnc) Chris@87: slobj[axis] = slice(None, index) Chris@87: new[slobj] = arr[slobj] Chris@87: slobj[axis] = slice(index, index+numnew) Chris@87: new[slobj] = values Chris@87: slobj[axis] = slice(index+numnew, None) Chris@87: slobj2 = [slice(None)] * ndim Chris@87: slobj2[axis] = slice(index, None) Chris@87: new[slobj] = arr[slobj2] Chris@87: if wrap: Chris@87: return wrap(new) Chris@87: return new Chris@87: elif indices.size == 0 and not isinstance(obj, np.ndarray): Chris@87: # Can safely cast the empty list to intp Chris@87: indices = indices.astype(intp) Chris@87: Chris@87: if not np.can_cast(indices, intp, 'same_kind'): Chris@87: warnings.warn( Chris@87: "using a non-integer array as obj in insert will result in an " Chris@87: "error in the future", DeprecationWarning) Chris@87: indices = indices.astype(intp) Chris@87: Chris@87: indices[indices < 0] += N Chris@87: Chris@87: numnew = len(indices) Chris@87: order = indices.argsort(kind='mergesort') # stable sort Chris@87: indices[order] += np.arange(numnew) Chris@87: Chris@87: newshape[axis] += numnew Chris@87: old_mask = ones(newshape[axis], dtype=bool) Chris@87: old_mask[indices] = False Chris@87: Chris@87: new = empty(newshape, arr.dtype, arr.flags.fnc) Chris@87: slobj2 = [slice(None)]*ndim Chris@87: slobj[axis] = indices Chris@87: slobj2[axis] = old_mask Chris@87: new[slobj] = values Chris@87: new[slobj2] = arr Chris@87: Chris@87: if wrap: Chris@87: return wrap(new) Chris@87: return new Chris@87: Chris@87: Chris@87: def append(arr, values, axis=None): Chris@87: """ Chris@87: Append values to the end of an array. Chris@87: Chris@87: Parameters Chris@87: ---------- Chris@87: arr : array_like Chris@87: Values are appended to a copy of this array. Chris@87: values : array_like Chris@87: These values are appended to a copy of `arr`. It must be of the Chris@87: correct shape (the same shape as `arr`, excluding `axis`). If Chris@87: `axis` is not specified, `values` can be any shape and will be Chris@87: flattened before use. Chris@87: axis : int, optional Chris@87: The axis along which `values` are appended. If `axis` is not Chris@87: given, both `arr` and `values` are flattened before use. Chris@87: Chris@87: Returns Chris@87: ------- Chris@87: append : ndarray Chris@87: A copy of `arr` with `values` appended to `axis`. Note that Chris@87: `append` does not occur in-place: a new array is allocated and Chris@87: filled. If `axis` is None, `out` is a flattened array. Chris@87: Chris@87: See Also Chris@87: -------- Chris@87: insert : Insert elements into an array. Chris@87: delete : Delete elements from an array. Chris@87: Chris@87: Examples Chris@87: -------- Chris@87: >>> np.append([1, 2, 3], [[4, 5, 6], [7, 8, 9]]) Chris@87: array([1, 2, 3, 4, 5, 6, 7, 8, 9]) Chris@87: Chris@87: When `axis` is specified, `values` must have the correct shape. Chris@87: Chris@87: >>> np.append([[1, 2, 3], [4, 5, 6]], [[7, 8, 9]], axis=0) Chris@87: array([[1, 2, 3], Chris@87: [4, 5, 6], Chris@87: [7, 8, 9]]) Chris@87: >>> np.append([[1, 2, 3], [4, 5, 6]], [7, 8, 9], axis=0) Chris@87: Traceback (most recent call last): Chris@87: ... Chris@87: ValueError: arrays must have same number of dimensions Chris@87: Chris@87: """ Chris@87: arr = asanyarray(arr) Chris@87: if axis is None: Chris@87: if arr.ndim != 1: Chris@87: arr = arr.ravel() Chris@87: values = ravel(values) Chris@87: axis = arr.ndim-1 Chris@87: return concatenate((arr, values), axis=axis)