Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1import collections.abc
2import functools
3import re
4import sys
5import warnings
7import numpy as np
8import numpy.core.numeric as _nx
9from numpy.core import transpose
10from numpy.core.numeric import (
11 ones, zeros_like, arange, concatenate, array, asarray, asanyarray, empty,
12 ndarray, take, dot, where, intp, integer, isscalar, absolute
13 )
14from numpy.core.umath import (
15 pi, add, arctan2, frompyfunc, cos, less_equal, sqrt, sin,
16 mod, exp, not_equal, subtract
17 )
18from numpy.core.fromnumeric import (
19 ravel, nonzero, partition, mean, any, sum
20 )
21from numpy.core.numerictypes import typecodes
22from numpy.core.overrides import set_module
23from numpy.core import overrides
24from numpy.core.function_base import add_newdoc
25from numpy.lib.twodim_base import diag
26from numpy.core.multiarray import (
27 _insert, add_docstring, bincount, normalize_axis_index, _monotonicity,
28 interp as compiled_interp, interp_complex as compiled_interp_complex
29 )
30from numpy.core.umath import _add_newdoc_ufunc as add_newdoc_ufunc
32import builtins
34# needed in this module for compatibility
35from numpy.lib.histograms import histogram, histogramdd # noqa: F401
38array_function_dispatch = functools.partial(
39 overrides.array_function_dispatch, module='numpy')
42__all__ = [
43 'select', 'piecewise', 'trim_zeros', 'copy', 'iterable', 'percentile',
44 'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'disp', 'flip',
45 'rot90', 'extract', 'place', 'vectorize', 'asarray_chkfinite', 'average',
46 'bincount', 'digitize', 'cov', 'corrcoef',
47 'msort', 'median', 'sinc', 'hamming', 'hanning', 'bartlett',
48 'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc', 'add_docstring',
49 'meshgrid', 'delete', 'insert', 'append', 'interp', 'add_newdoc_ufunc',
50 'quantile'
51 ]
53# _QuantileMethods is a dictionary listing all the supported methods to
54# compute quantile/percentile.
55#
56# Below virtual_index refer to the index of the element where the percentile
57# would be found in the sorted sample.
58# When the sample contains exactly the percentile wanted, the virtual_index is
59# an integer to the index of this element.
60# When the percentile wanted is in between two elements, the virtual_index
61# is made of a integer part (a.k.a 'i' or 'left') and a fractional part
62# (a.k.a 'g' or 'gamma')
63#
64# Each method in _QuantileMethods has two properties
65# get_virtual_index : Callable
66# The function used to compute the virtual_index.
67# fix_gamma : Callable
68# A function used for discret methods to force the index to a specific value.
69_QuantileMethods = dict(
70 # --- HYNDMAN and FAN METHODS
71 # Discrete methods
72 inverted_cdf=dict(
73 get_virtual_index=lambda n, quantiles: _inverted_cdf(n, quantiles),
74 fix_gamma=lambda gamma, _: gamma, # should never be called
75 ),
76 averaged_inverted_cdf=dict(
77 get_virtual_index=lambda n, quantiles: (n * quantiles) - 1,
78 fix_gamma=lambda gamma, _: _get_gamma_mask(
79 shape=gamma.shape,
80 default_value=1.,
81 conditioned_value=0.5,
82 where=gamma == 0),
83 ),
84 closest_observation=dict(
85 get_virtual_index=lambda n, quantiles: _closest_observation(n,
86 quantiles),
87 fix_gamma=lambda gamma, _: gamma, # should never be called
88 ),
89 # Continuous methods
90 interpolated_inverted_cdf=dict(
91 get_virtual_index=lambda n, quantiles:
92 _compute_virtual_index(n, quantiles, 0, 1),
93 fix_gamma=lambda gamma, _: gamma,
94 ),
95 hazen=dict(
96 get_virtual_index=lambda n, quantiles:
97 _compute_virtual_index(n, quantiles, 0.5, 0.5),
98 fix_gamma=lambda gamma, _: gamma,
99 ),
100 weibull=dict(
101 get_virtual_index=lambda n, quantiles:
102 _compute_virtual_index(n, quantiles, 0, 0),
103 fix_gamma=lambda gamma, _: gamma,
104 ),
105 # Default method.
106 # To avoid some rounding issues, `(n-1) * quantiles` is preferred to
107 # `_compute_virtual_index(n, quantiles, 1, 1)`.
108 # They are mathematically equivalent.
109 linear=dict(
110 get_virtual_index=lambda n, quantiles: (n - 1) * quantiles,
111 fix_gamma=lambda gamma, _: gamma,
112 ),
113 median_unbiased=dict(
114 get_virtual_index=lambda n, quantiles:
115 _compute_virtual_index(n, quantiles, 1 / 3.0, 1 / 3.0),
116 fix_gamma=lambda gamma, _: gamma,
117 ),
118 normal_unbiased=dict(
119 get_virtual_index=lambda n, quantiles:
120 _compute_virtual_index(n, quantiles, 3 / 8.0, 3 / 8.0),
121 fix_gamma=lambda gamma, _: gamma,
122 ),
123 # --- OTHER METHODS
124 lower=dict(
125 get_virtual_index=lambda n, quantiles: np.floor(
126 (n - 1) * quantiles).astype(np.intp),
127 fix_gamma=lambda gamma, _: gamma,
128 # should never be called, index dtype is int
129 ),
130 higher=dict(
131 get_virtual_index=lambda n, quantiles: np.ceil(
132 (n - 1) * quantiles).astype(np.intp),
133 fix_gamma=lambda gamma, _: gamma,
134 # should never be called, index dtype is int
135 ),
136 midpoint=dict(
137 get_virtual_index=lambda n, quantiles: 0.5 * (
138 np.floor((n - 1) * quantiles)
139 + np.ceil((n - 1) * quantiles)),
140 fix_gamma=lambda gamma, index: _get_gamma_mask(
141 shape=gamma.shape,
142 default_value=0.5,
143 conditioned_value=0.,
144 where=index % 1 == 0),
145 ),
146 nearest=dict(
147 get_virtual_index=lambda n, quantiles: np.around(
148 (n - 1) * quantiles).astype(np.intp),
149 fix_gamma=lambda gamma, _: gamma,
150 # should never be called, index dtype is int
151 ))
154def _rot90_dispatcher(m, k=None, axes=None):
155 return (m,)
158@array_function_dispatch(_rot90_dispatcher)
159def rot90(m, k=1, axes=(0, 1)):
160 """
161 Rotate an array by 90 degrees in the plane specified by axes.
163 Rotation direction is from the first towards the second axis.
165 Parameters
166 ----------
167 m : array_like
168 Array of two or more dimensions.
169 k : integer
170 Number of times the array is rotated by 90 degrees.
171 axes : (2,) array_like
172 The array is rotated in the plane defined by the axes.
173 Axes must be different.
175 .. versionadded:: 1.12.0
177 Returns
178 -------
179 y : ndarray
180 A rotated view of `m`.
182 See Also
183 --------
184 flip : Reverse the order of elements in an array along the given axis.
185 fliplr : Flip an array horizontally.
186 flipud : Flip an array vertically.
188 Notes
189 -----
190 ``rot90(m, k=1, axes=(1,0))`` is the reverse of
191 ``rot90(m, k=1, axes=(0,1))``
193 ``rot90(m, k=1, axes=(1,0))`` is equivalent to
194 ``rot90(m, k=-1, axes=(0,1))``
196 Examples
197 --------
198 >>> m = np.array([[1,2],[3,4]], int)
199 >>> m
200 array([[1, 2],
201 [3, 4]])
202 >>> np.rot90(m)
203 array([[2, 4],
204 [1, 3]])
205 >>> np.rot90(m, 2)
206 array([[4, 3],
207 [2, 1]])
208 >>> m = np.arange(8).reshape((2,2,2))
209 >>> np.rot90(m, 1, (1,2))
210 array([[[1, 3],
211 [0, 2]],
212 [[5, 7],
213 [4, 6]]])
215 """
216 axes = tuple(axes)
217 if len(axes) != 2:
218 raise ValueError("len(axes) must be 2.")
220 m = asanyarray(m)
222 if axes[0] == axes[1] or absolute(axes[0] - axes[1]) == m.ndim:
223 raise ValueError("Axes must be different.")
225 if (axes[0] >= m.ndim or axes[0] < -m.ndim
226 or axes[1] >= m.ndim or axes[1] < -m.ndim):
227 raise ValueError("Axes={} out of range for array of ndim={}."
228 .format(axes, m.ndim))
230 k %= 4
232 if k == 0:
233 return m[:]
234 if k == 2:
235 return flip(flip(m, axes[0]), axes[1])
237 axes_list = arange(0, m.ndim)
238 (axes_list[axes[0]], axes_list[axes[1]]) = (axes_list[axes[1]],
239 axes_list[axes[0]])
241 if k == 1:
242 return transpose(flip(m, axes[1]), axes_list)
243 else:
244 # k == 3
245 return flip(transpose(m, axes_list), axes[1])
248def _flip_dispatcher(m, axis=None):
249 return (m,)
252@array_function_dispatch(_flip_dispatcher)
253def flip(m, axis=None):
254 """
255 Reverse the order of elements in an array along the given axis.
257 The shape of the array is preserved, but the elements are reordered.
259 .. versionadded:: 1.12.0
261 Parameters
262 ----------
263 m : array_like
264 Input array.
265 axis : None or int or tuple of ints, optional
266 Axis or axes along which to flip over. The default,
267 axis=None, will flip over all of the axes of the input array.
268 If axis is negative it counts from the last to the first axis.
270 If axis is a tuple of ints, flipping is performed on all of the axes
271 specified in the tuple.
273 .. versionchanged:: 1.15.0
274 None and tuples of axes are supported
276 Returns
277 -------
278 out : array_like
279 A view of `m` with the entries of axis reversed. Since a view is
280 returned, this operation is done in constant time.
282 See Also
283 --------
284 flipud : Flip an array vertically (axis=0).
285 fliplr : Flip an array horizontally (axis=1).
287 Notes
288 -----
289 flip(m, 0) is equivalent to flipud(m).
291 flip(m, 1) is equivalent to fliplr(m).
293 flip(m, n) corresponds to ``m[...,::-1,...]`` with ``::-1`` at position n.
295 flip(m) corresponds to ``m[::-1,::-1,...,::-1]`` with ``::-1`` at all
296 positions.
298 flip(m, (0, 1)) corresponds to ``m[::-1,::-1,...]`` with ``::-1`` at
299 position 0 and position 1.
301 Examples
302 --------
303 >>> A = np.arange(8).reshape((2,2,2))
304 >>> A
305 array([[[0, 1],
306 [2, 3]],
307 [[4, 5],
308 [6, 7]]])
309 >>> np.flip(A, 0)
310 array([[[4, 5],
311 [6, 7]],
312 [[0, 1],
313 [2, 3]]])
314 >>> np.flip(A, 1)
315 array([[[2, 3],
316 [0, 1]],
317 [[6, 7],
318 [4, 5]]])
319 >>> np.flip(A)
320 array([[[7, 6],
321 [5, 4]],
322 [[3, 2],
323 [1, 0]]])
324 >>> np.flip(A, (0, 2))
325 array([[[5, 4],
326 [7, 6]],
327 [[1, 0],
328 [3, 2]]])
329 >>> A = np.random.randn(3,4,5)
330 >>> np.all(np.flip(A,2) == A[:,:,::-1,...])
331 True
332 """
333 if not hasattr(m, 'ndim'):
334 m = asarray(m)
335 if axis is None:
336 indexer = (np.s_[::-1],) * m.ndim
337 else:
338 axis = _nx.normalize_axis_tuple(axis, m.ndim)
339 indexer = [np.s_[:]] * m.ndim
340 for ax in axis:
341 indexer[ax] = np.s_[::-1]
342 indexer = tuple(indexer)
343 return m[indexer]
346@set_module('numpy')
347def iterable(y):
348 """
349 Check whether or not an object can be iterated over.
351 Parameters
352 ----------
353 y : object
354 Input object.
356 Returns
357 -------
358 b : bool
359 Return ``True`` if the object has an iterator method or is a
360 sequence and ``False`` otherwise.
363 Examples
364 --------
365 >>> np.iterable([1, 2, 3])
366 True
367 >>> np.iterable(2)
368 False
370 Notes
371 -----
372 In most cases, the results of ``np.iterable(obj)`` are consistent with
373 ``isinstance(obj, collections.abc.Iterable)``. One notable exception is
374 the treatment of 0-dimensional arrays::
376 >>> from collections.abc import Iterable
377 >>> a = np.array(1.0) # 0-dimensional numpy array
378 >>> isinstance(a, Iterable)
379 True
380 >>> np.iterable(a)
381 False
383 """
384 try:
385 iter(y)
386 except TypeError:
387 return False
388 return True
391def _average_dispatcher(a, axis=None, weights=None, returned=None, *,
392 keepdims=None):
393 return (a, weights)
396@array_function_dispatch(_average_dispatcher)
397def average(a, axis=None, weights=None, returned=False, *,
398 keepdims=np._NoValue):
399 """
400 Compute the weighted average along the specified axis.
402 Parameters
403 ----------
404 a : array_like
405 Array containing data to be averaged. If `a` is not an array, a
406 conversion is attempted.
407 axis : None or int or tuple of ints, optional
408 Axis or axes along which to average `a`. The default,
409 axis=None, will average over all of the elements of the input array.
410 If axis is negative it counts from the last to the first axis.
412 .. versionadded:: 1.7.0
414 If axis is a tuple of ints, averaging is performed on all of the axes
415 specified in the tuple instead of a single axis or all the axes as
416 before.
417 weights : array_like, optional
418 An array of weights associated with the values in `a`. Each value in
419 `a` contributes to the average according to its associated weight.
420 The weights array can either be 1-D (in which case its length must be
421 the size of `a` along the given axis) or of the same shape as `a`.
422 If `weights=None`, then all data in `a` are assumed to have a
423 weight equal to one. The 1-D calculation is::
425 avg = sum(a * weights) / sum(weights)
427 The only constraint on `weights` is that `sum(weights)` must not be 0.
428 returned : bool, optional
429 Default is `False`. If `True`, the tuple (`average`, `sum_of_weights`)
430 is returned, otherwise only the average is returned.
431 If `weights=None`, `sum_of_weights` is equivalent to the number of
432 elements over which the average is taken.
433 keepdims : bool, optional
434 If this is set to True, the axes which are reduced are left
435 in the result as dimensions with size one. With this option,
436 the result will broadcast correctly against the original `a`.
437 *Note:* `keepdims` will not work with instances of `numpy.matrix`
438 or other classes whose methods do not support `keepdims`.
440 .. versionadded:: 1.23.0
442 Returns
443 -------
444 retval, [sum_of_weights] : array_type or double
445 Return the average along the specified axis. When `returned` is `True`,
446 return a tuple with the average as the first element and the sum
447 of the weights as the second element. `sum_of_weights` is of the
448 same type as `retval`. The result dtype follows a genereal pattern.
449 If `weights` is None, the result dtype will be that of `a` , or ``float64``
450 if `a` is integral. Otherwise, if `weights` is not None and `a` is non-
451 integral, the result type will be the type of lowest precision capable of
452 representing values of both `a` and `weights`. If `a` happens to be
453 integral, the previous rules still applies but the result dtype will
454 at least be ``float64``.
456 Raises
457 ------
458 ZeroDivisionError
459 When all weights along axis are zero. See `numpy.ma.average` for a
460 version robust to this type of error.
461 TypeError
462 When the length of 1D `weights` is not the same as the shape of `a`
463 along axis.
465 See Also
466 --------
467 mean
469 ma.average : average for masked arrays -- useful if your data contains
470 "missing" values
471 numpy.result_type : Returns the type that results from applying the
472 numpy type promotion rules to the arguments.
474 Examples
475 --------
476 >>> data = np.arange(1, 5)
477 >>> data
478 array([1, 2, 3, 4])
479 >>> np.average(data)
480 2.5
481 >>> np.average(np.arange(1, 11), weights=np.arange(10, 0, -1))
482 4.0
484 >>> data = np.arange(6).reshape((3, 2))
485 >>> data
486 array([[0, 1],
487 [2, 3],
488 [4, 5]])
489 >>> np.average(data, axis=1, weights=[1./4, 3./4])
490 array([0.75, 2.75, 4.75])
491 >>> np.average(data, weights=[1./4, 3./4])
492 Traceback (most recent call last):
493 ...
494 TypeError: Axis must be specified when shapes of a and weights differ.
496 >>> a = np.ones(5, dtype=np.float128)
497 >>> w = np.ones(5, dtype=np.complex64)
498 >>> avg = np.average(a, weights=w)
499 >>> print(avg.dtype)
500 complex256
502 With ``keepdims=True``, the following result has shape (3, 1).
504 >>> np.average(data, axis=1, keepdims=True)
505 array([[0.5],
506 [2.5],
507 [4.5]])
508 """
509 a = np.asanyarray(a)
511 if keepdims is np._NoValue:
512 # Don't pass on the keepdims argument if one wasn't given.
513 keepdims_kw = {}
514 else:
515 keepdims_kw = {'keepdims': keepdims}
517 if weights is None:
518 avg = a.mean(axis, **keepdims_kw)
519 avg_as_array = np.asanyarray(avg)
520 scl = avg_as_array.dtype.type(a.size/avg_as_array.size)
521 else:
522 wgt = np.asanyarray(weights)
524 if issubclass(a.dtype.type, (np.integer, np.bool_)):
525 result_dtype = np.result_type(a.dtype, wgt.dtype, 'f8')
526 else:
527 result_dtype = np.result_type(a.dtype, wgt.dtype)
529 # Sanity checks
530 if a.shape != wgt.shape:
531 if axis is None:
532 raise TypeError(
533 "Axis must be specified when shapes of a and weights "
534 "differ.")
535 if wgt.ndim != 1:
536 raise TypeError(
537 "1D weights expected when shapes of a and weights differ.")
538 if wgt.shape[0] != a.shape[axis]:
539 raise ValueError(
540 "Length of weights not compatible with specified axis.")
542 # setup wgt to broadcast along axis
543 wgt = np.broadcast_to(wgt, (a.ndim-1)*(1,) + wgt.shape)
544 wgt = wgt.swapaxes(-1, axis)
546 scl = wgt.sum(axis=axis, dtype=result_dtype, **keepdims_kw)
547 if np.any(scl == 0.0):
548 raise ZeroDivisionError(
549 "Weights sum to zero, can't be normalized")
551 avg = avg_as_array = np.multiply(a, wgt,
552 dtype=result_dtype).sum(axis, **keepdims_kw) / scl
554 if returned:
555 if scl.shape != avg_as_array.shape:
556 scl = np.broadcast_to(scl, avg_as_array.shape).copy()
557 return avg, scl
558 else:
559 return avg
562@set_module('numpy')
563def asarray_chkfinite(a, dtype=None, order=None):
564 """Convert the input to an array, checking for NaNs or Infs.
566 Parameters
567 ----------
568 a : array_like
569 Input data, in any form that can be converted to an array. This
570 includes lists, lists of tuples, tuples, tuples of tuples, tuples
571 of lists and ndarrays. Success requires no NaNs or Infs.
572 dtype : data-type, optional
573 By default, the data-type is inferred from the input data.
574 order : {'C', 'F', 'A', 'K'}, optional
575 Memory layout. 'A' and 'K' depend on the order of input array a.
576 'C' row-major (C-style),
577 'F' column-major (Fortran-style) memory representation.
578 'A' (any) means 'F' if `a` is Fortran contiguous, 'C' otherwise
579 'K' (keep) preserve input order
580 Defaults to 'C'.
582 Returns
583 -------
584 out : ndarray
585 Array interpretation of `a`. No copy is performed if the input
586 is already an ndarray. If `a` is a subclass of ndarray, a base
587 class ndarray is returned.
589 Raises
590 ------
591 ValueError
592 Raises ValueError if `a` contains NaN (Not a Number) or Inf (Infinity).
594 See Also
595 --------
596 asarray : Create and array.
597 asanyarray : Similar function which passes through subclasses.
598 ascontiguousarray : Convert input to a contiguous array.
599 asfarray : Convert input to a floating point ndarray.
600 asfortranarray : Convert input to an ndarray with column-major
601 memory order.
602 fromiter : Create an array from an iterator.
603 fromfunction : Construct an array by executing a function on grid
604 positions.
606 Examples
607 --------
608 Convert a list into an array. If all elements are finite
609 ``asarray_chkfinite`` is identical to ``asarray``.
611 >>> a = [1, 2]
612 >>> np.asarray_chkfinite(a, dtype=float)
613 array([1., 2.])
615 Raises ValueError if array_like contains Nans or Infs.
617 >>> a = [1, 2, np.inf]
618 >>> try:
619 ... np.asarray_chkfinite(a)
620 ... except ValueError:
621 ... print('ValueError')
622 ...
623 ValueError
625 """
626 a = asarray(a, dtype=dtype, order=order)
627 if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
628 raise ValueError(
629 "array must not contain infs or NaNs")
630 return a
633def _piecewise_dispatcher(x, condlist, funclist, *args, **kw):
634 yield x
635 # support the undocumented behavior of allowing scalars
636 if np.iterable(condlist):
637 yield from condlist
640@array_function_dispatch(_piecewise_dispatcher)
641def piecewise(x, condlist, funclist, *args, **kw):
642 """
643 Evaluate a piecewise-defined function.
645 Given a set of conditions and corresponding functions, evaluate each
646 function on the input data wherever its condition is true.
648 Parameters
649 ----------
650 x : ndarray or scalar
651 The input domain.
652 condlist : list of bool arrays or bool scalars
653 Each boolean array corresponds to a function in `funclist`. Wherever
654 `condlist[i]` is True, `funclist[i](x)` is used as the output value.
656 Each boolean array in `condlist` selects a piece of `x`,
657 and should therefore be of the same shape as `x`.
659 The length of `condlist` must correspond to that of `funclist`.
660 If one extra function is given, i.e. if
661 ``len(funclist) == len(condlist) + 1``, then that extra function
662 is the default value, used wherever all conditions are false.
663 funclist : list of callables, f(x,*args,**kw), or scalars
664 Each function is evaluated over `x` wherever its corresponding
665 condition is True. It should take a 1d array as input and give an 1d
666 array or a scalar value as output. If, instead of a callable,
667 a scalar is provided then a constant function (``lambda x: scalar``) is
668 assumed.
669 args : tuple, optional
670 Any further arguments given to `piecewise` are passed to the functions
671 upon execution, i.e., if called ``piecewise(..., ..., 1, 'a')``, then
672 each function is called as ``f(x, 1, 'a')``.
673 kw : dict, optional
674 Keyword arguments used in calling `piecewise` are passed to the
675 functions upon execution, i.e., if called
676 ``piecewise(..., ..., alpha=1)``, then each function is called as
677 ``f(x, alpha=1)``.
679 Returns
680 -------
681 out : ndarray
682 The output is the same shape and type as x and is found by
683 calling the functions in `funclist` on the appropriate portions of `x`,
684 as defined by the boolean arrays in `condlist`. Portions not covered
685 by any condition have a default value of 0.
688 See Also
689 --------
690 choose, select, where
692 Notes
693 -----
694 This is similar to choose or select, except that functions are
695 evaluated on elements of `x` that satisfy the corresponding condition from
696 `condlist`.
698 The result is::
700 |--
701 |funclist[0](x[condlist[0]])
702 out = |funclist[1](x[condlist[1]])
703 |...
704 |funclist[n2](x[condlist[n2]])
705 |--
707 Examples
708 --------
709 Define the sigma function, which is -1 for ``x < 0`` and +1 for ``x >= 0``.
711 >>> x = np.linspace(-2.5, 2.5, 6)
712 >>> np.piecewise(x, [x < 0, x >= 0], [-1, 1])
713 array([-1., -1., -1., 1., 1., 1.])
715 Define the absolute value, which is ``-x`` for ``x <0`` and ``x`` for
716 ``x >= 0``.
718 >>> np.piecewise(x, [x < 0, x >= 0], [lambda x: -x, lambda x: x])
719 array([2.5, 1.5, 0.5, 0.5, 1.5, 2.5])
721 Apply the same function to a scalar value.
723 >>> y = -2
724 >>> np.piecewise(y, [y < 0, y >= 0], [lambda x: -x, lambda x: x])
725 array(2)
727 """
728 x = asanyarray(x)
729 n2 = len(funclist)
731 # undocumented: single condition is promoted to a list of one condition
732 if isscalar(condlist) or (
733 not isinstance(condlist[0], (list, ndarray)) and x.ndim != 0):
734 condlist = [condlist]
736 condlist = asarray(condlist, dtype=bool)
737 n = len(condlist)
739 if n == n2 - 1: # compute the "otherwise" condition.
740 condelse = ~np.any(condlist, axis=0, keepdims=True)
741 condlist = np.concatenate([condlist, condelse], axis=0)
742 n += 1
743 elif n != n2:
744 raise ValueError(
745 "with {} condition(s), either {} or {} functions are expected"
746 .format(n, n, n+1)
747 )
749 y = zeros_like(x)
750 for cond, func in zip(condlist, funclist):
751 if not isinstance(func, collections.abc.Callable):
752 y[cond] = func
753 else:
754 vals = x[cond]
755 if vals.size > 0:
756 y[cond] = func(vals, *args, **kw)
758 return y
761def _select_dispatcher(condlist, choicelist, default=None):
762 yield from condlist
763 yield from choicelist
766@array_function_dispatch(_select_dispatcher)
767def select(condlist, choicelist, default=0):
768 """
769 Return an array drawn from elements in choicelist, depending on conditions.
771 Parameters
772 ----------
773 condlist : list of bool ndarrays
774 The list of conditions which determine from which array in `choicelist`
775 the output elements are taken. When multiple conditions are satisfied,
776 the first one encountered in `condlist` is used.
777 choicelist : list of ndarrays
778 The list of arrays from which the output elements are taken. It has
779 to be of the same length as `condlist`.
780 default : scalar, optional
781 The element inserted in `output` when all conditions evaluate to False.
783 Returns
784 -------
785 output : ndarray
786 The output at position m is the m-th element of the array in
787 `choicelist` where the m-th element of the corresponding array in
788 `condlist` is True.
790 See Also
791 --------
792 where : Return elements from one of two arrays depending on condition.
793 take, choose, compress, diag, diagonal
795 Examples
796 --------
797 >>> x = np.arange(6)
798 >>> condlist = [x<3, x>3]
799 >>> choicelist = [x, x**2]
800 >>> np.select(condlist, choicelist, 42)
801 array([ 0, 1, 2, 42, 16, 25])
803 >>> condlist = [x<=4, x>3]
804 >>> choicelist = [x, x**2]
805 >>> np.select(condlist, choicelist, 55)
806 array([ 0, 1, 2, 3, 4, 25])
808 """
809 # Check the size of condlist and choicelist are the same, or abort.
810 if len(condlist) != len(choicelist):
811 raise ValueError(
812 'list of cases must be same length as list of conditions')
814 # Now that the dtype is known, handle the deprecated select([], []) case
815 if len(condlist) == 0:
816 raise ValueError("select with an empty condition list is not possible")
818 choicelist = [np.asarray(choice) for choice in choicelist]
820 try:
821 intermediate_dtype = np.result_type(*choicelist)
822 except TypeError as e:
823 msg = f'Choicelist elements do not have a common dtype: {e}'
824 raise TypeError(msg) from None
825 default_array = np.asarray(default)
826 choicelist.append(default_array)
828 # need to get the result type before broadcasting for correct scalar
829 # behaviour
830 try:
831 dtype = np.result_type(intermediate_dtype, default_array)
832 except TypeError as e:
833 msg = f'Choicelists and default value do not have a common dtype: {e}'
834 raise TypeError(msg) from None
836 # Convert conditions to arrays and broadcast conditions and choices
837 # as the shape is needed for the result. Doing it separately optimizes
838 # for example when all choices are scalars.
839 condlist = np.broadcast_arrays(*condlist)
840 choicelist = np.broadcast_arrays(*choicelist)
842 # If cond array is not an ndarray in boolean format or scalar bool, abort.
843 for i, cond in enumerate(condlist):
844 if cond.dtype.type is not np.bool_:
845 raise TypeError(
846 'invalid entry {} in condlist: should be boolean ndarray'.format(i))
848 if choicelist[0].ndim == 0:
849 # This may be common, so avoid the call.
850 result_shape = condlist[0].shape
851 else:
852 result_shape = np.broadcast_arrays(condlist[0], choicelist[0])[0].shape
854 result = np.full(result_shape, choicelist[-1], dtype)
856 # Use np.copyto to burn each choicelist array onto result, using the
857 # corresponding condlist as a boolean mask. This is done in reverse
858 # order since the first choice should take precedence.
859 choicelist = choicelist[-2::-1]
860 condlist = condlist[::-1]
861 for choice, cond in zip(choicelist, condlist):
862 np.copyto(result, choice, where=cond)
864 return result
867def _copy_dispatcher(a, order=None, subok=None):
868 return (a,)
871@array_function_dispatch(_copy_dispatcher)
872def copy(a, order='K', subok=False):
873 """
874 Return an array copy of the given object.
876 Parameters
877 ----------
878 a : array_like
879 Input data.
880 order : {'C', 'F', 'A', 'K'}, optional
881 Controls the memory layout of the copy. 'C' means C-order,
882 'F' means F-order, 'A' means 'F' if `a` is Fortran contiguous,
883 'C' otherwise. 'K' means match the layout of `a` as closely
884 as possible. (Note that this function and :meth:`ndarray.copy` are very
885 similar, but have different default values for their order=
886 arguments.)
887 subok : bool, optional
888 If True, then sub-classes will be passed-through, otherwise the
889 returned array will be forced to be a base-class array (defaults to False).
891 .. versionadded:: 1.19.0
893 Returns
894 -------
895 arr : ndarray
896 Array interpretation of `a`.
898 See Also
899 --------
900 ndarray.copy : Preferred method for creating an array copy
902 Notes
903 -----
904 This is equivalent to:
906 >>> np.array(a, copy=True) #doctest: +SKIP
908 Examples
909 --------
910 Create an array x, with a reference y and a copy z:
912 >>> x = np.array([1, 2, 3])
913 >>> y = x
914 >>> z = np.copy(x)
916 Note that, when we modify x, y changes, but not z:
918 >>> x[0] = 10
919 >>> x[0] == y[0]
920 True
921 >>> x[0] == z[0]
922 False
924 Note that, np.copy clears previously set WRITEABLE=False flag.
926 >>> a = np.array([1, 2, 3])
927 >>> a.flags["WRITEABLE"] = False
928 >>> b = np.copy(a)
929 >>> b.flags["WRITEABLE"]
930 True
931 >>> b[0] = 3
932 >>> b
933 array([3, 2, 3])
935 Note that np.copy is a shallow copy and will not copy object
936 elements within arrays. This is mainly important for arrays
937 containing Python objects. The new array will contain the
938 same object which may lead to surprises if that object can
939 be modified (is mutable):
941 >>> a = np.array([1, 'm', [2, 3, 4]], dtype=object)
942 >>> b = np.copy(a)
943 >>> b[2][0] = 10
944 >>> a
945 array([1, 'm', list([10, 3, 4])], dtype=object)
947 To ensure all elements within an ``object`` array are copied,
948 use `copy.deepcopy`:
950 >>> import copy
951 >>> a = np.array([1, 'm', [2, 3, 4]], dtype=object)
952 >>> c = copy.deepcopy(a)
953 >>> c[2][0] = 10
954 >>> c
955 array([1, 'm', list([10, 3, 4])], dtype=object)
956 >>> a
957 array([1, 'm', list([2, 3, 4])], dtype=object)
959 """
960 return array(a, order=order, subok=subok, copy=True)
962# Basic operations
965def _gradient_dispatcher(f, *varargs, axis=None, edge_order=None):
966 yield f
967 yield from varargs
970@array_function_dispatch(_gradient_dispatcher)
971def gradient(f, *varargs, axis=None, edge_order=1):
972 """
973 Return the gradient of an N-dimensional array.
975 The gradient is computed using second order accurate central differences
976 in the interior points and either first or second order accurate one-sides
977 (forward or backwards) differences at the boundaries.
978 The returned gradient hence has the same shape as the input array.
980 Parameters
981 ----------
982 f : array_like
983 An N-dimensional array containing samples of a scalar function.
984 varargs : list of scalar or array, optional
985 Spacing between f values. Default unitary spacing for all dimensions.
986 Spacing can be specified using:
988 1. single scalar to specify a sample distance for all dimensions.
989 2. N scalars to specify a constant sample distance for each dimension.
990 i.e. `dx`, `dy`, `dz`, ...
991 3. N arrays to specify the coordinates of the values along each
992 dimension of F. The length of the array must match the size of
993 the corresponding dimension
994 4. Any combination of N scalars/arrays with the meaning of 2. and 3.
996 If `axis` is given, the number of varargs must equal the number of axes.
997 Default: 1.
999 edge_order : {1, 2}, optional
1000 Gradient is calculated using N-th order accurate differences
1001 at the boundaries. Default: 1.
1003 .. versionadded:: 1.9.1
1005 axis : None or int or tuple of ints, optional
1006 Gradient is calculated only along the given axis or axes
1007 The default (axis = None) is to calculate the gradient for all the axes
1008 of the input array. axis may be negative, in which case it counts from
1009 the last to the first axis.
1011 .. versionadded:: 1.11.0
1013 Returns
1014 -------
1015 gradient : ndarray or list of ndarray
1016 A list of ndarrays (or a single ndarray if there is only one dimension)
1017 corresponding to the derivatives of f with respect to each dimension.
1018 Each derivative has the same shape as f.
1020 Examples
1021 --------
1022 >>> f = np.array([1, 2, 4, 7, 11, 16], dtype=float)
1023 >>> np.gradient(f)
1024 array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
1025 >>> np.gradient(f, 2)
1026 array([0.5 , 0.75, 1.25, 1.75, 2.25, 2.5 ])
1028 Spacing can be also specified with an array that represents the coordinates
1029 of the values F along the dimensions.
1030 For instance a uniform spacing:
1032 >>> x = np.arange(f.size)
1033 >>> np.gradient(f, x)
1034 array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
1036 Or a non uniform one:
1038 >>> x = np.array([0., 1., 1.5, 3.5, 4., 6.], dtype=float)
1039 >>> np.gradient(f, x)
1040 array([1. , 3. , 3.5, 6.7, 6.9, 2.5])
1042 For two dimensional arrays, the return will be two arrays ordered by
1043 axis. In this example the first array stands for the gradient in
1044 rows and the second one in columns direction:
1046 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float))
1047 [array([[ 2., 2., -1.],
1048 [ 2., 2., -1.]]), array([[1. , 2.5, 4. ],
1049 [1. , 1. , 1. ]])]
1051 In this example the spacing is also specified:
1052 uniform for axis=0 and non uniform for axis=1
1054 >>> dx = 2.
1055 >>> y = [1., 1.5, 3.5]
1056 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), dx, y)
1057 [array([[ 1. , 1. , -0.5],
1058 [ 1. , 1. , -0.5]]), array([[2. , 2. , 2. ],
1059 [2. , 1.7, 0.5]])]
1061 It is possible to specify how boundaries are treated using `edge_order`
1063 >>> x = np.array([0, 1, 2, 3, 4])
1064 >>> f = x**2
1065 >>> np.gradient(f, edge_order=1)
1066 array([1., 2., 4., 6., 7.])
1067 >>> np.gradient(f, edge_order=2)
1068 array([0., 2., 4., 6., 8.])
1070 The `axis` keyword can be used to specify a subset of axes of which the
1071 gradient is calculated
1073 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), axis=0)
1074 array([[ 2., 2., -1.],
1075 [ 2., 2., -1.]])
1077 Notes
1078 -----
1079 Assuming that :math:`f\\in C^{3}` (i.e., :math:`f` has at least 3 continuous
1080 derivatives) and let :math:`h_{*}` be a non-homogeneous stepsize, we
1081 minimize the "consistency error" :math:`\\eta_{i}` between the true gradient
1082 and its estimate from a linear combination of the neighboring grid-points:
1084 .. math::
1086 \\eta_{i} = f_{i}^{\\left(1\\right)} -
1087 \\left[ \\alpha f\\left(x_{i}\\right) +
1088 \\beta f\\left(x_{i} + h_{d}\\right) +
1089 \\gamma f\\left(x_{i}-h_{s}\\right)
1090 \\right]
1092 By substituting :math:`f(x_{i} + h_{d})` and :math:`f(x_{i} - h_{s})`
1093 with their Taylor series expansion, this translates into solving
1094 the following the linear system:
1096 .. math::
1098 \\left\\{
1099 \\begin{array}{r}
1100 \\alpha+\\beta+\\gamma=0 \\\\
1101 \\beta h_{d}-\\gamma h_{s}=1 \\\\
1102 \\beta h_{d}^{2}+\\gamma h_{s}^{2}=0
1103 \\end{array}
1104 \\right.
1106 The resulting approximation of :math:`f_{i}^{(1)}` is the following:
1108 .. math::
1110 \\hat f_{i}^{(1)} =
1111 \\frac{
1112 h_{s}^{2}f\\left(x_{i} + h_{d}\\right)
1113 + \\left(h_{d}^{2} - h_{s}^{2}\\right)f\\left(x_{i}\\right)
1114 - h_{d}^{2}f\\left(x_{i}-h_{s}\\right)}
1115 { h_{s}h_{d}\\left(h_{d} + h_{s}\\right)}
1116 + \\mathcal{O}\\left(\\frac{h_{d}h_{s}^{2}
1117 + h_{s}h_{d}^{2}}{h_{d}
1118 + h_{s}}\\right)
1120 It is worth noting that if :math:`h_{s}=h_{d}`
1121 (i.e., data are evenly spaced)
1122 we find the standard second order approximation:
1124 .. math::
1126 \\hat f_{i}^{(1)}=
1127 \\frac{f\\left(x_{i+1}\\right) - f\\left(x_{i-1}\\right)}{2h}
1128 + \\mathcal{O}\\left(h^{2}\\right)
1130 With a similar procedure the forward/backward approximations used for
1131 boundaries can be derived.
1133 References
1134 ----------
1135 .. [1] Quarteroni A., Sacco R., Saleri F. (2007) Numerical Mathematics
1136 (Texts in Applied Mathematics). New York: Springer.
1137 .. [2] Durran D. R. (1999) Numerical Methods for Wave Equations
1138 in Geophysical Fluid Dynamics. New York: Springer.
1139 .. [3] Fornberg B. (1988) Generation of Finite Difference Formulas on
1140 Arbitrarily Spaced Grids,
1141 Mathematics of Computation 51, no. 184 : 699-706.
1142 `PDF <http://www.ams.org/journals/mcom/1988-51-184/
1143 S0025-5718-1988-0935077-0/S0025-5718-1988-0935077-0.pdf>`_.
1144 """
1145 f = np.asanyarray(f)
1146 N = f.ndim # number of dimensions
1148 if axis is None:
1149 axes = tuple(range(N))
1150 else:
1151 axes = _nx.normalize_axis_tuple(axis, N)
1153 len_axes = len(axes)
1154 n = len(varargs)
1155 if n == 0:
1156 # no spacing argument - use 1 in all axes
1157 dx = [1.0] * len_axes
1158 elif n == 1 and np.ndim(varargs[0]) == 0:
1159 # single scalar for all axes
1160 dx = varargs * len_axes
1161 elif n == len_axes:
1162 # scalar or 1d array for each axis
1163 dx = list(varargs)
1164 for i, distances in enumerate(dx):
1165 distances = np.asanyarray(distances)
1166 if distances.ndim == 0:
1167 continue
1168 elif distances.ndim != 1:
1169 raise ValueError("distances must be either scalars or 1d")
1170 if len(distances) != f.shape[axes[i]]:
1171 raise ValueError("when 1d, distances must match "
1172 "the length of the corresponding dimension")
1173 if np.issubdtype(distances.dtype, np.integer):
1174 # Convert numpy integer types to float64 to avoid modular
1175 # arithmetic in np.diff(distances).
1176 distances = distances.astype(np.float64)
1177 diffx = np.diff(distances)
1178 # if distances are constant reduce to the scalar case
1179 # since it brings a consistent speedup
1180 if (diffx == diffx[0]).all():
1181 diffx = diffx[0]
1182 dx[i] = diffx
1183 else:
1184 raise TypeError("invalid number of arguments")
1186 if edge_order > 2:
1187 raise ValueError("'edge_order' greater than 2 not supported")
1189 # use central differences on interior and one-sided differences on the
1190 # endpoints. This preserves second order-accuracy over the full domain.
1192 outvals = []
1194 # create slice objects --- initially all are [:, :, ..., :]
1195 slice1 = [slice(None)]*N
1196 slice2 = [slice(None)]*N
1197 slice3 = [slice(None)]*N
1198 slice4 = [slice(None)]*N
1200 otype = f.dtype
1201 if otype.type is np.datetime64:
1202 # the timedelta dtype with the same unit information
1203 otype = np.dtype(otype.name.replace('datetime', 'timedelta'))
1204 # view as timedelta to allow addition
1205 f = f.view(otype)
1206 elif otype.type is np.timedelta64:
1207 pass
1208 elif np.issubdtype(otype, np.inexact):
1209 pass
1210 else:
1211 # All other types convert to floating point.
1212 # First check if f is a numpy integer type; if so, convert f to float64
1213 # to avoid modular arithmetic when computing the changes in f.
1214 if np.issubdtype(otype, np.integer):
1215 f = f.astype(np.float64)
1216 otype = np.float64
1218 for axis, ax_dx in zip(axes, dx):
1219 if f.shape[axis] < edge_order + 1:
1220 raise ValueError(
1221 "Shape of array too small to calculate a numerical gradient, "
1222 "at least (edge_order + 1) elements are required.")
1223 # result allocation
1224 out = np.empty_like(f, dtype=otype)
1226 # spacing for the current axis
1227 uniform_spacing = np.ndim(ax_dx) == 0
1229 # Numerical differentiation: 2nd order interior
1230 slice1[axis] = slice(1, -1)
1231 slice2[axis] = slice(None, -2)
1232 slice3[axis] = slice(1, -1)
1233 slice4[axis] = slice(2, None)
1235 if uniform_spacing:
1236 out[tuple(slice1)] = (f[tuple(slice4)] - f[tuple(slice2)]) / (2. * ax_dx)
1237 else:
1238 dx1 = ax_dx[0:-1]
1239 dx2 = ax_dx[1:]
1240 a = -(dx2)/(dx1 * (dx1 + dx2))
1241 b = (dx2 - dx1) / (dx1 * dx2)
1242 c = dx1 / (dx2 * (dx1 + dx2))
1243 # fix the shape for broadcasting
1244 shape = np.ones(N, dtype=int)
1245 shape[axis] = -1
1246 a.shape = b.shape = c.shape = shape
1247 # 1D equivalent -- out[1:-1] = a * f[:-2] + b * f[1:-1] + c * f[2:]
1248 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
1250 # Numerical differentiation: 1st order edges
1251 if edge_order == 1:
1252 slice1[axis] = 0
1253 slice2[axis] = 1
1254 slice3[axis] = 0
1255 dx_0 = ax_dx if uniform_spacing else ax_dx[0]
1256 # 1D equivalent -- out[0] = (f[1] - f[0]) / (x[1] - x[0])
1257 out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0
1259 slice1[axis] = -1
1260 slice2[axis] = -1
1261 slice3[axis] = -2
1262 dx_n = ax_dx if uniform_spacing else ax_dx[-1]
1263 # 1D equivalent -- out[-1] = (f[-1] - f[-2]) / (x[-1] - x[-2])
1264 out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_n
1266 # Numerical differentiation: 2nd order edges
1267 else:
1268 slice1[axis] = 0
1269 slice2[axis] = 0
1270 slice3[axis] = 1
1271 slice4[axis] = 2
1272 if uniform_spacing:
1273 a = -1.5 / ax_dx
1274 b = 2. / ax_dx
1275 c = -0.5 / ax_dx
1276 else:
1277 dx1 = ax_dx[0]
1278 dx2 = ax_dx[1]
1279 a = -(2. * dx1 + dx2)/(dx1 * (dx1 + dx2))
1280 b = (dx1 + dx2) / (dx1 * dx2)
1281 c = - dx1 / (dx2 * (dx1 + dx2))
1282 # 1D equivalent -- out[0] = a * f[0] + b * f[1] + c * f[2]
1283 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
1285 slice1[axis] = -1
1286 slice2[axis] = -3
1287 slice3[axis] = -2
1288 slice4[axis] = -1
1289 if uniform_spacing:
1290 a = 0.5 / ax_dx
1291 b = -2. / ax_dx
1292 c = 1.5 / ax_dx
1293 else:
1294 dx1 = ax_dx[-2]
1295 dx2 = ax_dx[-1]
1296 a = (dx2) / (dx1 * (dx1 + dx2))
1297 b = - (dx2 + dx1) / (dx1 * dx2)
1298 c = (2. * dx2 + dx1) / (dx2 * (dx1 + dx2))
1299 # 1D equivalent -- out[-1] = a * f[-3] + b * f[-2] + c * f[-1]
1300 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
1302 outvals.append(out)
1304 # reset the slice object in this dimension to ":"
1305 slice1[axis] = slice(None)
1306 slice2[axis] = slice(None)
1307 slice3[axis] = slice(None)
1308 slice4[axis] = slice(None)
1310 if len_axes == 1:
1311 return outvals[0]
1312 else:
1313 return outvals
1316def _diff_dispatcher(a, n=None, axis=None, prepend=None, append=None):
1317 return (a, prepend, append)
1320@array_function_dispatch(_diff_dispatcher)
1321def diff(a, n=1, axis=-1, prepend=np._NoValue, append=np._NoValue):
1322 """
1323 Calculate the n-th discrete difference along the given axis.
1325 The first difference is given by ``out[i] = a[i+1] - a[i]`` along
1326 the given axis, higher differences are calculated by using `diff`
1327 recursively.
1329 Parameters
1330 ----------
1331 a : array_like
1332 Input array
1333 n : int, optional
1334 The number of times values are differenced. If zero, the input
1335 is returned as-is.
1336 axis : int, optional
1337 The axis along which the difference is taken, default is the
1338 last axis.
1339 prepend, append : array_like, optional
1340 Values to prepend or append to `a` along axis prior to
1341 performing the difference. Scalar values are expanded to
1342 arrays with length 1 in the direction of axis and the shape
1343 of the input array in along all other axes. Otherwise the
1344 dimension and shape must match `a` except along axis.
1346 .. versionadded:: 1.16.0
1348 Returns
1349 -------
1350 diff : ndarray
1351 The n-th differences. The shape of the output is the same as `a`
1352 except along `axis` where the dimension is smaller by `n`. The
1353 type of the output is the same as the type of the difference
1354 between any two elements of `a`. This is the same as the type of
1355 `a` in most cases. A notable exception is `datetime64`, which
1356 results in a `timedelta64` output array.
1358 See Also
1359 --------
1360 gradient, ediff1d, cumsum
1362 Notes
1363 -----
1364 Type is preserved for boolean arrays, so the result will contain
1365 `False` when consecutive elements are the same and `True` when they
1366 differ.
1368 For unsigned integer arrays, the results will also be unsigned. This
1369 should not be surprising, as the result is consistent with
1370 calculating the difference directly:
1372 >>> u8_arr = np.array([1, 0], dtype=np.uint8)
1373 >>> np.diff(u8_arr)
1374 array([255], dtype=uint8)
1375 >>> u8_arr[1,...] - u8_arr[0,...]
1376 255
1378 If this is not desirable, then the array should be cast to a larger
1379 integer type first:
1381 >>> i16_arr = u8_arr.astype(np.int16)
1382 >>> np.diff(i16_arr)
1383 array([-1], dtype=int16)
1385 Examples
1386 --------
1387 >>> x = np.array([1, 2, 4, 7, 0])
1388 >>> np.diff(x)
1389 array([ 1, 2, 3, -7])
1390 >>> np.diff(x, n=2)
1391 array([ 1, 1, -10])
1393 >>> x = np.array([[1, 3, 6, 10], [0, 5, 6, 8]])
1394 >>> np.diff(x)
1395 array([[2, 3, 4],
1396 [5, 1, 2]])
1397 >>> np.diff(x, axis=0)
1398 array([[-1, 2, 0, -2]])
1400 >>> x = np.arange('1066-10-13', '1066-10-16', dtype=np.datetime64)
1401 >>> np.diff(x)
1402 array([1, 1], dtype='timedelta64[D]')
1404 """
1405 if n == 0:
1406 return a
1407 if n < 0:
1408 raise ValueError(
1409 "order must be non-negative but got " + repr(n))
1411 a = asanyarray(a)
1412 nd = a.ndim
1413 if nd == 0:
1414 raise ValueError("diff requires input that is at least one dimensional")
1415 axis = normalize_axis_index(axis, nd)
1417 combined = []
1418 if prepend is not np._NoValue:
1419 prepend = np.asanyarray(prepend)
1420 if prepend.ndim == 0:
1421 shape = list(a.shape)
1422 shape[axis] = 1
1423 prepend = np.broadcast_to(prepend, tuple(shape))
1424 combined.append(prepend)
1426 combined.append(a)
1428 if append is not np._NoValue:
1429 append = np.asanyarray(append)
1430 if append.ndim == 0:
1431 shape = list(a.shape)
1432 shape[axis] = 1
1433 append = np.broadcast_to(append, tuple(shape))
1434 combined.append(append)
1436 if len(combined) > 1:
1437 a = np.concatenate(combined, axis)
1439 slice1 = [slice(None)] * nd
1440 slice2 = [slice(None)] * nd
1441 slice1[axis] = slice(1, None)
1442 slice2[axis] = slice(None, -1)
1443 slice1 = tuple(slice1)
1444 slice2 = tuple(slice2)
1446 op = not_equal if a.dtype == np.bool_ else subtract
1447 for _ in range(n):
1448 a = op(a[slice1], a[slice2])
1450 return a
1453def _interp_dispatcher(x, xp, fp, left=None, right=None, period=None):
1454 return (x, xp, fp)
1457@array_function_dispatch(_interp_dispatcher)
1458def interp(x, xp, fp, left=None, right=None, period=None):
1459 """
1460 One-dimensional linear interpolation for monotonically increasing sample points.
1462 Returns the one-dimensional piecewise linear interpolant to a function
1463 with given discrete data points (`xp`, `fp`), evaluated at `x`.
1465 Parameters
1466 ----------
1467 x : array_like
1468 The x-coordinates at which to evaluate the interpolated values.
1470 xp : 1-D sequence of floats
1471 The x-coordinates of the data points, must be increasing if argument
1472 `period` is not specified. Otherwise, `xp` is internally sorted after
1473 normalizing the periodic boundaries with ``xp = xp % period``.
1475 fp : 1-D sequence of float or complex
1476 The y-coordinates of the data points, same length as `xp`.
1478 left : optional float or complex corresponding to fp
1479 Value to return for `x < xp[0]`, default is `fp[0]`.
1481 right : optional float or complex corresponding to fp
1482 Value to return for `x > xp[-1]`, default is `fp[-1]`.
1484 period : None or float, optional
1485 A period for the x-coordinates. This parameter allows the proper
1486 interpolation of angular x-coordinates. Parameters `left` and `right`
1487 are ignored if `period` is specified.
1489 .. versionadded:: 1.10.0
1491 Returns
1492 -------
1493 y : float or complex (corresponding to fp) or ndarray
1494 The interpolated values, same shape as `x`.
1496 Raises
1497 ------
1498 ValueError
1499 If `xp` and `fp` have different length
1500 If `xp` or `fp` are not 1-D sequences
1501 If `period == 0`
1503 See Also
1504 --------
1505 scipy.interpolate
1507 Warnings
1508 --------
1509 The x-coordinate sequence is expected to be increasing, but this is not
1510 explicitly enforced. However, if the sequence `xp` is non-increasing,
1511 interpolation results are meaningless.
1513 Note that, since NaN is unsortable, `xp` also cannot contain NaNs.
1515 A simple check for `xp` being strictly increasing is::
1517 np.all(np.diff(xp) > 0)
1519 Examples
1520 --------
1521 >>> xp = [1, 2, 3]
1522 >>> fp = [3, 2, 0]
1523 >>> np.interp(2.5, xp, fp)
1524 1.0
1525 >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
1526 array([3. , 3. , 2.5 , 0.56, 0. ])
1527 >>> UNDEF = -99.0
1528 >>> np.interp(3.14, xp, fp, right=UNDEF)
1529 -99.0
1531 Plot an interpolant to the sine function:
1533 >>> x = np.linspace(0, 2*np.pi, 10)
1534 >>> y = np.sin(x)
1535 >>> xvals = np.linspace(0, 2*np.pi, 50)
1536 >>> yinterp = np.interp(xvals, x, y)
1537 >>> import matplotlib.pyplot as plt
1538 >>> plt.plot(x, y, 'o')
1539 [<matplotlib.lines.Line2D object at 0x...>]
1540 >>> plt.plot(xvals, yinterp, '-x')
1541 [<matplotlib.lines.Line2D object at 0x...>]
1542 >>> plt.show()
1544 Interpolation with periodic x-coordinates:
1546 >>> x = [-180, -170, -185, 185, -10, -5, 0, 365]
1547 >>> xp = [190, -190, 350, -350]
1548 >>> fp = [5, 10, 3, 4]
1549 >>> np.interp(x, xp, fp, period=360)
1550 array([7.5 , 5. , 8.75, 6.25, 3. , 3.25, 3.5 , 3.75])
1552 Complex interpolation:
1554 >>> x = [1.5, 4.0]
1555 >>> xp = [2,3,5]
1556 >>> fp = [1.0j, 0, 2+3j]
1557 >>> np.interp(x, xp, fp)
1558 array([0.+1.j , 1.+1.5j])
1560 """
1562 fp = np.asarray(fp)
1564 if np.iscomplexobj(fp):
1565 interp_func = compiled_interp_complex
1566 input_dtype = np.complex128
1567 else:
1568 interp_func = compiled_interp
1569 input_dtype = np.float64
1571 if period is not None:
1572 if period == 0:
1573 raise ValueError("period must be a non-zero value")
1574 period = abs(period)
1575 left = None
1576 right = None
1578 x = np.asarray(x, dtype=np.float64)
1579 xp = np.asarray(xp, dtype=np.float64)
1580 fp = np.asarray(fp, dtype=input_dtype)
1582 if xp.ndim != 1 or fp.ndim != 1:
1583 raise ValueError("Data points must be 1-D sequences")
1584 if xp.shape[0] != fp.shape[0]:
1585 raise ValueError("fp and xp are not of the same length")
1586 # normalizing periodic boundaries
1587 x = x % period
1588 xp = xp % period
1589 asort_xp = np.argsort(xp)
1590 xp = xp[asort_xp]
1591 fp = fp[asort_xp]
1592 xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period))
1593 fp = np.concatenate((fp[-1:], fp, fp[0:1]))
1595 return interp_func(x, xp, fp, left, right)
1598def _angle_dispatcher(z, deg=None):
1599 return (z,)
1602@array_function_dispatch(_angle_dispatcher)
1603def angle(z, deg=False):
1604 """
1605 Return the angle of the complex argument.
1607 Parameters
1608 ----------
1609 z : array_like
1610 A complex number or sequence of complex numbers.
1611 deg : bool, optional
1612 Return angle in degrees if True, radians if False (default).
1614 Returns
1615 -------
1616 angle : ndarray or scalar
1617 The counterclockwise angle from the positive real axis on the complex
1618 plane in the range ``(-pi, pi]``, with dtype as numpy.float64.
1620 .. versionchanged:: 1.16.0
1621 This function works on subclasses of ndarray like `ma.array`.
1623 See Also
1624 --------
1625 arctan2
1626 absolute
1628 Notes
1629 -----
1630 Although the angle of the complex number 0 is undefined, ``numpy.angle(0)``
1631 returns the value 0.
1633 Examples
1634 --------
1635 >>> np.angle([1.0, 1.0j, 1+1j]) # in radians
1636 array([ 0. , 1.57079633, 0.78539816]) # may vary
1637 >>> np.angle(1+1j, deg=True) # in degrees
1638 45.0
1640 """
1641 z = asanyarray(z)
1642 if issubclass(z.dtype.type, _nx.complexfloating):
1643 zimag = z.imag
1644 zreal = z.real
1645 else:
1646 zimag = 0
1647 zreal = z
1649 a = arctan2(zimag, zreal)
1650 if deg:
1651 a *= 180/pi
1652 return a
1655def _unwrap_dispatcher(p, discont=None, axis=None, *, period=None):
1656 return (p,)
1659@array_function_dispatch(_unwrap_dispatcher)
1660def unwrap(p, discont=None, axis=-1, *, period=2*pi):
1661 r"""
1662 Unwrap by taking the complement of large deltas with respect to the period.
1664 This unwraps a signal `p` by changing elements which have an absolute
1665 difference from their predecessor of more than ``max(discont, period/2)``
1666 to their `period`-complementary values.
1668 For the default case where `period` is :math:`2\pi` and `discont` is
1669 :math:`\pi`, this unwraps a radian phase `p` such that adjacent differences
1670 are never greater than :math:`\pi` by adding :math:`2k\pi` for some
1671 integer :math:`k`.
1673 Parameters
1674 ----------
1675 p : array_like
1676 Input array.
1677 discont : float, optional
1678 Maximum discontinuity between values, default is ``period/2``.
1679 Values below ``period/2`` are treated as if they were ``period/2``.
1680 To have an effect different from the default, `discont` should be
1681 larger than ``period/2``.
1682 axis : int, optional
1683 Axis along which unwrap will operate, default is the last axis.
1684 period : float, optional
1685 Size of the range over which the input wraps. By default, it is
1686 ``2 pi``.
1688 .. versionadded:: 1.21.0
1690 Returns
1691 -------
1692 out : ndarray
1693 Output array.
1695 See Also
1696 --------
1697 rad2deg, deg2rad
1699 Notes
1700 -----
1701 If the discontinuity in `p` is smaller than ``period/2``,
1702 but larger than `discont`, no unwrapping is done because taking
1703 the complement would only make the discontinuity larger.
1705 Examples
1706 --------
1707 >>> phase = np.linspace(0, np.pi, num=5)
1708 >>> phase[3:] += np.pi
1709 >>> phase
1710 array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531]) # may vary
1711 >>> np.unwrap(phase)
1712 array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ]) # may vary
1713 >>> np.unwrap([0, 1, 2, -1, 0], period=4)
1714 array([0, 1, 2, 3, 4])
1715 >>> np.unwrap([ 1, 2, 3, 4, 5, 6, 1, 2, 3], period=6)
1716 array([1, 2, 3, 4, 5, 6, 7, 8, 9])
1717 >>> np.unwrap([2, 3, 4, 5, 2, 3, 4, 5], period=4)
1718 array([2, 3, 4, 5, 6, 7, 8, 9])
1719 >>> phase_deg = np.mod(np.linspace(0 ,720, 19), 360) - 180
1720 >>> np.unwrap(phase_deg, period=360)
1721 array([-180., -140., -100., -60., -20., 20., 60., 100., 140.,
1722 180., 220., 260., 300., 340., 380., 420., 460., 500.,
1723 540.])
1724 """
1725 p = asarray(p)
1726 nd = p.ndim
1727 dd = diff(p, axis=axis)
1728 if discont is None:
1729 discont = period/2
1730 slice1 = [slice(None, None)]*nd # full slices
1731 slice1[axis] = slice(1, None)
1732 slice1 = tuple(slice1)
1733 dtype = np.result_type(dd, period)
1734 if _nx.issubdtype(dtype, _nx.integer):
1735 interval_high, rem = divmod(period, 2)
1736 boundary_ambiguous = rem == 0
1737 else:
1738 interval_high = period / 2
1739 boundary_ambiguous = True
1740 interval_low = -interval_high
1741 ddmod = mod(dd - interval_low, period) + interval_low
1742 if boundary_ambiguous:
1743 # for `mask = (abs(dd) == period/2)`, the above line made
1744 # `ddmod[mask] == -period/2`. correct these such that
1745 # `ddmod[mask] == sign(dd[mask])*period/2`.
1746 _nx.copyto(ddmod, interval_high,
1747 where=(ddmod == interval_low) & (dd > 0))
1748 ph_correct = ddmod - dd
1749 _nx.copyto(ph_correct, 0, where=abs(dd) < discont)
1750 up = array(p, copy=True, dtype=dtype)
1751 up[slice1] = p[slice1] + ph_correct.cumsum(axis)
1752 return up
1755def _sort_complex(a):
1756 return (a,)
1759@array_function_dispatch(_sort_complex)
1760def sort_complex(a):
1761 """
1762 Sort a complex array using the real part first, then the imaginary part.
1764 Parameters
1765 ----------
1766 a : array_like
1767 Input array
1769 Returns
1770 -------
1771 out : complex ndarray
1772 Always returns a sorted complex array.
1774 Examples
1775 --------
1776 >>> np.sort_complex([5, 3, 6, 2, 1])
1777 array([1.+0.j, 2.+0.j, 3.+0.j, 5.+0.j, 6.+0.j])
1779 >>> np.sort_complex([1 + 2j, 2 - 1j, 3 - 2j, 3 - 3j, 3 + 5j])
1780 array([1.+2.j, 2.-1.j, 3.-3.j, 3.-2.j, 3.+5.j])
1782 """
1783 b = array(a, copy=True)
1784 b.sort()
1785 if not issubclass(b.dtype.type, _nx.complexfloating):
1786 if b.dtype.char in 'bhBH':
1787 return b.astype('F')
1788 elif b.dtype.char == 'g':
1789 return b.astype('G')
1790 else:
1791 return b.astype('D')
1792 else:
1793 return b
1796def _trim_zeros(filt, trim=None):
1797 return (filt,)
1800@array_function_dispatch(_trim_zeros)
1801def trim_zeros(filt, trim='fb'):
1802 """
1803 Trim the leading and/or trailing zeros from a 1-D array or sequence.
1805 Parameters
1806 ----------
1807 filt : 1-D array or sequence
1808 Input array.
1809 trim : str, optional
1810 A string with 'f' representing trim from front and 'b' to trim from
1811 back. Default is 'fb', trim zeros from both front and back of the
1812 array.
1814 Returns
1815 -------
1816 trimmed : 1-D array or sequence
1817 The result of trimming the input. The input data type is preserved.
1819 Examples
1820 --------
1821 >>> a = np.array((0, 0, 0, 1, 2, 3, 0, 2, 1, 0))
1822 >>> np.trim_zeros(a)
1823 array([1, 2, 3, 0, 2, 1])
1825 >>> np.trim_zeros(a, 'b')
1826 array([0, 0, 0, ..., 0, 2, 1])
1828 The input data type is preserved, list/tuple in means list/tuple out.
1830 >>> np.trim_zeros([0, 1, 2, 0])
1831 [1, 2]
1833 """
1835 first = 0
1836 trim = trim.upper()
1837 if 'F' in trim:
1838 for i in filt:
1839 if i != 0.:
1840 break
1841 else:
1842 first = first + 1
1843 last = len(filt)
1844 if 'B' in trim:
1845 for i in filt[::-1]:
1846 if i != 0.:
1847 break
1848 else:
1849 last = last - 1
1850 return filt[first:last]
1853def _extract_dispatcher(condition, arr):
1854 return (condition, arr)
1857@array_function_dispatch(_extract_dispatcher)
1858def extract(condition, arr):
1859 """
1860 Return the elements of an array that satisfy some condition.
1862 This is equivalent to ``np.compress(ravel(condition), ravel(arr))``. If
1863 `condition` is boolean ``np.extract`` is equivalent to ``arr[condition]``.
1865 Note that `place` does the exact opposite of `extract`.
1867 Parameters
1868 ----------
1869 condition : array_like
1870 An array whose nonzero or True entries indicate the elements of `arr`
1871 to extract.
1872 arr : array_like
1873 Input array of the same size as `condition`.
1875 Returns
1876 -------
1877 extract : ndarray
1878 Rank 1 array of values from `arr` where `condition` is True.
1880 See Also
1881 --------
1882 take, put, copyto, compress, place
1884 Examples
1885 --------
1886 >>> arr = np.arange(12).reshape((3, 4))
1887 >>> arr
1888 array([[ 0, 1, 2, 3],
1889 [ 4, 5, 6, 7],
1890 [ 8, 9, 10, 11]])
1891 >>> condition = np.mod(arr, 3)==0
1892 >>> condition
1893 array([[ True, False, False, True],
1894 [False, False, True, False],
1895 [False, True, False, False]])
1896 >>> np.extract(condition, arr)
1897 array([0, 3, 6, 9])
1900 If `condition` is boolean:
1902 >>> arr[condition]
1903 array([0, 3, 6, 9])
1905 """
1906 return _nx.take(ravel(arr), nonzero(ravel(condition))[0])
1909def _place_dispatcher(arr, mask, vals):
1910 return (arr, mask, vals)
1913@array_function_dispatch(_place_dispatcher)
1914def place(arr, mask, vals):
1915 """
1916 Change elements of an array based on conditional and input values.
1918 Similar to ``np.copyto(arr, vals, where=mask)``, the difference is that
1919 `place` uses the first N elements of `vals`, where N is the number of
1920 True values in `mask`, while `copyto` uses the elements where `mask`
1921 is True.
1923 Note that `extract` does the exact opposite of `place`.
1925 Parameters
1926 ----------
1927 arr : ndarray
1928 Array to put data into.
1929 mask : array_like
1930 Boolean mask array. Must have the same size as `a`.
1931 vals : 1-D sequence
1932 Values to put into `a`. Only the first N elements are used, where
1933 N is the number of True values in `mask`. If `vals` is smaller
1934 than N, it will be repeated, and if elements of `a` are to be masked,
1935 this sequence must be non-empty.
1937 See Also
1938 --------
1939 copyto, put, take, extract
1941 Examples
1942 --------
1943 >>> arr = np.arange(6).reshape(2, 3)
1944 >>> np.place(arr, arr>2, [44, 55])
1945 >>> arr
1946 array([[ 0, 1, 2],
1947 [44, 55, 44]])
1949 """
1950 if not isinstance(arr, np.ndarray):
1951 raise TypeError("argument 1 must be numpy.ndarray, "
1952 "not {name}".format(name=type(arr).__name__))
1954 return _insert(arr, mask, vals)
1957def disp(mesg, device=None, linefeed=True):
1958 """
1959 Display a message on a device.
1961 Parameters
1962 ----------
1963 mesg : str
1964 Message to display.
1965 device : object
1966 Device to write message. If None, defaults to ``sys.stdout`` which is
1967 very similar to ``print``. `device` needs to have ``write()`` and
1968 ``flush()`` methods.
1969 linefeed : bool, optional
1970 Option whether to print a line feed or not. Defaults to True.
1972 Raises
1973 ------
1974 AttributeError
1975 If `device` does not have a ``write()`` or ``flush()`` method.
1977 Examples
1978 --------
1979 Besides ``sys.stdout``, a file-like object can also be used as it has
1980 both required methods:
1982 >>> from io import StringIO
1983 >>> buf = StringIO()
1984 >>> np.disp(u'"Display" in a file', device=buf)
1985 >>> buf.getvalue()
1986 '"Display" in a file\\n'
1988 """
1989 if device is None:
1990 device = sys.stdout
1991 if linefeed:
1992 device.write('%s\n' % mesg)
1993 else:
1994 device.write('%s' % mesg)
1995 device.flush()
1996 return
1999# See https://docs.scipy.org/doc/numpy/reference/c-api.generalized-ufuncs.html
2000_DIMENSION_NAME = r'\w+'
2001_CORE_DIMENSION_LIST = '(?:{0:}(?:,{0:})*)?'.format(_DIMENSION_NAME)
2002_ARGUMENT = r'\({}\)'.format(_CORE_DIMENSION_LIST)
2003_ARGUMENT_LIST = '{0:}(?:,{0:})*'.format(_ARGUMENT)
2004_SIGNATURE = '^{0:}->{0:}$'.format(_ARGUMENT_LIST)
2007def _parse_gufunc_signature(signature):
2008 """
2009 Parse string signatures for a generalized universal function.
2011 Arguments
2012 ---------
2013 signature : string
2014 Generalized universal function signature, e.g., ``(m,n),(n,p)->(m,p)``
2015 for ``np.matmul``.
2017 Returns
2018 -------
2019 Tuple of input and output core dimensions parsed from the signature, each
2020 of the form List[Tuple[str, ...]].
2021 """
2022 signature = re.sub(r'\s+', '', signature)
2024 if not re.match(_SIGNATURE, signature):
2025 raise ValueError(
2026 'not a valid gufunc signature: {}'.format(signature))
2027 return tuple([tuple(re.findall(_DIMENSION_NAME, arg))
2028 for arg in re.findall(_ARGUMENT, arg_list)]
2029 for arg_list in signature.split('->'))
2032def _update_dim_sizes(dim_sizes, arg, core_dims):
2033 """
2034 Incrementally check and update core dimension sizes for a single argument.
2036 Arguments
2037 ---------
2038 dim_sizes : Dict[str, int]
2039 Sizes of existing core dimensions. Will be updated in-place.
2040 arg : ndarray
2041 Argument to examine.
2042 core_dims : Tuple[str, ...]
2043 Core dimensions for this argument.
2044 """
2045 if not core_dims:
2046 return
2048 num_core_dims = len(core_dims)
2049 if arg.ndim < num_core_dims:
2050 raise ValueError(
2051 '%d-dimensional argument does not have enough '
2052 'dimensions for all core dimensions %r'
2053 % (arg.ndim, core_dims))
2055 core_shape = arg.shape[-num_core_dims:]
2056 for dim, size in zip(core_dims, core_shape):
2057 if dim in dim_sizes:
2058 if size != dim_sizes[dim]:
2059 raise ValueError(
2060 'inconsistent size for core dimension %r: %r vs %r'
2061 % (dim, size, dim_sizes[dim]))
2062 else:
2063 dim_sizes[dim] = size
2066def _parse_input_dimensions(args, input_core_dims):
2067 """
2068 Parse broadcast and core dimensions for vectorize with a signature.
2070 Arguments
2071 ---------
2072 args : Tuple[ndarray, ...]
2073 Tuple of input arguments to examine.
2074 input_core_dims : List[Tuple[str, ...]]
2075 List of core dimensions corresponding to each input.
2077 Returns
2078 -------
2079 broadcast_shape : Tuple[int, ...]
2080 Common shape to broadcast all non-core dimensions to.
2081 dim_sizes : Dict[str, int]
2082 Common sizes for named core dimensions.
2083 """
2084 broadcast_args = []
2085 dim_sizes = {}
2086 for arg, core_dims in zip(args, input_core_dims):
2087 _update_dim_sizes(dim_sizes, arg, core_dims)
2088 ndim = arg.ndim - len(core_dims)
2089 dummy_array = np.lib.stride_tricks.as_strided(0, arg.shape[:ndim])
2090 broadcast_args.append(dummy_array)
2091 broadcast_shape = np.lib.stride_tricks._broadcast_shape(*broadcast_args)
2092 return broadcast_shape, dim_sizes
2095def _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims):
2096 """Helper for calculating broadcast shapes with core dimensions."""
2097 return [broadcast_shape + tuple(dim_sizes[dim] for dim in core_dims)
2098 for core_dims in list_of_core_dims]
2101def _create_arrays(broadcast_shape, dim_sizes, list_of_core_dims, dtypes,
2102 results=None):
2103 """Helper for creating output arrays in vectorize."""
2104 shapes = _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims)
2105 if dtypes is None:
2106 dtypes = [None] * len(shapes)
2107 if results is None:
2108 arrays = tuple(np.empty(shape=shape, dtype=dtype)
2109 for shape, dtype in zip(shapes, dtypes))
2110 else:
2111 arrays = tuple(np.empty_like(result, shape=shape, dtype=dtype)
2112 for result, shape, dtype
2113 in zip(results, shapes, dtypes))
2114 return arrays
2117@set_module('numpy')
2118class vectorize:
2119 """
2120 vectorize(pyfunc, otypes=None, doc=None, excluded=None, cache=False,
2121 signature=None)
2123 Generalized function class.
2125 Define a vectorized function which takes a nested sequence of objects or
2126 numpy arrays as inputs and returns a single numpy array or a tuple of numpy
2127 arrays. The vectorized function evaluates `pyfunc` over successive tuples
2128 of the input arrays like the python map function, except it uses the
2129 broadcasting rules of numpy.
2131 The data type of the output of `vectorized` is determined by calling
2132 the function with the first element of the input. This can be avoided
2133 by specifying the `otypes` argument.
2135 Parameters
2136 ----------
2137 pyfunc : callable
2138 A python function or method.
2139 otypes : str or list of dtypes, optional
2140 The output data type. It must be specified as either a string of
2141 typecode characters or a list of data type specifiers. There should
2142 be one data type specifier for each output.
2143 doc : str, optional
2144 The docstring for the function. If None, the docstring will be the
2145 ``pyfunc.__doc__``.
2146 excluded : set, optional
2147 Set of strings or integers representing the positional or keyword
2148 arguments for which the function will not be vectorized. These will be
2149 passed directly to `pyfunc` unmodified.
2151 .. versionadded:: 1.7.0
2153 cache : bool, optional
2154 If `True`, then cache the first function call that determines the number
2155 of outputs if `otypes` is not provided.
2157 .. versionadded:: 1.7.0
2159 signature : string, optional
2160 Generalized universal function signature, e.g., ``(m,n),(n)->(m)`` for
2161 vectorized matrix-vector multiplication. If provided, ``pyfunc`` will
2162 be called with (and expected to return) arrays with shapes given by the
2163 size of corresponding core dimensions. By default, ``pyfunc`` is
2164 assumed to take scalars as input and output.
2166 .. versionadded:: 1.12.0
2168 Returns
2169 -------
2170 vectorized : callable
2171 Vectorized function.
2173 See Also
2174 --------
2175 frompyfunc : Takes an arbitrary Python function and returns a ufunc
2177 Notes
2178 -----
2179 The `vectorize` function is provided primarily for convenience, not for
2180 performance. The implementation is essentially a for loop.
2182 If `otypes` is not specified, then a call to the function with the
2183 first argument will be used to determine the number of outputs. The
2184 results of this call will be cached if `cache` is `True` to prevent
2185 calling the function twice. However, to implement the cache, the
2186 original function must be wrapped which will slow down subsequent
2187 calls, so only do this if your function is expensive.
2189 The new keyword argument interface and `excluded` argument support
2190 further degrades performance.
2192 References
2193 ----------
2194 .. [1] :doc:`/reference/c-api/generalized-ufuncs`
2196 Examples
2197 --------
2198 >>> def myfunc(a, b):
2199 ... "Return a-b if a>b, otherwise return a+b"
2200 ... if a > b:
2201 ... return a - b
2202 ... else:
2203 ... return a + b
2205 >>> vfunc = np.vectorize(myfunc)
2206 >>> vfunc([1, 2, 3, 4], 2)
2207 array([3, 4, 1, 2])
2209 The docstring is taken from the input function to `vectorize` unless it
2210 is specified:
2212 >>> vfunc.__doc__
2213 'Return a-b if a>b, otherwise return a+b'
2214 >>> vfunc = np.vectorize(myfunc, doc='Vectorized `myfunc`')
2215 >>> vfunc.__doc__
2216 'Vectorized `myfunc`'
2218 The output type is determined by evaluating the first element of the input,
2219 unless it is specified:
2221 >>> out = vfunc([1, 2, 3, 4], 2)
2222 >>> type(out[0])
2223 <class 'numpy.int64'>
2224 >>> vfunc = np.vectorize(myfunc, otypes=[float])
2225 >>> out = vfunc([1, 2, 3, 4], 2)
2226 >>> type(out[0])
2227 <class 'numpy.float64'>
2229 The `excluded` argument can be used to prevent vectorizing over certain
2230 arguments. This can be useful for array-like arguments of a fixed length
2231 such as the coefficients for a polynomial as in `polyval`:
2233 >>> def mypolyval(p, x):
2234 ... _p = list(p)
2235 ... res = _p.pop(0)
2236 ... while _p:
2237 ... res = res*x + _p.pop(0)
2238 ... return res
2239 >>> vpolyval = np.vectorize(mypolyval, excluded=['p'])
2240 >>> vpolyval(p=[1, 2, 3], x=[0, 1])
2241 array([3, 6])
2243 Positional arguments may also be excluded by specifying their position:
2245 >>> vpolyval.excluded.add(0)
2246 >>> vpolyval([1, 2, 3], x=[0, 1])
2247 array([3, 6])
2249 The `signature` argument allows for vectorizing functions that act on
2250 non-scalar arrays of fixed length. For example, you can use it for a
2251 vectorized calculation of Pearson correlation coefficient and its p-value:
2253 >>> import scipy.stats
2254 >>> pearsonr = np.vectorize(scipy.stats.pearsonr,
2255 ... signature='(n),(n)->(),()')
2256 >>> pearsonr([[0, 1, 2, 3]], [[1, 2, 3, 4], [4, 3, 2, 1]])
2257 (array([ 1., -1.]), array([ 0., 0.]))
2259 Or for a vectorized convolution:
2261 >>> convolve = np.vectorize(np.convolve, signature='(n),(m)->(k)')
2262 >>> convolve(np.eye(4), [1, 2, 1])
2263 array([[1., 2., 1., 0., 0., 0.],
2264 [0., 1., 2., 1., 0., 0.],
2265 [0., 0., 1., 2., 1., 0.],
2266 [0., 0., 0., 1., 2., 1.]])
2268 """
2269 def __init__(self, pyfunc, otypes=None, doc=None, excluded=None,
2270 cache=False, signature=None):
2271 self.pyfunc = pyfunc
2272 self.cache = cache
2273 self.signature = signature
2274 self._ufunc = {} # Caching to improve default performance
2276 if doc is None:
2277 self.__doc__ = pyfunc.__doc__
2278 else:
2279 self.__doc__ = doc
2281 if isinstance(otypes, str):
2282 for char in otypes:
2283 if char not in typecodes['All']:
2284 raise ValueError("Invalid otype specified: %s" % (char,))
2285 elif iterable(otypes):
2286 otypes = ''.join([_nx.dtype(x).char for x in otypes])
2287 elif otypes is not None:
2288 raise ValueError("Invalid otype specification")
2289 self.otypes = otypes
2291 # Excluded variable support
2292 if excluded is None:
2293 excluded = set()
2294 self.excluded = set(excluded)
2296 if signature is not None:
2297 self._in_and_out_core_dims = _parse_gufunc_signature(signature)
2298 else:
2299 self._in_and_out_core_dims = None
2301 def __call__(self, *args, **kwargs):
2302 """
2303 Return arrays with the results of `pyfunc` broadcast (vectorized) over
2304 `args` and `kwargs` not in `excluded`.
2305 """
2306 excluded = self.excluded
2307 if not kwargs and not excluded:
2308 func = self.pyfunc
2309 vargs = args
2310 else:
2311 # The wrapper accepts only positional arguments: we use `names` and
2312 # `inds` to mutate `the_args` and `kwargs` to pass to the original
2313 # function.
2314 nargs = len(args)
2316 names = [_n for _n in kwargs if _n not in excluded]
2317 inds = [_i for _i in range(nargs) if _i not in excluded]
2318 the_args = list(args)
2320 def func(*vargs):
2321 for _n, _i in enumerate(inds):
2322 the_args[_i] = vargs[_n]
2323 kwargs.update(zip(names, vargs[len(inds):]))
2324 return self.pyfunc(*the_args, **kwargs)
2326 vargs = [args[_i] for _i in inds]
2327 vargs.extend([kwargs[_n] for _n in names])
2329 return self._vectorize_call(func=func, args=vargs)
2331 def _get_ufunc_and_otypes(self, func, args):
2332 """Return (ufunc, otypes)."""
2333 # frompyfunc will fail if args is empty
2334 if not args:
2335 raise ValueError('args can not be empty')
2337 if self.otypes is not None:
2338 otypes = self.otypes
2340 # self._ufunc is a dictionary whose keys are the number of
2341 # arguments (i.e. len(args)) and whose values are ufuncs created
2342 # by frompyfunc. len(args) can be different for different calls if
2343 # self.pyfunc has parameters with default values. We only use the
2344 # cache when func is self.pyfunc, which occurs when the call uses
2345 # only positional arguments and no arguments are excluded.
2347 nin = len(args)
2348 nout = len(self.otypes)
2349 if func is not self.pyfunc or nin not in self._ufunc:
2350 ufunc = frompyfunc(func, nin, nout)
2351 else:
2352 ufunc = None # We'll get it from self._ufunc
2353 if func is self.pyfunc:
2354 ufunc = self._ufunc.setdefault(nin, ufunc)
2355 else:
2356 # Get number of outputs and output types by calling the function on
2357 # the first entries of args. We also cache the result to prevent
2358 # the subsequent call when the ufunc is evaluated.
2359 # Assumes that ufunc first evaluates the 0th elements in the input
2360 # arrays (the input values are not checked to ensure this)
2361 args = [asarray(arg) for arg in args]
2362 if builtins.any(arg.size == 0 for arg in args):
2363 raise ValueError('cannot call `vectorize` on size 0 inputs '
2364 'unless `otypes` is set')
2366 inputs = [arg.flat[0] for arg in args]
2367 outputs = func(*inputs)
2369 # Performance note: profiling indicates that -- for simple
2370 # functions at least -- this wrapping can almost double the
2371 # execution time.
2372 # Hence we make it optional.
2373 if self.cache:
2374 _cache = [outputs]
2376 def _func(*vargs):
2377 if _cache:
2378 return _cache.pop()
2379 else:
2380 return func(*vargs)
2381 else:
2382 _func = func
2384 if isinstance(outputs, tuple):
2385 nout = len(outputs)
2386 else:
2387 nout = 1
2388 outputs = (outputs,)
2390 otypes = ''.join([asarray(outputs[_k]).dtype.char
2391 for _k in range(nout)])
2393 # Performance note: profiling indicates that creating the ufunc is
2394 # not a significant cost compared with wrapping so it seems not
2395 # worth trying to cache this.
2396 ufunc = frompyfunc(_func, len(args), nout)
2398 return ufunc, otypes
2400 def _vectorize_call(self, func, args):
2401 """Vectorized call to `func` over positional `args`."""
2402 if self.signature is not None:
2403 res = self._vectorize_call_with_signature(func, args)
2404 elif not args:
2405 res = func()
2406 else:
2407 ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args)
2409 # Convert args to object arrays first
2410 inputs = [asanyarray(a, dtype=object) for a in args]
2412 outputs = ufunc(*inputs)
2414 if ufunc.nout == 1:
2415 res = asanyarray(outputs, dtype=otypes[0])
2416 else:
2417 res = tuple([asanyarray(x, dtype=t)
2418 for x, t in zip(outputs, otypes)])
2419 return res
2421 def _vectorize_call_with_signature(self, func, args):
2422 """Vectorized call over positional arguments with a signature."""
2423 input_core_dims, output_core_dims = self._in_and_out_core_dims
2425 if len(args) != len(input_core_dims):
2426 raise TypeError('wrong number of positional arguments: '
2427 'expected %r, got %r'
2428 % (len(input_core_dims), len(args)))
2429 args = tuple(asanyarray(arg) for arg in args)
2431 broadcast_shape, dim_sizes = _parse_input_dimensions(
2432 args, input_core_dims)
2433 input_shapes = _calculate_shapes(broadcast_shape, dim_sizes,
2434 input_core_dims)
2435 args = [np.broadcast_to(arg, shape, subok=True)
2436 for arg, shape in zip(args, input_shapes)]
2438 outputs = None
2439 otypes = self.otypes
2440 nout = len(output_core_dims)
2442 for index in np.ndindex(*broadcast_shape):
2443 results = func(*(arg[index] for arg in args))
2445 n_results = len(results) if isinstance(results, tuple) else 1
2447 if nout != n_results:
2448 raise ValueError(
2449 'wrong number of outputs from pyfunc: expected %r, got %r'
2450 % (nout, n_results))
2452 if nout == 1:
2453 results = (results,)
2455 if outputs is None:
2456 for result, core_dims in zip(results, output_core_dims):
2457 _update_dim_sizes(dim_sizes, result, core_dims)
2459 outputs = _create_arrays(broadcast_shape, dim_sizes,
2460 output_core_dims, otypes, results)
2462 for output, result in zip(outputs, results):
2463 output[index] = result
2465 if outputs is None:
2466 # did not call the function even once
2467 if otypes is None:
2468 raise ValueError('cannot call `vectorize` on size 0 inputs '
2469 'unless `otypes` is set')
2470 if builtins.any(dim not in dim_sizes
2471 for dims in output_core_dims
2472 for dim in dims):
2473 raise ValueError('cannot call `vectorize` with a signature '
2474 'including new output dimensions on size 0 '
2475 'inputs')
2476 outputs = _create_arrays(broadcast_shape, dim_sizes,
2477 output_core_dims, otypes)
2479 return outputs[0] if nout == 1 else outputs
2482def _cov_dispatcher(m, y=None, rowvar=None, bias=None, ddof=None,
2483 fweights=None, aweights=None, *, dtype=None):
2484 return (m, y, fweights, aweights)
2487@array_function_dispatch(_cov_dispatcher)
2488def cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None,
2489 aweights=None, *, dtype=None):
2490 """
2491 Estimate a covariance matrix, given data and weights.
2493 Covariance indicates the level to which two variables vary together.
2494 If we examine N-dimensional samples, :math:`X = [x_1, x_2, ... x_N]^T`,
2495 then the covariance matrix element :math:`C_{ij}` is the covariance of
2496 :math:`x_i` and :math:`x_j`. The element :math:`C_{ii}` is the variance
2497 of :math:`x_i`.
2499 See the notes for an outline of the algorithm.
2501 Parameters
2502 ----------
2503 m : array_like
2504 A 1-D or 2-D array containing multiple variables and observations.
2505 Each row of `m` represents a variable, and each column a single
2506 observation of all those variables. Also see `rowvar` below.
2507 y : array_like, optional
2508 An additional set of variables and observations. `y` has the same form
2509 as that of `m`.
2510 rowvar : bool, optional
2511 If `rowvar` is True (default), then each row represents a
2512 variable, with observations in the columns. Otherwise, the relationship
2513 is transposed: each column represents a variable, while the rows
2514 contain observations.
2515 bias : bool, optional
2516 Default normalization (False) is by ``(N - 1)``, where ``N`` is the
2517 number of observations given (unbiased estimate). If `bias` is True,
2518 then normalization is by ``N``. These values can be overridden by using
2519 the keyword ``ddof`` in numpy versions >= 1.5.
2520 ddof : int, optional
2521 If not ``None`` the default value implied by `bias` is overridden.
2522 Note that ``ddof=1`` will return the unbiased estimate, even if both
2523 `fweights` and `aweights` are specified, and ``ddof=0`` will return
2524 the simple average. See the notes for the details. The default value
2525 is ``None``.
2527 .. versionadded:: 1.5
2528 fweights : array_like, int, optional
2529 1-D array of integer frequency weights; the number of times each
2530 observation vector should be repeated.
2532 .. versionadded:: 1.10
2533 aweights : array_like, optional
2534 1-D array of observation vector weights. These relative weights are
2535 typically large for observations considered "important" and smaller for
2536 observations considered less "important". If ``ddof=0`` the array of
2537 weights can be used to assign probabilities to observation vectors.
2539 .. versionadded:: 1.10
2540 dtype : data-type, optional
2541 Data-type of the result. By default, the return data-type will have
2542 at least `numpy.float64` precision.
2544 .. versionadded:: 1.20
2546 Returns
2547 -------
2548 out : ndarray
2549 The covariance matrix of the variables.
2551 See Also
2552 --------
2553 corrcoef : Normalized covariance matrix
2555 Notes
2556 -----
2557 Assume that the observations are in the columns of the observation
2558 array `m` and let ``f = fweights`` and ``a = aweights`` for brevity. The
2559 steps to compute the weighted covariance are as follows::
2561 >>> m = np.arange(10, dtype=np.float64)
2562 >>> f = np.arange(10) * 2
2563 >>> a = np.arange(10) ** 2.
2564 >>> ddof = 1
2565 >>> w = f * a
2566 >>> v1 = np.sum(w)
2567 >>> v2 = np.sum(w * a)
2568 >>> m -= np.sum(m * w, axis=None, keepdims=True) / v1
2569 >>> cov = np.dot(m * w, m.T) * v1 / (v1**2 - ddof * v2)
2571 Note that when ``a == 1``, the normalization factor
2572 ``v1 / (v1**2 - ddof * v2)`` goes over to ``1 / (np.sum(f) - ddof)``
2573 as it should.
2575 Examples
2576 --------
2577 Consider two variables, :math:`x_0` and :math:`x_1`, which
2578 correlate perfectly, but in opposite directions:
2580 >>> x = np.array([[0, 2], [1, 1], [2, 0]]).T
2581 >>> x
2582 array([[0, 1, 2],
2583 [2, 1, 0]])
2585 Note how :math:`x_0` increases while :math:`x_1` decreases. The covariance
2586 matrix shows this clearly:
2588 >>> np.cov(x)
2589 array([[ 1., -1.],
2590 [-1., 1.]])
2592 Note that element :math:`C_{0,1}`, which shows the correlation between
2593 :math:`x_0` and :math:`x_1`, is negative.
2595 Further, note how `x` and `y` are combined:
2597 >>> x = [-2.1, -1, 4.3]
2598 >>> y = [3, 1.1, 0.12]
2599 >>> X = np.stack((x, y), axis=0)
2600 >>> np.cov(X)
2601 array([[11.71 , -4.286 ], # may vary
2602 [-4.286 , 2.144133]])
2603 >>> np.cov(x, y)
2604 array([[11.71 , -4.286 ], # may vary
2605 [-4.286 , 2.144133]])
2606 >>> np.cov(x)
2607 array(11.71)
2609 """
2610 # Check inputs
2611 if ddof is not None and ddof != int(ddof):
2612 raise ValueError(
2613 "ddof must be integer")
2615 # Handles complex arrays too
2616 m = np.asarray(m)
2617 if m.ndim > 2:
2618 raise ValueError("m has more than 2 dimensions")
2620 if y is not None:
2621 y = np.asarray(y)
2622 if y.ndim > 2:
2623 raise ValueError("y has more than 2 dimensions")
2625 if dtype is None:
2626 if y is None:
2627 dtype = np.result_type(m, np.float64)
2628 else:
2629 dtype = np.result_type(m, y, np.float64)
2631 X = array(m, ndmin=2, dtype=dtype)
2632 if not rowvar and X.shape[0] != 1:
2633 X = X.T
2634 if X.shape[0] == 0:
2635 return np.array([]).reshape(0, 0)
2636 if y is not None:
2637 y = array(y, copy=False, ndmin=2, dtype=dtype)
2638 if not rowvar and y.shape[0] != 1:
2639 y = y.T
2640 X = np.concatenate((X, y), axis=0)
2642 if ddof is None:
2643 if bias == 0:
2644 ddof = 1
2645 else:
2646 ddof = 0
2648 # Get the product of frequencies and weights
2649 w = None
2650 if fweights is not None:
2651 fweights = np.asarray(fweights, dtype=float)
2652 if not np.all(fweights == np.around(fweights)):
2653 raise TypeError(
2654 "fweights must be integer")
2655 if fweights.ndim > 1:
2656 raise RuntimeError(
2657 "cannot handle multidimensional fweights")
2658 if fweights.shape[0] != X.shape[1]:
2659 raise RuntimeError(
2660 "incompatible numbers of samples and fweights")
2661 if any(fweights < 0):
2662 raise ValueError(
2663 "fweights cannot be negative")
2664 w = fweights
2665 if aweights is not None:
2666 aweights = np.asarray(aweights, dtype=float)
2667 if aweights.ndim > 1:
2668 raise RuntimeError(
2669 "cannot handle multidimensional aweights")
2670 if aweights.shape[0] != X.shape[1]:
2671 raise RuntimeError(
2672 "incompatible numbers of samples and aweights")
2673 if any(aweights < 0):
2674 raise ValueError(
2675 "aweights cannot be negative")
2676 if w is None:
2677 w = aweights
2678 else:
2679 w *= aweights
2681 avg, w_sum = average(X, axis=1, weights=w, returned=True)
2682 w_sum = w_sum[0]
2684 # Determine the normalization
2685 if w is None:
2686 fact = X.shape[1] - ddof
2687 elif ddof == 0:
2688 fact = w_sum
2689 elif aweights is None:
2690 fact = w_sum - ddof
2691 else:
2692 fact = w_sum - ddof*sum(w*aweights)/w_sum
2694 if fact <= 0:
2695 warnings.warn("Degrees of freedom <= 0 for slice",
2696 RuntimeWarning, stacklevel=3)
2697 fact = 0.0
2699 X -= avg[:, None]
2700 if w is None:
2701 X_T = X.T
2702 else:
2703 X_T = (X*w).T
2704 c = dot(X, X_T.conj())
2705 c *= np.true_divide(1, fact)
2706 return c.squeeze()
2709def _corrcoef_dispatcher(x, y=None, rowvar=None, bias=None, ddof=None, *,
2710 dtype=None):
2711 return (x, y)
2714@array_function_dispatch(_corrcoef_dispatcher)
2715def corrcoef(x, y=None, rowvar=True, bias=np._NoValue, ddof=np._NoValue, *,
2716 dtype=None):
2717 """
2718 Return Pearson product-moment correlation coefficients.
2720 Please refer to the documentation for `cov` for more detail. The
2721 relationship between the correlation coefficient matrix, `R`, and the
2722 covariance matrix, `C`, is
2724 .. math:: R_{ij} = \\frac{ C_{ij} } { \\sqrt{ C_{ii} C_{jj} } }
2726 The values of `R` are between -1 and 1, inclusive.
2728 Parameters
2729 ----------
2730 x : array_like
2731 A 1-D or 2-D array containing multiple variables and observations.
2732 Each row of `x` represents a variable, and each column a single
2733 observation of all those variables. Also see `rowvar` below.
2734 y : array_like, optional
2735 An additional set of variables and observations. `y` has the same
2736 shape as `x`.
2737 rowvar : bool, optional
2738 If `rowvar` is True (default), then each row represents a
2739 variable, with observations in the columns. Otherwise, the relationship
2740 is transposed: each column represents a variable, while the rows
2741 contain observations.
2742 bias : _NoValue, optional
2743 Has no effect, do not use.
2745 .. deprecated:: 1.10.0
2746 ddof : _NoValue, optional
2747 Has no effect, do not use.
2749 .. deprecated:: 1.10.0
2750 dtype : data-type, optional
2751 Data-type of the result. By default, the return data-type will have
2752 at least `numpy.float64` precision.
2754 .. versionadded:: 1.20
2756 Returns
2757 -------
2758 R : ndarray
2759 The correlation coefficient matrix of the variables.
2761 See Also
2762 --------
2763 cov : Covariance matrix
2765 Notes
2766 -----
2767 Due to floating point rounding the resulting array may not be Hermitian,
2768 the diagonal elements may not be 1, and the elements may not satisfy the
2769 inequality abs(a) <= 1. The real and imaginary parts are clipped to the
2770 interval [-1, 1] in an attempt to improve on that situation but is not
2771 much help in the complex case.
2773 This function accepts but discards arguments `bias` and `ddof`. This is
2774 for backwards compatibility with previous versions of this function. These
2775 arguments had no effect on the return values of the function and can be
2776 safely ignored in this and previous versions of numpy.
2778 Examples
2779 --------
2780 In this example we generate two random arrays, ``xarr`` and ``yarr``, and
2781 compute the row-wise and column-wise Pearson correlation coefficients,
2782 ``R``. Since ``rowvar`` is true by default, we first find the row-wise
2783 Pearson correlation coefficients between the variables of ``xarr``.
2785 >>> import numpy as np
2786 >>> rng = np.random.default_rng(seed=42)
2787 >>> xarr = rng.random((3, 3))
2788 >>> xarr
2789 array([[0.77395605, 0.43887844, 0.85859792],
2790 [0.69736803, 0.09417735, 0.97562235],
2791 [0.7611397 , 0.78606431, 0.12811363]])
2792 >>> R1 = np.corrcoef(xarr)
2793 >>> R1
2794 array([[ 1. , 0.99256089, -0.68080986],
2795 [ 0.99256089, 1. , -0.76492172],
2796 [-0.68080986, -0.76492172, 1. ]])
2798 If we add another set of variables and observations ``yarr``, we can
2799 compute the row-wise Pearson correlation coefficients between the
2800 variables in ``xarr`` and ``yarr``.
2802 >>> yarr = rng.random((3, 3))
2803 >>> yarr
2804 array([[0.45038594, 0.37079802, 0.92676499],
2805 [0.64386512, 0.82276161, 0.4434142 ],
2806 [0.22723872, 0.55458479, 0.06381726]])
2807 >>> R2 = np.corrcoef(xarr, yarr)
2808 >>> R2
2809 array([[ 1. , 0.99256089, -0.68080986, 0.75008178, -0.934284 ,
2810 -0.99004057],
2811 [ 0.99256089, 1. , -0.76492172, 0.82502011, -0.97074098,
2812 -0.99981569],
2813 [-0.68080986, -0.76492172, 1. , -0.99507202, 0.89721355,
2814 0.77714685],
2815 [ 0.75008178, 0.82502011, -0.99507202, 1. , -0.93657855,
2816 -0.83571711],
2817 [-0.934284 , -0.97074098, 0.89721355, -0.93657855, 1. ,
2818 0.97517215],
2819 [-0.99004057, -0.99981569, 0.77714685, -0.83571711, 0.97517215,
2820 1. ]])
2822 Finally if we use the option ``rowvar=False``, the columns are now
2823 being treated as the variables and we will find the column-wise Pearson
2824 correlation coefficients between variables in ``xarr`` and ``yarr``.
2826 >>> R3 = np.corrcoef(xarr, yarr, rowvar=False)
2827 >>> R3
2828 array([[ 1. , 0.77598074, -0.47458546, -0.75078643, -0.9665554 ,
2829 0.22423734],
2830 [ 0.77598074, 1. , -0.92346708, -0.99923895, -0.58826587,
2831 -0.44069024],
2832 [-0.47458546, -0.92346708, 1. , 0.93773029, 0.23297648,
2833 0.75137473],
2834 [-0.75078643, -0.99923895, 0.93773029, 1. , 0.55627469,
2835 0.47536961],
2836 [-0.9665554 , -0.58826587, 0.23297648, 0.55627469, 1. ,
2837 -0.46666491],
2838 [ 0.22423734, -0.44069024, 0.75137473, 0.47536961, -0.46666491,
2839 1. ]])
2841 """
2842 if bias is not np._NoValue or ddof is not np._NoValue:
2843 # 2015-03-15, 1.10
2844 warnings.warn('bias and ddof have no effect and are deprecated',
2845 DeprecationWarning, stacklevel=3)
2846 c = cov(x, y, rowvar, dtype=dtype)
2847 try:
2848 d = diag(c)
2849 except ValueError:
2850 # scalar covariance
2851 # nan if incorrect value (nan, inf, 0), 1 otherwise
2852 return c / c
2853 stddev = sqrt(d.real)
2854 c /= stddev[:, None]
2855 c /= stddev[None, :]
2857 # Clip real and imaginary parts to [-1, 1]. This does not guarantee
2858 # abs(a[i,j]) <= 1 for complex arrays, but is the best we can do without
2859 # excessive work.
2860 np.clip(c.real, -1, 1, out=c.real)
2861 if np.iscomplexobj(c):
2862 np.clip(c.imag, -1, 1, out=c.imag)
2864 return c
2867@set_module('numpy')
2868def blackman(M):
2869 """
2870 Return the Blackman window.
2872 The Blackman window is a taper formed by using the first three
2873 terms of a summation of cosines. It was designed to have close to the
2874 minimal leakage possible. It is close to optimal, only slightly worse
2875 than a Kaiser window.
2877 Parameters
2878 ----------
2879 M : int
2880 Number of points in the output window. If zero or less, an empty
2881 array is returned.
2883 Returns
2884 -------
2885 out : ndarray
2886 The window, with the maximum value normalized to one (the value one
2887 appears only if the number of samples is odd).
2889 See Also
2890 --------
2891 bartlett, hamming, hanning, kaiser
2893 Notes
2894 -----
2895 The Blackman window is defined as
2897 .. math:: w(n) = 0.42 - 0.5 \\cos(2\\pi n/M) + 0.08 \\cos(4\\pi n/M)
2899 Most references to the Blackman window come from the signal processing
2900 literature, where it is used as one of many windowing functions for
2901 smoothing values. It is also known as an apodization (which means
2902 "removing the foot", i.e. smoothing discontinuities at the beginning
2903 and end of the sampled signal) or tapering function. It is known as a
2904 "near optimal" tapering function, almost as good (by some measures)
2905 as the kaiser window.
2907 References
2908 ----------
2909 Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra,
2910 Dover Publications, New York.
2912 Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing.
2913 Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471.
2915 Examples
2916 --------
2917 >>> import matplotlib.pyplot as plt
2918 >>> np.blackman(12)
2919 array([-1.38777878e-17, 3.26064346e-02, 1.59903635e-01, # may vary
2920 4.14397981e-01, 7.36045180e-01, 9.67046769e-01,
2921 9.67046769e-01, 7.36045180e-01, 4.14397981e-01,
2922 1.59903635e-01, 3.26064346e-02, -1.38777878e-17])
2924 Plot the window and the frequency response:
2926 >>> from numpy.fft import fft, fftshift
2927 >>> window = np.blackman(51)
2928 >>> plt.plot(window)
2929 [<matplotlib.lines.Line2D object at 0x...>]
2930 >>> plt.title("Blackman window")
2931 Text(0.5, 1.0, 'Blackman window')
2932 >>> plt.ylabel("Amplitude")
2933 Text(0, 0.5, 'Amplitude')
2934 >>> plt.xlabel("Sample")
2935 Text(0.5, 0, 'Sample')
2936 >>> plt.show()
2938 >>> plt.figure()
2939 <Figure size 640x480 with 0 Axes>
2940 >>> A = fft(window, 2048) / 25.5
2941 >>> mag = np.abs(fftshift(A))
2942 >>> freq = np.linspace(-0.5, 0.5, len(A))
2943 >>> with np.errstate(divide='ignore', invalid='ignore'):
2944 ... response = 20 * np.log10(mag)
2945 ...
2946 >>> response = np.clip(response, -100, 100)
2947 >>> plt.plot(freq, response)
2948 [<matplotlib.lines.Line2D object at 0x...>]
2949 >>> plt.title("Frequency response of Blackman window")
2950 Text(0.5, 1.0, 'Frequency response of Blackman window')
2951 >>> plt.ylabel("Magnitude [dB]")
2952 Text(0, 0.5, 'Magnitude [dB]')
2953 >>> plt.xlabel("Normalized frequency [cycles per sample]")
2954 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
2955 >>> _ = plt.axis('tight')
2956 >>> plt.show()
2958 """
2959 if M < 1:
2960 return array([], dtype=np.result_type(M, 0.0))
2961 if M == 1:
2962 return ones(1, dtype=np.result_type(M, 0.0))
2963 n = arange(1-M, M, 2)
2964 return 0.42 + 0.5*cos(pi*n/(M-1)) + 0.08*cos(2.0*pi*n/(M-1))
2967@set_module('numpy')
2968def bartlett(M):
2969 """
2970 Return the Bartlett window.
2972 The Bartlett window is very similar to a triangular window, except
2973 that the end points are at zero. It is often used in signal
2974 processing for tapering a signal, without generating too much
2975 ripple in the frequency domain.
2977 Parameters
2978 ----------
2979 M : int
2980 Number of points in the output window. If zero or less, an
2981 empty array is returned.
2983 Returns
2984 -------
2985 out : array
2986 The triangular window, with the maximum value normalized to one
2987 (the value one appears only if the number of samples is odd), with
2988 the first and last samples equal to zero.
2990 See Also
2991 --------
2992 blackman, hamming, hanning, kaiser
2994 Notes
2995 -----
2996 The Bartlett window is defined as
2998 .. math:: w(n) = \\frac{2}{M-1} \\left(
2999 \\frac{M-1}{2} - \\left|n - \\frac{M-1}{2}\\right|
3000 \\right)
3002 Most references to the Bartlett window come from the signal processing
3003 literature, where it is used as one of many windowing functions for
3004 smoothing values. Note that convolution with this window produces linear
3005 interpolation. It is also known as an apodization (which means "removing
3006 the foot", i.e. smoothing discontinuities at the beginning and end of the
3007 sampled signal) or tapering function. The Fourier transform of the
3008 Bartlett window is the product of two sinc functions. Note the excellent
3009 discussion in Kanasewich [2]_.
3011 References
3012 ----------
3013 .. [1] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra",
3014 Biometrika 37, 1-16, 1950.
3015 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
3016 The University of Alberta Press, 1975, pp. 109-110.
3017 .. [3] A.V. Oppenheim and R.W. Schafer, "Discrete-Time Signal
3018 Processing", Prentice-Hall, 1999, pp. 468-471.
3019 .. [4] Wikipedia, "Window function",
3020 https://en.wikipedia.org/wiki/Window_function
3021 .. [5] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3022 "Numerical Recipes", Cambridge University Press, 1986, page 429.
3024 Examples
3025 --------
3026 >>> import matplotlib.pyplot as plt
3027 >>> np.bartlett(12)
3028 array([ 0. , 0.18181818, 0.36363636, 0.54545455, 0.72727273, # may vary
3029 0.90909091, 0.90909091, 0.72727273, 0.54545455, 0.36363636,
3030 0.18181818, 0. ])
3032 Plot the window and its frequency response (requires SciPy and matplotlib):
3034 >>> from numpy.fft import fft, fftshift
3035 >>> window = np.bartlett(51)
3036 >>> plt.plot(window)
3037 [<matplotlib.lines.Line2D object at 0x...>]
3038 >>> plt.title("Bartlett window")
3039 Text(0.5, 1.0, 'Bartlett window')
3040 >>> plt.ylabel("Amplitude")
3041 Text(0, 0.5, 'Amplitude')
3042 >>> plt.xlabel("Sample")
3043 Text(0.5, 0, 'Sample')
3044 >>> plt.show()
3046 >>> plt.figure()
3047 <Figure size 640x480 with 0 Axes>
3048 >>> A = fft(window, 2048) / 25.5
3049 >>> mag = np.abs(fftshift(A))
3050 >>> freq = np.linspace(-0.5, 0.5, len(A))
3051 >>> with np.errstate(divide='ignore', invalid='ignore'):
3052 ... response = 20 * np.log10(mag)
3053 ...
3054 >>> response = np.clip(response, -100, 100)
3055 >>> plt.plot(freq, response)
3056 [<matplotlib.lines.Line2D object at 0x...>]
3057 >>> plt.title("Frequency response of Bartlett window")
3058 Text(0.5, 1.0, 'Frequency response of Bartlett window')
3059 >>> plt.ylabel("Magnitude [dB]")
3060 Text(0, 0.5, 'Magnitude [dB]')
3061 >>> plt.xlabel("Normalized frequency [cycles per sample]")
3062 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
3063 >>> _ = plt.axis('tight')
3064 >>> plt.show()
3066 """
3067 if M < 1:
3068 return array([], dtype=np.result_type(M, 0.0))
3069 if M == 1:
3070 return ones(1, dtype=np.result_type(M, 0.0))
3071 n = arange(1-M, M, 2)
3072 return where(less_equal(n, 0), 1 + n/(M-1), 1 - n/(M-1))
3075@set_module('numpy')
3076def hanning(M):
3077 """
3078 Return the Hanning window.
3080 The Hanning window is a taper formed by using a weighted cosine.
3082 Parameters
3083 ----------
3084 M : int
3085 Number of points in the output window. If zero or less, an
3086 empty array is returned.
3088 Returns
3089 -------
3090 out : ndarray, shape(M,)
3091 The window, with the maximum value normalized to one (the value
3092 one appears only if `M` is odd).
3094 See Also
3095 --------
3096 bartlett, blackman, hamming, kaiser
3098 Notes
3099 -----
3100 The Hanning window is defined as
3102 .. math:: w(n) = 0.5 - 0.5\\cos\\left(\\frac{2\\pi{n}}{M-1}\\right)
3103 \\qquad 0 \\leq n \\leq M-1
3105 The Hanning was named for Julius von Hann, an Austrian meteorologist.
3106 It is also known as the Cosine Bell. Some authors prefer that it be
3107 called a Hann window, to help avoid confusion with the very similar
3108 Hamming window.
3110 Most references to the Hanning window come from the signal processing
3111 literature, where it is used as one of many windowing functions for
3112 smoothing values. It is also known as an apodization (which means
3113 "removing the foot", i.e. smoothing discontinuities at the beginning
3114 and end of the sampled signal) or tapering function.
3116 References
3117 ----------
3118 .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
3119 spectra, Dover Publications, New York.
3120 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
3121 The University of Alberta Press, 1975, pp. 106-108.
3122 .. [3] Wikipedia, "Window function",
3123 https://en.wikipedia.org/wiki/Window_function
3124 .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3125 "Numerical Recipes", Cambridge University Press, 1986, page 425.
3127 Examples
3128 --------
3129 >>> np.hanning(12)
3130 array([0. , 0.07937323, 0.29229249, 0.57115742, 0.82743037,
3131 0.97974649, 0.97974649, 0.82743037, 0.57115742, 0.29229249,
3132 0.07937323, 0. ])
3134 Plot the window and its frequency response:
3136 >>> import matplotlib.pyplot as plt
3137 >>> from numpy.fft import fft, fftshift
3138 >>> window = np.hanning(51)
3139 >>> plt.plot(window)
3140 [<matplotlib.lines.Line2D object at 0x...>]
3141 >>> plt.title("Hann window")
3142 Text(0.5, 1.0, 'Hann window')
3143 >>> plt.ylabel("Amplitude")
3144 Text(0, 0.5, 'Amplitude')
3145 >>> plt.xlabel("Sample")
3146 Text(0.5, 0, 'Sample')
3147 >>> plt.show()
3149 >>> plt.figure()
3150 <Figure size 640x480 with 0 Axes>
3151 >>> A = fft(window, 2048) / 25.5
3152 >>> mag = np.abs(fftshift(A))
3153 >>> freq = np.linspace(-0.5, 0.5, len(A))
3154 >>> with np.errstate(divide='ignore', invalid='ignore'):
3155 ... response = 20 * np.log10(mag)
3156 ...
3157 >>> response = np.clip(response, -100, 100)
3158 >>> plt.plot(freq, response)
3159 [<matplotlib.lines.Line2D object at 0x...>]
3160 >>> plt.title("Frequency response of the Hann window")
3161 Text(0.5, 1.0, 'Frequency response of the Hann window')
3162 >>> plt.ylabel("Magnitude [dB]")
3163 Text(0, 0.5, 'Magnitude [dB]')
3164 >>> plt.xlabel("Normalized frequency [cycles per sample]")
3165 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
3166 >>> plt.axis('tight')
3167 ...
3168 >>> plt.show()
3170 """
3171 if M < 1:
3172 return array([], dtype=np.result_type(M, 0.0))
3173 if M == 1:
3174 return ones(1, dtype=np.result_type(M, 0.0))
3175 n = arange(1-M, M, 2)
3176 return 0.5 + 0.5*cos(pi*n/(M-1))
3179@set_module('numpy')
3180def hamming(M):
3181 """
3182 Return the Hamming window.
3184 The Hamming window is a taper formed by using a weighted cosine.
3186 Parameters
3187 ----------
3188 M : int
3189 Number of points in the output window. If zero or less, an
3190 empty array is returned.
3192 Returns
3193 -------
3194 out : ndarray
3195 The window, with the maximum value normalized to one (the value
3196 one appears only if the number of samples is odd).
3198 See Also
3199 --------
3200 bartlett, blackman, hanning, kaiser
3202 Notes
3203 -----
3204 The Hamming window is defined as
3206 .. math:: w(n) = 0.54 - 0.46\\cos\\left(\\frac{2\\pi{n}}{M-1}\\right)
3207 \\qquad 0 \\leq n \\leq M-1
3209 The Hamming was named for R. W. Hamming, an associate of J. W. Tukey
3210 and is described in Blackman and Tukey. It was recommended for
3211 smoothing the truncated autocovariance function in the time domain.
3212 Most references to the Hamming window come from the signal processing
3213 literature, where it is used as one of many windowing functions for
3214 smoothing values. It is also known as an apodization (which means
3215 "removing the foot", i.e. smoothing discontinuities at the beginning
3216 and end of the sampled signal) or tapering function.
3218 References
3219 ----------
3220 .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
3221 spectra, Dover Publications, New York.
3222 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
3223 University of Alberta Press, 1975, pp. 109-110.
3224 .. [3] Wikipedia, "Window function",
3225 https://en.wikipedia.org/wiki/Window_function
3226 .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3227 "Numerical Recipes", Cambridge University Press, 1986, page 425.
3229 Examples
3230 --------
3231 >>> np.hamming(12)
3232 array([ 0.08 , 0.15302337, 0.34890909, 0.60546483, 0.84123594, # may vary
3233 0.98136677, 0.98136677, 0.84123594, 0.60546483, 0.34890909,
3234 0.15302337, 0.08 ])
3236 Plot the window and the frequency response:
3238 >>> import matplotlib.pyplot as plt
3239 >>> from numpy.fft import fft, fftshift
3240 >>> window = np.hamming(51)
3241 >>> plt.plot(window)
3242 [<matplotlib.lines.Line2D object at 0x...>]
3243 >>> plt.title("Hamming window")
3244 Text(0.5, 1.0, 'Hamming window')
3245 >>> plt.ylabel("Amplitude")
3246 Text(0, 0.5, 'Amplitude')
3247 >>> plt.xlabel("Sample")
3248 Text(0.5, 0, 'Sample')
3249 >>> plt.show()
3251 >>> plt.figure()
3252 <Figure size 640x480 with 0 Axes>
3253 >>> A = fft(window, 2048) / 25.5
3254 >>> mag = np.abs(fftshift(A))
3255 >>> freq = np.linspace(-0.5, 0.5, len(A))
3256 >>> response = 20 * np.log10(mag)
3257 >>> response = np.clip(response, -100, 100)
3258 >>> plt.plot(freq, response)
3259 [<matplotlib.lines.Line2D object at 0x...>]
3260 >>> plt.title("Frequency response of Hamming window")
3261 Text(0.5, 1.0, 'Frequency response of Hamming window')
3262 >>> plt.ylabel("Magnitude [dB]")
3263 Text(0, 0.5, 'Magnitude [dB]')
3264 >>> plt.xlabel("Normalized frequency [cycles per sample]")
3265 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
3266 >>> plt.axis('tight')
3267 ...
3268 >>> plt.show()
3270 """
3271 if M < 1:
3272 return array([], dtype=np.result_type(M, 0.0))
3273 if M == 1:
3274 return ones(1, dtype=np.result_type(M, 0.0))
3275 n = arange(1-M, M, 2)
3276 return 0.54 + 0.46*cos(pi*n/(M-1))
3279## Code from cephes for i0
3281_i0A = [
3282 -4.41534164647933937950E-18,
3283 3.33079451882223809783E-17,
3284 -2.43127984654795469359E-16,
3285 1.71539128555513303061E-15,
3286 -1.16853328779934516808E-14,
3287 7.67618549860493561688E-14,
3288 -4.85644678311192946090E-13,
3289 2.95505266312963983461E-12,
3290 -1.72682629144155570723E-11,
3291 9.67580903537323691224E-11,
3292 -5.18979560163526290666E-10,
3293 2.65982372468238665035E-9,
3294 -1.30002500998624804212E-8,
3295 6.04699502254191894932E-8,
3296 -2.67079385394061173391E-7,
3297 1.11738753912010371815E-6,
3298 -4.41673835845875056359E-6,
3299 1.64484480707288970893E-5,
3300 -5.75419501008210370398E-5,
3301 1.88502885095841655729E-4,
3302 -5.76375574538582365885E-4,
3303 1.63947561694133579842E-3,
3304 -4.32430999505057594430E-3,
3305 1.05464603945949983183E-2,
3306 -2.37374148058994688156E-2,
3307 4.93052842396707084878E-2,
3308 -9.49010970480476444210E-2,
3309 1.71620901522208775349E-1,
3310 -3.04682672343198398683E-1,
3311 6.76795274409476084995E-1
3312 ]
3314_i0B = [
3315 -7.23318048787475395456E-18,
3316 -4.83050448594418207126E-18,
3317 4.46562142029675999901E-17,
3318 3.46122286769746109310E-17,
3319 -2.82762398051658348494E-16,
3320 -3.42548561967721913462E-16,
3321 1.77256013305652638360E-15,
3322 3.81168066935262242075E-15,
3323 -9.55484669882830764870E-15,
3324 -4.15056934728722208663E-14,
3325 1.54008621752140982691E-14,
3326 3.85277838274214270114E-13,
3327 7.18012445138366623367E-13,
3328 -1.79417853150680611778E-12,
3329 -1.32158118404477131188E-11,
3330 -3.14991652796324136454E-11,
3331 1.18891471078464383424E-11,
3332 4.94060238822496958910E-10,
3333 3.39623202570838634515E-9,
3334 2.26666899049817806459E-8,
3335 2.04891858946906374183E-7,
3336 2.89137052083475648297E-6,
3337 6.88975834691682398426E-5,
3338 3.36911647825569408990E-3,
3339 8.04490411014108831608E-1
3340 ]
3343def _chbevl(x, vals):
3344 b0 = vals[0]
3345 b1 = 0.0
3347 for i in range(1, len(vals)):
3348 b2 = b1
3349 b1 = b0
3350 b0 = x*b1 - b2 + vals[i]
3352 return 0.5*(b0 - b2)
3355def _i0_1(x):
3356 return exp(x) * _chbevl(x/2.0-2, _i0A)
3359def _i0_2(x):
3360 return exp(x) * _chbevl(32.0/x - 2.0, _i0B) / sqrt(x)
3363def _i0_dispatcher(x):
3364 return (x,)
3367@array_function_dispatch(_i0_dispatcher)
3368def i0(x):
3369 """
3370 Modified Bessel function of the first kind, order 0.
3372 Usually denoted :math:`I_0`.
3374 Parameters
3375 ----------
3376 x : array_like of float
3377 Argument of the Bessel function.
3379 Returns
3380 -------
3381 out : ndarray, shape = x.shape, dtype = float
3382 The modified Bessel function evaluated at each of the elements of `x`.
3384 See Also
3385 --------
3386 scipy.special.i0, scipy.special.iv, scipy.special.ive
3388 Notes
3389 -----
3390 The scipy implementation is recommended over this function: it is a
3391 proper ufunc written in C, and more than an order of magnitude faster.
3393 We use the algorithm published by Clenshaw [1]_ and referenced by
3394 Abramowitz and Stegun [2]_, for which the function domain is
3395 partitioned into the two intervals [0,8] and (8,inf), and Chebyshev
3396 polynomial expansions are employed in each interval. Relative error on
3397 the domain [0,30] using IEEE arithmetic is documented [3]_ as having a
3398 peak of 5.8e-16 with an rms of 1.4e-16 (n = 30000).
3400 References
3401 ----------
3402 .. [1] C. W. Clenshaw, "Chebyshev series for mathematical functions", in
3403 *National Physical Laboratory Mathematical Tables*, vol. 5, London:
3404 Her Majesty's Stationery Office, 1962.
3405 .. [2] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical
3406 Functions*, 10th printing, New York: Dover, 1964, pp. 379.
3407 https://personal.math.ubc.ca/~cbm/aands/page_379.htm
3408 .. [3] https://metacpan.org/pod/distribution/Math-Cephes/lib/Math/Cephes.pod#i0:-Modified-Bessel-function-of-order-zero
3410 Examples
3411 --------
3412 >>> np.i0(0.)
3413 array(1.0)
3414 >>> np.i0([0, 1, 2, 3])
3415 array([1. , 1.26606588, 2.2795853 , 4.88079259])
3417 """
3418 x = np.asanyarray(x)
3419 if x.dtype.kind == 'c':
3420 raise TypeError("i0 not supported for complex values")
3421 if x.dtype.kind != 'f':
3422 x = x.astype(float)
3423 x = np.abs(x)
3424 return piecewise(x, [x <= 8.0], [_i0_1, _i0_2])
3426## End of cephes code for i0
3429@set_module('numpy')
3430def kaiser(M, beta):
3431 """
3432 Return the Kaiser window.
3434 The Kaiser window is a taper formed by using a Bessel function.
3436 Parameters
3437 ----------
3438 M : int
3439 Number of points in the output window. If zero or less, an
3440 empty array is returned.
3441 beta : float
3442 Shape parameter for window.
3444 Returns
3445 -------
3446 out : array
3447 The window, with the maximum value normalized to one (the value
3448 one appears only if the number of samples is odd).
3450 See Also
3451 --------
3452 bartlett, blackman, hamming, hanning
3454 Notes
3455 -----
3456 The Kaiser window is defined as
3458 .. math:: w(n) = I_0\\left( \\beta \\sqrt{1-\\frac{4n^2}{(M-1)^2}}
3459 \\right)/I_0(\\beta)
3461 with
3463 .. math:: \\quad -\\frac{M-1}{2} \\leq n \\leq \\frac{M-1}{2},
3465 where :math:`I_0` is the modified zeroth-order Bessel function.
3467 The Kaiser was named for Jim Kaiser, who discovered a simple
3468 approximation to the DPSS window based on Bessel functions. The Kaiser
3469 window is a very good approximation to the Digital Prolate Spheroidal
3470 Sequence, or Slepian window, which is the transform which maximizes the
3471 energy in the main lobe of the window relative to total energy.
3473 The Kaiser can approximate many other windows by varying the beta
3474 parameter.
3476 ==== =======================
3477 beta Window shape
3478 ==== =======================
3479 0 Rectangular
3480 5 Similar to a Hamming
3481 6 Similar to a Hanning
3482 8.6 Similar to a Blackman
3483 ==== =======================
3485 A beta value of 14 is probably a good starting point. Note that as beta
3486 gets large, the window narrows, and so the number of samples needs to be
3487 large enough to sample the increasingly narrow spike, otherwise NaNs will
3488 get returned.
3490 Most references to the Kaiser window come from the signal processing
3491 literature, where it is used as one of many windowing functions for
3492 smoothing values. It is also known as an apodization (which means
3493 "removing the foot", i.e. smoothing discontinuities at the beginning
3494 and end of the sampled signal) or tapering function.
3496 References
3497 ----------
3498 .. [1] J. F. Kaiser, "Digital Filters" - Ch 7 in "Systems analysis by
3499 digital computer", Editors: F.F. Kuo and J.F. Kaiser, p 218-285.
3500 John Wiley and Sons, New York, (1966).
3501 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
3502 University of Alberta Press, 1975, pp. 177-178.
3503 .. [3] Wikipedia, "Window function",
3504 https://en.wikipedia.org/wiki/Window_function
3506 Examples
3507 --------
3508 >>> import matplotlib.pyplot as plt
3509 >>> np.kaiser(12, 14)
3510 array([7.72686684e-06, 3.46009194e-03, 4.65200189e-02, # may vary
3511 2.29737120e-01, 5.99885316e-01, 9.45674898e-01,
3512 9.45674898e-01, 5.99885316e-01, 2.29737120e-01,
3513 4.65200189e-02, 3.46009194e-03, 7.72686684e-06])
3516 Plot the window and the frequency response:
3518 >>> from numpy.fft import fft, fftshift
3519 >>> window = np.kaiser(51, 14)
3520 >>> plt.plot(window)
3521 [<matplotlib.lines.Line2D object at 0x...>]
3522 >>> plt.title("Kaiser window")
3523 Text(0.5, 1.0, 'Kaiser window')
3524 >>> plt.ylabel("Amplitude")
3525 Text(0, 0.5, 'Amplitude')
3526 >>> plt.xlabel("Sample")
3527 Text(0.5, 0, 'Sample')
3528 >>> plt.show()
3530 >>> plt.figure()
3531 <Figure size 640x480 with 0 Axes>
3532 >>> A = fft(window, 2048) / 25.5
3533 >>> mag = np.abs(fftshift(A))
3534 >>> freq = np.linspace(-0.5, 0.5, len(A))
3535 >>> response = 20 * np.log10(mag)
3536 >>> response = np.clip(response, -100, 100)
3537 >>> plt.plot(freq, response)
3538 [<matplotlib.lines.Line2D object at 0x...>]
3539 >>> plt.title("Frequency response of Kaiser window")
3540 Text(0.5, 1.0, 'Frequency response of Kaiser window')
3541 >>> plt.ylabel("Magnitude [dB]")
3542 Text(0, 0.5, 'Magnitude [dB]')
3543 >>> plt.xlabel("Normalized frequency [cycles per sample]")
3544 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
3545 >>> plt.axis('tight')
3546 (-0.5, 0.5, -100.0, ...) # may vary
3547 >>> plt.show()
3549 """
3550 if M == 1:
3551 return np.ones(1, dtype=np.result_type(M, 0.0))
3552 n = arange(0, M)
3553 alpha = (M-1)/2.0
3554 return i0(beta * sqrt(1-((n-alpha)/alpha)**2.0))/i0(float(beta))
3557def _sinc_dispatcher(x):
3558 return (x,)
3561@array_function_dispatch(_sinc_dispatcher)
3562def sinc(x):
3563 r"""
3564 Return the normalized sinc function.
3566 The sinc function is equal to :math:`\sin(\pi x)/(\pi x)` for any argument
3567 :math:`x\ne 0`. ``sinc(0)`` takes the limit value 1, making ``sinc`` not
3568 only everywhere continuous but also infinitely differentiable.
3570 .. note::
3572 Note the normalization factor of ``pi`` used in the definition.
3573 This is the most commonly used definition in signal processing.
3574 Use ``sinc(x / np.pi)`` to obtain the unnormalized sinc function
3575 :math:`\sin(x)/x` that is more common in mathematics.
3577 Parameters
3578 ----------
3579 x : ndarray
3580 Array (possibly multi-dimensional) of values for which to calculate
3581 ``sinc(x)``.
3583 Returns
3584 -------
3585 out : ndarray
3586 ``sinc(x)``, which has the same shape as the input.
3588 Notes
3589 -----
3590 The name sinc is short for "sine cardinal" or "sinus cardinalis".
3592 The sinc function is used in various signal processing applications,
3593 including in anti-aliasing, in the construction of a Lanczos resampling
3594 filter, and in interpolation.
3596 For bandlimited interpolation of discrete-time signals, the ideal
3597 interpolation kernel is proportional to the sinc function.
3599 References
3600 ----------
3601 .. [1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram Web
3602 Resource. http://mathworld.wolfram.com/SincFunction.html
3603 .. [2] Wikipedia, "Sinc function",
3604 https://en.wikipedia.org/wiki/Sinc_function
3606 Examples
3607 --------
3608 >>> import matplotlib.pyplot as plt
3609 >>> x = np.linspace(-4, 4, 41)
3610 >>> np.sinc(x)
3611 array([-3.89804309e-17, -4.92362781e-02, -8.40918587e-02, # may vary
3612 -8.90384387e-02, -5.84680802e-02, 3.89804309e-17,
3613 6.68206631e-02, 1.16434881e-01, 1.26137788e-01,
3614 8.50444803e-02, -3.89804309e-17, -1.03943254e-01,
3615 -1.89206682e-01, -2.16236208e-01, -1.55914881e-01,
3616 3.89804309e-17, 2.33872321e-01, 5.04551152e-01,
3617 7.56826729e-01, 9.35489284e-01, 1.00000000e+00,
3618 9.35489284e-01, 7.56826729e-01, 5.04551152e-01,
3619 2.33872321e-01, 3.89804309e-17, -1.55914881e-01,
3620 -2.16236208e-01, -1.89206682e-01, -1.03943254e-01,
3621 -3.89804309e-17, 8.50444803e-02, 1.26137788e-01,
3622 1.16434881e-01, 6.68206631e-02, 3.89804309e-17,
3623 -5.84680802e-02, -8.90384387e-02, -8.40918587e-02,
3624 -4.92362781e-02, -3.89804309e-17])
3626 >>> plt.plot(x, np.sinc(x))
3627 [<matplotlib.lines.Line2D object at 0x...>]
3628 >>> plt.title("Sinc Function")
3629 Text(0.5, 1.0, 'Sinc Function')
3630 >>> plt.ylabel("Amplitude")
3631 Text(0, 0.5, 'Amplitude')
3632 >>> plt.xlabel("X")
3633 Text(0.5, 0, 'X')
3634 >>> plt.show()
3636 """
3637 x = np.asanyarray(x)
3638 y = pi * where(x == 0, 1.0e-20, x)
3639 return sin(y)/y
3642def _msort_dispatcher(a):
3643 return (a,)
3646@array_function_dispatch(_msort_dispatcher)
3647def msort(a):
3648 """
3649 Return a copy of an array sorted along the first axis.
3651 .. deprecated:: 1.24
3653 msort is deprecated, use ``np.sort(a, axis=0)`` instead.
3655 Parameters
3656 ----------
3657 a : array_like
3658 Array to be sorted.
3660 Returns
3661 -------
3662 sorted_array : ndarray
3663 Array of the same type and shape as `a`.
3665 See Also
3666 --------
3667 sort
3669 Notes
3670 -----
3671 ``np.msort(a)`` is equivalent to ``np.sort(a, axis=0)``.
3673 Examples
3674 --------
3675 >>> a = np.array([[1, 4], [3, 1]])
3676 >>> np.msort(a) # sort along the first axis
3677 array([[1, 1],
3678 [3, 4]])
3680 """
3681 # 2022-10-20 1.24
3682 warnings.warn(
3683 "msort is deprecated, use np.sort(a, axis=0) instead",
3684 DeprecationWarning,
3685 stacklevel=3,
3686 )
3687 b = array(a, subok=True, copy=True)
3688 b.sort(0)
3689 return b
3692def _ureduce(a, func, keepdims=False, **kwargs):
3693 """
3694 Internal Function.
3695 Call `func` with `a` as first argument swapping the axes to use extended
3696 axis on functions that don't support it natively.
3698 Returns result and a.shape with axis dims set to 1.
3700 Parameters
3701 ----------
3702 a : array_like
3703 Input array or object that can be converted to an array.
3704 func : callable
3705 Reduction function capable of receiving a single axis argument.
3706 It is called with `a` as first argument followed by `kwargs`.
3707 kwargs : keyword arguments
3708 additional keyword arguments to pass to `func`.
3710 Returns
3711 -------
3712 result : tuple
3713 Result of func(a, **kwargs) and a.shape with axis dims set to 1
3714 which can be used to reshape the result to the same shape a ufunc with
3715 keepdims=True would produce.
3717 """
3718 a = np.asanyarray(a)
3719 axis = kwargs.get('axis', None)
3720 out = kwargs.get('out', None)
3722 if keepdims is np._NoValue:
3723 keepdims = False
3725 nd = a.ndim
3726 if axis is not None:
3727 axis = _nx.normalize_axis_tuple(axis, nd)
3729 if keepdims:
3730 if out is not None:
3731 index_out = tuple(
3732 0 if i in axis else slice(None) for i in range(nd))
3733 kwargs['out'] = out[(Ellipsis, ) + index_out]
3735 if len(axis) == 1:
3736 kwargs['axis'] = axis[0]
3737 else:
3738 keep = set(range(nd)) - set(axis)
3739 nkeep = len(keep)
3740 # swap axis that should not be reduced to front
3741 for i, s in enumerate(sorted(keep)):
3742 a = a.swapaxes(i, s)
3743 # merge reduced axis
3744 a = a.reshape(a.shape[:nkeep] + (-1,))
3745 kwargs['axis'] = -1
3746 else:
3747 if keepdims:
3748 if out is not None:
3749 index_out = (0, ) * nd
3750 kwargs['out'] = out[(Ellipsis, ) + index_out]
3752 r = func(a, **kwargs)
3754 if out is not None:
3755 return out
3757 if keepdims:
3758 if axis is None:
3759 index_r = (np.newaxis, ) * nd
3760 else:
3761 index_r = tuple(
3762 np.newaxis if i in axis else slice(None)
3763 for i in range(nd))
3764 r = r[(Ellipsis, ) + index_r]
3766 return r
3769def _median_dispatcher(
3770 a, axis=None, out=None, overwrite_input=None, keepdims=None):
3771 return (a, out)
3774@array_function_dispatch(_median_dispatcher)
3775def median(a, axis=None, out=None, overwrite_input=False, keepdims=False):
3776 """
3777 Compute the median along the specified axis.
3779 Returns the median of the array elements.
3781 Parameters
3782 ----------
3783 a : array_like
3784 Input array or object that can be converted to an array.
3785 axis : {int, sequence of int, None}, optional
3786 Axis or axes along which the medians are computed. The default
3787 is to compute the median along a flattened version of the array.
3788 A sequence of axes is supported since version 1.9.0.
3789 out : ndarray, optional
3790 Alternative output array in which to place the result. It must
3791 have the same shape and buffer length as the expected output,
3792 but the type (of the output) will be cast if necessary.
3793 overwrite_input : bool, optional
3794 If True, then allow use of memory of input array `a` for
3795 calculations. The input array will be modified by the call to
3796 `median`. This will save memory when you do not need to preserve
3797 the contents of the input array. Treat the input as undefined,
3798 but it will probably be fully or partially sorted. Default is
3799 False. If `overwrite_input` is ``True`` and `a` is not already an
3800 `ndarray`, an error will be raised.
3801 keepdims : bool, optional
3802 If this is set to True, the axes which are reduced are left
3803 in the result as dimensions with size one. With this option,
3804 the result will broadcast correctly against the original `arr`.
3806 .. versionadded:: 1.9.0
3808 Returns
3809 -------
3810 median : ndarray
3811 A new array holding the result. If the input contains integers
3812 or floats smaller than ``float64``, then the output data-type is
3813 ``np.float64``. Otherwise, the data-type of the output is the
3814 same as that of the input. If `out` is specified, that array is
3815 returned instead.
3817 See Also
3818 --------
3819 mean, percentile
3821 Notes
3822 -----
3823 Given a vector ``V`` of length ``N``, the median of ``V`` is the
3824 middle value of a sorted copy of ``V``, ``V_sorted`` - i
3825 e., ``V_sorted[(N-1)/2]``, when ``N`` is odd, and the average of the
3826 two middle values of ``V_sorted`` when ``N`` is even.
3828 Examples
3829 --------
3830 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
3831 >>> a
3832 array([[10, 7, 4],
3833 [ 3, 2, 1]])
3834 >>> np.median(a)
3835 3.5
3836 >>> np.median(a, axis=0)
3837 array([6.5, 4.5, 2.5])
3838 >>> np.median(a, axis=1)
3839 array([7., 2.])
3840 >>> m = np.median(a, axis=0)
3841 >>> out = np.zeros_like(m)
3842 >>> np.median(a, axis=0, out=m)
3843 array([6.5, 4.5, 2.5])
3844 >>> m
3845 array([6.5, 4.5, 2.5])
3846 >>> b = a.copy()
3847 >>> np.median(b, axis=1, overwrite_input=True)
3848 array([7., 2.])
3849 >>> assert not np.all(a==b)
3850 >>> b = a.copy()
3851 >>> np.median(b, axis=None, overwrite_input=True)
3852 3.5
3853 >>> assert not np.all(a==b)
3855 """
3856 return _ureduce(a, func=_median, keepdims=keepdims, axis=axis, out=out,
3857 overwrite_input=overwrite_input)
3860def _median(a, axis=None, out=None, overwrite_input=False):
3861 # can't be reasonably be implemented in terms of percentile as we have to
3862 # call mean to not break astropy
3863 a = np.asanyarray(a)
3865 # Set the partition indexes
3866 if axis is None:
3867 sz = a.size
3868 else:
3869 sz = a.shape[axis]
3870 if sz % 2 == 0:
3871 szh = sz // 2
3872 kth = [szh - 1, szh]
3873 else:
3874 kth = [(sz - 1) // 2]
3875 # Check if the array contains any nan's
3876 if np.issubdtype(a.dtype, np.inexact):
3877 kth.append(-1)
3879 if overwrite_input:
3880 if axis is None:
3881 part = a.ravel()
3882 part.partition(kth)
3883 else:
3884 a.partition(kth, axis=axis)
3885 part = a
3886 else:
3887 part = partition(a, kth, axis=axis)
3889 if part.shape == ():
3890 # make 0-D arrays work
3891 return part.item()
3892 if axis is None:
3893 axis = 0
3895 indexer = [slice(None)] * part.ndim
3896 index = part.shape[axis] // 2
3897 if part.shape[axis] % 2 == 1:
3898 # index with slice to allow mean (below) to work
3899 indexer[axis] = slice(index, index+1)
3900 else:
3901 indexer[axis] = slice(index-1, index+1)
3902 indexer = tuple(indexer)
3904 # Use mean in both odd and even case to coerce data type,
3905 # using out array if needed.
3906 rout = mean(part[indexer], axis=axis, out=out)
3907 # Check if the array contains any nan's
3908 if np.issubdtype(a.dtype, np.inexact) and sz > 0:
3909 # If nans are possible, warn and replace by nans like mean would.
3910 rout = np.lib.utils._median_nancheck(part, rout, axis)
3912 return rout
3915def _percentile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
3916 method=None, keepdims=None, *, interpolation=None):
3917 return (a, q, out)
3920@array_function_dispatch(_percentile_dispatcher)
3921def percentile(a,
3922 q,
3923 axis=None,
3924 out=None,
3925 overwrite_input=False,
3926 method="linear",
3927 keepdims=False,
3928 *,
3929 interpolation=None):
3930 """
3931 Compute the q-th percentile of the data along the specified axis.
3933 Returns the q-th percentile(s) of the array elements.
3935 Parameters
3936 ----------
3937 a : array_like
3938 Input array or object that can be converted to an array.
3939 q : array_like of float
3940 Percentile or sequence of percentiles to compute, which must be between
3941 0 and 100 inclusive.
3942 axis : {int, tuple of int, None}, optional
3943 Axis or axes along which the percentiles are computed. The
3944 default is to compute the percentile(s) along a flattened
3945 version of the array.
3947 .. versionchanged:: 1.9.0
3948 A tuple of axes is supported
3949 out : ndarray, optional
3950 Alternative output array in which to place the result. It must
3951 have the same shape and buffer length as the expected output,
3952 but the type (of the output) will be cast if necessary.
3953 overwrite_input : bool, optional
3954 If True, then allow the input array `a` to be modified by intermediate
3955 calculations, to save memory. In this case, the contents of the input
3956 `a` after this function completes is undefined.
3957 method : str, optional
3958 This parameter specifies the method to use for estimating the
3959 percentile. There are many different methods, some unique to NumPy.
3960 See the notes for explanation. The options sorted by their R type
3961 as summarized in the H&F paper [1]_ are:
3963 1. 'inverted_cdf'
3964 2. 'averaged_inverted_cdf'
3965 3. 'closest_observation'
3966 4. 'interpolated_inverted_cdf'
3967 5. 'hazen'
3968 6. 'weibull'
3969 7. 'linear' (default)
3970 8. 'median_unbiased'
3971 9. 'normal_unbiased'
3973 The first three methods are discontinuous. NumPy further defines the
3974 following discontinuous variations of the default 'linear' (7.) option:
3976 * 'lower'
3977 * 'higher',
3978 * 'midpoint'
3979 * 'nearest'
3981 .. versionchanged:: 1.22.0
3982 This argument was previously called "interpolation" and only
3983 offered the "linear" default and last four options.
3985 keepdims : bool, optional
3986 If this is set to True, the axes which are reduced are left in
3987 the result as dimensions with size one. With this option, the
3988 result will broadcast correctly against the original array `a`.
3990 .. versionadded:: 1.9.0
3992 interpolation : str, optional
3993 Deprecated name for the method keyword argument.
3995 .. deprecated:: 1.22.0
3997 Returns
3998 -------
3999 percentile : scalar or ndarray
4000 If `q` is a single percentile and `axis=None`, then the result
4001 is a scalar. If multiple percentiles are given, first axis of
4002 the result corresponds to the percentiles. The other axes are
4003 the axes that remain after the reduction of `a`. If the input
4004 contains integers or floats smaller than ``float64``, the output
4005 data-type is ``float64``. Otherwise, the output data-type is the
4006 same as that of the input. If `out` is specified, that array is
4007 returned instead.
4009 See Also
4010 --------
4011 mean
4012 median : equivalent to ``percentile(..., 50)``
4013 nanpercentile
4014 quantile : equivalent to percentile, except q in the range [0, 1].
4016 Notes
4017 -----
4018 Given a vector ``V`` of length ``n``, the q-th percentile of ``V`` is
4019 the value ``q/100`` of the way from the minimum to the maximum in a
4020 sorted copy of ``V``. The values and distances of the two nearest
4021 neighbors as well as the `method` parameter will determine the
4022 percentile if the normalized ranking does not match the location of
4023 ``q`` exactly. This function is the same as the median if ``q=50``, the
4024 same as the minimum if ``q=0`` and the same as the maximum if
4025 ``q=100``.
4027 The optional `method` parameter specifies the method to use when the
4028 desired percentile lies between two indexes ``i`` and ``j = i + 1``.
4029 In that case, we first determine ``i + g``, a virtual index that lies
4030 between ``i`` and ``j``, where ``i`` is the floor and ``g`` is the
4031 fractional part of the index. The final result is, then, an interpolation
4032 of ``a[i]`` and ``a[j]`` based on ``g``. During the computation of ``g``,
4033 ``i`` and ``j`` are modified using correction constants ``alpha`` and
4034 ``beta`` whose choices depend on the ``method`` used. Finally, note that
4035 since Python uses 0-based indexing, the code subtracts another 1 from the
4036 index internally.
4038 The following formula determines the virtual index ``i + g``, the location
4039 of the percentile in the sorted sample:
4041 .. math::
4042 i + g = (q / 100) * ( n - alpha - beta + 1 ) + alpha
4044 The different methods then work as follows
4046 inverted_cdf:
4047 method 1 of H&F [1]_.
4048 This method gives discontinuous results:
4050 * if g > 0 ; then take j
4051 * if g = 0 ; then take i
4053 averaged_inverted_cdf:
4054 method 2 of H&F [1]_.
4055 This method give discontinuous results:
4057 * if g > 0 ; then take j
4058 * if g = 0 ; then average between bounds
4060 closest_observation:
4061 method 3 of H&F [1]_.
4062 This method give discontinuous results:
4064 * if g > 0 ; then take j
4065 * if g = 0 and index is odd ; then take j
4066 * if g = 0 and index is even ; then take i
4068 interpolated_inverted_cdf:
4069 method 4 of H&F [1]_.
4070 This method give continuous results using:
4072 * alpha = 0
4073 * beta = 1
4075 hazen:
4076 method 5 of H&F [1]_.
4077 This method give continuous results using:
4079 * alpha = 1/2
4080 * beta = 1/2
4082 weibull:
4083 method 6 of H&F [1]_.
4084 This method give continuous results using:
4086 * alpha = 0
4087 * beta = 0
4089 linear:
4090 method 7 of H&F [1]_.
4091 This method give continuous results using:
4093 * alpha = 1
4094 * beta = 1
4096 median_unbiased:
4097 method 8 of H&F [1]_.
4098 This method is probably the best method if the sample
4099 distribution function is unknown (see reference).
4100 This method give continuous results using:
4102 * alpha = 1/3
4103 * beta = 1/3
4105 normal_unbiased:
4106 method 9 of H&F [1]_.
4107 This method is probably the best method if the sample
4108 distribution function is known to be normal.
4109 This method give continuous results using:
4111 * alpha = 3/8
4112 * beta = 3/8
4114 lower:
4115 NumPy method kept for backwards compatibility.
4116 Takes ``i`` as the interpolation point.
4118 higher:
4119 NumPy method kept for backwards compatibility.
4120 Takes ``j`` as the interpolation point.
4122 nearest:
4123 NumPy method kept for backwards compatibility.
4124 Takes ``i`` or ``j``, whichever is nearest.
4126 midpoint:
4127 NumPy method kept for backwards compatibility.
4128 Uses ``(i + j) / 2``.
4130 Examples
4131 --------
4132 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
4133 >>> a
4134 array([[10, 7, 4],
4135 [ 3, 2, 1]])
4136 >>> np.percentile(a, 50)
4137 3.5
4138 >>> np.percentile(a, 50, axis=0)
4139 array([6.5, 4.5, 2.5])
4140 >>> np.percentile(a, 50, axis=1)
4141 array([7., 2.])
4142 >>> np.percentile(a, 50, axis=1, keepdims=True)
4143 array([[7.],
4144 [2.]])
4146 >>> m = np.percentile(a, 50, axis=0)
4147 >>> out = np.zeros_like(m)
4148 >>> np.percentile(a, 50, axis=0, out=out)
4149 array([6.5, 4.5, 2.5])
4150 >>> m
4151 array([6.5, 4.5, 2.5])
4153 >>> b = a.copy()
4154 >>> np.percentile(b, 50, axis=1, overwrite_input=True)
4155 array([7., 2.])
4156 >>> assert not np.all(a == b)
4158 The different methods can be visualized graphically:
4160 .. plot::
4162 import matplotlib.pyplot as plt
4164 a = np.arange(4)
4165 p = np.linspace(0, 100, 6001)
4166 ax = plt.gca()
4167 lines = [
4168 ('linear', '-', 'C0'),
4169 ('inverted_cdf', ':', 'C1'),
4170 # Almost the same as `inverted_cdf`:
4171 ('averaged_inverted_cdf', '-.', 'C1'),
4172 ('closest_observation', ':', 'C2'),
4173 ('interpolated_inverted_cdf', '--', 'C1'),
4174 ('hazen', '--', 'C3'),
4175 ('weibull', '-.', 'C4'),
4176 ('median_unbiased', '--', 'C5'),
4177 ('normal_unbiased', '-.', 'C6'),
4178 ]
4179 for method, style, color in lines:
4180 ax.plot(
4181 p, np.percentile(a, p, method=method),
4182 label=method, linestyle=style, color=color)
4183 ax.set(
4184 title='Percentiles for different methods and data: ' + str(a),
4185 xlabel='Percentile',
4186 ylabel='Estimated percentile value',
4187 yticks=a)
4188 ax.legend()
4189 plt.show()
4191 References
4192 ----------
4193 .. [1] R. J. Hyndman and Y. Fan,
4194 "Sample quantiles in statistical packages,"
4195 The American Statistician, 50(4), pp. 361-365, 1996
4197 """
4198 if interpolation is not None:
4199 method = _check_interpolation_as_method(
4200 method, interpolation, "percentile")
4201 q = np.true_divide(q, 100)
4202 q = asanyarray(q) # undo any decay that the ufunc performed (see gh-13105)
4203 if not _quantile_is_valid(q):
4204 raise ValueError("Percentiles must be in the range [0, 100]")
4205 return _quantile_unchecked(
4206 a, q, axis, out, overwrite_input, method, keepdims)
4209def _quantile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
4210 method=None, keepdims=None, *, interpolation=None):
4211 return (a, q, out)
4214@array_function_dispatch(_quantile_dispatcher)
4215def quantile(a,
4216 q,
4217 axis=None,
4218 out=None,
4219 overwrite_input=False,
4220 method="linear",
4221 keepdims=False,
4222 *,
4223 interpolation=None):
4224 """
4225 Compute the q-th quantile of the data along the specified axis.
4227 .. versionadded:: 1.15.0
4229 Parameters
4230 ----------
4231 a : array_like
4232 Input array or object that can be converted to an array.
4233 q : array_like of float
4234 Quantile or sequence of quantiles to compute, which must be between
4235 0 and 1 inclusive.
4236 axis : {int, tuple of int, None}, optional
4237 Axis or axes along which the quantiles are computed. The default is
4238 to compute the quantile(s) along a flattened version of the array.
4239 out : ndarray, optional
4240 Alternative output array in which to place the result. It must have
4241 the same shape and buffer length as the expected output, but the
4242 type (of the output) will be cast if necessary.
4243 overwrite_input : bool, optional
4244 If True, then allow the input array `a` to be modified by
4245 intermediate calculations, to save memory. In this case, the
4246 contents of the input `a` after this function completes is
4247 undefined.
4248 method : str, optional
4249 This parameter specifies the method to use for estimating the
4250 quantile. There are many different methods, some unique to NumPy.
4251 See the notes for explanation. The options sorted by their R type
4252 as summarized in the H&F paper [1]_ are:
4254 1. 'inverted_cdf'
4255 2. 'averaged_inverted_cdf'
4256 3. 'closest_observation'
4257 4. 'interpolated_inverted_cdf'
4258 5. 'hazen'
4259 6. 'weibull'
4260 7. 'linear' (default)
4261 8. 'median_unbiased'
4262 9. 'normal_unbiased'
4264 The first three methods are discontinuous. NumPy further defines the
4265 following discontinuous variations of the default 'linear' (7.) option:
4267 * 'lower'
4268 * 'higher',
4269 * 'midpoint'
4270 * 'nearest'
4272 .. versionchanged:: 1.22.0
4273 This argument was previously called "interpolation" and only
4274 offered the "linear" default and last four options.
4276 keepdims : bool, optional
4277 If this is set to True, the axes which are reduced are left in
4278 the result as dimensions with size one. With this option, the
4279 result will broadcast correctly against the original array `a`.
4281 interpolation : str, optional
4282 Deprecated name for the method keyword argument.
4284 .. deprecated:: 1.22.0
4286 Returns
4287 -------
4288 quantile : scalar or ndarray
4289 If `q` is a single quantile and `axis=None`, then the result
4290 is a scalar. If multiple quantiles are given, first axis of
4291 the result corresponds to the quantiles. The other axes are
4292 the axes that remain after the reduction of `a`. If the input
4293 contains integers or floats smaller than ``float64``, the output
4294 data-type is ``float64``. Otherwise, the output data-type is the
4295 same as that of the input. If `out` is specified, that array is
4296 returned instead.
4298 See Also
4299 --------
4300 mean
4301 percentile : equivalent to quantile, but with q in the range [0, 100].
4302 median : equivalent to ``quantile(..., 0.5)``
4303 nanquantile
4305 Notes
4306 -----
4307 Given a vector ``V`` of length ``n``, the q-th quantile of ``V`` is
4308 the value ``q`` of the way from the minimum to the maximum in a
4309 sorted copy of ``V``. The values and distances of the two nearest
4310 neighbors as well as the `method` parameter will determine the
4311 quantile if the normalized ranking does not match the location of
4312 ``q`` exactly. This function is the same as the median if ``q=0.5``, the
4313 same as the minimum if ``q=0.0`` and the same as the maximum if
4314 ``q=1.0``.
4316 The optional `method` parameter specifies the method to use when the
4317 desired quantile lies between two indexes ``i`` and ``j = i + 1``.
4318 In that case, we first determine ``i + g``, a virtual index that lies
4319 between ``i`` and ``j``, where ``i`` is the floor and ``g`` is the
4320 fractional part of the index. The final result is, then, an interpolation
4321 of ``a[i]`` and ``a[j]`` based on ``g``. During the computation of ``g``,
4322 ``i`` and ``j`` are modified using correction constants ``alpha`` and
4323 ``beta`` whose choices depend on the ``method`` used. Finally, note that
4324 since Python uses 0-based indexing, the code subtracts another 1 from the
4325 index internally.
4327 The following formula determines the virtual index ``i + g``, the location
4328 of the quantile in the sorted sample:
4330 .. math::
4331 i + g = q * ( n - alpha - beta + 1 ) + alpha
4333 The different methods then work as follows
4335 inverted_cdf:
4336 method 1 of H&F [1]_.
4337 This method gives discontinuous results:
4339 * if g > 0 ; then take j
4340 * if g = 0 ; then take i
4342 averaged_inverted_cdf:
4343 method 2 of H&F [1]_.
4344 This method gives discontinuous results:
4346 * if g > 0 ; then take j
4347 * if g = 0 ; then average between bounds
4349 closest_observation:
4350 method 3 of H&F [1]_.
4351 This method gives discontinuous results:
4353 * if g > 0 ; then take j
4354 * if g = 0 and index is odd ; then take j
4355 * if g = 0 and index is even ; then take i
4357 interpolated_inverted_cdf:
4358 method 4 of H&F [1]_.
4359 This method gives continuous results using:
4361 * alpha = 0
4362 * beta = 1
4364 hazen:
4365 method 5 of H&F [1]_.
4366 This method gives continuous results using:
4368 * alpha = 1/2
4369 * beta = 1/2
4371 weibull:
4372 method 6 of H&F [1]_.
4373 This method gives continuous results using:
4375 * alpha = 0
4376 * beta = 0
4378 linear:
4379 method 7 of H&F [1]_.
4380 This method gives continuous results using:
4382 * alpha = 1
4383 * beta = 1
4385 median_unbiased:
4386 method 8 of H&F [1]_.
4387 This method is probably the best method if the sample
4388 distribution function is unknown (see reference).
4389 This method gives continuous results using:
4391 * alpha = 1/3
4392 * beta = 1/3
4394 normal_unbiased:
4395 method 9 of H&F [1]_.
4396 This method is probably the best method if the sample
4397 distribution function is known to be normal.
4398 This method gives continuous results using:
4400 * alpha = 3/8
4401 * beta = 3/8
4403 lower:
4404 NumPy method kept for backwards compatibility.
4405 Takes ``i`` as the interpolation point.
4407 higher:
4408 NumPy method kept for backwards compatibility.
4409 Takes ``j`` as the interpolation point.
4411 nearest:
4412 NumPy method kept for backwards compatibility.
4413 Takes ``i`` or ``j``, whichever is nearest.
4415 midpoint:
4416 NumPy method kept for backwards compatibility.
4417 Uses ``(i + j) / 2``.
4419 Examples
4420 --------
4421 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
4422 >>> a
4423 array([[10, 7, 4],
4424 [ 3, 2, 1]])
4425 >>> np.quantile(a, 0.5)
4426 3.5
4427 >>> np.quantile(a, 0.5, axis=0)
4428 array([6.5, 4.5, 2.5])
4429 >>> np.quantile(a, 0.5, axis=1)
4430 array([7., 2.])
4431 >>> np.quantile(a, 0.5, axis=1, keepdims=True)
4432 array([[7.],
4433 [2.]])
4434 >>> m = np.quantile(a, 0.5, axis=0)
4435 >>> out = np.zeros_like(m)
4436 >>> np.quantile(a, 0.5, axis=0, out=out)
4437 array([6.5, 4.5, 2.5])
4438 >>> m
4439 array([6.5, 4.5, 2.5])
4440 >>> b = a.copy()
4441 >>> np.quantile(b, 0.5, axis=1, overwrite_input=True)
4442 array([7., 2.])
4443 >>> assert not np.all(a == b)
4445 See also `numpy.percentile` for a visualization of most methods.
4447 References
4448 ----------
4449 .. [1] R. J. Hyndman and Y. Fan,
4450 "Sample quantiles in statistical packages,"
4451 The American Statistician, 50(4), pp. 361-365, 1996
4453 """
4454 if interpolation is not None:
4455 method = _check_interpolation_as_method(
4456 method, interpolation, "quantile")
4458 q = np.asanyarray(q)
4459 if not _quantile_is_valid(q):
4460 raise ValueError("Quantiles must be in the range [0, 1]")
4461 return _quantile_unchecked(
4462 a, q, axis, out, overwrite_input, method, keepdims)
4465def _quantile_unchecked(a,
4466 q,
4467 axis=None,
4468 out=None,
4469 overwrite_input=False,
4470 method="linear",
4471 keepdims=False):
4472 """Assumes that q is in [0, 1], and is an ndarray"""
4473 return _ureduce(a,
4474 func=_quantile_ureduce_func,
4475 q=q,
4476 keepdims=keepdims,
4477 axis=axis,
4478 out=out,
4479 overwrite_input=overwrite_input,
4480 method=method)
4483def _quantile_is_valid(q):
4484 # avoid expensive reductions, relevant for arrays with < O(1000) elements
4485 if q.ndim == 1 and q.size < 10:
4486 for i in range(q.size):
4487 if not (0.0 <= q[i] <= 1.0):
4488 return False
4489 else:
4490 if not (np.all(0 <= q) and np.all(q <= 1)):
4491 return False
4492 return True
4495def _check_interpolation_as_method(method, interpolation, fname):
4496 # Deprecated NumPy 1.22, 2021-11-08
4497 warnings.warn(
4498 f"the `interpolation=` argument to {fname} was renamed to "
4499 "`method=`, which has additional options.\n"
4500 "Users of the modes 'nearest', 'lower', 'higher', or "
4501 "'midpoint' are encouraged to review the method they used. "
4502 "(Deprecated NumPy 1.22)",
4503 DeprecationWarning, stacklevel=4)
4504 if method != "linear":
4505 # sanity check, we assume this basically never happens
4506 raise TypeError(
4507 "You shall not pass both `method` and `interpolation`!\n"
4508 "(`interpolation` is Deprecated in favor of `method`)")
4509 return interpolation
4512def _compute_virtual_index(n, quantiles, alpha: float, beta: float):
4513 """
4514 Compute the floating point indexes of an array for the linear
4515 interpolation of quantiles.
4516 n : array_like
4517 The sample sizes.
4518 quantiles : array_like
4519 The quantiles values.
4520 alpha : float
4521 A constant used to correct the index computed.
4522 beta : float
4523 A constant used to correct the index computed.
4525 alpha and beta values depend on the chosen method
4526 (see quantile documentation)
4528 Reference:
4529 Hyndman&Fan paper "Sample Quantiles in Statistical Packages",
4530 DOI: 10.1080/00031305.1996.10473566
4531 """
4532 return n * quantiles + (
4533 alpha + quantiles * (1 - alpha - beta)
4534 ) - 1
4537def _get_gamma(virtual_indexes, previous_indexes, method):
4538 """
4539 Compute gamma (a.k.a 'm' or 'weight') for the linear interpolation
4540 of quantiles.
4542 virtual_indexes : array_like
4543 The indexes where the percentile is supposed to be found in the sorted
4544 sample.
4545 previous_indexes : array_like
4546 The floor values of virtual_indexes.
4547 interpolation : dict
4548 The interpolation method chosen, which may have a specific rule
4549 modifying gamma.
4551 gamma is usually the fractional part of virtual_indexes but can be modified
4552 by the interpolation method.
4553 """
4554 gamma = np.asanyarray(virtual_indexes - previous_indexes)
4555 gamma = method["fix_gamma"](gamma, virtual_indexes)
4556 return np.asanyarray(gamma)
4559def _lerp(a, b, t, out=None):
4560 """
4561 Compute the linear interpolation weighted by gamma on each point of
4562 two same shape array.
4564 a : array_like
4565 Left bound.
4566 b : array_like
4567 Right bound.
4568 t : array_like
4569 The interpolation weight.
4570 out : array_like
4571 Output array.
4572 """
4573 diff_b_a = subtract(b, a)
4574 # asanyarray is a stop-gap until gh-13105
4575 lerp_interpolation = asanyarray(add(a, diff_b_a * t, out=out))
4576 subtract(b, diff_b_a * (1 - t), out=lerp_interpolation, where=t >= 0.5)
4577 if lerp_interpolation.ndim == 0 and out is None:
4578 lerp_interpolation = lerp_interpolation[()] # unpack 0d arrays
4579 return lerp_interpolation
4582def _get_gamma_mask(shape, default_value, conditioned_value, where):
4583 out = np.full(shape, default_value)
4584 np.copyto(out, conditioned_value, where=where, casting="unsafe")
4585 return out
4588def _discret_interpolation_to_boundaries(index, gamma_condition_fun):
4589 previous = np.floor(index)
4590 next = previous + 1
4591 gamma = index - previous
4592 res = _get_gamma_mask(shape=index.shape,
4593 default_value=next,
4594 conditioned_value=previous,
4595 where=gamma_condition_fun(gamma, index)
4596 ).astype(np.intp)
4597 # Some methods can lead to out-of-bound integers, clip them:
4598 res[res < 0] = 0
4599 return res
4602def _closest_observation(n, quantiles):
4603 gamma_fun = lambda gamma, index: (gamma == 0) & (np.floor(index) % 2 == 0)
4604 return _discret_interpolation_to_boundaries((n * quantiles) - 1 - 0.5,
4605 gamma_fun)
4608def _inverted_cdf(n, quantiles):
4609 gamma_fun = lambda gamma, _: (gamma == 0)
4610 return _discret_interpolation_to_boundaries((n * quantiles) - 1,
4611 gamma_fun)
4614def _quantile_ureduce_func(
4615 a: np.array,
4616 q: np.array,
4617 axis: int = None,
4618 out=None,
4619 overwrite_input: bool = False,
4620 method="linear",
4621) -> np.array:
4622 if q.ndim > 2:
4623 # The code below works fine for nd, but it might not have useful
4624 # semantics. For now, keep the supported dimensions the same as it was
4625 # before.
4626 raise ValueError("q must be a scalar or 1d")
4627 if overwrite_input:
4628 if axis is None:
4629 axis = 0
4630 arr = a.ravel()
4631 else:
4632 arr = a
4633 else:
4634 if axis is None:
4635 axis = 0
4636 arr = a.flatten()
4637 else:
4638 arr = a.copy()
4639 result = _quantile(arr,
4640 quantiles=q,
4641 axis=axis,
4642 method=method,
4643 out=out)
4644 return result
4647def _get_indexes(arr, virtual_indexes, valid_values_count):
4648 """
4649 Get the valid indexes of arr neighbouring virtual_indexes.
4650 Note
4651 This is a companion function to linear interpolation of
4652 Quantiles
4654 Returns
4655 -------
4656 (previous_indexes, next_indexes): Tuple
4657 A Tuple of virtual_indexes neighbouring indexes
4658 """
4659 previous_indexes = np.asanyarray(np.floor(virtual_indexes))
4660 next_indexes = np.asanyarray(previous_indexes + 1)
4661 indexes_above_bounds = virtual_indexes >= valid_values_count - 1
4662 # When indexes is above max index, take the max value of the array
4663 if indexes_above_bounds.any():
4664 previous_indexes[indexes_above_bounds] = -1
4665 next_indexes[indexes_above_bounds] = -1
4666 # When indexes is below min index, take the min value of the array
4667 indexes_below_bounds = virtual_indexes < 0
4668 if indexes_below_bounds.any():
4669 previous_indexes[indexes_below_bounds] = 0
4670 next_indexes[indexes_below_bounds] = 0
4671 if np.issubdtype(arr.dtype, np.inexact):
4672 # After the sort, slices having NaNs will have for last element a NaN
4673 virtual_indexes_nans = np.isnan(virtual_indexes)
4674 if virtual_indexes_nans.any():
4675 previous_indexes[virtual_indexes_nans] = -1
4676 next_indexes[virtual_indexes_nans] = -1
4677 previous_indexes = previous_indexes.astype(np.intp)
4678 next_indexes = next_indexes.astype(np.intp)
4679 return previous_indexes, next_indexes
4682def _quantile(
4683 arr: np.array,
4684 quantiles: np.array,
4685 axis: int = -1,
4686 method="linear",
4687 out=None,
4688):
4689 """
4690 Private function that doesn't support extended axis or keepdims.
4691 These methods are extended to this function using _ureduce
4692 See nanpercentile for parameter usage
4693 It computes the quantiles of the array for the given axis.
4694 A linear interpolation is performed based on the `interpolation`.
4696 By default, the method is "linear" where alpha == beta == 1 which
4697 performs the 7th method of Hyndman&Fan.
4698 With "median_unbiased" we get alpha == beta == 1/3
4699 thus the 8th method of Hyndman&Fan.
4700 """
4701 # --- Setup
4702 arr = np.asanyarray(arr)
4703 values_count = arr.shape[axis]
4704 # The dimensions of `q` are prepended to the output shape, so we need the
4705 # axis being sampled from `arr` to be last.
4706 DATA_AXIS = 0
4707 if axis != DATA_AXIS: # But moveaxis is slow, so only call it if axis!=0.
4708 arr = np.moveaxis(arr, axis, destination=DATA_AXIS)
4709 # --- Computation of indexes
4710 # Index where to find the value in the sorted array.
4711 # Virtual because it is a floating point value, not an valid index.
4712 # The nearest neighbours are used for interpolation
4713 try:
4714 method = _QuantileMethods[method]
4715 except KeyError:
4716 raise ValueError(
4717 f"{method!r} is not a valid method. Use one of: "
4718 f"{_QuantileMethods.keys()}") from None
4719 virtual_indexes = method["get_virtual_index"](values_count, quantiles)
4720 virtual_indexes = np.asanyarray(virtual_indexes)
4721 if np.issubdtype(virtual_indexes.dtype, np.integer):
4722 # No interpolation needed, take the points along axis
4723 if np.issubdtype(arr.dtype, np.inexact):
4724 # may contain nan, which would sort to the end
4725 arr.partition(concatenate((virtual_indexes.ravel(), [-1])), axis=0)
4726 slices_having_nans = np.isnan(arr[-1])
4727 else:
4728 # cannot contain nan
4729 arr.partition(virtual_indexes.ravel(), axis=0)
4730 slices_having_nans = np.array(False, dtype=bool)
4731 result = take(arr, virtual_indexes, axis=0, out=out)
4732 else:
4733 previous_indexes, next_indexes = _get_indexes(arr,
4734 virtual_indexes,
4735 values_count)
4736 # --- Sorting
4737 arr.partition(
4738 np.unique(np.concatenate(([0, -1],
4739 previous_indexes.ravel(),
4740 next_indexes.ravel(),
4741 ))),
4742 axis=DATA_AXIS)
4743 if np.issubdtype(arr.dtype, np.inexact):
4744 slices_having_nans = np.isnan(
4745 take(arr, indices=-1, axis=DATA_AXIS)
4746 )
4747 else:
4748 slices_having_nans = None
4749 # --- Get values from indexes
4750 previous = np.take(arr, previous_indexes, axis=DATA_AXIS)
4751 next = np.take(arr, next_indexes, axis=DATA_AXIS)
4752 # --- Linear interpolation
4753 gamma = _get_gamma(virtual_indexes, previous_indexes, method)
4754 result_shape = virtual_indexes.shape + (1,) * (arr.ndim - 1)
4755 gamma = gamma.reshape(result_shape)
4756 result = _lerp(previous,
4757 next,
4758 gamma,
4759 out=out)
4760 if np.any(slices_having_nans):
4761 if result.ndim == 0 and out is None:
4762 # can't write to a scalar
4763 result = arr.dtype.type(np.nan)
4764 else:
4765 result[..., slices_having_nans] = np.nan
4766 return result
4769def _trapz_dispatcher(y, x=None, dx=None, axis=None):
4770 return (y, x)
4773@array_function_dispatch(_trapz_dispatcher)
4774def trapz(y, x=None, dx=1.0, axis=-1):
4775 r"""
4776 Integrate along the given axis using the composite trapezoidal rule.
4778 If `x` is provided, the integration happens in sequence along its
4779 elements - they are not sorted.
4781 Integrate `y` (`x`) along each 1d slice on the given axis, compute
4782 :math:`\int y(x) dx`.
4783 When `x` is specified, this integrates along the parametric curve,
4784 computing :math:`\int_t y(t) dt =
4785 \int_t y(t) \left.\frac{dx}{dt}\right|_{x=x(t)} dt`.
4787 Parameters
4788 ----------
4789 y : array_like
4790 Input array to integrate.
4791 x : array_like, optional
4792 The sample points corresponding to the `y` values. If `x` is None,
4793 the sample points are assumed to be evenly spaced `dx` apart. The
4794 default is None.
4795 dx : scalar, optional
4796 The spacing between sample points when `x` is None. The default is 1.
4797 axis : int, optional
4798 The axis along which to integrate.
4800 Returns
4801 -------
4802 trapz : float or ndarray
4803 Definite integral of `y` = n-dimensional array as approximated along
4804 a single axis by the trapezoidal rule. If `y` is a 1-dimensional array,
4805 then the result is a float. If `n` is greater than 1, then the result
4806 is an `n`-1 dimensional array.
4808 See Also
4809 --------
4810 sum, cumsum
4812 Notes
4813 -----
4814 Image [2]_ illustrates trapezoidal rule -- y-axis locations of points
4815 will be taken from `y` array, by default x-axis distances between
4816 points will be 1.0, alternatively they can be provided with `x` array
4817 or with `dx` scalar. Return value will be equal to combined area under
4818 the red lines.
4821 References
4822 ----------
4823 .. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule
4825 .. [2] Illustration image:
4826 https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
4828 Examples
4829 --------
4830 >>> np.trapz([1,2,3])
4831 4.0
4832 >>> np.trapz([1,2,3], x=[4,6,8])
4833 8.0
4834 >>> np.trapz([1,2,3], dx=2)
4835 8.0
4837 Using a decreasing `x` corresponds to integrating in reverse:
4839 >>> np.trapz([1,2,3], x=[8,6,4])
4840 -8.0
4842 More generally `x` is used to integrate along a parametric curve.
4843 This finds the area of a circle, noting we repeat the sample which closes
4844 the curve:
4846 >>> theta = np.linspace(0, 2 * np.pi, num=1000, endpoint=True)
4847 >>> np.trapz(np.cos(theta), x=np.sin(theta))
4848 3.141571941375841
4850 >>> a = np.arange(6).reshape(2, 3)
4851 >>> a
4852 array([[0, 1, 2],
4853 [3, 4, 5]])
4854 >>> np.trapz(a, axis=0)
4855 array([1.5, 2.5, 3.5])
4856 >>> np.trapz(a, axis=1)
4857 array([2., 8.])
4858 """
4859 y = asanyarray(y)
4860 if x is None:
4861 d = dx
4862 else:
4863 x = asanyarray(x)
4864 if x.ndim == 1:
4865 d = diff(x)
4866 # reshape to correct shape
4867 shape = [1]*y.ndim
4868 shape[axis] = d.shape[0]
4869 d = d.reshape(shape)
4870 else:
4871 d = diff(x, axis=axis)
4872 nd = y.ndim
4873 slice1 = [slice(None)]*nd
4874 slice2 = [slice(None)]*nd
4875 slice1[axis] = slice(1, None)
4876 slice2[axis] = slice(None, -1)
4877 try:
4878 ret = (d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0).sum(axis)
4879 except ValueError:
4880 # Operations didn't work, cast to ndarray
4881 d = np.asarray(d)
4882 y = np.asarray(y)
4883 ret = add.reduce(d * (y[tuple(slice1)]+y[tuple(slice2)])/2.0, axis)
4884 return ret
4887def _meshgrid_dispatcher(*xi, copy=None, sparse=None, indexing=None):
4888 return xi
4891# Based on scitools meshgrid
4892@array_function_dispatch(_meshgrid_dispatcher)
4893def meshgrid(*xi, copy=True, sparse=False, indexing='xy'):
4894 """
4895 Return coordinate matrices from coordinate vectors.
4897 Make N-D coordinate arrays for vectorized evaluations of
4898 N-D scalar/vector fields over N-D grids, given
4899 one-dimensional coordinate arrays x1, x2,..., xn.
4901 .. versionchanged:: 1.9
4902 1-D and 0-D cases are allowed.
4904 Parameters
4905 ----------
4906 x1, x2,..., xn : array_like
4907 1-D arrays representing the coordinates of a grid.
4908 indexing : {'xy', 'ij'}, optional
4909 Cartesian ('xy', default) or matrix ('ij') indexing of output.
4910 See Notes for more details.
4912 .. versionadded:: 1.7.0
4913 sparse : bool, optional
4914 If True the shape of the returned coordinate array for dimension *i*
4915 is reduced from ``(N1, ..., Ni, ... Nn)`` to
4916 ``(1, ..., 1, Ni, 1, ..., 1)``. These sparse coordinate grids are
4917 intended to be use with :ref:`basics.broadcasting`. When all
4918 coordinates are used in an expression, broadcasting still leads to a
4919 fully-dimensonal result array.
4921 Default is False.
4923 .. versionadded:: 1.7.0
4924 copy : bool, optional
4925 If False, a view into the original arrays are returned in order to
4926 conserve memory. Default is True. Please note that
4927 ``sparse=False, copy=False`` will likely return non-contiguous
4928 arrays. Furthermore, more than one element of a broadcast array
4929 may refer to a single memory location. If you need to write to the
4930 arrays, make copies first.
4932 .. versionadded:: 1.7.0
4934 Returns
4935 -------
4936 X1, X2,..., XN : ndarray
4937 For vectors `x1`, `x2`,..., `xn` with lengths ``Ni=len(xi)``,
4938 returns ``(N1, N2, N3,..., Nn)`` shaped arrays if indexing='ij'
4939 or ``(N2, N1, N3,..., Nn)`` shaped arrays if indexing='xy'
4940 with the elements of `xi` repeated to fill the matrix along
4941 the first dimension for `x1`, the second for `x2` and so on.
4943 Notes
4944 -----
4945 This function supports both indexing conventions through the indexing
4946 keyword argument. Giving the string 'ij' returns a meshgrid with
4947 matrix indexing, while 'xy' returns a meshgrid with Cartesian indexing.
4948 In the 2-D case with inputs of length M and N, the outputs are of shape
4949 (N, M) for 'xy' indexing and (M, N) for 'ij' indexing. In the 3-D case
4950 with inputs of length M, N and P, outputs are of shape (N, M, P) for
4951 'xy' indexing and (M, N, P) for 'ij' indexing. The difference is
4952 illustrated by the following code snippet::
4954 xv, yv = np.meshgrid(x, y, indexing='ij')
4955 for i in range(nx):
4956 for j in range(ny):
4957 # treat xv[i,j], yv[i,j]
4959 xv, yv = np.meshgrid(x, y, indexing='xy')
4960 for i in range(nx):
4961 for j in range(ny):
4962 # treat xv[j,i], yv[j,i]
4964 In the 1-D and 0-D case, the indexing and sparse keywords have no effect.
4966 See Also
4967 --------
4968 mgrid : Construct a multi-dimensional "meshgrid" using indexing notation.
4969 ogrid : Construct an open multi-dimensional "meshgrid" using indexing
4970 notation.
4971 how-to-index
4973 Examples
4974 --------
4975 >>> nx, ny = (3, 2)
4976 >>> x = np.linspace(0, 1, nx)
4977 >>> y = np.linspace(0, 1, ny)
4978 >>> xv, yv = np.meshgrid(x, y)
4979 >>> xv
4980 array([[0. , 0.5, 1. ],
4981 [0. , 0.5, 1. ]])
4982 >>> yv
4983 array([[0., 0., 0.],
4984 [1., 1., 1.]])
4986 The result of `meshgrid` is a coordinate grid:
4988 >>> import matplotlib.pyplot as plt
4989 >>> plt.plot(xv, yv, marker='o', color='k', linestyle='none')
4990 >>> plt.show()
4992 You can create sparse output arrays to save memory and computation time.
4994 >>> xv, yv = np.meshgrid(x, y, sparse=True)
4995 >>> xv
4996 array([[0. , 0.5, 1. ]])
4997 >>> yv
4998 array([[0.],
4999 [1.]])
5001 `meshgrid` is very useful to evaluate functions on a grid. If the
5002 function depends on all coordinates, both dense and sparse outputs can be
5003 used.
5005 >>> x = np.linspace(-5, 5, 101)
5006 >>> y = np.linspace(-5, 5, 101)
5007 >>> # full coordinate arrays
5008 >>> xx, yy = np.meshgrid(x, y)
5009 >>> zz = np.sqrt(xx**2 + yy**2)
5010 >>> xx.shape, yy.shape, zz.shape
5011 ((101, 101), (101, 101), (101, 101))
5012 >>> # sparse coordinate arrays
5013 >>> xs, ys = np.meshgrid(x, y, sparse=True)
5014 >>> zs = np.sqrt(xs**2 + ys**2)
5015 >>> xs.shape, ys.shape, zs.shape
5016 ((1, 101), (101, 1), (101, 101))
5017 >>> np.array_equal(zz, zs)
5018 True
5020 >>> h = plt.contourf(x, y, zs)
5021 >>> plt.axis('scaled')
5022 >>> plt.colorbar()
5023 >>> plt.show()
5024 """
5025 ndim = len(xi)
5027 if indexing not in ['xy', 'ij']:
5028 raise ValueError(
5029 "Valid values for `indexing` are 'xy' and 'ij'.")
5031 s0 = (1,) * ndim
5032 output = [np.asanyarray(x).reshape(s0[:i] + (-1,) + s0[i + 1:])
5033 for i, x in enumerate(xi)]
5035 if indexing == 'xy' and ndim > 1:
5036 # switch first and second axis
5037 output[0].shape = (1, -1) + s0[2:]
5038 output[1].shape = (-1, 1) + s0[2:]
5040 if not sparse:
5041 # Return the full N-D matrix (not only the 1-D vector)
5042 output = np.broadcast_arrays(*output, subok=True)
5044 if copy:
5045 output = [x.copy() for x in output]
5047 return output
5050def _delete_dispatcher(arr, obj, axis=None):
5051 return (arr, obj)
5054@array_function_dispatch(_delete_dispatcher)
5055def delete(arr, obj, axis=None):
5056 """
5057 Return a new array with sub-arrays along an axis deleted. For a one
5058 dimensional array, this returns those entries not returned by
5059 `arr[obj]`.
5061 Parameters
5062 ----------
5063 arr : array_like
5064 Input array.
5065 obj : slice, int or array of ints
5066 Indicate indices of sub-arrays to remove along the specified axis.
5068 .. versionchanged:: 1.19.0
5069 Boolean indices are now treated as a mask of elements to remove,
5070 rather than being cast to the integers 0 and 1.
5072 axis : int, optional
5073 The axis along which to delete the subarray defined by `obj`.
5074 If `axis` is None, `obj` is applied to the flattened array.
5076 Returns
5077 -------
5078 out : ndarray
5079 A copy of `arr` with the elements specified by `obj` removed. Note
5080 that `delete` does not occur in-place. If `axis` is None, `out` is
5081 a flattened array.
5083 See Also
5084 --------
5085 insert : Insert elements into an array.
5086 append : Append elements at the end of an array.
5088 Notes
5089 -----
5090 Often it is preferable to use a boolean mask. For example:
5092 >>> arr = np.arange(12) + 1
5093 >>> mask = np.ones(len(arr), dtype=bool)
5094 >>> mask[[0,2,4]] = False
5095 >>> result = arr[mask,...]
5097 Is equivalent to ``np.delete(arr, [0,2,4], axis=0)``, but allows further
5098 use of `mask`.
5100 Examples
5101 --------
5102 >>> arr = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]])
5103 >>> arr
5104 array([[ 1, 2, 3, 4],
5105 [ 5, 6, 7, 8],
5106 [ 9, 10, 11, 12]])
5107 >>> np.delete(arr, 1, 0)
5108 array([[ 1, 2, 3, 4],
5109 [ 9, 10, 11, 12]])
5111 >>> np.delete(arr, np.s_[::2], 1)
5112 array([[ 2, 4],
5113 [ 6, 8],
5114 [10, 12]])
5115 >>> np.delete(arr, [1,3,5], None)
5116 array([ 1, 3, 5, 7, 8, 9, 10, 11, 12])
5118 """
5119 wrap = None
5120 if type(arr) is not ndarray:
5121 try:
5122 wrap = arr.__array_wrap__
5123 except AttributeError:
5124 pass
5126 arr = asarray(arr)
5127 ndim = arr.ndim
5128 arrorder = 'F' if arr.flags.fnc else 'C'
5129 if axis is None:
5130 if ndim != 1:
5131 arr = arr.ravel()
5132 # needed for np.matrix, which is still not 1d after being ravelled
5133 ndim = arr.ndim
5134 axis = ndim - 1
5135 else:
5136 axis = normalize_axis_index(axis, ndim)
5138 slobj = [slice(None)]*ndim
5139 N = arr.shape[axis]
5140 newshape = list(arr.shape)
5142 if isinstance(obj, slice):
5143 start, stop, step = obj.indices(N)
5144 xr = range(start, stop, step)
5145 numtodel = len(xr)
5147 if numtodel <= 0:
5148 if wrap:
5149 return wrap(arr.copy(order=arrorder))
5150 else:
5151 return arr.copy(order=arrorder)
5153 # Invert if step is negative:
5154 if step < 0:
5155 step = -step
5156 start = xr[-1]
5157 stop = xr[0] + 1
5159 newshape[axis] -= numtodel
5160 new = empty(newshape, arr.dtype, arrorder)
5161 # copy initial chunk
5162 if start == 0:
5163 pass
5164 else:
5165 slobj[axis] = slice(None, start)
5166 new[tuple(slobj)] = arr[tuple(slobj)]
5167 # copy end chunk
5168 if stop == N:
5169 pass
5170 else:
5171 slobj[axis] = slice(stop-numtodel, None)
5172 slobj2 = [slice(None)]*ndim
5173 slobj2[axis] = slice(stop, None)
5174 new[tuple(slobj)] = arr[tuple(slobj2)]
5175 # copy middle pieces
5176 if step == 1:
5177 pass
5178 else: # use array indexing.
5179 keep = ones(stop-start, dtype=bool)
5180 keep[:stop-start:step] = False
5181 slobj[axis] = slice(start, stop-numtodel)
5182 slobj2 = [slice(None)]*ndim
5183 slobj2[axis] = slice(start, stop)
5184 arr = arr[tuple(slobj2)]
5185 slobj2[axis] = keep
5186 new[tuple(slobj)] = arr[tuple(slobj2)]
5187 if wrap:
5188 return wrap(new)
5189 else:
5190 return new
5192 if isinstance(obj, (int, integer)) and not isinstance(obj, bool):
5193 single_value = True
5194 else:
5195 single_value = False
5196 _obj = obj
5197 obj = np.asarray(obj)
5198 # `size == 0` to allow empty lists similar to indexing, but (as there)
5199 # is really too generic:
5200 if obj.size == 0 and not isinstance(_obj, np.ndarray):
5201 obj = obj.astype(intp)
5202 elif obj.size == 1 and obj.dtype.kind in "ui":
5203 # For a size 1 integer array we can use the single-value path
5204 # (most dtypes, except boolean, should just fail later).
5205 obj = obj.item()
5206 single_value = True
5208 if single_value:
5209 # optimization for a single value
5210 if (obj < -N or obj >= N):
5211 raise IndexError(
5212 "index %i is out of bounds for axis %i with "
5213 "size %i" % (obj, axis, N))
5214 if (obj < 0):
5215 obj += N
5216 newshape[axis] -= 1
5217 new = empty(newshape, arr.dtype, arrorder)
5218 slobj[axis] = slice(None, obj)
5219 new[tuple(slobj)] = arr[tuple(slobj)]
5220 slobj[axis] = slice(obj, None)
5221 slobj2 = [slice(None)]*ndim
5222 slobj2[axis] = slice(obj+1, None)
5223 new[tuple(slobj)] = arr[tuple(slobj2)]
5224 else:
5225 if obj.dtype == bool:
5226 if obj.shape != (N,):
5227 raise ValueError('boolean array argument obj to delete '
5228 'must be one dimensional and match the axis '
5229 'length of {}'.format(N))
5231 # optimization, the other branch is slower
5232 keep = ~obj
5233 else:
5234 keep = ones(N, dtype=bool)
5235 keep[obj,] = False
5237 slobj[axis] = keep
5238 new = arr[tuple(slobj)]
5240 if wrap:
5241 return wrap(new)
5242 else:
5243 return new
5246def _insert_dispatcher(arr, obj, values, axis=None):
5247 return (arr, obj, values)
5250@array_function_dispatch(_insert_dispatcher)
5251def insert(arr, obj, values, axis=None):
5252 """
5253 Insert values along the given axis before the given indices.
5255 Parameters
5256 ----------
5257 arr : array_like
5258 Input array.
5259 obj : int, slice or sequence of ints
5260 Object that defines the index or indices before which `values` is
5261 inserted.
5263 .. versionadded:: 1.8.0
5265 Support for multiple insertions when `obj` is a single scalar or a
5266 sequence with one element (similar to calling insert multiple
5267 times).
5268 values : array_like
5269 Values to insert into `arr`. If the type of `values` is different
5270 from that of `arr`, `values` is converted to the type of `arr`.
5271 `values` should be shaped so that ``arr[...,obj,...] = values``
5272 is legal.
5273 axis : int, optional
5274 Axis along which to insert `values`. If `axis` is None then `arr`
5275 is flattened first.
5277 Returns
5278 -------
5279 out : ndarray
5280 A copy of `arr` with `values` inserted. Note that `insert`
5281 does not occur in-place: a new array is returned. If
5282 `axis` is None, `out` is a flattened array.
5284 See Also
5285 --------
5286 append : Append elements at the end of an array.
5287 concatenate : Join a sequence of arrays along an existing axis.
5288 delete : Delete elements from an array.
5290 Notes
5291 -----
5292 Note that for higher dimensional inserts ``obj=0`` behaves very different
5293 from ``obj=[0]`` just like ``arr[:,0,:] = values`` is different from
5294 ``arr[:,[0],:] = values``.
5296 Examples
5297 --------
5298 >>> a = np.array([[1, 1], [2, 2], [3, 3]])
5299 >>> a
5300 array([[1, 1],
5301 [2, 2],
5302 [3, 3]])
5303 >>> np.insert(a, 1, 5)
5304 array([1, 5, 1, ..., 2, 3, 3])
5305 >>> np.insert(a, 1, 5, axis=1)
5306 array([[1, 5, 1],
5307 [2, 5, 2],
5308 [3, 5, 3]])
5310 Difference between sequence and scalars:
5312 >>> np.insert(a, [1], [[1],[2],[3]], axis=1)
5313 array([[1, 1, 1],
5314 [2, 2, 2],
5315 [3, 3, 3]])
5316 >>> np.array_equal(np.insert(a, 1, [1, 2, 3], axis=1),
5317 ... np.insert(a, [1], [[1],[2],[3]], axis=1))
5318 True
5320 >>> b = a.flatten()
5321 >>> b
5322 array([1, 1, 2, 2, 3, 3])
5323 >>> np.insert(b, [2, 2], [5, 6])
5324 array([1, 1, 5, ..., 2, 3, 3])
5326 >>> np.insert(b, slice(2, 4), [5, 6])
5327 array([1, 1, 5, ..., 2, 3, 3])
5329 >>> np.insert(b, [2, 2], [7.13, False]) # type casting
5330 array([1, 1, 7, ..., 2, 3, 3])
5332 >>> x = np.arange(8).reshape(2, 4)
5333 >>> idx = (1, 3)
5334 >>> np.insert(x, idx, 999, axis=1)
5335 array([[ 0, 999, 1, 2, 999, 3],
5336 [ 4, 999, 5, 6, 999, 7]])
5338 """
5339 wrap = None
5340 if type(arr) is not ndarray:
5341 try:
5342 wrap = arr.__array_wrap__
5343 except AttributeError:
5344 pass
5346 arr = asarray(arr)
5347 ndim = arr.ndim
5348 arrorder = 'F' if arr.flags.fnc else 'C'
5349 if axis is None:
5350 if ndim != 1:
5351 arr = arr.ravel()
5352 # needed for np.matrix, which is still not 1d after being ravelled
5353 ndim = arr.ndim
5354 axis = ndim - 1
5355 else:
5356 axis = normalize_axis_index(axis, ndim)
5357 slobj = [slice(None)]*ndim
5358 N = arr.shape[axis]
5359 newshape = list(arr.shape)
5361 if isinstance(obj, slice):
5362 # turn it into a range object
5363 indices = arange(*obj.indices(N), dtype=intp)
5364 else:
5365 # need to copy obj, because indices will be changed in-place
5366 indices = np.array(obj)
5367 if indices.dtype == bool:
5368 # See also delete
5369 # 2012-10-11, NumPy 1.8
5370 warnings.warn(
5371 "in the future insert will treat boolean arrays and "
5372 "array-likes as a boolean index instead of casting it to "
5373 "integer", FutureWarning, stacklevel=3)
5374 indices = indices.astype(intp)
5375 # Code after warning period:
5376 #if obj.ndim != 1:
5377 # raise ValueError('boolean array argument obj to insert '
5378 # 'must be one dimensional')
5379 #indices = np.flatnonzero(obj)
5380 elif indices.ndim > 1:
5381 raise ValueError(
5382 "index array argument obj to insert must be one dimensional "
5383 "or scalar")
5384 if indices.size == 1:
5385 index = indices.item()
5386 if index < -N or index > N:
5387 raise IndexError(f"index {obj} is out of bounds for axis {axis} "
5388 f"with size {N}")
5389 if (index < 0):
5390 index += N
5392 # There are some object array corner cases here, but we cannot avoid
5393 # that:
5394 values = array(values, copy=False, ndmin=arr.ndim, dtype=arr.dtype)
5395 if indices.ndim == 0:
5396 # broadcasting is very different here, since a[:,0,:] = ... behaves
5397 # very different from a[:,[0],:] = ...! This changes values so that
5398 # it works likes the second case. (here a[:,0:1,:])
5399 values = np.moveaxis(values, 0, axis)
5400 numnew = values.shape[axis]
5401 newshape[axis] += numnew
5402 new = empty(newshape, arr.dtype, arrorder)
5403 slobj[axis] = slice(None, index)
5404 new[tuple(slobj)] = arr[tuple(slobj)]
5405 slobj[axis] = slice(index, index+numnew)
5406 new[tuple(slobj)] = values
5407 slobj[axis] = slice(index+numnew, None)
5408 slobj2 = [slice(None)] * ndim
5409 slobj2[axis] = slice(index, None)
5410 new[tuple(slobj)] = arr[tuple(slobj2)]
5411 if wrap:
5412 return wrap(new)
5413 return new
5414 elif indices.size == 0 and not isinstance(obj, np.ndarray):
5415 # Can safely cast the empty list to intp
5416 indices = indices.astype(intp)
5418 indices[indices < 0] += N
5420 numnew = len(indices)
5421 order = indices.argsort(kind='mergesort') # stable sort
5422 indices[order] += np.arange(numnew)
5424 newshape[axis] += numnew
5425 old_mask = ones(newshape[axis], dtype=bool)
5426 old_mask[indices] = False
5428 new = empty(newshape, arr.dtype, arrorder)
5429 slobj2 = [slice(None)]*ndim
5430 slobj[axis] = indices
5431 slobj2[axis] = old_mask
5432 new[tuple(slobj)] = values
5433 new[tuple(slobj2)] = arr
5435 if wrap:
5436 return wrap(new)
5437 return new
5440def _append_dispatcher(arr, values, axis=None):
5441 return (arr, values)
5444@array_function_dispatch(_append_dispatcher)
5445def append(arr, values, axis=None):
5446 """
5447 Append values to the end of an array.
5449 Parameters
5450 ----------
5451 arr : array_like
5452 Values are appended to a copy of this array.
5453 values : array_like
5454 These values are appended to a copy of `arr`. It must be of the
5455 correct shape (the same shape as `arr`, excluding `axis`). If
5456 `axis` is not specified, `values` can be any shape and will be
5457 flattened before use.
5458 axis : int, optional
5459 The axis along which `values` are appended. If `axis` is not
5460 given, both `arr` and `values` are flattened before use.
5462 Returns
5463 -------
5464 append : ndarray
5465 A copy of `arr` with `values` appended to `axis`. Note that
5466 `append` does not occur in-place: a new array is allocated and
5467 filled. If `axis` is None, `out` is a flattened array.
5469 See Also
5470 --------
5471 insert : Insert elements into an array.
5472 delete : Delete elements from an array.
5474 Examples
5475 --------
5476 >>> np.append([1, 2, 3], [[4, 5, 6], [7, 8, 9]])
5477 array([1, 2, 3, ..., 7, 8, 9])
5479 When `axis` is specified, `values` must have the correct shape.
5481 >>> np.append([[1, 2, 3], [4, 5, 6]], [[7, 8, 9]], axis=0)
5482 array([[1, 2, 3],
5483 [4, 5, 6],
5484 [7, 8, 9]])
5485 >>> np.append([[1, 2, 3], [4, 5, 6]], [7, 8, 9], axis=0)
5486 Traceback (most recent call last):
5487 ...
5488 ValueError: all the input arrays must have same number of dimensions, but
5489 the array at index 0 has 2 dimension(s) and the array at index 1 has 1
5490 dimension(s)
5492 """
5493 arr = asanyarray(arr)
5494 if axis is None:
5495 if arr.ndim != 1:
5496 arr = arr.ravel()
5497 values = ravel(values)
5498 axis = arr.ndim-1
5499 return concatenate((arr, values), axis=axis)
5502def _digitize_dispatcher(x, bins, right=None):
5503 return (x, bins)
5506@array_function_dispatch(_digitize_dispatcher)
5507def digitize(x, bins, right=False):
5508 """
5509 Return the indices of the bins to which each value in input array belongs.
5511 ========= ============= ============================
5512 `right` order of bins returned index `i` satisfies
5513 ========= ============= ============================
5514 ``False`` increasing ``bins[i-1] <= x < bins[i]``
5515 ``True`` increasing ``bins[i-1] < x <= bins[i]``
5516 ``False`` decreasing ``bins[i-1] > x >= bins[i]``
5517 ``True`` decreasing ``bins[i-1] >= x > bins[i]``
5518 ========= ============= ============================
5520 If values in `x` are beyond the bounds of `bins`, 0 or ``len(bins)`` is
5521 returned as appropriate.
5523 Parameters
5524 ----------
5525 x : array_like
5526 Input array to be binned. Prior to NumPy 1.10.0, this array had to
5527 be 1-dimensional, but can now have any shape.
5528 bins : array_like
5529 Array of bins. It has to be 1-dimensional and monotonic.
5530 right : bool, optional
5531 Indicating whether the intervals include the right or the left bin
5532 edge. Default behavior is (right==False) indicating that the interval
5533 does not include the right edge. The left bin end is open in this
5534 case, i.e., bins[i-1] <= x < bins[i] is the default behavior for
5535 monotonically increasing bins.
5537 Returns
5538 -------
5539 indices : ndarray of ints
5540 Output array of indices, of same shape as `x`.
5542 Raises
5543 ------
5544 ValueError
5545 If `bins` is not monotonic.
5546 TypeError
5547 If the type of the input is complex.
5549 See Also
5550 --------
5551 bincount, histogram, unique, searchsorted
5553 Notes
5554 -----
5555 If values in `x` are such that they fall outside the bin range,
5556 attempting to index `bins` with the indices that `digitize` returns
5557 will result in an IndexError.
5559 .. versionadded:: 1.10.0
5561 `np.digitize` is implemented in terms of `np.searchsorted`. This means
5562 that a binary search is used to bin the values, which scales much better
5563 for larger number of bins than the previous linear search. It also removes
5564 the requirement for the input array to be 1-dimensional.
5566 For monotonically _increasing_ `bins`, the following are equivalent::
5568 np.digitize(x, bins, right=True)
5569 np.searchsorted(bins, x, side='left')
5571 Note that as the order of the arguments are reversed, the side must be too.
5572 The `searchsorted` call is marginally faster, as it does not do any
5573 monotonicity checks. Perhaps more importantly, it supports all dtypes.
5575 Examples
5576 --------
5577 >>> x = np.array([0.2, 6.4, 3.0, 1.6])
5578 >>> bins = np.array([0.0, 1.0, 2.5, 4.0, 10.0])
5579 >>> inds = np.digitize(x, bins)
5580 >>> inds
5581 array([1, 4, 3, 2])
5582 >>> for n in range(x.size):
5583 ... print(bins[inds[n]-1], "<=", x[n], "<", bins[inds[n]])
5584 ...
5585 0.0 <= 0.2 < 1.0
5586 4.0 <= 6.4 < 10.0
5587 2.5 <= 3.0 < 4.0
5588 1.0 <= 1.6 < 2.5
5590 >>> x = np.array([1.2, 10.0, 12.4, 15.5, 20.])
5591 >>> bins = np.array([0, 5, 10, 15, 20])
5592 >>> np.digitize(x,bins,right=True)
5593 array([1, 2, 3, 4, 4])
5594 >>> np.digitize(x,bins,right=False)
5595 array([1, 3, 3, 4, 5])
5596 """
5597 x = _nx.asarray(x)
5598 bins = _nx.asarray(bins)
5600 # here for compatibility, searchsorted below is happy to take this
5601 if np.issubdtype(x.dtype, _nx.complexfloating):
5602 raise TypeError("x may not be complex")
5604 mono = _monotonicity(bins)
5605 if mono == 0:
5606 raise ValueError("bins must be monotonically increasing or decreasing")
5608 # this is backwards because the arguments below are swapped
5609 side = 'left' if right else 'right'
5610 if mono == -1:
5611 # reverse the bins, and invert the results
5612 return len(bins) - _nx.searchsorted(bins[::-1], x, side=side)
5613 else:
5614 return _nx.searchsorted(bins, x, side=side)