diff options
Diffstat (limited to 'numpy/lib/function_base.py')
| -rw-r--r-- | numpy/lib/function_base.py | 149 |
1 files changed, 108 insertions, 41 deletions
diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py index 0a1d05f77..135053e43 100644 --- a/numpy/lib/function_base.py +++ b/numpy/lib/function_base.py @@ -263,7 +263,7 @@ def histogramdd(sample, bins=10, range=None, normed=False, weights=None): normed : bool, optional If False, returns the number of samples in each bin. If True, returns the bin density ``bin_count / sample_count / bin_volume``. - weights : array_like (N,), optional + weights : (N,) array_like, optional An array of values `w_i` weighing each sample `(x_i, y_i, z_i, ...)`. Weights are normalized to 1 if normed is True. If normed is False, the values of the returned histogram are equal to the sum of the @@ -337,6 +337,11 @@ def histogramdd(sample, bins=10, range=None, normed=False, weights=None): smin[i] = smin[i] - .5 smax[i] = smax[i] + .5 + # avoid rounding issues for comparisons when dealing with inexact types + if np.issubdtype(sample.dtype, np.inexact): + edge_dt = sample.dtype + else: + edge_dt = float # Create edge arrays for i in arange(D): if isscalar(bins[i]): @@ -345,9 +350,9 @@ def histogramdd(sample, bins=10, range=None, normed=False, weights=None): "Element at index %s in `bins` should be a positive " "integer." % i) nbin[i] = bins[i] + 2 # +2 for outlier bins - edges[i] = linspace(smin[i], smax[i], nbin[i]-1) + edges[i] = linspace(smin[i], smax[i], nbin[i]-1, dtype=edge_dt) else: - edges[i] = asarray(bins[i], float) + edges[i] = asarray(bins[i], edge_dt) nbin[i] = len(edges[i]) + 1 # +1 for outlier bins dedges[i] = diff(edges[i]) if np.any(np.asarray(dedges[i]) <= 0): @@ -456,7 +461,7 @@ def average(a, axis=None, weights=None, returned=False): Returns ------- - average, [sum_of_weights] : {array_type, double} + average, [sum_of_weights] : array_type or double Return the average along the specified axis. When returned is `True`, return a tuple with the average as the first element and the sum of the weights as the second element. The return type is `Float` @@ -510,8 +515,7 @@ def average(a, axis=None, weights=None, returned=False): scl = avg.dtype.type(a.size/avg.size) else: a = a + 0.0 - wgt = np.array(weights, dtype=a.dtype, copy=0) - + wgt = np.asarray(weights) # Sanity checks if a.shape != wgt.shape: if axis is None: @@ -528,7 +532,7 @@ def average(a, axis=None, weights=None, returned=False): # setup wgt to broadcast along axis wgt = np.array(wgt, copy=0, ndmin=a.ndim).swapaxes(-1, axis) - scl = wgt.sum(axis=axis) + scl = wgt.sum(axis=axis, dtype=np.result_type(a.dtype, wgt.dtype)) if (scl == 0.0).any(): raise ZeroDivisionError( "Weights sum to zero, can't be normalized") @@ -878,28 +882,33 @@ def copy(a, order='K'): # Basic operations -def gradient(f, *varargs): +def gradient(f, *varargs, **kwargs): """ Return the gradient of an N-dimensional array. The gradient is computed using second order accurate central differences - in the interior and second order accurate one-sides (forward or backwards) - differences at the boundaries. The returned gradient hence has the same - shape as the input array. + in the interior and either first differences or second order accurate + one-sides (forward or backwards) differences at the boundaries. The + returned gradient hence has the same shape as the input array. Parameters ---------- f : array_like - An N-dimensional array containing samples of a scalar function. - `*varargs` : scalars - 0, 1, or N scalars specifying the sample distances in each direction, - that is: `dx`, `dy`, `dz`, ... The default distance is 1. + An N-dimensional array containing samples of a scalar function. + varargs : list of scalar, optional + N scalars specifying the sample distances for each dimension, + i.e. `dx`, `dy`, `dz`, ... Default distance: 1. + edge_order : {1, 2}, optional + Gradient is calculated using N\ :sup:`th` order accurate differences + at the boundaries. Default: 1. + + .. versionadded:: 1.9.1 Returns ------- gradient : ndarray - N arrays of the same shape as `f` giving the derivative of `f` with - respect to each dimension. + N arrays of the same shape as `f` giving the derivative of `f` with + respect to each dimension. Examples -------- @@ -911,15 +920,14 @@ def gradient(f, *varargs): >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=np.float)) [array([[ 2., 2., -1.], - [ 2., 2., -1.]]), - array([[ 1. , 2.5, 4. ], - [ 1. , 1. , 1. ]])] + [ 2., 2., -1.]]), array([[ 1. , 2.5, 4. ], + [ 1. , 1. , 1. ]])] - >>> x = np.array([0,1,2,3,4]) - >>> dx = gradient(x) + >>> x = np.array([0, 1, 2, 3, 4]) + >>> dx = np.gradient(x) >>> y = x**2 - >>> gradient(y,dx) - array([0., 2., 4., 6., 8.]) + >>> np.gradient(y, dx, edge_order=2) + array([-0., 2., 4., 6., 8.]) """ f = np.asanyarray(f) N = len(f.shape) # number of dimensions @@ -934,6 +942,13 @@ def gradient(f, *varargs): raise SyntaxError( "invalid number of arguments") + edge_order = kwargs.pop('edge_order', 1) + if kwargs: + raise TypeError('"{}" are not valid keyword arguments.'.format( + '", "'.join(kwargs.keys()))) + if edge_order > 2: + raise ValueError("'edge_order' greater than 2 not supported") + # use central differences on interior and one-sided differences on the # endpoints. This preserves second order-accuracy over the full domain. @@ -973,7 +988,7 @@ def gradient(f, *varargs): "at least two elements are required.") # Numerical differentiation: 1st order edges, 2nd order interior - if y.shape[axis] == 2: + if y.shape[axis] == 2 or edge_order == 1: # Use first order differences for time data out = np.empty_like(y, dtype=otype) @@ -1021,7 +1036,8 @@ def gradient(f, *varargs): out[slice1] = (3.0*y[slice2] - 4.0*y[slice3] + y[slice4])/2.0 # divide by step size - outvals.append(out / dx[axis]) + out /= dx[axis] + outvals.append(out) # reset the slice object in this dimension to ":" slice1[axis] = slice(None) @@ -1097,7 +1113,7 @@ def diff(a, n=1, axis=-1): return a[slice1]-a[slice2] -def interp(x, xp, fp, left=None, right=None): +def interp(x, xp, fp, left=None, right=None, period=None): """ One-dimensional linear interpolation. @@ -1110,7 +1126,9 @@ def interp(x, xp, fp, left=None, right=None): The x-coordinates of the interpolated values. xp : 1-D sequence of floats - The x-coordinates of the data points, must be increasing. + The x-coordinates of the data points, must be increasing if argument + `period` is not specified. Otherwise, `xp` is internally sorted after + normalizing the periodic boundaries with ``xp = xp % period``. fp : 1-D sequence of floats The y-coordinates of the data points, same length as `xp`. @@ -1121,15 +1139,23 @@ def interp(x, xp, fp, left=None, right=None): right : float, optional Value to return for `x > xp[-1]`, default is `fp[-1]`. + period : None or float, optional + .. versionadded:: 1.10.0 + A period for the x-coordinates. This parameter allows the proper + interpolation of angular x-coordinates. Parameters `left` and `right` + are ignored if `period` is specified. + Returns ------- - y : {float, ndarray} + y : float or ndarray The interpolated values, same shape as `x`. Raises ------ ValueError If `xp` and `fp` have different length + If `xp` or `fp` are not 1-D sequences + If `period == 0` Notes ----- @@ -1139,7 +1165,6 @@ def interp(x, xp, fp, left=None, right=None): np.all(np.diff(xp) > 0) - Examples -------- >>> xp = [1, 2, 3] @@ -1165,13 +1190,51 @@ def interp(x, xp, fp, left=None, right=None): [<matplotlib.lines.Line2D object at 0x...>] >>> plt.show() + Interpolation with periodic x-coordinates: + + >>> x = [-180, -170, -185, 185, -10, -5, 0, 365] + >>> xp = [190, -190, 350, -350] + >>> fp = [5, 10, 3, 4] + >>> np.interp(x, xp, fp, period=360) + array([7.5, 5., 8.75, 6.25, 3., 3.25, 3.5, 3.75]) + """ - if isinstance(x, (float, int, number)): - return compiled_interp([x], xp, fp, left, right).item() - elif isinstance(x, np.ndarray) and x.ndim == 0: - return compiled_interp([x], xp, fp, left, right).item() + if period is None: + if isinstance(x, (float, int, number)): + return compiled_interp([x], xp, fp, left, right).item() + elif isinstance(x, np.ndarray) and x.ndim == 0: + return compiled_interp([x], xp, fp, left, right).item() + else: + return compiled_interp(x, xp, fp, left, right) else: - return compiled_interp(x, xp, fp, left, right) + if period == 0: + raise ValueError("period must be a non-zero value") + period = abs(period) + left = None + right = None + return_array = True + if isinstance(x, (float, int, number)): + return_array = False + x = [x] + x = np.asarray(x, dtype=np.float64) + xp = np.asarray(xp, dtype=np.float64) + fp = np.asarray(fp, dtype=np.float64) + if xp.ndim != 1 or fp.ndim != 1: + raise ValueError("Data points must be 1-D sequences") + if xp.shape[0] != fp.shape[0]: + raise ValueError("fp and xp are not of the same length") + # normalizing periodic boundaries + x = x % period + xp = xp % period + asort_xp = np.argsort(xp) + xp = xp[asort_xp] + fp = fp[asort_xp] + xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period)) + fp = np.concatenate((fp[-1:], fp, fp[0:1])) + if return_array: + return compiled_interp(x, xp, fp, left, right) + else: + return compiled_interp(x, xp, fp, left, right).item() def angle(z, deg=0): @@ -1187,7 +1250,7 @@ def angle(z, deg=0): Returns ------- - angle : {ndarray, scalar} + angle : ndarray or scalar The counterclockwise angle from the positive real axis on the complex plane, with dtype as numpy.float64. @@ -1387,6 +1450,8 @@ def extract(condition, arr): This is equivalent to ``np.compress(ravel(condition), ravel(arr))``. If `condition` is boolean ``np.extract`` is equivalent to ``arr[condition]``. + Note that `place` does the exact opposite of `extract`. + Parameters ---------- condition : array_like @@ -1402,7 +1467,7 @@ def extract(condition, arr): See Also -------- - take, put, copyto, compress + take, put, copyto, compress, place Examples -------- @@ -1915,7 +1980,7 @@ def corrcoef(x, y=None, rowvar=1, bias=0, ddof=None): observations (unbiased estimate). If `bias` is 1, then normalization is by ``N``. These values can be overridden by using the keyword ``ddof`` in numpy versions >= 1.5. - ddof : {None, int}, optional + ddof : int, optional .. versionadded:: 1.5 If not ``None`` normalization is by ``(N - ddof)``, where ``N`` is the number of observations; this overrides the value implied by @@ -2998,7 +3063,7 @@ def percentile(a, q, axis=None, out=None, nearest neighbors as well as the `interpolation` parameter will determine the percentile if the normalized ranking does not match q exactly. This function is the same as the median if ``q=50``, the same - as the minimum if ``q=0``and the same as the maximum if ``q=100``. + as the minimum if ``q=0`` and the same as the maximum if ``q=100``. Examples -------- @@ -3031,7 +3096,7 @@ def percentile(a, q, axis=None, out=None, array([ 3.5]) """ - q = asarray(q, dtype=np.float64) + q = array(q, dtype=np.float64, copy=True) r, k = _ureduce(a, func=_percentile, q=q, axis=axis, out=out, overwrite_input=overwrite_input, interpolation=interpolation) @@ -3758,7 +3823,9 @@ def insert(arr, obj, values, axis=None): if (index < 0): index += N - values = array(values, copy=False, ndmin=arr.ndim) + # There are some object array corner cases here, but we cannot avoid + # that: + values = array(values, copy=False, ndmin=arr.ndim, dtype=arr.dtype) if indices.ndim == 0: # broadcasting is very different here, since a[:,0,:] = ... behaves # very different from a[:,[0],:] = ...! This changes values so that |