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
7from .._utils import set_module
8import numpy as np
9import numpy.core.numeric as _nx
10from numpy.core import transpose
11from numpy.core.numeric import (
12 ones, zeros_like, arange, concatenate, array, asarray, asanyarray, empty,
13 ndarray, take, dot, where, intp, integer, isscalar, absolute
14 )
15from numpy.core.umath import (
16 pi, add, arctan2, frompyfunc, cos, less_equal, sqrt, sin,
17 mod, exp, not_equal, subtract
18 )
19from numpy.core.fromnumeric import (
20 ravel, nonzero, partition, mean, any, sum
21 )
22from numpy.core.numerictypes import typecodes
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 _place, 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.
164 This means for a 2D array with the default `k` and `axes`, the
165 rotation will be counterclockwise.
167 Parameters
168 ----------
169 m : array_like
170 Array of two or more dimensions.
171 k : integer
172 Number of times the array is rotated by 90 degrees.
173 axes : (2,) array_like
174 The array is rotated in the plane defined by the axes.
175 Axes must be different.
177 .. versionadded:: 1.12.0
179 Returns
180 -------
181 y : ndarray
182 A rotated view of `m`.
184 See Also
185 --------
186 flip : Reverse the order of elements in an array along the given axis.
187 fliplr : Flip an array horizontally.
188 flipud : Flip an array vertically.
190 Notes
191 -----
192 ``rot90(m, k=1, axes=(1,0))`` is the reverse of
193 ``rot90(m, k=1, axes=(0,1))``
195 ``rot90(m, k=1, axes=(1,0))`` is equivalent to
196 ``rot90(m, k=-1, axes=(0,1))``
198 Examples
199 --------
200 >>> m = np.array([[1,2],[3,4]], int)
201 >>> m
202 array([[1, 2],
203 [3, 4]])
204 >>> np.rot90(m)
205 array([[2, 4],
206 [1, 3]])
207 >>> np.rot90(m, 2)
208 array([[4, 3],
209 [2, 1]])
210 >>> m = np.arange(8).reshape((2,2,2))
211 >>> np.rot90(m, 1, (1,2))
212 array([[[1, 3],
213 [0, 2]],
214 [[5, 7],
215 [4, 6]]])
217 """
218 axes = tuple(axes)
219 if len(axes) != 2:
220 raise ValueError("len(axes) must be 2.")
222 m = asanyarray(m)
224 if axes[0] == axes[1] or absolute(axes[0] - axes[1]) == m.ndim:
225 raise ValueError("Axes must be different.")
227 if (axes[0] >= m.ndim or axes[0] < -m.ndim
228 or axes[1] >= m.ndim or axes[1] < -m.ndim):
229 raise ValueError("Axes={} out of range for array of ndim={}."
230 .format(axes, m.ndim))
232 k %= 4
234 if k == 0:
235 return m[:]
236 if k == 2:
237 return flip(flip(m, axes[0]), axes[1])
239 axes_list = arange(0, m.ndim)
240 (axes_list[axes[0]], axes_list[axes[1]]) = (axes_list[axes[1]],
241 axes_list[axes[0]])
243 if k == 1:
244 return transpose(flip(m, axes[1]), axes_list)
245 else:
246 # k == 3
247 return flip(transpose(m, axes_list), axes[1])
250def _flip_dispatcher(m, axis=None):
251 return (m,)
254@array_function_dispatch(_flip_dispatcher)
255def flip(m, axis=None):
256 """
257 Reverse the order of elements in an array along the given axis.
259 The shape of the array is preserved, but the elements are reordered.
261 .. versionadded:: 1.12.0
263 Parameters
264 ----------
265 m : array_like
266 Input array.
267 axis : None or int or tuple of ints, optional
268 Axis or axes along which to flip over. The default,
269 axis=None, will flip over all of the axes of the input array.
270 If axis is negative it counts from the last to the first axis.
272 If axis is a tuple of ints, flipping is performed on all of the axes
273 specified in the tuple.
275 .. versionchanged:: 1.15.0
276 None and tuples of axes are supported
278 Returns
279 -------
280 out : array_like
281 A view of `m` with the entries of axis reversed. Since a view is
282 returned, this operation is done in constant time.
284 See Also
285 --------
286 flipud : Flip an array vertically (axis=0).
287 fliplr : Flip an array horizontally (axis=1).
289 Notes
290 -----
291 flip(m, 0) is equivalent to flipud(m).
293 flip(m, 1) is equivalent to fliplr(m).
295 flip(m, n) corresponds to ``m[...,::-1,...]`` with ``::-1`` at position n.
297 flip(m) corresponds to ``m[::-1,::-1,...,::-1]`` with ``::-1`` at all
298 positions.
300 flip(m, (0, 1)) corresponds to ``m[::-1,::-1,...]`` with ``::-1`` at
301 position 0 and position 1.
303 Examples
304 --------
305 >>> A = np.arange(8).reshape((2,2,2))
306 >>> A
307 array([[[0, 1],
308 [2, 3]],
309 [[4, 5],
310 [6, 7]]])
311 >>> np.flip(A, 0)
312 array([[[4, 5],
313 [6, 7]],
314 [[0, 1],
315 [2, 3]]])
316 >>> np.flip(A, 1)
317 array([[[2, 3],
318 [0, 1]],
319 [[6, 7],
320 [4, 5]]])
321 >>> np.flip(A)
322 array([[[7, 6],
323 [5, 4]],
324 [[3, 2],
325 [1, 0]]])
326 >>> np.flip(A, (0, 2))
327 array([[[5, 4],
328 [7, 6]],
329 [[1, 0],
330 [3, 2]]])
331 >>> A = np.random.randn(3,4,5)
332 >>> np.all(np.flip(A,2) == A[:,:,::-1,...])
333 True
334 """
335 if not hasattr(m, 'ndim'):
336 m = asarray(m)
337 if axis is None:
338 indexer = (np.s_[::-1],) * m.ndim
339 else:
340 axis = _nx.normalize_axis_tuple(axis, m.ndim)
341 indexer = [np.s_[:]] * m.ndim
342 for ax in axis:
343 indexer[ax] = np.s_[::-1]
344 indexer = tuple(indexer)
345 return m[indexer]
348@set_module('numpy')
349def iterable(y):
350 """
351 Check whether or not an object can be iterated over.
353 Parameters
354 ----------
355 y : object
356 Input object.
358 Returns
359 -------
360 b : bool
361 Return ``True`` if the object has an iterator method or is a
362 sequence and ``False`` otherwise.
365 Examples
366 --------
367 >>> np.iterable([1, 2, 3])
368 True
369 >>> np.iterable(2)
370 False
372 Notes
373 -----
374 In most cases, the results of ``np.iterable(obj)`` are consistent with
375 ``isinstance(obj, collections.abc.Iterable)``. One notable exception is
376 the treatment of 0-dimensional arrays::
378 >>> from collections.abc import Iterable
379 >>> a = np.array(1.0) # 0-dimensional numpy array
380 >>> isinstance(a, Iterable)
381 True
382 >>> np.iterable(a)
383 False
385 """
386 try:
387 iter(y)
388 except TypeError:
389 return False
390 return True
393def _average_dispatcher(a, axis=None, weights=None, returned=None, *,
394 keepdims=None):
395 return (a, weights)
398@array_function_dispatch(_average_dispatcher)
399def average(a, axis=None, weights=None, returned=False, *,
400 keepdims=np._NoValue):
401 """
402 Compute the weighted average along the specified axis.
404 Parameters
405 ----------
406 a : array_like
407 Array containing data to be averaged. If `a` is not an array, a
408 conversion is attempted.
409 axis : None or int or tuple of ints, optional
410 Axis or axes along which to average `a`. The default,
411 axis=None, will average over all of the elements of the input array.
412 If axis is negative it counts from the last to the first axis.
414 .. versionadded:: 1.7.0
416 If axis is a tuple of ints, averaging is performed on all of the axes
417 specified in the tuple instead of a single axis or all the axes as
418 before.
419 weights : array_like, optional
420 An array of weights associated with the values in `a`. Each value in
421 `a` contributes to the average according to its associated weight.
422 The weights array can either be 1-D (in which case its length must be
423 the size of `a` along the given axis) or of the same shape as `a`.
424 If `weights=None`, then all data in `a` are assumed to have a
425 weight equal to one. The 1-D calculation is::
427 avg = sum(a * weights) / sum(weights)
429 The only constraint on `weights` is that `sum(weights)` must not be 0.
430 returned : bool, optional
431 Default is `False`. If `True`, the tuple (`average`, `sum_of_weights`)
432 is returned, otherwise only the average is returned.
433 If `weights=None`, `sum_of_weights` is equivalent to the number of
434 elements over which the average is taken.
435 keepdims : bool, optional
436 If this is set to True, the axes which are reduced are left
437 in the result as dimensions with size one. With this option,
438 the result will broadcast correctly against the original `a`.
439 *Note:* `keepdims` will not work with instances of `numpy.matrix`
440 or other classes whose methods do not support `keepdims`.
442 .. versionadded:: 1.23.0
444 Returns
445 -------
446 retval, [sum_of_weights] : array_type or double
447 Return the average along the specified axis. When `returned` is `True`,
448 return a tuple with the average as the first element and the sum
449 of the weights as the second element. `sum_of_weights` is of the
450 same type as `retval`. The result dtype follows a genereal pattern.
451 If `weights` is None, the result dtype will be that of `a` , or ``float64``
452 if `a` is integral. Otherwise, if `weights` is not None and `a` is non-
453 integral, the result type will be the type of lowest precision capable of
454 representing values of both `a` and `weights`. If `a` happens to be
455 integral, the previous rules still applies but the result dtype will
456 at least be ``float64``.
458 Raises
459 ------
460 ZeroDivisionError
461 When all weights along axis are zero. See `numpy.ma.average` for a
462 version robust to this type of error.
463 TypeError
464 When the length of 1D `weights` is not the same as the shape of `a`
465 along axis.
467 See Also
468 --------
469 mean
471 ma.average : average for masked arrays -- useful if your data contains
472 "missing" values
473 numpy.result_type : Returns the type that results from applying the
474 numpy type promotion rules to the arguments.
476 Examples
477 --------
478 >>> data = np.arange(1, 5)
479 >>> data
480 array([1, 2, 3, 4])
481 >>> np.average(data)
482 2.5
483 >>> np.average(np.arange(1, 11), weights=np.arange(10, 0, -1))
484 4.0
486 >>> data = np.arange(6).reshape((3, 2))
487 >>> data
488 array([[0, 1],
489 [2, 3],
490 [4, 5]])
491 >>> np.average(data, axis=1, weights=[1./4, 3./4])
492 array([0.75, 2.75, 4.75])
493 >>> np.average(data, weights=[1./4, 3./4])
494 Traceback (most recent call last):
495 ...
496 TypeError: Axis must be specified when shapes of a and weights differ.
498 >>> a = np.ones(5, dtype=np.float128)
499 >>> w = np.ones(5, dtype=np.complex64)
500 >>> avg = np.average(a, weights=w)
501 >>> print(avg.dtype)
502 complex256
504 With ``keepdims=True``, the following result has shape (3, 1).
506 >>> np.average(data, axis=1, keepdims=True)
507 array([[0.5],
508 [2.5],
509 [4.5]])
510 """
511 a = np.asanyarray(a)
513 if keepdims is np._NoValue:
514 # Don't pass on the keepdims argument if one wasn't given.
515 keepdims_kw = {}
516 else:
517 keepdims_kw = {'keepdims': keepdims}
519 if weights is None:
520 avg = a.mean(axis, **keepdims_kw)
521 avg_as_array = np.asanyarray(avg)
522 scl = avg_as_array.dtype.type(a.size/avg_as_array.size)
523 else:
524 wgt = np.asanyarray(weights)
526 if issubclass(a.dtype.type, (np.integer, np.bool_)):
527 result_dtype = np.result_type(a.dtype, wgt.dtype, 'f8')
528 else:
529 result_dtype = np.result_type(a.dtype, wgt.dtype)
531 # Sanity checks
532 if a.shape != wgt.shape:
533 if axis is None:
534 raise TypeError(
535 "Axis must be specified when shapes of a and weights "
536 "differ.")
537 if wgt.ndim != 1:
538 raise TypeError(
539 "1D weights expected when shapes of a and weights differ.")
540 if wgt.shape[0] != a.shape[axis]:
541 raise ValueError(
542 "Length of weights not compatible with specified axis.")
544 # setup wgt to broadcast along axis
545 wgt = np.broadcast_to(wgt, (a.ndim-1)*(1,) + wgt.shape)
546 wgt = wgt.swapaxes(-1, axis)
548 scl = wgt.sum(axis=axis, dtype=result_dtype, **keepdims_kw)
549 if np.any(scl == 0.0):
550 raise ZeroDivisionError(
551 "Weights sum to zero, can't be normalized")
553 avg = avg_as_array = np.multiply(a, wgt,
554 dtype=result_dtype).sum(axis, **keepdims_kw) / scl
556 if returned:
557 if scl.shape != avg_as_array.shape:
558 scl = np.broadcast_to(scl, avg_as_array.shape).copy()
559 return avg, scl
560 else:
561 return avg
564@set_module('numpy')
565def asarray_chkfinite(a, dtype=None, order=None):
566 """Convert the input to an array, checking for NaNs or Infs.
568 Parameters
569 ----------
570 a : array_like
571 Input data, in any form that can be converted to an array. This
572 includes lists, lists of tuples, tuples, tuples of tuples, tuples
573 of lists and ndarrays. Success requires no NaNs or Infs.
574 dtype : data-type, optional
575 By default, the data-type is inferred from the input data.
576 order : {'C', 'F', 'A', 'K'}, optional
577 Memory layout. 'A' and 'K' depend on the order of input array a.
578 'C' row-major (C-style),
579 'F' column-major (Fortran-style) memory representation.
580 'A' (any) means 'F' if `a` is Fortran contiguous, 'C' otherwise
581 'K' (keep) preserve input order
582 Defaults to 'C'.
584 Returns
585 -------
586 out : ndarray
587 Array interpretation of `a`. No copy is performed if the input
588 is already an ndarray. If `a` is a subclass of ndarray, a base
589 class ndarray is returned.
591 Raises
592 ------
593 ValueError
594 Raises ValueError if `a` contains NaN (Not a Number) or Inf (Infinity).
596 See Also
597 --------
598 asarray : Create and array.
599 asanyarray : Similar function which passes through subclasses.
600 ascontiguousarray : Convert input to a contiguous array.
601 asfarray : Convert input to a floating point ndarray.
602 asfortranarray : Convert input to an ndarray with column-major
603 memory order.
604 fromiter : Create an array from an iterator.
605 fromfunction : Construct an array by executing a function on grid
606 positions.
608 Examples
609 --------
610 Convert a list into an array. If all elements are finite
611 ``asarray_chkfinite`` is identical to ``asarray``.
613 >>> a = [1, 2]
614 >>> np.asarray_chkfinite(a, dtype=float)
615 array([1., 2.])
617 Raises ValueError if array_like contains Nans or Infs.
619 >>> a = [1, 2, np.inf]
620 >>> try:
621 ... np.asarray_chkfinite(a)
622 ... except ValueError:
623 ... print('ValueError')
624 ...
625 ValueError
627 """
628 a = asarray(a, dtype=dtype, order=order)
629 if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
630 raise ValueError(
631 "array must not contain infs or NaNs")
632 return a
635def _piecewise_dispatcher(x, condlist, funclist, *args, **kw):
636 yield x
637 # support the undocumented behavior of allowing scalars
638 if np.iterable(condlist):
639 yield from condlist
642@array_function_dispatch(_piecewise_dispatcher)
643def piecewise(x, condlist, funclist, *args, **kw):
644 """
645 Evaluate a piecewise-defined function.
647 Given a set of conditions and corresponding functions, evaluate each
648 function on the input data wherever its condition is true.
650 Parameters
651 ----------
652 x : ndarray or scalar
653 The input domain.
654 condlist : list of bool arrays or bool scalars
655 Each boolean array corresponds to a function in `funclist`. Wherever
656 `condlist[i]` is True, `funclist[i](x)` is used as the output value.
658 Each boolean array in `condlist` selects a piece of `x`,
659 and should therefore be of the same shape as `x`.
661 The length of `condlist` must correspond to that of `funclist`.
662 If one extra function is given, i.e. if
663 ``len(funclist) == len(condlist) + 1``, then that extra function
664 is the default value, used wherever all conditions are false.
665 funclist : list of callables, f(x,*args,**kw), or scalars
666 Each function is evaluated over `x` wherever its corresponding
667 condition is True. It should take a 1d array as input and give an 1d
668 array or a scalar value as output. If, instead of a callable,
669 a scalar is provided then a constant function (``lambda x: scalar``) is
670 assumed.
671 args : tuple, optional
672 Any further arguments given to `piecewise` are passed to the functions
673 upon execution, i.e., if called ``piecewise(..., ..., 1, 'a')``, then
674 each function is called as ``f(x, 1, 'a')``.
675 kw : dict, optional
676 Keyword arguments used in calling `piecewise` are passed to the
677 functions upon execution, i.e., if called
678 ``piecewise(..., ..., alpha=1)``, then each function is called as
679 ``f(x, alpha=1)``.
681 Returns
682 -------
683 out : ndarray
684 The output is the same shape and type as x and is found by
685 calling the functions in `funclist` on the appropriate portions of `x`,
686 as defined by the boolean arrays in `condlist`. Portions not covered
687 by any condition have a default value of 0.
690 See Also
691 --------
692 choose, select, where
694 Notes
695 -----
696 This is similar to choose or select, except that functions are
697 evaluated on elements of `x` that satisfy the corresponding condition from
698 `condlist`.
700 The result is::
702 |--
703 |funclist[0](x[condlist[0]])
704 out = |funclist[1](x[condlist[1]])
705 |...
706 |funclist[n2](x[condlist[n2]])
707 |--
709 Examples
710 --------
711 Define the sigma function, which is -1 for ``x < 0`` and +1 for ``x >= 0``.
713 >>> x = np.linspace(-2.5, 2.5, 6)
714 >>> np.piecewise(x, [x < 0, x >= 0], [-1, 1])
715 array([-1., -1., -1., 1., 1., 1.])
717 Define the absolute value, which is ``-x`` for ``x <0`` and ``x`` for
718 ``x >= 0``.
720 >>> np.piecewise(x, [x < 0, x >= 0], [lambda x: -x, lambda x: x])
721 array([2.5, 1.5, 0.5, 0.5, 1.5, 2.5])
723 Apply the same function to a scalar value.
725 >>> y = -2
726 >>> np.piecewise(y, [y < 0, y >= 0], [lambda x: -x, lambda x: x])
727 array(2)
729 """
730 x = asanyarray(x)
731 n2 = len(funclist)
733 # undocumented: single condition is promoted to a list of one condition
734 if isscalar(condlist) or (
735 not isinstance(condlist[0], (list, ndarray)) and x.ndim != 0):
736 condlist = [condlist]
738 condlist = asarray(condlist, dtype=bool)
739 n = len(condlist)
741 if n == n2 - 1: # compute the "otherwise" condition.
742 condelse = ~np.any(condlist, axis=0, keepdims=True)
743 condlist = np.concatenate([condlist, condelse], axis=0)
744 n += 1
745 elif n != n2:
746 raise ValueError(
747 "with {} condition(s), either {} or {} functions are expected"
748 .format(n, n, n+1)
749 )
751 y = zeros_like(x)
752 for cond, func in zip(condlist, funclist):
753 if not isinstance(func, collections.abc.Callable):
754 y[cond] = func
755 else:
756 vals = x[cond]
757 if vals.size > 0:
758 y[cond] = func(vals, *args, **kw)
760 return y
763def _select_dispatcher(condlist, choicelist, default=None):
764 yield from condlist
765 yield from choicelist
768@array_function_dispatch(_select_dispatcher)
769def select(condlist, choicelist, default=0):
770 """
771 Return an array drawn from elements in choicelist, depending on conditions.
773 Parameters
774 ----------
775 condlist : list of bool ndarrays
776 The list of conditions which determine from which array in `choicelist`
777 the output elements are taken. When multiple conditions are satisfied,
778 the first one encountered in `condlist` is used.
779 choicelist : list of ndarrays
780 The list of arrays from which the output elements are taken. It has
781 to be of the same length as `condlist`.
782 default : scalar, optional
783 The element inserted in `output` when all conditions evaluate to False.
785 Returns
786 -------
787 output : ndarray
788 The output at position m is the m-th element of the array in
789 `choicelist` where the m-th element of the corresponding array in
790 `condlist` is True.
792 See Also
793 --------
794 where : Return elements from one of two arrays depending on condition.
795 take, choose, compress, diag, diagonal
797 Examples
798 --------
799 >>> x = np.arange(6)
800 >>> condlist = [x<3, x>3]
801 >>> choicelist = [x, x**2]
802 >>> np.select(condlist, choicelist, 42)
803 array([ 0, 1, 2, 42, 16, 25])
805 >>> condlist = [x<=4, x>3]
806 >>> choicelist = [x, x**2]
807 >>> np.select(condlist, choicelist, 55)
808 array([ 0, 1, 2, 3, 4, 25])
810 """
811 # Check the size of condlist and choicelist are the same, or abort.
812 if len(condlist) != len(choicelist):
813 raise ValueError(
814 'list of cases must be same length as list of conditions')
816 # Now that the dtype is known, handle the deprecated select([], []) case
817 if len(condlist) == 0:
818 raise ValueError("select with an empty condition list is not possible")
820 choicelist = [np.asarray(choice) for choice in choicelist]
822 try:
823 intermediate_dtype = np.result_type(*choicelist)
824 except TypeError as e:
825 msg = f'Choicelist elements do not have a common dtype: {e}'
826 raise TypeError(msg) from None
827 default_array = np.asarray(default)
828 choicelist.append(default_array)
830 # need to get the result type before broadcasting for correct scalar
831 # behaviour
832 try:
833 dtype = np.result_type(intermediate_dtype, default_array)
834 except TypeError as e:
835 msg = f'Choicelists and default value do not have a common dtype: {e}'
836 raise TypeError(msg) from None
838 # Convert conditions to arrays and broadcast conditions and choices
839 # as the shape is needed for the result. Doing it separately optimizes
840 # for example when all choices are scalars.
841 condlist = np.broadcast_arrays(*condlist)
842 choicelist = np.broadcast_arrays(*choicelist)
844 # If cond array is not an ndarray in boolean format or scalar bool, abort.
845 for i, cond in enumerate(condlist):
846 if cond.dtype.type is not np.bool_:
847 raise TypeError(
848 'invalid entry {} in condlist: should be boolean ndarray'.format(i))
850 if choicelist[0].ndim == 0:
851 # This may be common, so avoid the call.
852 result_shape = condlist[0].shape
853 else:
854 result_shape = np.broadcast_arrays(condlist[0], choicelist[0])[0].shape
856 result = np.full(result_shape, choicelist[-1], dtype)
858 # Use np.copyto to burn each choicelist array onto result, using the
859 # corresponding condlist as a boolean mask. This is done in reverse
860 # order since the first choice should take precedence.
861 choicelist = choicelist[-2::-1]
862 condlist = condlist[::-1]
863 for choice, cond in zip(choicelist, condlist):
864 np.copyto(result, choice, where=cond)
866 return result
869def _copy_dispatcher(a, order=None, subok=None):
870 return (a,)
873@array_function_dispatch(_copy_dispatcher)
874def copy(a, order='K', subok=False):
875 """
876 Return an array copy of the given object.
878 Parameters
879 ----------
880 a : array_like
881 Input data.
882 order : {'C', 'F', 'A', 'K'}, optional
883 Controls the memory layout of the copy. 'C' means C-order,
884 'F' means F-order, 'A' means 'F' if `a` is Fortran contiguous,
885 'C' otherwise. 'K' means match the layout of `a` as closely
886 as possible. (Note that this function and :meth:`ndarray.copy` are very
887 similar, but have different default values for their order=
888 arguments.)
889 subok : bool, optional
890 If True, then sub-classes will be passed-through, otherwise the
891 returned array will be forced to be a base-class array (defaults to False).
893 .. versionadded:: 1.19.0
895 Returns
896 -------
897 arr : ndarray
898 Array interpretation of `a`.
900 See Also
901 --------
902 ndarray.copy : Preferred method for creating an array copy
904 Notes
905 -----
906 This is equivalent to:
908 >>> np.array(a, copy=True) #doctest: +SKIP
910 Examples
911 --------
912 Create an array x, with a reference y and a copy z:
914 >>> x = np.array([1, 2, 3])
915 >>> y = x
916 >>> z = np.copy(x)
918 Note that, when we modify x, y changes, but not z:
920 >>> x[0] = 10
921 >>> x[0] == y[0]
922 True
923 >>> x[0] == z[0]
924 False
926 Note that, np.copy clears previously set WRITEABLE=False flag.
928 >>> a = np.array([1, 2, 3])
929 >>> a.flags["WRITEABLE"] = False
930 >>> b = np.copy(a)
931 >>> b.flags["WRITEABLE"]
932 True
933 >>> b[0] = 3
934 >>> b
935 array([3, 2, 3])
937 Note that np.copy is a shallow copy and will not copy object
938 elements within arrays. This is mainly important for arrays
939 containing Python objects. The new array will contain the
940 same object which may lead to surprises if that object can
941 be modified (is mutable):
943 >>> a = np.array([1, 'm', [2, 3, 4]], dtype=object)
944 >>> b = np.copy(a)
945 >>> b[2][0] = 10
946 >>> a
947 array([1, 'm', list([10, 3, 4])], dtype=object)
949 To ensure all elements within an ``object`` array are copied,
950 use `copy.deepcopy`:
952 >>> import copy
953 >>> a = np.array([1, 'm', [2, 3, 4]], dtype=object)
954 >>> c = copy.deepcopy(a)
955 >>> c[2][0] = 10
956 >>> c
957 array([1, 'm', list([10, 3, 4])], dtype=object)
958 >>> a
959 array([1, 'm', list([2, 3, 4])], dtype=object)
961 """
962 return array(a, order=order, subok=subok, copy=True)
964# Basic operations
967def _gradient_dispatcher(f, *varargs, axis=None, edge_order=None):
968 yield f
969 yield from varargs
972@array_function_dispatch(_gradient_dispatcher)
973def gradient(f, *varargs, axis=None, edge_order=1):
974 """
975 Return the gradient of an N-dimensional array.
977 The gradient is computed using second order accurate central differences
978 in the interior points and either first or second order accurate one-sides
979 (forward or backwards) differences at the boundaries.
980 The returned gradient hence has the same shape as the input array.
982 Parameters
983 ----------
984 f : array_like
985 An N-dimensional array containing samples of a scalar function.
986 varargs : list of scalar or array, optional
987 Spacing between f values. Default unitary spacing for all dimensions.
988 Spacing can be specified using:
990 1. single scalar to specify a sample distance for all dimensions.
991 2. N scalars to specify a constant sample distance for each dimension.
992 i.e. `dx`, `dy`, `dz`, ...
993 3. N arrays to specify the coordinates of the values along each
994 dimension of F. The length of the array must match the size of
995 the corresponding dimension
996 4. Any combination of N scalars/arrays with the meaning of 2. and 3.
998 If `axis` is given, the number of varargs must equal the number of axes.
999 Default: 1.
1001 edge_order : {1, 2}, optional
1002 Gradient is calculated using N-th order accurate differences
1003 at the boundaries. Default: 1.
1005 .. versionadded:: 1.9.1
1007 axis : None or int or tuple of ints, optional
1008 Gradient is calculated only along the given axis or axes
1009 The default (axis = None) is to calculate the gradient for all the axes
1010 of the input array. axis may be negative, in which case it counts from
1011 the last to the first axis.
1013 .. versionadded:: 1.11.0
1015 Returns
1016 -------
1017 gradient : ndarray or list of ndarray
1018 A list of ndarrays (or a single ndarray if there is only one dimension)
1019 corresponding to the derivatives of f with respect to each dimension.
1020 Each derivative has the same shape as f.
1022 Examples
1023 --------
1024 >>> f = np.array([1, 2, 4, 7, 11, 16], dtype=float)
1025 >>> np.gradient(f)
1026 array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
1027 >>> np.gradient(f, 2)
1028 array([0.5 , 0.75, 1.25, 1.75, 2.25, 2.5 ])
1030 Spacing can be also specified with an array that represents the coordinates
1031 of the values F along the dimensions.
1032 For instance a uniform spacing:
1034 >>> x = np.arange(f.size)
1035 >>> np.gradient(f, x)
1036 array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
1038 Or a non uniform one:
1040 >>> x = np.array([0., 1., 1.5, 3.5, 4., 6.], dtype=float)
1041 >>> np.gradient(f, x)
1042 array([1. , 3. , 3.5, 6.7, 6.9, 2.5])
1044 For two dimensional arrays, the return will be two arrays ordered by
1045 axis. In this example the first array stands for the gradient in
1046 rows and the second one in columns direction:
1048 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float))
1049 [array([[ 2., 2., -1.],
1050 [ 2., 2., -1.]]), array([[1. , 2.5, 4. ],
1051 [1. , 1. , 1. ]])]
1053 In this example the spacing is also specified:
1054 uniform for axis=0 and non uniform for axis=1
1056 >>> dx = 2.
1057 >>> y = [1., 1.5, 3.5]
1058 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), dx, y)
1059 [array([[ 1. , 1. , -0.5],
1060 [ 1. , 1. , -0.5]]), array([[2. , 2. , 2. ],
1061 [2. , 1.7, 0.5]])]
1063 It is possible to specify how boundaries are treated using `edge_order`
1065 >>> x = np.array([0, 1, 2, 3, 4])
1066 >>> f = x**2
1067 >>> np.gradient(f, edge_order=1)
1068 array([1., 2., 4., 6., 7.])
1069 >>> np.gradient(f, edge_order=2)
1070 array([0., 2., 4., 6., 8.])
1072 The `axis` keyword can be used to specify a subset of axes of which the
1073 gradient is calculated
1075 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), axis=0)
1076 array([[ 2., 2., -1.],
1077 [ 2., 2., -1.]])
1079 Notes
1080 -----
1081 Assuming that :math:`f\\in C^{3}` (i.e., :math:`f` has at least 3 continuous
1082 derivatives) and let :math:`h_{*}` be a non-homogeneous stepsize, we
1083 minimize the "consistency error" :math:`\\eta_{i}` between the true gradient
1084 and its estimate from a linear combination of the neighboring grid-points:
1086 .. math::
1088 \\eta_{i} = f_{i}^{\\left(1\\right)} -
1089 \\left[ \\alpha f\\left(x_{i}\\right) +
1090 \\beta f\\left(x_{i} + h_{d}\\right) +
1091 \\gamma f\\left(x_{i}-h_{s}\\right)
1092 \\right]
1094 By substituting :math:`f(x_{i} + h_{d})` and :math:`f(x_{i} - h_{s})`
1095 with their Taylor series expansion, this translates into solving
1096 the following the linear system:
1098 .. math::
1100 \\left\\{
1101 \\begin{array}{r}
1102 \\alpha+\\beta+\\gamma=0 \\\\
1103 \\beta h_{d}-\\gamma h_{s}=1 \\\\
1104 \\beta h_{d}^{2}+\\gamma h_{s}^{2}=0
1105 \\end{array}
1106 \\right.
1108 The resulting approximation of :math:`f_{i}^{(1)}` is the following:
1110 .. math::
1112 \\hat f_{i}^{(1)} =
1113 \\frac{
1114 h_{s}^{2}f\\left(x_{i} + h_{d}\\right)
1115 + \\left(h_{d}^{2} - h_{s}^{2}\\right)f\\left(x_{i}\\right)
1116 - h_{d}^{2}f\\left(x_{i}-h_{s}\\right)}
1117 { h_{s}h_{d}\\left(h_{d} + h_{s}\\right)}
1118 + \\mathcal{O}\\left(\\frac{h_{d}h_{s}^{2}
1119 + h_{s}h_{d}^{2}}{h_{d}
1120 + h_{s}}\\right)
1122 It is worth noting that if :math:`h_{s}=h_{d}`
1123 (i.e., data are evenly spaced)
1124 we find the standard second order approximation:
1126 .. math::
1128 \\hat f_{i}^{(1)}=
1129 \\frac{f\\left(x_{i+1}\\right) - f\\left(x_{i-1}\\right)}{2h}
1130 + \\mathcal{O}\\left(h^{2}\\right)
1132 With a similar procedure the forward/backward approximations used for
1133 boundaries can be derived.
1135 References
1136 ----------
1137 .. [1] Quarteroni A., Sacco R., Saleri F. (2007) Numerical Mathematics
1138 (Texts in Applied Mathematics). New York: Springer.
1139 .. [2] Durran D. R. (1999) Numerical Methods for Wave Equations
1140 in Geophysical Fluid Dynamics. New York: Springer.
1141 .. [3] Fornberg B. (1988) Generation of Finite Difference Formulas on
1142 Arbitrarily Spaced Grids,
1143 Mathematics of Computation 51, no. 184 : 699-706.
1144 `PDF <http://www.ams.org/journals/mcom/1988-51-184/
1145 S0025-5718-1988-0935077-0/S0025-5718-1988-0935077-0.pdf>`_.
1146 """
1147 f = np.asanyarray(f)
1148 N = f.ndim # number of dimensions
1150 if axis is None:
1151 axes = tuple(range(N))
1152 else:
1153 axes = _nx.normalize_axis_tuple(axis, N)
1155 len_axes = len(axes)
1156 n = len(varargs)
1157 if n == 0:
1158 # no spacing argument - use 1 in all axes
1159 dx = [1.0] * len_axes
1160 elif n == 1 and np.ndim(varargs[0]) == 0:
1161 # single scalar for all axes
1162 dx = varargs * len_axes
1163 elif n == len_axes:
1164 # scalar or 1d array for each axis
1165 dx = list(varargs)
1166 for i, distances in enumerate(dx):
1167 distances = np.asanyarray(distances)
1168 if distances.ndim == 0:
1169 continue
1170 elif distances.ndim != 1:
1171 raise ValueError("distances must be either scalars or 1d")
1172 if len(distances) != f.shape[axes[i]]:
1173 raise ValueError("when 1d, distances must match "
1174 "the length of the corresponding dimension")
1175 if np.issubdtype(distances.dtype, np.integer):
1176 # Convert numpy integer types to float64 to avoid modular
1177 # arithmetic in np.diff(distances).
1178 distances = distances.astype(np.float64)
1179 diffx = np.diff(distances)
1180 # if distances are constant reduce to the scalar case
1181 # since it brings a consistent speedup
1182 if (diffx == diffx[0]).all():
1183 diffx = diffx[0]
1184 dx[i] = diffx
1185 else:
1186 raise TypeError("invalid number of arguments")
1188 if edge_order > 2:
1189 raise ValueError("'edge_order' greater than 2 not supported")
1191 # use central differences on interior and one-sided differences on the
1192 # endpoints. This preserves second order-accuracy over the full domain.
1194 outvals = []
1196 # create slice objects --- initially all are [:, :, ..., :]
1197 slice1 = [slice(None)]*N
1198 slice2 = [slice(None)]*N
1199 slice3 = [slice(None)]*N
1200 slice4 = [slice(None)]*N
1202 otype = f.dtype
1203 if otype.type is np.datetime64:
1204 # the timedelta dtype with the same unit information
1205 otype = np.dtype(otype.name.replace('datetime', 'timedelta'))
1206 # view as timedelta to allow addition
1207 f = f.view(otype)
1208 elif otype.type is np.timedelta64:
1209 pass
1210 elif np.issubdtype(otype, np.inexact):
1211 pass
1212 else:
1213 # All other types convert to floating point.
1214 # First check if f is a numpy integer type; if so, convert f to float64
1215 # to avoid modular arithmetic when computing the changes in f.
1216 if np.issubdtype(otype, np.integer):
1217 f = f.astype(np.float64)
1218 otype = np.float64
1220 for axis, ax_dx in zip(axes, dx):
1221 if f.shape[axis] < edge_order + 1:
1222 raise ValueError(
1223 "Shape of array too small to calculate a numerical gradient, "
1224 "at least (edge_order + 1) elements are required.")
1225 # result allocation
1226 out = np.empty_like(f, dtype=otype)
1228 # spacing for the current axis
1229 uniform_spacing = np.ndim(ax_dx) == 0
1231 # Numerical differentiation: 2nd order interior
1232 slice1[axis] = slice(1, -1)
1233 slice2[axis] = slice(None, -2)
1234 slice3[axis] = slice(1, -1)
1235 slice4[axis] = slice(2, None)
1237 if uniform_spacing:
1238 out[tuple(slice1)] = (f[tuple(slice4)] - f[tuple(slice2)]) / (2. * ax_dx)
1239 else:
1240 dx1 = ax_dx[0:-1]
1241 dx2 = ax_dx[1:]
1242 a = -(dx2)/(dx1 * (dx1 + dx2))
1243 b = (dx2 - dx1) / (dx1 * dx2)
1244 c = dx1 / (dx2 * (dx1 + dx2))
1245 # fix the shape for broadcasting
1246 shape = np.ones(N, dtype=int)
1247 shape[axis] = -1
1248 a.shape = b.shape = c.shape = shape
1249 # 1D equivalent -- out[1:-1] = a * f[:-2] + b * f[1:-1] + c * f[2:]
1250 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
1252 # Numerical differentiation: 1st order edges
1253 if edge_order == 1:
1254 slice1[axis] = 0
1255 slice2[axis] = 1
1256 slice3[axis] = 0
1257 dx_0 = ax_dx if uniform_spacing else ax_dx[0]
1258 # 1D equivalent -- out[0] = (f[1] - f[0]) / (x[1] - x[0])
1259 out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0
1261 slice1[axis] = -1
1262 slice2[axis] = -1
1263 slice3[axis] = -2
1264 dx_n = ax_dx if uniform_spacing else ax_dx[-1]
1265 # 1D equivalent -- out[-1] = (f[-1] - f[-2]) / (x[-1] - x[-2])
1266 out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_n
1268 # Numerical differentiation: 2nd order edges
1269 else:
1270 slice1[axis] = 0
1271 slice2[axis] = 0
1272 slice3[axis] = 1
1273 slice4[axis] = 2
1274 if uniform_spacing:
1275 a = -1.5 / ax_dx
1276 b = 2. / ax_dx
1277 c = -0.5 / ax_dx
1278 else:
1279 dx1 = ax_dx[0]
1280 dx2 = ax_dx[1]
1281 a = -(2. * dx1 + dx2)/(dx1 * (dx1 + dx2))
1282 b = (dx1 + dx2) / (dx1 * dx2)
1283 c = - dx1 / (dx2 * (dx1 + dx2))
1284 # 1D equivalent -- out[0] = a * f[0] + b * f[1] + c * f[2]
1285 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
1287 slice1[axis] = -1
1288 slice2[axis] = -3
1289 slice3[axis] = -2
1290 slice4[axis] = -1
1291 if uniform_spacing:
1292 a = 0.5 / ax_dx
1293 b = -2. / ax_dx
1294 c = 1.5 / ax_dx
1295 else:
1296 dx1 = ax_dx[-2]
1297 dx2 = ax_dx[-1]
1298 a = (dx2) / (dx1 * (dx1 + dx2))
1299 b = - (dx2 + dx1) / (dx1 * dx2)
1300 c = (2. * dx2 + dx1) / (dx2 * (dx1 + dx2))
1301 # 1D equivalent -- out[-1] = a * f[-3] + b * f[-2] + c * f[-1]
1302 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
1304 outvals.append(out)
1306 # reset the slice object in this dimension to ":"
1307 slice1[axis] = slice(None)
1308 slice2[axis] = slice(None)
1309 slice3[axis] = slice(None)
1310 slice4[axis] = slice(None)
1312 if len_axes == 1:
1313 return outvals[0]
1314 elif np._using_numpy2_behavior():
1315 return tuple(outvals)
1316 else:
1317 return outvals
1320def _diff_dispatcher(a, n=None, axis=None, prepend=None, append=None):
1321 return (a, prepend, append)
1324@array_function_dispatch(_diff_dispatcher)
1325def diff(a, n=1, axis=-1, prepend=np._NoValue, append=np._NoValue):
1326 """
1327 Calculate the n-th discrete difference along the given axis.
1329 The first difference is given by ``out[i] = a[i+1] - a[i]`` along
1330 the given axis, higher differences are calculated by using `diff`
1331 recursively.
1333 Parameters
1334 ----------
1335 a : array_like
1336 Input array
1337 n : int, optional
1338 The number of times values are differenced. If zero, the input
1339 is returned as-is.
1340 axis : int, optional
1341 The axis along which the difference is taken, default is the
1342 last axis.
1343 prepend, append : array_like, optional
1344 Values to prepend or append to `a` along axis prior to
1345 performing the difference. Scalar values are expanded to
1346 arrays with length 1 in the direction of axis and the shape
1347 of the input array in along all other axes. Otherwise the
1348 dimension and shape must match `a` except along axis.
1350 .. versionadded:: 1.16.0
1352 Returns
1353 -------
1354 diff : ndarray
1355 The n-th differences. The shape of the output is the same as `a`
1356 except along `axis` where the dimension is smaller by `n`. The
1357 type of the output is the same as the type of the difference
1358 between any two elements of `a`. This is the same as the type of
1359 `a` in most cases. A notable exception is `datetime64`, which
1360 results in a `timedelta64` output array.
1362 See Also
1363 --------
1364 gradient, ediff1d, cumsum
1366 Notes
1367 -----
1368 Type is preserved for boolean arrays, so the result will contain
1369 `False` when consecutive elements are the same and `True` when they
1370 differ.
1372 For unsigned integer arrays, the results will also be unsigned. This
1373 should not be surprising, as the result is consistent with
1374 calculating the difference directly:
1376 >>> u8_arr = np.array([1, 0], dtype=np.uint8)
1377 >>> np.diff(u8_arr)
1378 array([255], dtype=uint8)
1379 >>> u8_arr[1,...] - u8_arr[0,...]
1380 255
1382 If this is not desirable, then the array should be cast to a larger
1383 integer type first:
1385 >>> i16_arr = u8_arr.astype(np.int16)
1386 >>> np.diff(i16_arr)
1387 array([-1], dtype=int16)
1389 Examples
1390 --------
1391 >>> x = np.array([1, 2, 4, 7, 0])
1392 >>> np.diff(x)
1393 array([ 1, 2, 3, -7])
1394 >>> np.diff(x, n=2)
1395 array([ 1, 1, -10])
1397 >>> x = np.array([[1, 3, 6, 10], [0, 5, 6, 8]])
1398 >>> np.diff(x)
1399 array([[2, 3, 4],
1400 [5, 1, 2]])
1401 >>> np.diff(x, axis=0)
1402 array([[-1, 2, 0, -2]])
1404 >>> x = np.arange('1066-10-13', '1066-10-16', dtype=np.datetime64)
1405 >>> np.diff(x)
1406 array([1, 1], dtype='timedelta64[D]')
1408 """
1409 if n == 0:
1410 return a
1411 if n < 0:
1412 raise ValueError(
1413 "order must be non-negative but got " + repr(n))
1415 a = asanyarray(a)
1416 nd = a.ndim
1417 if nd == 0:
1418 raise ValueError("diff requires input that is at least one dimensional")
1419 axis = normalize_axis_index(axis, nd)
1421 combined = []
1422 if prepend is not np._NoValue:
1423 prepend = np.asanyarray(prepend)
1424 if prepend.ndim == 0:
1425 shape = list(a.shape)
1426 shape[axis] = 1
1427 prepend = np.broadcast_to(prepend, tuple(shape))
1428 combined.append(prepend)
1430 combined.append(a)
1432 if append is not np._NoValue:
1433 append = np.asanyarray(append)
1434 if append.ndim == 0:
1435 shape = list(a.shape)
1436 shape[axis] = 1
1437 append = np.broadcast_to(append, tuple(shape))
1438 combined.append(append)
1440 if len(combined) > 1:
1441 a = np.concatenate(combined, axis)
1443 slice1 = [slice(None)] * nd
1444 slice2 = [slice(None)] * nd
1445 slice1[axis] = slice(1, None)
1446 slice2[axis] = slice(None, -1)
1447 slice1 = tuple(slice1)
1448 slice2 = tuple(slice2)
1450 op = not_equal if a.dtype == np.bool_ else subtract
1451 for _ in range(n):
1452 a = op(a[slice1], a[slice2])
1454 return a
1457def _interp_dispatcher(x, xp, fp, left=None, right=None, period=None):
1458 return (x, xp, fp)
1461@array_function_dispatch(_interp_dispatcher)
1462def interp(x, xp, fp, left=None, right=None, period=None):
1463 """
1464 One-dimensional linear interpolation for monotonically increasing sample points.
1466 Returns the one-dimensional piecewise linear interpolant to a function
1467 with given discrete data points (`xp`, `fp`), evaluated at `x`.
1469 Parameters
1470 ----------
1471 x : array_like
1472 The x-coordinates at which to evaluate the interpolated values.
1474 xp : 1-D sequence of floats
1475 The x-coordinates of the data points, must be increasing if argument
1476 `period` is not specified. Otherwise, `xp` is internally sorted after
1477 normalizing the periodic boundaries with ``xp = xp % period``.
1479 fp : 1-D sequence of float or complex
1480 The y-coordinates of the data points, same length as `xp`.
1482 left : optional float or complex corresponding to fp
1483 Value to return for `x < xp[0]`, default is `fp[0]`.
1485 right : optional float or complex corresponding to fp
1486 Value to return for `x > xp[-1]`, default is `fp[-1]`.
1488 period : None or float, optional
1489 A period for the x-coordinates. This parameter allows the proper
1490 interpolation of angular x-coordinates. Parameters `left` and `right`
1491 are ignored if `period` is specified.
1493 .. versionadded:: 1.10.0
1495 Returns
1496 -------
1497 y : float or complex (corresponding to fp) or ndarray
1498 The interpolated values, same shape as `x`.
1500 Raises
1501 ------
1502 ValueError
1503 If `xp` and `fp` have different length
1504 If `xp` or `fp` are not 1-D sequences
1505 If `period == 0`
1507 See Also
1508 --------
1509 scipy.interpolate
1511 Warnings
1512 --------
1513 The x-coordinate sequence is expected to be increasing, but this is not
1514 explicitly enforced. However, if the sequence `xp` is non-increasing,
1515 interpolation results are meaningless.
1517 Note that, since NaN is unsortable, `xp` also cannot contain NaNs.
1519 A simple check for `xp` being strictly increasing is::
1521 np.all(np.diff(xp) > 0)
1523 Examples
1524 --------
1525 >>> xp = [1, 2, 3]
1526 >>> fp = [3, 2, 0]
1527 >>> np.interp(2.5, xp, fp)
1528 1.0
1529 >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
1530 array([3. , 3. , 2.5 , 0.56, 0. ])
1531 >>> UNDEF = -99.0
1532 >>> np.interp(3.14, xp, fp, right=UNDEF)
1533 -99.0
1535 Plot an interpolant to the sine function:
1537 >>> x = np.linspace(0, 2*np.pi, 10)
1538 >>> y = np.sin(x)
1539 >>> xvals = np.linspace(0, 2*np.pi, 50)
1540 >>> yinterp = np.interp(xvals, x, y)
1541 >>> import matplotlib.pyplot as plt
1542 >>> plt.plot(x, y, 'o')
1543 [<matplotlib.lines.Line2D object at 0x...>]
1544 >>> plt.plot(xvals, yinterp, '-x')
1545 [<matplotlib.lines.Line2D object at 0x...>]
1546 >>> plt.show()
1548 Interpolation with periodic x-coordinates:
1550 >>> x = [-180, -170, -185, 185, -10, -5, 0, 365]
1551 >>> xp = [190, -190, 350, -350]
1552 >>> fp = [5, 10, 3, 4]
1553 >>> np.interp(x, xp, fp, period=360)
1554 array([7.5 , 5. , 8.75, 6.25, 3. , 3.25, 3.5 , 3.75])
1556 Complex interpolation:
1558 >>> x = [1.5, 4.0]
1559 >>> xp = [2,3,5]
1560 >>> fp = [1.0j, 0, 2+3j]
1561 >>> np.interp(x, xp, fp)
1562 array([0.+1.j , 1.+1.5j])
1564 """
1566 fp = np.asarray(fp)
1568 if np.iscomplexobj(fp):
1569 interp_func = compiled_interp_complex
1570 input_dtype = np.complex128
1571 else:
1572 interp_func = compiled_interp
1573 input_dtype = np.float64
1575 if period is not None:
1576 if period == 0:
1577 raise ValueError("period must be a non-zero value")
1578 period = abs(period)
1579 left = None
1580 right = None
1582 x = np.asarray(x, dtype=np.float64)
1583 xp = np.asarray(xp, dtype=np.float64)
1584 fp = np.asarray(fp, dtype=input_dtype)
1586 if xp.ndim != 1 or fp.ndim != 1:
1587 raise ValueError("Data points must be 1-D sequences")
1588 if xp.shape[0] != fp.shape[0]:
1589 raise ValueError("fp and xp are not of the same length")
1590 # normalizing periodic boundaries
1591 x = x % period
1592 xp = xp % period
1593 asort_xp = np.argsort(xp)
1594 xp = xp[asort_xp]
1595 fp = fp[asort_xp]
1596 xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period))
1597 fp = np.concatenate((fp[-1:], fp, fp[0:1]))
1599 return interp_func(x, xp, fp, left, right)
1602def _angle_dispatcher(z, deg=None):
1603 return (z,)
1606@array_function_dispatch(_angle_dispatcher)
1607def angle(z, deg=False):
1608 """
1609 Return the angle of the complex argument.
1611 Parameters
1612 ----------
1613 z : array_like
1614 A complex number or sequence of complex numbers.
1615 deg : bool, optional
1616 Return angle in degrees if True, radians if False (default).
1618 Returns
1619 -------
1620 angle : ndarray or scalar
1621 The counterclockwise angle from the positive real axis on the complex
1622 plane in the range ``(-pi, pi]``, with dtype as numpy.float64.
1624 .. versionchanged:: 1.16.0
1625 This function works on subclasses of ndarray like `ma.array`.
1627 See Also
1628 --------
1629 arctan2
1630 absolute
1632 Notes
1633 -----
1634 Although the angle of the complex number 0 is undefined, ``numpy.angle(0)``
1635 returns the value 0.
1637 Examples
1638 --------
1639 >>> np.angle([1.0, 1.0j, 1+1j]) # in radians
1640 array([ 0. , 1.57079633, 0.78539816]) # may vary
1641 >>> np.angle(1+1j, deg=True) # in degrees
1642 45.0
1644 """
1645 z = asanyarray(z)
1646 if issubclass(z.dtype.type, _nx.complexfloating):
1647 zimag = z.imag
1648 zreal = z.real
1649 else:
1650 zimag = 0
1651 zreal = z
1653 a = arctan2(zimag, zreal)
1654 if deg:
1655 a *= 180/pi
1656 return a
1659def _unwrap_dispatcher(p, discont=None, axis=None, *, period=None):
1660 return (p,)
1663@array_function_dispatch(_unwrap_dispatcher)
1664def unwrap(p, discont=None, axis=-1, *, period=2*pi):
1665 r"""
1666 Unwrap by taking the complement of large deltas with respect to the period.
1668 This unwraps a signal `p` by changing elements which have an absolute
1669 difference from their predecessor of more than ``max(discont, period/2)``
1670 to their `period`-complementary values.
1672 For the default case where `period` is :math:`2\pi` and `discont` is
1673 :math:`\pi`, this unwraps a radian phase `p` such that adjacent differences
1674 are never greater than :math:`\pi` by adding :math:`2k\pi` for some
1675 integer :math:`k`.
1677 Parameters
1678 ----------
1679 p : array_like
1680 Input array.
1681 discont : float, optional
1682 Maximum discontinuity between values, default is ``period/2``.
1683 Values below ``period/2`` are treated as if they were ``period/2``.
1684 To have an effect different from the default, `discont` should be
1685 larger than ``period/2``.
1686 axis : int, optional
1687 Axis along which unwrap will operate, default is the last axis.
1688 period : float, optional
1689 Size of the range over which the input wraps. By default, it is
1690 ``2 pi``.
1692 .. versionadded:: 1.21.0
1694 Returns
1695 -------
1696 out : ndarray
1697 Output array.
1699 See Also
1700 --------
1701 rad2deg, deg2rad
1703 Notes
1704 -----
1705 If the discontinuity in `p` is smaller than ``period/2``,
1706 but larger than `discont`, no unwrapping is done because taking
1707 the complement would only make the discontinuity larger.
1709 Examples
1710 --------
1711 >>> phase = np.linspace(0, np.pi, num=5)
1712 >>> phase[3:] += np.pi
1713 >>> phase
1714 array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531]) # may vary
1715 >>> np.unwrap(phase)
1716 array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ]) # may vary
1717 >>> np.unwrap([0, 1, 2, -1, 0], period=4)
1718 array([0, 1, 2, 3, 4])
1719 >>> np.unwrap([ 1, 2, 3, 4, 5, 6, 1, 2, 3], period=6)
1720 array([1, 2, 3, 4, 5, 6, 7, 8, 9])
1721 >>> np.unwrap([2, 3, 4, 5, 2, 3, 4, 5], period=4)
1722 array([2, 3, 4, 5, 6, 7, 8, 9])
1723 >>> phase_deg = np.mod(np.linspace(0 ,720, 19), 360) - 180
1724 >>> np.unwrap(phase_deg, period=360)
1725 array([-180., -140., -100., -60., -20., 20., 60., 100., 140.,
1726 180., 220., 260., 300., 340., 380., 420., 460., 500.,
1727 540.])
1728 """
1729 p = asarray(p)
1730 nd = p.ndim
1731 dd = diff(p, axis=axis)
1732 if discont is None:
1733 discont = period/2
1734 slice1 = [slice(None, None)]*nd # full slices
1735 slice1[axis] = slice(1, None)
1736 slice1 = tuple(slice1)
1737 dtype = np.result_type(dd, period)
1738 if _nx.issubdtype(dtype, _nx.integer):
1739 interval_high, rem = divmod(period, 2)
1740 boundary_ambiguous = rem == 0
1741 else:
1742 interval_high = period / 2
1743 boundary_ambiguous = True
1744 interval_low = -interval_high
1745 ddmod = mod(dd - interval_low, period) + interval_low
1746 if boundary_ambiguous:
1747 # for `mask = (abs(dd) == period/2)`, the above line made
1748 # `ddmod[mask] == -period/2`. correct these such that
1749 # `ddmod[mask] == sign(dd[mask])*period/2`.
1750 _nx.copyto(ddmod, interval_high,
1751 where=(ddmod == interval_low) & (dd > 0))
1752 ph_correct = ddmod - dd
1753 _nx.copyto(ph_correct, 0, where=abs(dd) < discont)
1754 up = array(p, copy=True, dtype=dtype)
1755 up[slice1] = p[slice1] + ph_correct.cumsum(axis)
1756 return up
1759def _sort_complex(a):
1760 return (a,)
1763@array_function_dispatch(_sort_complex)
1764def sort_complex(a):
1765 """
1766 Sort a complex array using the real part first, then the imaginary part.
1768 Parameters
1769 ----------
1770 a : array_like
1771 Input array
1773 Returns
1774 -------
1775 out : complex ndarray
1776 Always returns a sorted complex array.
1778 Examples
1779 --------
1780 >>> np.sort_complex([5, 3, 6, 2, 1])
1781 array([1.+0.j, 2.+0.j, 3.+0.j, 5.+0.j, 6.+0.j])
1783 >>> np.sort_complex([1 + 2j, 2 - 1j, 3 - 2j, 3 - 3j, 3 + 5j])
1784 array([1.+2.j, 2.-1.j, 3.-3.j, 3.-2.j, 3.+5.j])
1786 """
1787 b = array(a, copy=True)
1788 b.sort()
1789 if not issubclass(b.dtype.type, _nx.complexfloating):
1790 if b.dtype.char in 'bhBH':
1791 return b.astype('F')
1792 elif b.dtype.char == 'g':
1793 return b.astype('G')
1794 else:
1795 return b.astype('D')
1796 else:
1797 return b
1800def _trim_zeros(filt, trim=None):
1801 return (filt,)
1804@array_function_dispatch(_trim_zeros)
1805def trim_zeros(filt, trim='fb'):
1806 """
1807 Trim the leading and/or trailing zeros from a 1-D array or sequence.
1809 Parameters
1810 ----------
1811 filt : 1-D array or sequence
1812 Input array.
1813 trim : str, optional
1814 A string with 'f' representing trim from front and 'b' to trim from
1815 back. Default is 'fb', trim zeros from both front and back of the
1816 array.
1818 Returns
1819 -------
1820 trimmed : 1-D array or sequence
1821 The result of trimming the input. The input data type is preserved.
1823 Examples
1824 --------
1825 >>> a = np.array((0, 0, 0, 1, 2, 3, 0, 2, 1, 0))
1826 >>> np.trim_zeros(a)
1827 array([1, 2, 3, 0, 2, 1])
1829 >>> np.trim_zeros(a, 'b')
1830 array([0, 0, 0, ..., 0, 2, 1])
1832 The input data type is preserved, list/tuple in means list/tuple out.
1834 >>> np.trim_zeros([0, 1, 2, 0])
1835 [1, 2]
1837 """
1839 first = 0
1840 trim = trim.upper()
1841 if 'F' in trim:
1842 for i in filt:
1843 if i != 0.:
1844 break
1845 else:
1846 first = first + 1
1847 last = len(filt)
1848 if 'B' in trim:
1849 for i in filt[::-1]:
1850 if i != 0.:
1851 break
1852 else:
1853 last = last - 1
1854 return filt[first:last]
1857def _extract_dispatcher(condition, arr):
1858 return (condition, arr)
1861@array_function_dispatch(_extract_dispatcher)
1862def extract(condition, arr):
1863 """
1864 Return the elements of an array that satisfy some condition.
1866 This is equivalent to ``np.compress(ravel(condition), ravel(arr))``. If
1867 `condition` is boolean ``np.extract`` is equivalent to ``arr[condition]``.
1869 Note that `place` does the exact opposite of `extract`.
1871 Parameters
1872 ----------
1873 condition : array_like
1874 An array whose nonzero or True entries indicate the elements of `arr`
1875 to extract.
1876 arr : array_like
1877 Input array of the same size as `condition`.
1879 Returns
1880 -------
1881 extract : ndarray
1882 Rank 1 array of values from `arr` where `condition` is True.
1884 See Also
1885 --------
1886 take, put, copyto, compress, place
1888 Examples
1889 --------
1890 >>> arr = np.arange(12).reshape((3, 4))
1891 >>> arr
1892 array([[ 0, 1, 2, 3],
1893 [ 4, 5, 6, 7],
1894 [ 8, 9, 10, 11]])
1895 >>> condition = np.mod(arr, 3)==0
1896 >>> condition
1897 array([[ True, False, False, True],
1898 [False, False, True, False],
1899 [False, True, False, False]])
1900 >>> np.extract(condition, arr)
1901 array([0, 3, 6, 9])
1904 If `condition` is boolean:
1906 >>> arr[condition]
1907 array([0, 3, 6, 9])
1909 """
1910 return _nx.take(ravel(arr), nonzero(ravel(condition))[0])
1913def _place_dispatcher(arr, mask, vals):
1914 return (arr, mask, vals)
1917@array_function_dispatch(_place_dispatcher)
1918def place(arr, mask, vals):
1919 """
1920 Change elements of an array based on conditional and input values.
1922 Similar to ``np.copyto(arr, vals, where=mask)``, the difference is that
1923 `place` uses the first N elements of `vals`, where N is the number of
1924 True values in `mask`, while `copyto` uses the elements where `mask`
1925 is True.
1927 Note that `extract` does the exact opposite of `place`.
1929 Parameters
1930 ----------
1931 arr : ndarray
1932 Array to put data into.
1933 mask : array_like
1934 Boolean mask array. Must have the same size as `a`.
1935 vals : 1-D sequence
1936 Values to put into `a`. Only the first N elements are used, where
1937 N is the number of True values in `mask`. If `vals` is smaller
1938 than N, it will be repeated, and if elements of `a` are to be masked,
1939 this sequence must be non-empty.
1941 See Also
1942 --------
1943 copyto, put, take, extract
1945 Examples
1946 --------
1947 >>> arr = np.arange(6).reshape(2, 3)
1948 >>> np.place(arr, arr>2, [44, 55])
1949 >>> arr
1950 array([[ 0, 1, 2],
1951 [44, 55, 44]])
1953 """
1954 return _place(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=np._NoValue, otypes=None, doc=None, excluded=None,
2121 cache=False, signature=None)
2123 Returns an object that acts like pyfunc, but takes arrays as input.
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, optional
2138 A python function or method.
2139 Can be omitted to produce a decorator with keyword arguments.
2140 otypes : str or list of dtypes, optional
2141 The output data type. It must be specified as either a string of
2142 typecode characters or a list of data type specifiers. There should
2143 be one data type specifier for each output.
2144 doc : str, optional
2145 The docstring for the function. If None, the docstring will be the
2146 ``pyfunc.__doc__``.
2147 excluded : set, optional
2148 Set of strings or integers representing the positional or keyword
2149 arguments for which the function will not be vectorized. These will be
2150 passed directly to `pyfunc` unmodified.
2152 .. versionadded:: 1.7.0
2154 cache : bool, optional
2155 If `True`, then cache the first function call that determines the number
2156 of outputs if `otypes` is not provided.
2158 .. versionadded:: 1.7.0
2160 signature : string, optional
2161 Generalized universal function signature, e.g., ``(m,n),(n)->(m)`` for
2162 vectorized matrix-vector multiplication. If provided, ``pyfunc`` will
2163 be called with (and expected to return) arrays with shapes given by the
2164 size of corresponding core dimensions. By default, ``pyfunc`` is
2165 assumed to take scalars as input and output.
2167 .. versionadded:: 1.12.0
2169 Returns
2170 -------
2171 out : callable
2172 A vectorized function if ``pyfunc`` was provided,
2173 a decorator otherwise.
2175 See Also
2176 --------
2177 frompyfunc : Takes an arbitrary Python function and returns a ufunc
2179 Notes
2180 -----
2181 The `vectorize` function is provided primarily for convenience, not for
2182 performance. The implementation is essentially a for loop.
2184 If `otypes` is not specified, then a call to the function with the
2185 first argument will be used to determine the number of outputs. The
2186 results of this call will be cached if `cache` is `True` to prevent
2187 calling the function twice. However, to implement the cache, the
2188 original function must be wrapped which will slow down subsequent
2189 calls, so only do this if your function is expensive.
2191 The new keyword argument interface and `excluded` argument support
2192 further degrades performance.
2194 References
2195 ----------
2196 .. [1] :doc:`/reference/c-api/generalized-ufuncs`
2198 Examples
2199 --------
2200 >>> def myfunc(a, b):
2201 ... "Return a-b if a>b, otherwise return a+b"
2202 ... if a > b:
2203 ... return a - b
2204 ... else:
2205 ... return a + b
2207 >>> vfunc = np.vectorize(myfunc)
2208 >>> vfunc([1, 2, 3, 4], 2)
2209 array([3, 4, 1, 2])
2211 The docstring is taken from the input function to `vectorize` unless it
2212 is specified:
2214 >>> vfunc.__doc__
2215 'Return a-b if a>b, otherwise return a+b'
2216 >>> vfunc = np.vectorize(myfunc, doc='Vectorized `myfunc`')
2217 >>> vfunc.__doc__
2218 'Vectorized `myfunc`'
2220 The output type is determined by evaluating the first element of the input,
2221 unless it is specified:
2223 >>> out = vfunc([1, 2, 3, 4], 2)
2224 >>> type(out[0])
2225 <class 'numpy.int64'>
2226 >>> vfunc = np.vectorize(myfunc, otypes=[float])
2227 >>> out = vfunc([1, 2, 3, 4], 2)
2228 >>> type(out[0])
2229 <class 'numpy.float64'>
2231 The `excluded` argument can be used to prevent vectorizing over certain
2232 arguments. This can be useful for array-like arguments of a fixed length
2233 such as the coefficients for a polynomial as in `polyval`:
2235 >>> def mypolyval(p, x):
2236 ... _p = list(p)
2237 ... res = _p.pop(0)
2238 ... while _p:
2239 ... res = res*x + _p.pop(0)
2240 ... return res
2241 >>> vpolyval = np.vectorize(mypolyval, excluded=['p'])
2242 >>> vpolyval(p=[1, 2, 3], x=[0, 1])
2243 array([3, 6])
2245 Positional arguments may also be excluded by specifying their position:
2247 >>> vpolyval.excluded.add(0)
2248 >>> vpolyval([1, 2, 3], x=[0, 1])
2249 array([3, 6])
2251 The `signature` argument allows for vectorizing functions that act on
2252 non-scalar arrays of fixed length. For example, you can use it for a
2253 vectorized calculation of Pearson correlation coefficient and its p-value:
2255 >>> import scipy.stats
2256 >>> pearsonr = np.vectorize(scipy.stats.pearsonr,
2257 ... signature='(n),(n)->(),()')
2258 >>> pearsonr([[0, 1, 2, 3]], [[1, 2, 3, 4], [4, 3, 2, 1]])
2259 (array([ 1., -1.]), array([ 0., 0.]))
2261 Or for a vectorized convolution:
2263 >>> convolve = np.vectorize(np.convolve, signature='(n),(m)->(k)')
2264 >>> convolve(np.eye(4), [1, 2, 1])
2265 array([[1., 2., 1., 0., 0., 0.],
2266 [0., 1., 2., 1., 0., 0.],
2267 [0., 0., 1., 2., 1., 0.],
2268 [0., 0., 0., 1., 2., 1.]])
2270 Decorator syntax is supported. The decorator can be called as
2271 a function to provide keyword arguments.
2272 >>>@np.vectorize
2273 ...def identity(x):
2274 ... return x
2275 ...
2276 >>>identity([0, 1, 2])
2277 array([0, 1, 2])
2278 >>>@np.vectorize(otypes=[float])
2279 ...def as_float(x):
2280 ... return x
2281 ...
2282 >>>as_float([0, 1, 2])
2283 array([0., 1., 2.])
2284 """
2285 def __init__(self, pyfunc=np._NoValue, otypes=None, doc=None,
2286 excluded=None, cache=False, signature=None):
2288 if (pyfunc != np._NoValue) and (not callable(pyfunc)):
2289 #Splitting the error message to keep
2290 #the length below 79 characters.
2291 part1 = "When used as a decorator, "
2292 part2 = "only accepts keyword arguments."
2293 raise TypeError(part1 + part2)
2295 self.pyfunc = pyfunc
2296 self.cache = cache
2297 self.signature = signature
2298 if pyfunc != np._NoValue and hasattr(pyfunc, '__name__'):
2299 self.__name__ = pyfunc.__name__
2301 self._ufunc = {} # Caching to improve default performance
2302 self._doc = None
2303 self.__doc__ = doc
2304 if doc is None and hasattr(pyfunc, '__doc__'):
2305 self.__doc__ = pyfunc.__doc__
2306 else:
2307 self._doc = doc
2309 if isinstance(otypes, str):
2310 for char in otypes:
2311 if char not in typecodes['All']:
2312 raise ValueError("Invalid otype specified: %s" % (char,))
2313 elif iterable(otypes):
2314 otypes = ''.join([_nx.dtype(x).char for x in otypes])
2315 elif otypes is not None:
2316 raise ValueError("Invalid otype specification")
2317 self.otypes = otypes
2319 # Excluded variable support
2320 if excluded is None:
2321 excluded = set()
2322 self.excluded = set(excluded)
2324 if signature is not None:
2325 self._in_and_out_core_dims = _parse_gufunc_signature(signature)
2326 else:
2327 self._in_and_out_core_dims = None
2329 def _init_stage_2(self, pyfunc, *args, **kwargs):
2330 self.__name__ = pyfunc.__name__
2331 self.pyfunc = pyfunc
2332 if self._doc is None:
2333 self.__doc__ = pyfunc.__doc__
2334 else:
2335 self.__doc__ = self._doc
2337 def _call_as_normal(self, *args, **kwargs):
2338 """
2339 Return arrays with the results of `pyfunc` broadcast (vectorized) over
2340 `args` and `kwargs` not in `excluded`.
2341 """
2342 excluded = self.excluded
2343 if not kwargs and not excluded:
2344 func = self.pyfunc
2345 vargs = args
2346 else:
2347 # The wrapper accepts only positional arguments: we use `names` and
2348 # `inds` to mutate `the_args` and `kwargs` to pass to the original
2349 # function.
2350 nargs = len(args)
2352 names = [_n for _n in kwargs if _n not in excluded]
2353 inds = [_i for _i in range(nargs) if _i not in excluded]
2354 the_args = list(args)
2356 def func(*vargs):
2357 for _n, _i in enumerate(inds):
2358 the_args[_i] = vargs[_n]
2359 kwargs.update(zip(names, vargs[len(inds):]))
2360 return self.pyfunc(*the_args, **kwargs)
2362 vargs = [args[_i] for _i in inds]
2363 vargs.extend([kwargs[_n] for _n in names])
2365 return self._vectorize_call(func=func, args=vargs)
2367 def __call__(self, *args, **kwargs):
2368 if self.pyfunc is np._NoValue:
2369 self._init_stage_2(*args, **kwargs)
2370 return self
2372 return self._call_as_normal(*args, **kwargs)
2374 def _get_ufunc_and_otypes(self, func, args):
2375 """Return (ufunc, otypes)."""
2376 # frompyfunc will fail if args is empty
2377 if not args:
2378 raise ValueError('args can not be empty')
2380 if self.otypes is not None:
2381 otypes = self.otypes
2383 # self._ufunc is a dictionary whose keys are the number of
2384 # arguments (i.e. len(args)) and whose values are ufuncs created
2385 # by frompyfunc. len(args) can be different for different calls if
2386 # self.pyfunc has parameters with default values. We only use the
2387 # cache when func is self.pyfunc, which occurs when the call uses
2388 # only positional arguments and no arguments are excluded.
2390 nin = len(args)
2391 nout = len(self.otypes)
2392 if func is not self.pyfunc or nin not in self._ufunc:
2393 ufunc = frompyfunc(func, nin, nout)
2394 else:
2395 ufunc = None # We'll get it from self._ufunc
2396 if func is self.pyfunc:
2397 ufunc = self._ufunc.setdefault(nin, ufunc)
2398 else:
2399 # Get number of outputs and output types by calling the function on
2400 # the first entries of args. We also cache the result to prevent
2401 # the subsequent call when the ufunc is evaluated.
2402 # Assumes that ufunc first evaluates the 0th elements in the input
2403 # arrays (the input values are not checked to ensure this)
2404 args = [asarray(arg) for arg in args]
2405 if builtins.any(arg.size == 0 for arg in args):
2406 raise ValueError('cannot call `vectorize` on size 0 inputs '
2407 'unless `otypes` is set')
2409 inputs = [arg.flat[0] for arg in args]
2410 outputs = func(*inputs)
2412 # Performance note: profiling indicates that -- for simple
2413 # functions at least -- this wrapping can almost double the
2414 # execution time.
2415 # Hence we make it optional.
2416 if self.cache:
2417 _cache = [outputs]
2419 def _func(*vargs):
2420 if _cache:
2421 return _cache.pop()
2422 else:
2423 return func(*vargs)
2424 else:
2425 _func = func
2427 if isinstance(outputs, tuple):
2428 nout = len(outputs)
2429 else:
2430 nout = 1
2431 outputs = (outputs,)
2433 otypes = ''.join([asarray(outputs[_k]).dtype.char
2434 for _k in range(nout)])
2436 # Performance note: profiling indicates that creating the ufunc is
2437 # not a significant cost compared with wrapping so it seems not
2438 # worth trying to cache this.
2439 ufunc = frompyfunc(_func, len(args), nout)
2441 return ufunc, otypes
2443 def _vectorize_call(self, func, args):
2444 """Vectorized call to `func` over positional `args`."""
2445 if self.signature is not None:
2446 res = self._vectorize_call_with_signature(func, args)
2447 elif not args:
2448 res = func()
2449 else:
2450 ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args)
2452 # Convert args to object arrays first
2453 inputs = [asanyarray(a, dtype=object) for a in args]
2455 outputs = ufunc(*inputs)
2457 if ufunc.nout == 1:
2458 res = asanyarray(outputs, dtype=otypes[0])
2459 else:
2460 res = tuple([asanyarray(x, dtype=t)
2461 for x, t in zip(outputs, otypes)])
2462 return res
2464 def _vectorize_call_with_signature(self, func, args):
2465 """Vectorized call over positional arguments with a signature."""
2466 input_core_dims, output_core_dims = self._in_and_out_core_dims
2468 if len(args) != len(input_core_dims):
2469 raise TypeError('wrong number of positional arguments: '
2470 'expected %r, got %r'
2471 % (len(input_core_dims), len(args)))
2472 args = tuple(asanyarray(arg) for arg in args)
2474 broadcast_shape, dim_sizes = _parse_input_dimensions(
2475 args, input_core_dims)
2476 input_shapes = _calculate_shapes(broadcast_shape, dim_sizes,
2477 input_core_dims)
2478 args = [np.broadcast_to(arg, shape, subok=True)
2479 for arg, shape in zip(args, input_shapes)]
2481 outputs = None
2482 otypes = self.otypes
2483 nout = len(output_core_dims)
2485 for index in np.ndindex(*broadcast_shape):
2486 results = func(*(arg[index] for arg in args))
2488 n_results = len(results) if isinstance(results, tuple) else 1
2490 if nout != n_results:
2491 raise ValueError(
2492 'wrong number of outputs from pyfunc: expected %r, got %r'
2493 % (nout, n_results))
2495 if nout == 1:
2496 results = (results,)
2498 if outputs is None:
2499 for result, core_dims in zip(results, output_core_dims):
2500 _update_dim_sizes(dim_sizes, result, core_dims)
2502 outputs = _create_arrays(broadcast_shape, dim_sizes,
2503 output_core_dims, otypes, results)
2505 for output, result in zip(outputs, results):
2506 output[index] = result
2508 if outputs is None:
2509 # did not call the function even once
2510 if otypes is None:
2511 raise ValueError('cannot call `vectorize` on size 0 inputs '
2512 'unless `otypes` is set')
2513 if builtins.any(dim not in dim_sizes
2514 for dims in output_core_dims
2515 for dim in dims):
2516 raise ValueError('cannot call `vectorize` with a signature '
2517 'including new output dimensions on size 0 '
2518 'inputs')
2519 outputs = _create_arrays(broadcast_shape, dim_sizes,
2520 output_core_dims, otypes)
2522 return outputs[0] if nout == 1 else outputs
2525def _cov_dispatcher(m, y=None, rowvar=None, bias=None, ddof=None,
2526 fweights=None, aweights=None, *, dtype=None):
2527 return (m, y, fweights, aweights)
2530@array_function_dispatch(_cov_dispatcher)
2531def cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None,
2532 aweights=None, *, dtype=None):
2533 """
2534 Estimate a covariance matrix, given data and weights.
2536 Covariance indicates the level to which two variables vary together.
2537 If we examine N-dimensional samples, :math:`X = [x_1, x_2, ... x_N]^T`,
2538 then the covariance matrix element :math:`C_{ij}` is the covariance of
2539 :math:`x_i` and :math:`x_j`. The element :math:`C_{ii}` is the variance
2540 of :math:`x_i`.
2542 See the notes for an outline of the algorithm.
2544 Parameters
2545 ----------
2546 m : array_like
2547 A 1-D or 2-D array containing multiple variables and observations.
2548 Each row of `m` represents a variable, and each column a single
2549 observation of all those variables. Also see `rowvar` below.
2550 y : array_like, optional
2551 An additional set of variables and observations. `y` has the same form
2552 as that of `m`.
2553 rowvar : bool, optional
2554 If `rowvar` is True (default), then each row represents a
2555 variable, with observations in the columns. Otherwise, the relationship
2556 is transposed: each column represents a variable, while the rows
2557 contain observations.
2558 bias : bool, optional
2559 Default normalization (False) is by ``(N - 1)``, where ``N`` is the
2560 number of observations given (unbiased estimate). If `bias` is True,
2561 then normalization is by ``N``. These values can be overridden by using
2562 the keyword ``ddof`` in numpy versions >= 1.5.
2563 ddof : int, optional
2564 If not ``None`` the default value implied by `bias` is overridden.
2565 Note that ``ddof=1`` will return the unbiased estimate, even if both
2566 `fweights` and `aweights` are specified, and ``ddof=0`` will return
2567 the simple average. See the notes for the details. The default value
2568 is ``None``.
2570 .. versionadded:: 1.5
2571 fweights : array_like, int, optional
2572 1-D array of integer frequency weights; the number of times each
2573 observation vector should be repeated.
2575 .. versionadded:: 1.10
2576 aweights : array_like, optional
2577 1-D array of observation vector weights. These relative weights are
2578 typically large for observations considered "important" and smaller for
2579 observations considered less "important". If ``ddof=0`` the array of
2580 weights can be used to assign probabilities to observation vectors.
2582 .. versionadded:: 1.10
2583 dtype : data-type, optional
2584 Data-type of the result. By default, the return data-type will have
2585 at least `numpy.float64` precision.
2587 .. versionadded:: 1.20
2589 Returns
2590 -------
2591 out : ndarray
2592 The covariance matrix of the variables.
2594 See Also
2595 --------
2596 corrcoef : Normalized covariance matrix
2598 Notes
2599 -----
2600 Assume that the observations are in the columns of the observation
2601 array `m` and let ``f = fweights`` and ``a = aweights`` for brevity. The
2602 steps to compute the weighted covariance are as follows::
2604 >>> m = np.arange(10, dtype=np.float64)
2605 >>> f = np.arange(10) * 2
2606 >>> a = np.arange(10) ** 2.
2607 >>> ddof = 1
2608 >>> w = f * a
2609 >>> v1 = np.sum(w)
2610 >>> v2 = np.sum(w * a)
2611 >>> m -= np.sum(m * w, axis=None, keepdims=True) / v1
2612 >>> cov = np.dot(m * w, m.T) * v1 / (v1**2 - ddof * v2)
2614 Note that when ``a == 1``, the normalization factor
2615 ``v1 / (v1**2 - ddof * v2)`` goes over to ``1 / (np.sum(f) - ddof)``
2616 as it should.
2618 Examples
2619 --------
2620 Consider two variables, :math:`x_0` and :math:`x_1`, which
2621 correlate perfectly, but in opposite directions:
2623 >>> x = np.array([[0, 2], [1, 1], [2, 0]]).T
2624 >>> x
2625 array([[0, 1, 2],
2626 [2, 1, 0]])
2628 Note how :math:`x_0` increases while :math:`x_1` decreases. The covariance
2629 matrix shows this clearly:
2631 >>> np.cov(x)
2632 array([[ 1., -1.],
2633 [-1., 1.]])
2635 Note that element :math:`C_{0,1}`, which shows the correlation between
2636 :math:`x_0` and :math:`x_1`, is negative.
2638 Further, note how `x` and `y` are combined:
2640 >>> x = [-2.1, -1, 4.3]
2641 >>> y = [3, 1.1, 0.12]
2642 >>> X = np.stack((x, y), axis=0)
2643 >>> np.cov(X)
2644 array([[11.71 , -4.286 ], # may vary
2645 [-4.286 , 2.144133]])
2646 >>> np.cov(x, y)
2647 array([[11.71 , -4.286 ], # may vary
2648 [-4.286 , 2.144133]])
2649 >>> np.cov(x)
2650 array(11.71)
2652 """
2653 # Check inputs
2654 if ddof is not None and ddof != int(ddof):
2655 raise ValueError(
2656 "ddof must be integer")
2658 # Handles complex arrays too
2659 m = np.asarray(m)
2660 if m.ndim > 2:
2661 raise ValueError("m has more than 2 dimensions")
2663 if y is not None:
2664 y = np.asarray(y)
2665 if y.ndim > 2:
2666 raise ValueError("y has more than 2 dimensions")
2668 if dtype is None:
2669 if y is None:
2670 dtype = np.result_type(m, np.float64)
2671 else:
2672 dtype = np.result_type(m, y, np.float64)
2674 X = array(m, ndmin=2, dtype=dtype)
2675 if not rowvar and X.shape[0] != 1:
2676 X = X.T
2677 if X.shape[0] == 0:
2678 return np.array([]).reshape(0, 0)
2679 if y is not None:
2680 y = array(y, copy=False, ndmin=2, dtype=dtype)
2681 if not rowvar and y.shape[0] != 1:
2682 y = y.T
2683 X = np.concatenate((X, y), axis=0)
2685 if ddof is None:
2686 if bias == 0:
2687 ddof = 1
2688 else:
2689 ddof = 0
2691 # Get the product of frequencies and weights
2692 w = None
2693 if fweights is not None:
2694 fweights = np.asarray(fweights, dtype=float)
2695 if not np.all(fweights == np.around(fweights)):
2696 raise TypeError(
2697 "fweights must be integer")
2698 if fweights.ndim > 1:
2699 raise RuntimeError(
2700 "cannot handle multidimensional fweights")
2701 if fweights.shape[0] != X.shape[1]:
2702 raise RuntimeError(
2703 "incompatible numbers of samples and fweights")
2704 if any(fweights < 0):
2705 raise ValueError(
2706 "fweights cannot be negative")
2707 w = fweights
2708 if aweights is not None:
2709 aweights = np.asarray(aweights, dtype=float)
2710 if aweights.ndim > 1:
2711 raise RuntimeError(
2712 "cannot handle multidimensional aweights")
2713 if aweights.shape[0] != X.shape[1]:
2714 raise RuntimeError(
2715 "incompatible numbers of samples and aweights")
2716 if any(aweights < 0):
2717 raise ValueError(
2718 "aweights cannot be negative")
2719 if w is None:
2720 w = aweights
2721 else:
2722 w *= aweights
2724 avg, w_sum = average(X, axis=1, weights=w, returned=True)
2725 w_sum = w_sum[0]
2727 # Determine the normalization
2728 if w is None:
2729 fact = X.shape[1] - ddof
2730 elif ddof == 0:
2731 fact = w_sum
2732 elif aweights is None:
2733 fact = w_sum - ddof
2734 else:
2735 fact = w_sum - ddof*sum(w*aweights)/w_sum
2737 if fact <= 0:
2738 warnings.warn("Degrees of freedom <= 0 for slice",
2739 RuntimeWarning, stacklevel=2)
2740 fact = 0.0
2742 X -= avg[:, None]
2743 if w is None:
2744 X_T = X.T
2745 else:
2746 X_T = (X*w).T
2747 c = dot(X, X_T.conj())
2748 c *= np.true_divide(1, fact)
2749 return c.squeeze()
2752def _corrcoef_dispatcher(x, y=None, rowvar=None, bias=None, ddof=None, *,
2753 dtype=None):
2754 return (x, y)
2757@array_function_dispatch(_corrcoef_dispatcher)
2758def corrcoef(x, y=None, rowvar=True, bias=np._NoValue, ddof=np._NoValue, *,
2759 dtype=None):
2760 """
2761 Return Pearson product-moment correlation coefficients.
2763 Please refer to the documentation for `cov` for more detail. The
2764 relationship between the correlation coefficient matrix, `R`, and the
2765 covariance matrix, `C`, is
2767 .. math:: R_{ij} = \\frac{ C_{ij} } { \\sqrt{ C_{ii} C_{jj} } }
2769 The values of `R` are between -1 and 1, inclusive.
2771 Parameters
2772 ----------
2773 x : array_like
2774 A 1-D or 2-D array containing multiple variables and observations.
2775 Each row of `x` represents a variable, and each column a single
2776 observation of all those variables. Also see `rowvar` below.
2777 y : array_like, optional
2778 An additional set of variables and observations. `y` has the same
2779 shape as `x`.
2780 rowvar : bool, optional
2781 If `rowvar` is True (default), then each row represents a
2782 variable, with observations in the columns. Otherwise, the relationship
2783 is transposed: each column represents a variable, while the rows
2784 contain observations.
2785 bias : _NoValue, optional
2786 Has no effect, do not use.
2788 .. deprecated:: 1.10.0
2789 ddof : _NoValue, optional
2790 Has no effect, do not use.
2792 .. deprecated:: 1.10.0
2793 dtype : data-type, optional
2794 Data-type of the result. By default, the return data-type will have
2795 at least `numpy.float64` precision.
2797 .. versionadded:: 1.20
2799 Returns
2800 -------
2801 R : ndarray
2802 The correlation coefficient matrix of the variables.
2804 See Also
2805 --------
2806 cov : Covariance matrix
2808 Notes
2809 -----
2810 Due to floating point rounding the resulting array may not be Hermitian,
2811 the diagonal elements may not be 1, and the elements may not satisfy the
2812 inequality abs(a) <= 1. The real and imaginary parts are clipped to the
2813 interval [-1, 1] in an attempt to improve on that situation but is not
2814 much help in the complex case.
2816 This function accepts but discards arguments `bias` and `ddof`. This is
2817 for backwards compatibility with previous versions of this function. These
2818 arguments had no effect on the return values of the function and can be
2819 safely ignored in this and previous versions of numpy.
2821 Examples
2822 --------
2823 In this example we generate two random arrays, ``xarr`` and ``yarr``, and
2824 compute the row-wise and column-wise Pearson correlation coefficients,
2825 ``R``. Since ``rowvar`` is true by default, we first find the row-wise
2826 Pearson correlation coefficients between the variables of ``xarr``.
2828 >>> import numpy as np
2829 >>> rng = np.random.default_rng(seed=42)
2830 >>> xarr = rng.random((3, 3))
2831 >>> xarr
2832 array([[0.77395605, 0.43887844, 0.85859792],
2833 [0.69736803, 0.09417735, 0.97562235],
2834 [0.7611397 , 0.78606431, 0.12811363]])
2835 >>> R1 = np.corrcoef(xarr)
2836 >>> R1
2837 array([[ 1. , 0.99256089, -0.68080986],
2838 [ 0.99256089, 1. , -0.76492172],
2839 [-0.68080986, -0.76492172, 1. ]])
2841 If we add another set of variables and observations ``yarr``, we can
2842 compute the row-wise Pearson correlation coefficients between the
2843 variables in ``xarr`` and ``yarr``.
2845 >>> yarr = rng.random((3, 3))
2846 >>> yarr
2847 array([[0.45038594, 0.37079802, 0.92676499],
2848 [0.64386512, 0.82276161, 0.4434142 ],
2849 [0.22723872, 0.55458479, 0.06381726]])
2850 >>> R2 = np.corrcoef(xarr, yarr)
2851 >>> R2
2852 array([[ 1. , 0.99256089, -0.68080986, 0.75008178, -0.934284 ,
2853 -0.99004057],
2854 [ 0.99256089, 1. , -0.76492172, 0.82502011, -0.97074098,
2855 -0.99981569],
2856 [-0.68080986, -0.76492172, 1. , -0.99507202, 0.89721355,
2857 0.77714685],
2858 [ 0.75008178, 0.82502011, -0.99507202, 1. , -0.93657855,
2859 -0.83571711],
2860 [-0.934284 , -0.97074098, 0.89721355, -0.93657855, 1. ,
2861 0.97517215],
2862 [-0.99004057, -0.99981569, 0.77714685, -0.83571711, 0.97517215,
2863 1. ]])
2865 Finally if we use the option ``rowvar=False``, the columns are now
2866 being treated as the variables and we will find the column-wise Pearson
2867 correlation coefficients between variables in ``xarr`` and ``yarr``.
2869 >>> R3 = np.corrcoef(xarr, yarr, rowvar=False)
2870 >>> R3
2871 array([[ 1. , 0.77598074, -0.47458546, -0.75078643, -0.9665554 ,
2872 0.22423734],
2873 [ 0.77598074, 1. , -0.92346708, -0.99923895, -0.58826587,
2874 -0.44069024],
2875 [-0.47458546, -0.92346708, 1. , 0.93773029, 0.23297648,
2876 0.75137473],
2877 [-0.75078643, -0.99923895, 0.93773029, 1. , 0.55627469,
2878 0.47536961],
2879 [-0.9665554 , -0.58826587, 0.23297648, 0.55627469, 1. ,
2880 -0.46666491],
2881 [ 0.22423734, -0.44069024, 0.75137473, 0.47536961, -0.46666491,
2882 1. ]])
2884 """
2885 if bias is not np._NoValue or ddof is not np._NoValue:
2886 # 2015-03-15, 1.10
2887 warnings.warn('bias and ddof have no effect and are deprecated',
2888 DeprecationWarning, stacklevel=2)
2889 c = cov(x, y, rowvar, dtype=dtype)
2890 try:
2891 d = diag(c)
2892 except ValueError:
2893 # scalar covariance
2894 # nan if incorrect value (nan, inf, 0), 1 otherwise
2895 return c / c
2896 stddev = sqrt(d.real)
2897 c /= stddev[:, None]
2898 c /= stddev[None, :]
2900 # Clip real and imaginary parts to [-1, 1]. This does not guarantee
2901 # abs(a[i,j]) <= 1 for complex arrays, but is the best we can do without
2902 # excessive work.
2903 np.clip(c.real, -1, 1, out=c.real)
2904 if np.iscomplexobj(c):
2905 np.clip(c.imag, -1, 1, out=c.imag)
2907 return c
2910@set_module('numpy')
2911def blackman(M):
2912 """
2913 Return the Blackman window.
2915 The Blackman window is a taper formed by using the first three
2916 terms of a summation of cosines. It was designed to have close to the
2917 minimal leakage possible. It is close to optimal, only slightly worse
2918 than a Kaiser window.
2920 Parameters
2921 ----------
2922 M : int
2923 Number of points in the output window. If zero or less, an empty
2924 array is returned.
2926 Returns
2927 -------
2928 out : ndarray
2929 The window, with the maximum value normalized to one (the value one
2930 appears only if the number of samples is odd).
2932 See Also
2933 --------
2934 bartlett, hamming, hanning, kaiser
2936 Notes
2937 -----
2938 The Blackman window is defined as
2940 .. math:: w(n) = 0.42 - 0.5 \\cos(2\\pi n/M) + 0.08 \\cos(4\\pi n/M)
2942 Most references to the Blackman window come from the signal processing
2943 literature, where it is used as one of many windowing functions for
2944 smoothing values. It is also known as an apodization (which means
2945 "removing the foot", i.e. smoothing discontinuities at the beginning
2946 and end of the sampled signal) or tapering function. It is known as a
2947 "near optimal" tapering function, almost as good (by some measures)
2948 as the kaiser window.
2950 References
2951 ----------
2952 Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra,
2953 Dover Publications, New York.
2955 Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing.
2956 Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471.
2958 Examples
2959 --------
2960 >>> import matplotlib.pyplot as plt
2961 >>> np.blackman(12)
2962 array([-1.38777878e-17, 3.26064346e-02, 1.59903635e-01, # may vary
2963 4.14397981e-01, 7.36045180e-01, 9.67046769e-01,
2964 9.67046769e-01, 7.36045180e-01, 4.14397981e-01,
2965 1.59903635e-01, 3.26064346e-02, -1.38777878e-17])
2967 Plot the window and the frequency response:
2969 >>> from numpy.fft import fft, fftshift
2970 >>> window = np.blackman(51)
2971 >>> plt.plot(window)
2972 [<matplotlib.lines.Line2D object at 0x...>]
2973 >>> plt.title("Blackman window")
2974 Text(0.5, 1.0, 'Blackman window')
2975 >>> plt.ylabel("Amplitude")
2976 Text(0, 0.5, 'Amplitude')
2977 >>> plt.xlabel("Sample")
2978 Text(0.5, 0, 'Sample')
2979 >>> plt.show()
2981 >>> plt.figure()
2982 <Figure size 640x480 with 0 Axes>
2983 >>> A = fft(window, 2048) / 25.5
2984 >>> mag = np.abs(fftshift(A))
2985 >>> freq = np.linspace(-0.5, 0.5, len(A))
2986 >>> with np.errstate(divide='ignore', invalid='ignore'):
2987 ... response = 20 * np.log10(mag)
2988 ...
2989 >>> response = np.clip(response, -100, 100)
2990 >>> plt.plot(freq, response)
2991 [<matplotlib.lines.Line2D object at 0x...>]
2992 >>> plt.title("Frequency response of Blackman window")
2993 Text(0.5, 1.0, 'Frequency response of Blackman window')
2994 >>> plt.ylabel("Magnitude [dB]")
2995 Text(0, 0.5, 'Magnitude [dB]')
2996 >>> plt.xlabel("Normalized frequency [cycles per sample]")
2997 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
2998 >>> _ = plt.axis('tight')
2999 >>> plt.show()
3001 """
3002 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3003 # to double is safe for a range.
3004 values = np.array([0.0, M])
3005 M = values[1]
3007 if M < 1:
3008 return array([], dtype=values.dtype)
3009 if M == 1:
3010 return ones(1, dtype=values.dtype)
3011 n = arange(1-M, M, 2)
3012 return 0.42 + 0.5*cos(pi*n/(M-1)) + 0.08*cos(2.0*pi*n/(M-1))
3015@set_module('numpy')
3016def bartlett(M):
3017 """
3018 Return the Bartlett window.
3020 The Bartlett window is very similar to a triangular window, except
3021 that the end points are at zero. It is often used in signal
3022 processing for tapering a signal, without generating too much
3023 ripple in the frequency domain.
3025 Parameters
3026 ----------
3027 M : int
3028 Number of points in the output window. If zero or less, an
3029 empty array is returned.
3031 Returns
3032 -------
3033 out : array
3034 The triangular window, with the maximum value normalized to one
3035 (the value one appears only if the number of samples is odd), with
3036 the first and last samples equal to zero.
3038 See Also
3039 --------
3040 blackman, hamming, hanning, kaiser
3042 Notes
3043 -----
3044 The Bartlett window is defined as
3046 .. math:: w(n) = \\frac{2}{M-1} \\left(
3047 \\frac{M-1}{2} - \\left|n - \\frac{M-1}{2}\\right|
3048 \\right)
3050 Most references to the Bartlett window come from the signal processing
3051 literature, where it is used as one of many windowing functions for
3052 smoothing values. Note that convolution with this window produces linear
3053 interpolation. It is also known as an apodization (which means "removing
3054 the foot", i.e. smoothing discontinuities at the beginning and end of the
3055 sampled signal) or tapering function. The Fourier transform of the
3056 Bartlett window is the product of two sinc functions. Note the excellent
3057 discussion in Kanasewich [2]_.
3059 References
3060 ----------
3061 .. [1] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra",
3062 Biometrika 37, 1-16, 1950.
3063 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
3064 The University of Alberta Press, 1975, pp. 109-110.
3065 .. [3] A.V. Oppenheim and R.W. Schafer, "Discrete-Time Signal
3066 Processing", Prentice-Hall, 1999, pp. 468-471.
3067 .. [4] Wikipedia, "Window function",
3068 https://en.wikipedia.org/wiki/Window_function
3069 .. [5] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3070 "Numerical Recipes", Cambridge University Press, 1986, page 429.
3072 Examples
3073 --------
3074 >>> import matplotlib.pyplot as plt
3075 >>> np.bartlett(12)
3076 array([ 0. , 0.18181818, 0.36363636, 0.54545455, 0.72727273, # may vary
3077 0.90909091, 0.90909091, 0.72727273, 0.54545455, 0.36363636,
3078 0.18181818, 0. ])
3080 Plot the window and its frequency response (requires SciPy and matplotlib):
3082 >>> from numpy.fft import fft, fftshift
3083 >>> window = np.bartlett(51)
3084 >>> plt.plot(window)
3085 [<matplotlib.lines.Line2D object at 0x...>]
3086 >>> plt.title("Bartlett window")
3087 Text(0.5, 1.0, 'Bartlett window')
3088 >>> plt.ylabel("Amplitude")
3089 Text(0, 0.5, 'Amplitude')
3090 >>> plt.xlabel("Sample")
3091 Text(0.5, 0, 'Sample')
3092 >>> plt.show()
3094 >>> plt.figure()
3095 <Figure size 640x480 with 0 Axes>
3096 >>> A = fft(window, 2048) / 25.5
3097 >>> mag = np.abs(fftshift(A))
3098 >>> freq = np.linspace(-0.5, 0.5, len(A))
3099 >>> with np.errstate(divide='ignore', invalid='ignore'):
3100 ... response = 20 * np.log10(mag)
3101 ...
3102 >>> response = np.clip(response, -100, 100)
3103 >>> plt.plot(freq, response)
3104 [<matplotlib.lines.Line2D object at 0x...>]
3105 >>> plt.title("Frequency response of Bartlett window")
3106 Text(0.5, 1.0, 'Frequency response of Bartlett window')
3107 >>> plt.ylabel("Magnitude [dB]")
3108 Text(0, 0.5, 'Magnitude [dB]')
3109 >>> plt.xlabel("Normalized frequency [cycles per sample]")
3110 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
3111 >>> _ = plt.axis('tight')
3112 >>> plt.show()
3114 """
3115 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3116 # to double is safe for a range.
3117 values = np.array([0.0, M])
3118 M = values[1]
3120 if M < 1:
3121 return array([], dtype=values.dtype)
3122 if M == 1:
3123 return ones(1, dtype=values.dtype)
3124 n = arange(1-M, M, 2)
3125 return where(less_equal(n, 0), 1 + n/(M-1), 1 - n/(M-1))
3128@set_module('numpy')
3129def hanning(M):
3130 """
3131 Return the Hanning window.
3133 The Hanning window is a taper formed by using a weighted cosine.
3135 Parameters
3136 ----------
3137 M : int
3138 Number of points in the output window. If zero or less, an
3139 empty array is returned.
3141 Returns
3142 -------
3143 out : ndarray, shape(M,)
3144 The window, with the maximum value normalized to one (the value
3145 one appears only if `M` is odd).
3147 See Also
3148 --------
3149 bartlett, blackman, hamming, kaiser
3151 Notes
3152 -----
3153 The Hanning window is defined as
3155 .. math:: w(n) = 0.5 - 0.5\\cos\\left(\\frac{2\\pi{n}}{M-1}\\right)
3156 \\qquad 0 \\leq n \\leq M-1
3158 The Hanning was named for Julius von Hann, an Austrian meteorologist.
3159 It is also known as the Cosine Bell. Some authors prefer that it be
3160 called a Hann window, to help avoid confusion with the very similar
3161 Hamming window.
3163 Most references to the Hanning window come from the signal processing
3164 literature, where it is used as one of many windowing functions for
3165 smoothing values. It is also known as an apodization (which means
3166 "removing the foot", i.e. smoothing discontinuities at the beginning
3167 and end of the sampled signal) or tapering function.
3169 References
3170 ----------
3171 .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
3172 spectra, Dover Publications, New York.
3173 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
3174 The University of Alberta Press, 1975, pp. 106-108.
3175 .. [3] Wikipedia, "Window function",
3176 https://en.wikipedia.org/wiki/Window_function
3177 .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3178 "Numerical Recipes", Cambridge University Press, 1986, page 425.
3180 Examples
3181 --------
3182 >>> np.hanning(12)
3183 array([0. , 0.07937323, 0.29229249, 0.57115742, 0.82743037,
3184 0.97974649, 0.97974649, 0.82743037, 0.57115742, 0.29229249,
3185 0.07937323, 0. ])
3187 Plot the window and its frequency response:
3189 >>> import matplotlib.pyplot as plt
3190 >>> from numpy.fft import fft, fftshift
3191 >>> window = np.hanning(51)
3192 >>> plt.plot(window)
3193 [<matplotlib.lines.Line2D object at 0x...>]
3194 >>> plt.title("Hann window")
3195 Text(0.5, 1.0, 'Hann window')
3196 >>> plt.ylabel("Amplitude")
3197 Text(0, 0.5, 'Amplitude')
3198 >>> plt.xlabel("Sample")
3199 Text(0.5, 0, 'Sample')
3200 >>> plt.show()
3202 >>> plt.figure()
3203 <Figure size 640x480 with 0 Axes>
3204 >>> A = fft(window, 2048) / 25.5
3205 >>> mag = np.abs(fftshift(A))
3206 >>> freq = np.linspace(-0.5, 0.5, len(A))
3207 >>> with np.errstate(divide='ignore', invalid='ignore'):
3208 ... response = 20 * np.log10(mag)
3209 ...
3210 >>> response = np.clip(response, -100, 100)
3211 >>> plt.plot(freq, response)
3212 [<matplotlib.lines.Line2D object at 0x...>]
3213 >>> plt.title("Frequency response of the Hann window")
3214 Text(0.5, 1.0, 'Frequency response of the Hann window')
3215 >>> plt.ylabel("Magnitude [dB]")
3216 Text(0, 0.5, 'Magnitude [dB]')
3217 >>> plt.xlabel("Normalized frequency [cycles per sample]")
3218 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
3219 >>> plt.axis('tight')
3220 ...
3221 >>> plt.show()
3223 """
3224 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3225 # to double is safe for a range.
3226 values = np.array([0.0, M])
3227 M = values[1]
3229 if M < 1:
3230 return array([], dtype=values.dtype)
3231 if M == 1:
3232 return ones(1, dtype=values.dtype)
3233 n = arange(1-M, M, 2)
3234 return 0.5 + 0.5*cos(pi*n/(M-1))
3237@set_module('numpy')
3238def hamming(M):
3239 """
3240 Return the Hamming window.
3242 The Hamming window is a taper formed by using a weighted cosine.
3244 Parameters
3245 ----------
3246 M : int
3247 Number of points in the output window. If zero or less, an
3248 empty array is returned.
3250 Returns
3251 -------
3252 out : ndarray
3253 The window, with the maximum value normalized to one (the value
3254 one appears only if the number of samples is odd).
3256 See Also
3257 --------
3258 bartlett, blackman, hanning, kaiser
3260 Notes
3261 -----
3262 The Hamming window is defined as
3264 .. math:: w(n) = 0.54 - 0.46\\cos\\left(\\frac{2\\pi{n}}{M-1}\\right)
3265 \\qquad 0 \\leq n \\leq M-1
3267 The Hamming was named for R. W. Hamming, an associate of J. W. Tukey
3268 and is described in Blackman and Tukey. It was recommended for
3269 smoothing the truncated autocovariance function in the time domain.
3270 Most references to the Hamming window come from the signal processing
3271 literature, where it is used as one of many windowing functions for
3272 smoothing values. It is also known as an apodization (which means
3273 "removing the foot", i.e. smoothing discontinuities at the beginning
3274 and end of the sampled signal) or tapering function.
3276 References
3277 ----------
3278 .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
3279 spectra, Dover Publications, New York.
3280 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
3281 University of Alberta Press, 1975, pp. 109-110.
3282 .. [3] Wikipedia, "Window function",
3283 https://en.wikipedia.org/wiki/Window_function
3284 .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3285 "Numerical Recipes", Cambridge University Press, 1986, page 425.
3287 Examples
3288 --------
3289 >>> np.hamming(12)
3290 array([ 0.08 , 0.15302337, 0.34890909, 0.60546483, 0.84123594, # may vary
3291 0.98136677, 0.98136677, 0.84123594, 0.60546483, 0.34890909,
3292 0.15302337, 0.08 ])
3294 Plot the window and the frequency response:
3296 >>> import matplotlib.pyplot as plt
3297 >>> from numpy.fft import fft, fftshift
3298 >>> window = np.hamming(51)
3299 >>> plt.plot(window)
3300 [<matplotlib.lines.Line2D object at 0x...>]
3301 >>> plt.title("Hamming window")
3302 Text(0.5, 1.0, 'Hamming window')
3303 >>> plt.ylabel("Amplitude")
3304 Text(0, 0.5, 'Amplitude')
3305 >>> plt.xlabel("Sample")
3306 Text(0.5, 0, 'Sample')
3307 >>> plt.show()
3309 >>> plt.figure()
3310 <Figure size 640x480 with 0 Axes>
3311 >>> A = fft(window, 2048) / 25.5
3312 >>> mag = np.abs(fftshift(A))
3313 >>> freq = np.linspace(-0.5, 0.5, len(A))
3314 >>> response = 20 * np.log10(mag)
3315 >>> response = np.clip(response, -100, 100)
3316 >>> plt.plot(freq, response)
3317 [<matplotlib.lines.Line2D object at 0x...>]
3318 >>> plt.title("Frequency response of Hamming window")
3319 Text(0.5, 1.0, 'Frequency response of Hamming window')
3320 >>> plt.ylabel("Magnitude [dB]")
3321 Text(0, 0.5, 'Magnitude [dB]')
3322 >>> plt.xlabel("Normalized frequency [cycles per sample]")
3323 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
3324 >>> plt.axis('tight')
3325 ...
3326 >>> plt.show()
3328 """
3329 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3330 # to double is safe for a range.
3331 values = np.array([0.0, M])
3332 M = values[1]
3334 if M < 1:
3335 return array([], dtype=values.dtype)
3336 if M == 1:
3337 return ones(1, dtype=values.dtype)
3338 n = arange(1-M, M, 2)
3339 return 0.54 + 0.46*cos(pi*n/(M-1))
3342## Code from cephes for i0
3344_i0A = [
3345 -4.41534164647933937950E-18,
3346 3.33079451882223809783E-17,
3347 -2.43127984654795469359E-16,
3348 1.71539128555513303061E-15,
3349 -1.16853328779934516808E-14,
3350 7.67618549860493561688E-14,
3351 -4.85644678311192946090E-13,
3352 2.95505266312963983461E-12,
3353 -1.72682629144155570723E-11,
3354 9.67580903537323691224E-11,
3355 -5.18979560163526290666E-10,
3356 2.65982372468238665035E-9,
3357 -1.30002500998624804212E-8,
3358 6.04699502254191894932E-8,
3359 -2.67079385394061173391E-7,
3360 1.11738753912010371815E-6,
3361 -4.41673835845875056359E-6,
3362 1.64484480707288970893E-5,
3363 -5.75419501008210370398E-5,
3364 1.88502885095841655729E-4,
3365 -5.76375574538582365885E-4,
3366 1.63947561694133579842E-3,
3367 -4.32430999505057594430E-3,
3368 1.05464603945949983183E-2,
3369 -2.37374148058994688156E-2,
3370 4.93052842396707084878E-2,
3371 -9.49010970480476444210E-2,
3372 1.71620901522208775349E-1,
3373 -3.04682672343198398683E-1,
3374 6.76795274409476084995E-1
3375 ]
3377_i0B = [
3378 -7.23318048787475395456E-18,
3379 -4.83050448594418207126E-18,
3380 4.46562142029675999901E-17,
3381 3.46122286769746109310E-17,
3382 -2.82762398051658348494E-16,
3383 -3.42548561967721913462E-16,
3384 1.77256013305652638360E-15,
3385 3.81168066935262242075E-15,
3386 -9.55484669882830764870E-15,
3387 -4.15056934728722208663E-14,
3388 1.54008621752140982691E-14,
3389 3.85277838274214270114E-13,
3390 7.18012445138366623367E-13,
3391 -1.79417853150680611778E-12,
3392 -1.32158118404477131188E-11,
3393 -3.14991652796324136454E-11,
3394 1.18891471078464383424E-11,
3395 4.94060238822496958910E-10,
3396 3.39623202570838634515E-9,
3397 2.26666899049817806459E-8,
3398 2.04891858946906374183E-7,
3399 2.89137052083475648297E-6,
3400 6.88975834691682398426E-5,
3401 3.36911647825569408990E-3,
3402 8.04490411014108831608E-1
3403 ]
3406def _chbevl(x, vals):
3407 b0 = vals[0]
3408 b1 = 0.0
3410 for i in range(1, len(vals)):
3411 b2 = b1
3412 b1 = b0
3413 b0 = x*b1 - b2 + vals[i]
3415 return 0.5*(b0 - b2)
3418def _i0_1(x):
3419 return exp(x) * _chbevl(x/2.0-2, _i0A)
3422def _i0_2(x):
3423 return exp(x) * _chbevl(32.0/x - 2.0, _i0B) / sqrt(x)
3426def _i0_dispatcher(x):
3427 return (x,)
3430@array_function_dispatch(_i0_dispatcher)
3431def i0(x):
3432 """
3433 Modified Bessel function of the first kind, order 0.
3435 Usually denoted :math:`I_0`.
3437 Parameters
3438 ----------
3439 x : array_like of float
3440 Argument of the Bessel function.
3442 Returns
3443 -------
3444 out : ndarray, shape = x.shape, dtype = float
3445 The modified Bessel function evaluated at each of the elements of `x`.
3447 See Also
3448 --------
3449 scipy.special.i0, scipy.special.iv, scipy.special.ive
3451 Notes
3452 -----
3453 The scipy implementation is recommended over this function: it is a
3454 proper ufunc written in C, and more than an order of magnitude faster.
3456 We use the algorithm published by Clenshaw [1]_ and referenced by
3457 Abramowitz and Stegun [2]_, for which the function domain is
3458 partitioned into the two intervals [0,8] and (8,inf), and Chebyshev
3459 polynomial expansions are employed in each interval. Relative error on
3460 the domain [0,30] using IEEE arithmetic is documented [3]_ as having a
3461 peak of 5.8e-16 with an rms of 1.4e-16 (n = 30000).
3463 References
3464 ----------
3465 .. [1] C. W. Clenshaw, "Chebyshev series for mathematical functions", in
3466 *National Physical Laboratory Mathematical Tables*, vol. 5, London:
3467 Her Majesty's Stationery Office, 1962.
3468 .. [2] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical
3469 Functions*, 10th printing, New York: Dover, 1964, pp. 379.
3470 https://personal.math.ubc.ca/~cbm/aands/page_379.htm
3471 .. [3] https://metacpan.org/pod/distribution/Math-Cephes/lib/Math/Cephes.pod#i0:-Modified-Bessel-function-of-order-zero
3473 Examples
3474 --------
3475 >>> np.i0(0.)
3476 array(1.0)
3477 >>> np.i0([0, 1, 2, 3])
3478 array([1. , 1.26606588, 2.2795853 , 4.88079259])
3480 """
3481 x = np.asanyarray(x)
3482 if x.dtype.kind == 'c':
3483 raise TypeError("i0 not supported for complex values")
3484 if x.dtype.kind != 'f':
3485 x = x.astype(float)
3486 x = np.abs(x)
3487 return piecewise(x, [x <= 8.0], [_i0_1, _i0_2])
3489## End of cephes code for i0
3492@set_module('numpy')
3493def kaiser(M, beta):
3494 """
3495 Return the Kaiser window.
3497 The Kaiser window is a taper formed by using a Bessel function.
3499 Parameters
3500 ----------
3501 M : int
3502 Number of points in the output window. If zero or less, an
3503 empty array is returned.
3504 beta : float
3505 Shape parameter for window.
3507 Returns
3508 -------
3509 out : array
3510 The window, with the maximum value normalized to one (the value
3511 one appears only if the number of samples is odd).
3513 See Also
3514 --------
3515 bartlett, blackman, hamming, hanning
3517 Notes
3518 -----
3519 The Kaiser window is defined as
3521 .. math:: w(n) = I_0\\left( \\beta \\sqrt{1-\\frac{4n^2}{(M-1)^2}}
3522 \\right)/I_0(\\beta)
3524 with
3526 .. math:: \\quad -\\frac{M-1}{2} \\leq n \\leq \\frac{M-1}{2},
3528 where :math:`I_0` is the modified zeroth-order Bessel function.
3530 The Kaiser was named for Jim Kaiser, who discovered a simple
3531 approximation to the DPSS window based on Bessel functions. The Kaiser
3532 window is a very good approximation to the Digital Prolate Spheroidal
3533 Sequence, or Slepian window, which is the transform which maximizes the
3534 energy in the main lobe of the window relative to total energy.
3536 The Kaiser can approximate many other windows by varying the beta
3537 parameter.
3539 ==== =======================
3540 beta Window shape
3541 ==== =======================
3542 0 Rectangular
3543 5 Similar to a Hamming
3544 6 Similar to a Hanning
3545 8.6 Similar to a Blackman
3546 ==== =======================
3548 A beta value of 14 is probably a good starting point. Note that as beta
3549 gets large, the window narrows, and so the number of samples needs to be
3550 large enough to sample the increasingly narrow spike, otherwise NaNs will
3551 get returned.
3553 Most references to the Kaiser window come from the signal processing
3554 literature, where it is used as one of many windowing functions for
3555 smoothing values. It is also known as an apodization (which means
3556 "removing the foot", i.e. smoothing discontinuities at the beginning
3557 and end of the sampled signal) or tapering function.
3559 References
3560 ----------
3561 .. [1] J. F. Kaiser, "Digital Filters" - Ch 7 in "Systems analysis by
3562 digital computer", Editors: F.F. Kuo and J.F. Kaiser, p 218-285.
3563 John Wiley and Sons, New York, (1966).
3564 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
3565 University of Alberta Press, 1975, pp. 177-178.
3566 .. [3] Wikipedia, "Window function",
3567 https://en.wikipedia.org/wiki/Window_function
3569 Examples
3570 --------
3571 >>> import matplotlib.pyplot as plt
3572 >>> np.kaiser(12, 14)
3573 array([7.72686684e-06, 3.46009194e-03, 4.65200189e-02, # may vary
3574 2.29737120e-01, 5.99885316e-01, 9.45674898e-01,
3575 9.45674898e-01, 5.99885316e-01, 2.29737120e-01,
3576 4.65200189e-02, 3.46009194e-03, 7.72686684e-06])
3579 Plot the window and the frequency response:
3581 >>> from numpy.fft import fft, fftshift
3582 >>> window = np.kaiser(51, 14)
3583 >>> plt.plot(window)
3584 [<matplotlib.lines.Line2D object at 0x...>]
3585 >>> plt.title("Kaiser window")
3586 Text(0.5, 1.0, 'Kaiser window')
3587 >>> plt.ylabel("Amplitude")
3588 Text(0, 0.5, 'Amplitude')
3589 >>> plt.xlabel("Sample")
3590 Text(0.5, 0, 'Sample')
3591 >>> plt.show()
3593 >>> plt.figure()
3594 <Figure size 640x480 with 0 Axes>
3595 >>> A = fft(window, 2048) / 25.5
3596 >>> mag = np.abs(fftshift(A))
3597 >>> freq = np.linspace(-0.5, 0.5, len(A))
3598 >>> response = 20 * np.log10(mag)
3599 >>> response = np.clip(response, -100, 100)
3600 >>> plt.plot(freq, response)
3601 [<matplotlib.lines.Line2D object at 0x...>]
3602 >>> plt.title("Frequency response of Kaiser window")
3603 Text(0.5, 1.0, 'Frequency response of Kaiser window')
3604 >>> plt.ylabel("Magnitude [dB]")
3605 Text(0, 0.5, 'Magnitude [dB]')
3606 >>> plt.xlabel("Normalized frequency [cycles per sample]")
3607 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
3608 >>> plt.axis('tight')
3609 (-0.5, 0.5, -100.0, ...) # may vary
3610 >>> plt.show()
3612 """
3613 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3614 # to double is safe for a range. (Simplified result_type with 0.0
3615 # strongly typed. result-type is not/less order sensitive, but that mainly
3616 # matters for integers anyway.)
3617 values = np.array([0.0, M, beta])
3618 M = values[1]
3619 beta = values[2]
3621 if M == 1:
3622 return np.ones(1, dtype=values.dtype)
3623 n = arange(0, M)
3624 alpha = (M-1)/2.0
3625 return i0(beta * sqrt(1-((n-alpha)/alpha)**2.0))/i0(beta)
3628def _sinc_dispatcher(x):
3629 return (x,)
3632@array_function_dispatch(_sinc_dispatcher)
3633def sinc(x):
3634 r"""
3635 Return the normalized sinc function.
3637 The sinc function is equal to :math:`\sin(\pi x)/(\pi x)` for any argument
3638 :math:`x\ne 0`. ``sinc(0)`` takes the limit value 1, making ``sinc`` not
3639 only everywhere continuous but also infinitely differentiable.
3641 .. note::
3643 Note the normalization factor of ``pi`` used in the definition.
3644 This is the most commonly used definition in signal processing.
3645 Use ``sinc(x / np.pi)`` to obtain the unnormalized sinc function
3646 :math:`\sin(x)/x` that is more common in mathematics.
3648 Parameters
3649 ----------
3650 x : ndarray
3651 Array (possibly multi-dimensional) of values for which to calculate
3652 ``sinc(x)``.
3654 Returns
3655 -------
3656 out : ndarray
3657 ``sinc(x)``, which has the same shape as the input.
3659 Notes
3660 -----
3661 The name sinc is short for "sine cardinal" or "sinus cardinalis".
3663 The sinc function is used in various signal processing applications,
3664 including in anti-aliasing, in the construction of a Lanczos resampling
3665 filter, and in interpolation.
3667 For bandlimited interpolation of discrete-time signals, the ideal
3668 interpolation kernel is proportional to the sinc function.
3670 References
3671 ----------
3672 .. [1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram Web
3673 Resource. http://mathworld.wolfram.com/SincFunction.html
3674 .. [2] Wikipedia, "Sinc function",
3675 https://en.wikipedia.org/wiki/Sinc_function
3677 Examples
3678 --------
3679 >>> import matplotlib.pyplot as plt
3680 >>> x = np.linspace(-4, 4, 41)
3681 >>> np.sinc(x)
3682 array([-3.89804309e-17, -4.92362781e-02, -8.40918587e-02, # may vary
3683 -8.90384387e-02, -5.84680802e-02, 3.89804309e-17,
3684 6.68206631e-02, 1.16434881e-01, 1.26137788e-01,
3685 8.50444803e-02, -3.89804309e-17, -1.03943254e-01,
3686 -1.89206682e-01, -2.16236208e-01, -1.55914881e-01,
3687 3.89804309e-17, 2.33872321e-01, 5.04551152e-01,
3688 7.56826729e-01, 9.35489284e-01, 1.00000000e+00,
3689 9.35489284e-01, 7.56826729e-01, 5.04551152e-01,
3690 2.33872321e-01, 3.89804309e-17, -1.55914881e-01,
3691 -2.16236208e-01, -1.89206682e-01, -1.03943254e-01,
3692 -3.89804309e-17, 8.50444803e-02, 1.26137788e-01,
3693 1.16434881e-01, 6.68206631e-02, 3.89804309e-17,
3694 -5.84680802e-02, -8.90384387e-02, -8.40918587e-02,
3695 -4.92362781e-02, -3.89804309e-17])
3697 >>> plt.plot(x, np.sinc(x))
3698 [<matplotlib.lines.Line2D object at 0x...>]
3699 >>> plt.title("Sinc Function")
3700 Text(0.5, 1.0, 'Sinc Function')
3701 >>> plt.ylabel("Amplitude")
3702 Text(0, 0.5, 'Amplitude')
3703 >>> plt.xlabel("X")
3704 Text(0.5, 0, 'X')
3705 >>> plt.show()
3707 """
3708 x = np.asanyarray(x)
3709 y = pi * where(x == 0, 1.0e-20, x)
3710 return sin(y)/y
3713def _msort_dispatcher(a):
3714 return (a,)
3717@array_function_dispatch(_msort_dispatcher)
3718def msort(a):
3719 """
3720 Return a copy of an array sorted along the first axis.
3722 .. deprecated:: 1.24
3724 msort is deprecated, use ``np.sort(a, axis=0)`` instead.
3726 Parameters
3727 ----------
3728 a : array_like
3729 Array to be sorted.
3731 Returns
3732 -------
3733 sorted_array : ndarray
3734 Array of the same type and shape as `a`.
3736 See Also
3737 --------
3738 sort
3740 Notes
3741 -----
3742 ``np.msort(a)`` is equivalent to ``np.sort(a, axis=0)``.
3744 Examples
3745 --------
3746 >>> a = np.array([[1, 4], [3, 1]])
3747 >>> np.msort(a) # sort along the first axis
3748 array([[1, 1],
3749 [3, 4]])
3751 """
3752 # 2022-10-20 1.24
3753 warnings.warn(
3754 "msort is deprecated, use np.sort(a, axis=0) instead",
3755 DeprecationWarning,
3756 stacklevel=2,
3757 )
3758 b = array(a, subok=True, copy=True)
3759 b.sort(0)
3760 return b
3763def _ureduce(a, func, keepdims=False, **kwargs):
3764 """
3765 Internal Function.
3766 Call `func` with `a` as first argument swapping the axes to use extended
3767 axis on functions that don't support it natively.
3769 Returns result and a.shape with axis dims set to 1.
3771 Parameters
3772 ----------
3773 a : array_like
3774 Input array or object that can be converted to an array.
3775 func : callable
3776 Reduction function capable of receiving a single axis argument.
3777 It is called with `a` as first argument followed by `kwargs`.
3778 kwargs : keyword arguments
3779 additional keyword arguments to pass to `func`.
3781 Returns
3782 -------
3783 result : tuple
3784 Result of func(a, **kwargs) and a.shape with axis dims set to 1
3785 which can be used to reshape the result to the same shape a ufunc with
3786 keepdims=True would produce.
3788 """
3789 a = np.asanyarray(a)
3790 axis = kwargs.get('axis', None)
3791 out = kwargs.get('out', None)
3793 if keepdims is np._NoValue:
3794 keepdims = False
3796 nd = a.ndim
3797 if axis is not None:
3798 axis = _nx.normalize_axis_tuple(axis, nd)
3800 if keepdims:
3801 if out is not None:
3802 index_out = tuple(
3803 0 if i in axis else slice(None) for i in range(nd))
3804 kwargs['out'] = out[(Ellipsis, ) + index_out]
3806 if len(axis) == 1:
3807 kwargs['axis'] = axis[0]
3808 else:
3809 keep = set(range(nd)) - set(axis)
3810 nkeep = len(keep)
3811 # swap axis that should not be reduced to front
3812 for i, s in enumerate(sorted(keep)):
3813 a = a.swapaxes(i, s)
3814 # merge reduced axis
3815 a = a.reshape(a.shape[:nkeep] + (-1,))
3816 kwargs['axis'] = -1
3817 else:
3818 if keepdims:
3819 if out is not None:
3820 index_out = (0, ) * nd
3821 kwargs['out'] = out[(Ellipsis, ) + index_out]
3823 r = func(a, **kwargs)
3825 if out is not None:
3826 return out
3828 if keepdims:
3829 if axis is None:
3830 index_r = (np.newaxis, ) * nd
3831 else:
3832 index_r = tuple(
3833 np.newaxis if i in axis else slice(None)
3834 for i in range(nd))
3835 r = r[(Ellipsis, ) + index_r]
3837 return r
3840def _median_dispatcher(
3841 a, axis=None, out=None, overwrite_input=None, keepdims=None):
3842 return (a, out)
3845@array_function_dispatch(_median_dispatcher)
3846def median(a, axis=None, out=None, overwrite_input=False, keepdims=False):
3847 """
3848 Compute the median along the specified axis.
3850 Returns the median of the array elements.
3852 Parameters
3853 ----------
3854 a : array_like
3855 Input array or object that can be converted to an array.
3856 axis : {int, sequence of int, None}, optional
3857 Axis or axes along which the medians are computed. The default
3858 is to compute the median along a flattened version of the array.
3859 A sequence of axes is supported since version 1.9.0.
3860 out : ndarray, optional
3861 Alternative output array in which to place the result. It must
3862 have the same shape and buffer length as the expected output,
3863 but the type (of the output) will be cast if necessary.
3864 overwrite_input : bool, optional
3865 If True, then allow use of memory of input array `a` for
3866 calculations. The input array will be modified by the call to
3867 `median`. This will save memory when you do not need to preserve
3868 the contents of the input array. Treat the input as undefined,
3869 but it will probably be fully or partially sorted. Default is
3870 False. If `overwrite_input` is ``True`` and `a` is not already an
3871 `ndarray`, an error will be raised.
3872 keepdims : bool, optional
3873 If this is set to True, the axes which are reduced are left
3874 in the result as dimensions with size one. With this option,
3875 the result will broadcast correctly against the original `arr`.
3877 .. versionadded:: 1.9.0
3879 Returns
3880 -------
3881 median : ndarray
3882 A new array holding the result. If the input contains integers
3883 or floats smaller than ``float64``, then the output data-type is
3884 ``np.float64``. Otherwise, the data-type of the output is the
3885 same as that of the input. If `out` is specified, that array is
3886 returned instead.
3888 See Also
3889 --------
3890 mean, percentile
3892 Notes
3893 -----
3894 Given a vector ``V`` of length ``N``, the median of ``V`` is the
3895 middle value of a sorted copy of ``V``, ``V_sorted`` - i
3896 e., ``V_sorted[(N-1)/2]``, when ``N`` is odd, and the average of the
3897 two middle values of ``V_sorted`` when ``N`` is even.
3899 Examples
3900 --------
3901 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
3902 >>> a
3903 array([[10, 7, 4],
3904 [ 3, 2, 1]])
3905 >>> np.median(a)
3906 3.5
3907 >>> np.median(a, axis=0)
3908 array([6.5, 4.5, 2.5])
3909 >>> np.median(a, axis=1)
3910 array([7., 2.])
3911 >>> m = np.median(a, axis=0)
3912 >>> out = np.zeros_like(m)
3913 >>> np.median(a, axis=0, out=m)
3914 array([6.5, 4.5, 2.5])
3915 >>> m
3916 array([6.5, 4.5, 2.5])
3917 >>> b = a.copy()
3918 >>> np.median(b, axis=1, overwrite_input=True)
3919 array([7., 2.])
3920 >>> assert not np.all(a==b)
3921 >>> b = a.copy()
3922 >>> np.median(b, axis=None, overwrite_input=True)
3923 3.5
3924 >>> assert not np.all(a==b)
3926 """
3927 return _ureduce(a, func=_median, keepdims=keepdims, axis=axis, out=out,
3928 overwrite_input=overwrite_input)
3931def _median(a, axis=None, out=None, overwrite_input=False):
3932 # can't be reasonably be implemented in terms of percentile as we have to
3933 # call mean to not break astropy
3934 a = np.asanyarray(a)
3936 # Set the partition indexes
3937 if axis is None:
3938 sz = a.size
3939 else:
3940 sz = a.shape[axis]
3941 if sz % 2 == 0:
3942 szh = sz // 2
3943 kth = [szh - 1, szh]
3944 else:
3945 kth = [(sz - 1) // 2]
3947 # We have to check for NaNs (as of writing 'M' doesn't actually work).
3948 supports_nans = np.issubdtype(a.dtype, np.inexact) or a.dtype.kind in 'Mm'
3949 if supports_nans:
3950 kth.append(-1)
3952 if overwrite_input:
3953 if axis is None:
3954 part = a.ravel()
3955 part.partition(kth)
3956 else:
3957 a.partition(kth, axis=axis)
3958 part = a
3959 else:
3960 part = partition(a, kth, axis=axis)
3962 if part.shape == ():
3963 # make 0-D arrays work
3964 return part.item()
3965 if axis is None:
3966 axis = 0
3968 indexer = [slice(None)] * part.ndim
3969 index = part.shape[axis] // 2
3970 if part.shape[axis] % 2 == 1:
3971 # index with slice to allow mean (below) to work
3972 indexer[axis] = slice(index, index+1)
3973 else:
3974 indexer[axis] = slice(index-1, index+1)
3975 indexer = tuple(indexer)
3977 # Use mean in both odd and even case to coerce data type,
3978 # using out array if needed.
3979 rout = mean(part[indexer], axis=axis, out=out)
3980 if supports_nans and sz > 0:
3981 # If nans are possible, warn and replace by nans like mean would.
3982 rout = np.lib.utils._median_nancheck(part, rout, axis)
3984 return rout
3987def _percentile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
3988 method=None, keepdims=None, *, interpolation=None):
3989 return (a, q, out)
3992@array_function_dispatch(_percentile_dispatcher)
3993def percentile(a,
3994 q,
3995 axis=None,
3996 out=None,
3997 overwrite_input=False,
3998 method="linear",
3999 keepdims=False,
4000 *,
4001 interpolation=None):
4002 """
4003 Compute the q-th percentile of the data along the specified axis.
4005 Returns the q-th percentile(s) of the array elements.
4007 Parameters
4008 ----------
4009 a : array_like of real numbers
4010 Input array or object that can be converted to an array.
4011 q : array_like of float
4012 Percentage or sequence of percentages for the percentiles to compute.
4013 Values must be between 0 and 100 inclusive.
4014 axis : {int, tuple of int, None}, optional
4015 Axis or axes along which the percentiles are computed. The
4016 default is to compute the percentile(s) along a flattened
4017 version of the array.
4019 .. versionchanged:: 1.9.0
4020 A tuple of axes is supported
4021 out : ndarray, optional
4022 Alternative output array in which to place the result. It must
4023 have the same shape and buffer length as the expected output,
4024 but the type (of the output) will be cast if necessary.
4025 overwrite_input : bool, optional
4026 If True, then allow the input array `a` to be modified by intermediate
4027 calculations, to save memory. In this case, the contents of the input
4028 `a` after this function completes is undefined.
4029 method : str, optional
4030 This parameter specifies the method to use for estimating the
4031 percentile. There are many different methods, some unique to NumPy.
4032 See the notes for explanation. The options sorted by their R type
4033 as summarized in the H&F paper [1]_ are:
4035 1. 'inverted_cdf'
4036 2. 'averaged_inverted_cdf'
4037 3. 'closest_observation'
4038 4. 'interpolated_inverted_cdf'
4039 5. 'hazen'
4040 6. 'weibull'
4041 7. 'linear' (default)
4042 8. 'median_unbiased'
4043 9. 'normal_unbiased'
4045 The first three methods are discontinuous. NumPy further defines the
4046 following discontinuous variations of the default 'linear' (7.) option:
4048 * 'lower'
4049 * 'higher',
4050 * 'midpoint'
4051 * 'nearest'
4053 .. versionchanged:: 1.22.0
4054 This argument was previously called "interpolation" and only
4055 offered the "linear" default and last four options.
4057 keepdims : bool, optional
4058 If this is set to True, the axes which are reduced are left in
4059 the result as dimensions with size one. With this option, the
4060 result will broadcast correctly against the original array `a`.
4062 .. versionadded:: 1.9.0
4064 interpolation : str, optional
4065 Deprecated name for the method keyword argument.
4067 .. deprecated:: 1.22.0
4069 Returns
4070 -------
4071 percentile : scalar or ndarray
4072 If `q` is a single percentile and `axis=None`, then the result
4073 is a scalar. If multiple percentiles are given, first axis of
4074 the result corresponds to the percentiles. The other axes are
4075 the axes that remain after the reduction of `a`. If the input
4076 contains integers or floats smaller than ``float64``, the output
4077 data-type is ``float64``. Otherwise, the output data-type is the
4078 same as that of the input. If `out` is specified, that array is
4079 returned instead.
4081 See Also
4082 --------
4083 mean
4084 median : equivalent to ``percentile(..., 50)``
4085 nanpercentile
4086 quantile : equivalent to percentile, except q in the range [0, 1].
4088 Notes
4089 -----
4090 Given a vector ``V`` of length ``n``, the q-th percentile of ``V`` is
4091 the value ``q/100`` of the way from the minimum to the maximum in a
4092 sorted copy of ``V``. The values and distances of the two nearest
4093 neighbors as well as the `method` parameter will determine the
4094 percentile if the normalized ranking does not match the location of
4095 ``q`` exactly. This function is the same as the median if ``q=50``, the
4096 same as the minimum if ``q=0`` and the same as the maximum if
4097 ``q=100``.
4099 The optional `method` parameter specifies the method to use when the
4100 desired percentile lies between two indexes ``i`` and ``j = i + 1``.
4101 In that case, we first determine ``i + g``, a virtual index that lies
4102 between ``i`` and ``j``, where ``i`` is the floor and ``g`` is the
4103 fractional part of the index. The final result is, then, an interpolation
4104 of ``a[i]`` and ``a[j]`` based on ``g``. During the computation of ``g``,
4105 ``i`` and ``j`` are modified using correction constants ``alpha`` and
4106 ``beta`` whose choices depend on the ``method`` used. Finally, note that
4107 since Python uses 0-based indexing, the code subtracts another 1 from the
4108 index internally.
4110 The following formula determines the virtual index ``i + g``, the location
4111 of the percentile in the sorted sample:
4113 .. math::
4114 i + g = (q / 100) * ( n - alpha - beta + 1 ) + alpha
4116 The different methods then work as follows
4118 inverted_cdf:
4119 method 1 of H&F [1]_.
4120 This method gives discontinuous results:
4122 * if g > 0 ; then take j
4123 * if g = 0 ; then take i
4125 averaged_inverted_cdf:
4126 method 2 of H&F [1]_.
4127 This method give discontinuous results:
4129 * if g > 0 ; then take j
4130 * if g = 0 ; then average between bounds
4132 closest_observation:
4133 method 3 of H&F [1]_.
4134 This method give discontinuous results:
4136 * if g > 0 ; then take j
4137 * if g = 0 and index is odd ; then take j
4138 * if g = 0 and index is even ; then take i
4140 interpolated_inverted_cdf:
4141 method 4 of H&F [1]_.
4142 This method give continuous results using:
4144 * alpha = 0
4145 * beta = 1
4147 hazen:
4148 method 5 of H&F [1]_.
4149 This method give continuous results using:
4151 * alpha = 1/2
4152 * beta = 1/2
4154 weibull:
4155 method 6 of H&F [1]_.
4156 This method give continuous results using:
4158 * alpha = 0
4159 * beta = 0
4161 linear:
4162 method 7 of H&F [1]_.
4163 This method give continuous results using:
4165 * alpha = 1
4166 * beta = 1
4168 median_unbiased:
4169 method 8 of H&F [1]_.
4170 This method is probably the best method if the sample
4171 distribution function is unknown (see reference).
4172 This method give continuous results using:
4174 * alpha = 1/3
4175 * beta = 1/3
4177 normal_unbiased:
4178 method 9 of H&F [1]_.
4179 This method is probably the best method if the sample
4180 distribution function is known to be normal.
4181 This method give continuous results using:
4183 * alpha = 3/8
4184 * beta = 3/8
4186 lower:
4187 NumPy method kept for backwards compatibility.
4188 Takes ``i`` as the interpolation point.
4190 higher:
4191 NumPy method kept for backwards compatibility.
4192 Takes ``j`` as the interpolation point.
4194 nearest:
4195 NumPy method kept for backwards compatibility.
4196 Takes ``i`` or ``j``, whichever is nearest.
4198 midpoint:
4199 NumPy method kept for backwards compatibility.
4200 Uses ``(i + j) / 2``.
4202 Examples
4203 --------
4204 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
4205 >>> a
4206 array([[10, 7, 4],
4207 [ 3, 2, 1]])
4208 >>> np.percentile(a, 50)
4209 3.5
4210 >>> np.percentile(a, 50, axis=0)
4211 array([6.5, 4.5, 2.5])
4212 >>> np.percentile(a, 50, axis=1)
4213 array([7., 2.])
4214 >>> np.percentile(a, 50, axis=1, keepdims=True)
4215 array([[7.],
4216 [2.]])
4218 >>> m = np.percentile(a, 50, axis=0)
4219 >>> out = np.zeros_like(m)
4220 >>> np.percentile(a, 50, axis=0, out=out)
4221 array([6.5, 4.5, 2.5])
4222 >>> m
4223 array([6.5, 4.5, 2.5])
4225 >>> b = a.copy()
4226 >>> np.percentile(b, 50, axis=1, overwrite_input=True)
4227 array([7., 2.])
4228 >>> assert not np.all(a == b)
4230 The different methods can be visualized graphically:
4232 .. plot::
4234 import matplotlib.pyplot as plt
4236 a = np.arange(4)
4237 p = np.linspace(0, 100, 6001)
4238 ax = plt.gca()
4239 lines = [
4240 ('linear', '-', 'C0'),
4241 ('inverted_cdf', ':', 'C1'),
4242 # Almost the same as `inverted_cdf`:
4243 ('averaged_inverted_cdf', '-.', 'C1'),
4244 ('closest_observation', ':', 'C2'),
4245 ('interpolated_inverted_cdf', '--', 'C1'),
4246 ('hazen', '--', 'C3'),
4247 ('weibull', '-.', 'C4'),
4248 ('median_unbiased', '--', 'C5'),
4249 ('normal_unbiased', '-.', 'C6'),
4250 ]
4251 for method, style, color in lines:
4252 ax.plot(
4253 p, np.percentile(a, p, method=method),
4254 label=method, linestyle=style, color=color)
4255 ax.set(
4256 title='Percentiles for different methods and data: ' + str(a),
4257 xlabel='Percentile',
4258 ylabel='Estimated percentile value',
4259 yticks=a)
4260 ax.legend(bbox_to_anchor=(1.03, 1))
4261 plt.tight_layout()
4262 plt.show()
4264 References
4265 ----------
4266 .. [1] R. J. Hyndman and Y. Fan,
4267 "Sample quantiles in statistical packages,"
4268 The American Statistician, 50(4), pp. 361-365, 1996
4270 """
4271 if interpolation is not None:
4272 method = _check_interpolation_as_method(
4273 method, interpolation, "percentile")
4275 a = np.asanyarray(a)
4276 if a.dtype.kind == "c":
4277 raise TypeError("a must be an array of real numbers")
4279 q = np.true_divide(q, 100)
4280 q = asanyarray(q) # undo any decay that the ufunc performed (see gh-13105)
4281 if not _quantile_is_valid(q):
4282 raise ValueError("Percentiles must be in the range [0, 100]")
4283 return _quantile_unchecked(
4284 a, q, axis, out, overwrite_input, method, keepdims)
4287def _quantile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
4288 method=None, keepdims=None, *, interpolation=None):
4289 return (a, q, out)
4292@array_function_dispatch(_quantile_dispatcher)
4293def quantile(a,
4294 q,
4295 axis=None,
4296 out=None,
4297 overwrite_input=False,
4298 method="linear",
4299 keepdims=False,
4300 *,
4301 interpolation=None):
4302 """
4303 Compute the q-th quantile of the data along the specified axis.
4305 .. versionadded:: 1.15.0
4307 Parameters
4308 ----------
4309 a : array_like of real numbers
4310 Input array or object that can be converted to an array.
4311 q : array_like of float
4312 Probability or sequence of probabilities for the quantiles to compute.
4313 Values must be between 0 and 1 inclusive.
4314 axis : {int, tuple of int, None}, optional
4315 Axis or axes along which the quantiles are computed. The default is
4316 to compute the quantile(s) along a flattened version of the array.
4317 out : ndarray, optional
4318 Alternative output array in which to place the result. It must have
4319 the same shape and buffer length as the expected output, but the
4320 type (of the output) will be cast if necessary.
4321 overwrite_input : bool, optional
4322 If True, then allow the input array `a` to be modified by
4323 intermediate calculations, to save memory. In this case, the
4324 contents of the input `a` after this function completes is
4325 undefined.
4326 method : str, optional
4327 This parameter specifies the method to use for estimating the
4328 quantile. There are many different methods, some unique to NumPy.
4329 See the notes for explanation. The options sorted by their R type
4330 as summarized in the H&F paper [1]_ are:
4332 1. 'inverted_cdf'
4333 2. 'averaged_inverted_cdf'
4334 3. 'closest_observation'
4335 4. 'interpolated_inverted_cdf'
4336 5. 'hazen'
4337 6. 'weibull'
4338 7. 'linear' (default)
4339 8. 'median_unbiased'
4340 9. 'normal_unbiased'
4342 The first three methods are discontinuous. NumPy further defines the
4343 following discontinuous variations of the default 'linear' (7.) option:
4345 * 'lower'
4346 * 'higher',
4347 * 'midpoint'
4348 * 'nearest'
4350 .. versionchanged:: 1.22.0
4351 This argument was previously called "interpolation" and only
4352 offered the "linear" default and last four options.
4354 keepdims : bool, optional
4355 If this is set to True, the axes which are reduced are left in
4356 the result as dimensions with size one. With this option, the
4357 result will broadcast correctly against the original array `a`.
4359 interpolation : str, optional
4360 Deprecated name for the method keyword argument.
4362 .. deprecated:: 1.22.0
4364 Returns
4365 -------
4366 quantile : scalar or ndarray
4367 If `q` is a single probability and `axis=None`, then the result
4368 is a scalar. If multiple probabilies levels are given, first axis of
4369 the result corresponds to the quantiles. The other axes are
4370 the axes that remain after the reduction of `a`. If the input
4371 contains integers or floats smaller than ``float64``, the output
4372 data-type is ``float64``. Otherwise, the output data-type is the
4373 same as that of the input. If `out` is specified, that array is
4374 returned instead.
4376 See Also
4377 --------
4378 mean
4379 percentile : equivalent to quantile, but with q in the range [0, 100].
4380 median : equivalent to ``quantile(..., 0.5)``
4381 nanquantile
4383 Notes
4384 -----
4385 Given a vector ``V`` of length ``n``, the q-th quantile of ``V`` is
4386 the value ``q`` of the way from the minimum to the maximum in a
4387 sorted copy of ``V``. The values and distances of the two nearest
4388 neighbors as well as the `method` parameter will determine the
4389 quantile if the normalized ranking does not match the location of
4390 ``q`` exactly. This function is the same as the median if ``q=0.5``, the
4391 same as the minimum if ``q=0.0`` and the same as the maximum if
4392 ``q=1.0``.
4394 The optional `method` parameter specifies the method to use when the
4395 desired quantile lies between two indexes ``i`` and ``j = i + 1``.
4396 In that case, we first determine ``i + g``, a virtual index that lies
4397 between ``i`` and ``j``, where ``i`` is the floor and ``g`` is the
4398 fractional part of the index. The final result is, then, an interpolation
4399 of ``a[i]`` and ``a[j]`` based on ``g``. During the computation of ``g``,
4400 ``i`` and ``j`` are modified using correction constants ``alpha`` and
4401 ``beta`` whose choices depend on the ``method`` used. Finally, note that
4402 since Python uses 0-based indexing, the code subtracts another 1 from the
4403 index internally.
4405 The following formula determines the virtual index ``i + g``, the location
4406 of the quantile in the sorted sample:
4408 .. math::
4409 i + g = q * ( n - alpha - beta + 1 ) + alpha
4411 The different methods then work as follows
4413 inverted_cdf:
4414 method 1 of H&F [1]_.
4415 This method gives discontinuous results:
4417 * if g > 0 ; then take j
4418 * if g = 0 ; then take i
4420 averaged_inverted_cdf:
4421 method 2 of H&F [1]_.
4422 This method gives discontinuous results:
4424 * if g > 0 ; then take j
4425 * if g = 0 ; then average between bounds
4427 closest_observation:
4428 method 3 of H&F [1]_.
4429 This method gives discontinuous results:
4431 * if g > 0 ; then take j
4432 * if g = 0 and index is odd ; then take j
4433 * if g = 0 and index is even ; then take i
4435 interpolated_inverted_cdf:
4436 method 4 of H&F [1]_.
4437 This method gives continuous results using:
4439 * alpha = 0
4440 * beta = 1
4442 hazen:
4443 method 5 of H&F [1]_.
4444 This method gives continuous results using:
4446 * alpha = 1/2
4447 * beta = 1/2
4449 weibull:
4450 method 6 of H&F [1]_.
4451 This method gives continuous results using:
4453 * alpha = 0
4454 * beta = 0
4456 linear:
4457 method 7 of H&F [1]_.
4458 This method gives continuous results using:
4460 * alpha = 1
4461 * beta = 1
4463 median_unbiased:
4464 method 8 of H&F [1]_.
4465 This method is probably the best method if the sample
4466 distribution function is unknown (see reference).
4467 This method gives continuous results using:
4469 * alpha = 1/3
4470 * beta = 1/3
4472 normal_unbiased:
4473 method 9 of H&F [1]_.
4474 This method is probably the best method if the sample
4475 distribution function is known to be normal.
4476 This method gives continuous results using:
4478 * alpha = 3/8
4479 * beta = 3/8
4481 lower:
4482 NumPy method kept for backwards compatibility.
4483 Takes ``i`` as the interpolation point.
4485 higher:
4486 NumPy method kept for backwards compatibility.
4487 Takes ``j`` as the interpolation point.
4489 nearest:
4490 NumPy method kept for backwards compatibility.
4491 Takes ``i`` or ``j``, whichever is nearest.
4493 midpoint:
4494 NumPy method kept for backwards compatibility.
4495 Uses ``(i + j) / 2``.
4497 Examples
4498 --------
4499 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
4500 >>> a
4501 array([[10, 7, 4],
4502 [ 3, 2, 1]])
4503 >>> np.quantile(a, 0.5)
4504 3.5
4505 >>> np.quantile(a, 0.5, axis=0)
4506 array([6.5, 4.5, 2.5])
4507 >>> np.quantile(a, 0.5, axis=1)
4508 array([7., 2.])
4509 >>> np.quantile(a, 0.5, axis=1, keepdims=True)
4510 array([[7.],
4511 [2.]])
4512 >>> m = np.quantile(a, 0.5, axis=0)
4513 >>> out = np.zeros_like(m)
4514 >>> np.quantile(a, 0.5, axis=0, out=out)
4515 array([6.5, 4.5, 2.5])
4516 >>> m
4517 array([6.5, 4.5, 2.5])
4518 >>> b = a.copy()
4519 >>> np.quantile(b, 0.5, axis=1, overwrite_input=True)
4520 array([7., 2.])
4521 >>> assert not np.all(a == b)
4523 See also `numpy.percentile` for a visualization of most methods.
4525 References
4526 ----------
4527 .. [1] R. J. Hyndman and Y. Fan,
4528 "Sample quantiles in statistical packages,"
4529 The American Statistician, 50(4), pp. 361-365, 1996
4531 """
4532 if interpolation is not None:
4533 method = _check_interpolation_as_method(
4534 method, interpolation, "quantile")
4536 a = np.asanyarray(a)
4537 if a.dtype.kind == "c":
4538 raise TypeError("a must be an array of real numbers")
4540 q = np.asanyarray(q)
4541 if not _quantile_is_valid(q):
4542 raise ValueError("Quantiles must be in the range [0, 1]")
4543 return _quantile_unchecked(
4544 a, q, axis, out, overwrite_input, method, keepdims)
4547def _quantile_unchecked(a,
4548 q,
4549 axis=None,
4550 out=None,
4551 overwrite_input=False,
4552 method="linear",
4553 keepdims=False):
4554 """Assumes that q is in [0, 1], and is an ndarray"""
4555 return _ureduce(a,
4556 func=_quantile_ureduce_func,
4557 q=q,
4558 keepdims=keepdims,
4559 axis=axis,
4560 out=out,
4561 overwrite_input=overwrite_input,
4562 method=method)
4565def _quantile_is_valid(q):
4566 # avoid expensive reductions, relevant for arrays with < O(1000) elements
4567 if q.ndim == 1 and q.size < 10:
4568 for i in range(q.size):
4569 if not (0.0 <= q[i] <= 1.0):
4570 return False
4571 else:
4572 if not (np.all(0 <= q) and np.all(q <= 1)):
4573 return False
4574 return True
4577def _check_interpolation_as_method(method, interpolation, fname):
4578 # Deprecated NumPy 1.22, 2021-11-08
4579 warnings.warn(
4580 f"the `interpolation=` argument to {fname} was renamed to "
4581 "`method=`, which has additional options.\n"
4582 "Users of the modes 'nearest', 'lower', 'higher', or "
4583 "'midpoint' are encouraged to review the method they used. "
4584 "(Deprecated NumPy 1.22)",
4585 DeprecationWarning, stacklevel=4)
4586 if method != "linear":
4587 # sanity check, we assume this basically never happens
4588 raise TypeError(
4589 "You shall not pass both `method` and `interpolation`!\n"
4590 "(`interpolation` is Deprecated in favor of `method`)")
4591 return interpolation
4594def _compute_virtual_index(n, quantiles, alpha: float, beta: float):
4595 """
4596 Compute the floating point indexes of an array for the linear
4597 interpolation of quantiles.
4598 n : array_like
4599 The sample sizes.
4600 quantiles : array_like
4601 The quantiles values.
4602 alpha : float
4603 A constant used to correct the index computed.
4604 beta : float
4605 A constant used to correct the index computed.
4607 alpha and beta values depend on the chosen method
4608 (see quantile documentation)
4610 Reference:
4611 Hyndman&Fan paper "Sample Quantiles in Statistical Packages",
4612 DOI: 10.1080/00031305.1996.10473566
4613 """
4614 return n * quantiles + (
4615 alpha + quantiles * (1 - alpha - beta)
4616 ) - 1
4619def _get_gamma(virtual_indexes, previous_indexes, method):
4620 """
4621 Compute gamma (a.k.a 'm' or 'weight') for the linear interpolation
4622 of quantiles.
4624 virtual_indexes : array_like
4625 The indexes where the percentile is supposed to be found in the sorted
4626 sample.
4627 previous_indexes : array_like
4628 The floor values of virtual_indexes.
4629 interpolation : dict
4630 The interpolation method chosen, which may have a specific rule
4631 modifying gamma.
4633 gamma is usually the fractional part of virtual_indexes but can be modified
4634 by the interpolation method.
4635 """
4636 gamma = np.asanyarray(virtual_indexes - previous_indexes)
4637 gamma = method["fix_gamma"](gamma, virtual_indexes)
4638 return np.asanyarray(gamma)
4641def _lerp(a, b, t, out=None):
4642 """
4643 Compute the linear interpolation weighted by gamma on each point of
4644 two same shape array.
4646 a : array_like
4647 Left bound.
4648 b : array_like
4649 Right bound.
4650 t : array_like
4651 The interpolation weight.
4652 out : array_like
4653 Output array.
4654 """
4655 diff_b_a = subtract(b, a)
4656 # asanyarray is a stop-gap until gh-13105
4657 lerp_interpolation = asanyarray(add(a, diff_b_a * t, out=out))
4658 subtract(b, diff_b_a * (1 - t), out=lerp_interpolation, where=t >= 0.5,
4659 casting='unsafe', dtype=type(lerp_interpolation.dtype))
4660 if lerp_interpolation.ndim == 0 and out is None:
4661 lerp_interpolation = lerp_interpolation[()] # unpack 0d arrays
4662 return lerp_interpolation
4665def _get_gamma_mask(shape, default_value, conditioned_value, where):
4666 out = np.full(shape, default_value)
4667 np.copyto(out, conditioned_value, where=where, casting="unsafe")
4668 return out
4671def _discret_interpolation_to_boundaries(index, gamma_condition_fun):
4672 previous = np.floor(index)
4673 next = previous + 1
4674 gamma = index - previous
4675 res = _get_gamma_mask(shape=index.shape,
4676 default_value=next,
4677 conditioned_value=previous,
4678 where=gamma_condition_fun(gamma, index)
4679 ).astype(np.intp)
4680 # Some methods can lead to out-of-bound integers, clip them:
4681 res[res < 0] = 0
4682 return res
4685def _closest_observation(n, quantiles):
4686 gamma_fun = lambda gamma, index: (gamma == 0) & (np.floor(index) % 2 == 0)
4687 return _discret_interpolation_to_boundaries((n * quantiles) - 1 - 0.5,
4688 gamma_fun)
4691def _inverted_cdf(n, quantiles):
4692 gamma_fun = lambda gamma, _: (gamma == 0)
4693 return _discret_interpolation_to_boundaries((n * quantiles) - 1,
4694 gamma_fun)
4697def _quantile_ureduce_func(
4698 a: np.array,
4699 q: np.array,
4700 axis: int = None,
4701 out=None,
4702 overwrite_input: bool = False,
4703 method="linear",
4704) -> np.array:
4705 if q.ndim > 2:
4706 # The code below works fine for nd, but it might not have useful
4707 # semantics. For now, keep the supported dimensions the same as it was
4708 # before.
4709 raise ValueError("q must be a scalar or 1d")
4710 if overwrite_input:
4711 if axis is None:
4712 axis = 0
4713 arr = a.ravel()
4714 else:
4715 arr = a
4716 else:
4717 if axis is None:
4718 axis = 0
4719 arr = a.flatten()
4720 else:
4721 arr = a.copy()
4722 result = _quantile(arr,
4723 quantiles=q,
4724 axis=axis,
4725 method=method,
4726 out=out)
4727 return result
4730def _get_indexes(arr, virtual_indexes, valid_values_count):
4731 """
4732 Get the valid indexes of arr neighbouring virtual_indexes.
4733 Note
4734 This is a companion function to linear interpolation of
4735 Quantiles
4737 Returns
4738 -------
4739 (previous_indexes, next_indexes): Tuple
4740 A Tuple of virtual_indexes neighbouring indexes
4741 """
4742 previous_indexes = np.asanyarray(np.floor(virtual_indexes))
4743 next_indexes = np.asanyarray(previous_indexes + 1)
4744 indexes_above_bounds = virtual_indexes >= valid_values_count - 1
4745 # When indexes is above max index, take the max value of the array
4746 if indexes_above_bounds.any():
4747 previous_indexes[indexes_above_bounds] = -1
4748 next_indexes[indexes_above_bounds] = -1
4749 # When indexes is below min index, take the min value of the array
4750 indexes_below_bounds = virtual_indexes < 0
4751 if indexes_below_bounds.any():
4752 previous_indexes[indexes_below_bounds] = 0
4753 next_indexes[indexes_below_bounds] = 0
4754 if np.issubdtype(arr.dtype, np.inexact):
4755 # After the sort, slices having NaNs will have for last element a NaN
4756 virtual_indexes_nans = np.isnan(virtual_indexes)
4757 if virtual_indexes_nans.any():
4758 previous_indexes[virtual_indexes_nans] = -1
4759 next_indexes[virtual_indexes_nans] = -1
4760 previous_indexes = previous_indexes.astype(np.intp)
4761 next_indexes = next_indexes.astype(np.intp)
4762 return previous_indexes, next_indexes
4765def _quantile(
4766 arr: np.array,
4767 quantiles: np.array,
4768 axis: int = -1,
4769 method="linear",
4770 out=None,
4771):
4772 """
4773 Private function that doesn't support extended axis or keepdims.
4774 These methods are extended to this function using _ureduce
4775 See nanpercentile for parameter usage
4776 It computes the quantiles of the array for the given axis.
4777 A linear interpolation is performed based on the `interpolation`.
4779 By default, the method is "linear" where alpha == beta == 1 which
4780 performs the 7th method of Hyndman&Fan.
4781 With "median_unbiased" we get alpha == beta == 1/3
4782 thus the 8th method of Hyndman&Fan.
4783 """
4784 # --- Setup
4785 arr = np.asanyarray(arr)
4786 values_count = arr.shape[axis]
4787 # The dimensions of `q` are prepended to the output shape, so we need the
4788 # axis being sampled from `arr` to be last.
4790 if axis != 0: # But moveaxis is slow, so only call it if necessary.
4791 arr = np.moveaxis(arr, axis, destination=0)
4792 # --- Computation of indexes
4793 # Index where to find the value in the sorted array.
4794 # Virtual because it is a floating point value, not an valid index.
4795 # The nearest neighbours are used for interpolation
4796 try:
4797 method = _QuantileMethods[method]
4798 except KeyError:
4799 raise ValueError(
4800 f"{method!r} is not a valid method. Use one of: "
4801 f"{_QuantileMethods.keys()}") from None
4802 virtual_indexes = method["get_virtual_index"](values_count, quantiles)
4803 virtual_indexes = np.asanyarray(virtual_indexes)
4805 supports_nans = (
4806 np.issubdtype(arr.dtype, np.inexact) or arr.dtype.kind in 'Mm')
4808 if np.issubdtype(virtual_indexes.dtype, np.integer):
4809 # No interpolation needed, take the points along axis
4810 if supports_nans:
4811 # may contain nan, which would sort to the end
4812 arr.partition(concatenate((virtual_indexes.ravel(), [-1])), axis=0)
4813 slices_having_nans = np.isnan(arr[-1, ...])
4814 else:
4815 # cannot contain nan
4816 arr.partition(virtual_indexes.ravel(), axis=0)
4817 slices_having_nans = np.array(False, dtype=bool)
4818 result = take(arr, virtual_indexes, axis=0, out=out)
4819 else:
4820 previous_indexes, next_indexes = _get_indexes(arr,
4821 virtual_indexes,
4822 values_count)
4823 # --- Sorting
4824 arr.partition(
4825 np.unique(np.concatenate(([0, -1],
4826 previous_indexes.ravel(),
4827 next_indexes.ravel(),
4828 ))),
4829 axis=0)
4830 if supports_nans:
4831 slices_having_nans = np.isnan(arr[-1, ...])
4832 else:
4833 slices_having_nans = None
4834 # --- Get values from indexes
4835 previous = arr[previous_indexes]
4836 next = arr[next_indexes]
4837 # --- Linear interpolation
4838 gamma = _get_gamma(virtual_indexes, previous_indexes, method)
4839 result_shape = virtual_indexes.shape + (1,) * (arr.ndim - 1)
4840 gamma = gamma.reshape(result_shape)
4841 result = _lerp(previous,
4842 next,
4843 gamma,
4844 out=out)
4845 if np.any(slices_having_nans):
4846 if result.ndim == 0 and out is None:
4847 # can't write to a scalar, but indexing will be correct
4848 result = arr[-1]
4849 else:
4850 np.copyto(result, arr[-1, ...], where=slices_having_nans)
4851 return result
4854def _trapz_dispatcher(y, x=None, dx=None, axis=None):
4855 return (y, x)
4858@array_function_dispatch(_trapz_dispatcher)
4859def trapz(y, x=None, dx=1.0, axis=-1):
4860 r"""
4861 Integrate along the given axis using the composite trapezoidal rule.
4863 If `x` is provided, the integration happens in sequence along its
4864 elements - they are not sorted.
4866 Integrate `y` (`x`) along each 1d slice on the given axis, compute
4867 :math:`\int y(x) dx`.
4868 When `x` is specified, this integrates along the parametric curve,
4869 computing :math:`\int_t y(t) dt =
4870 \int_t y(t) \left.\frac{dx}{dt}\right|_{x=x(t)} dt`.
4872 Parameters
4873 ----------
4874 y : array_like
4875 Input array to integrate.
4876 x : array_like, optional
4877 The sample points corresponding to the `y` values. If `x` is None,
4878 the sample points are assumed to be evenly spaced `dx` apart. The
4879 default is None.
4880 dx : scalar, optional
4881 The spacing between sample points when `x` is None. The default is 1.
4882 axis : int, optional
4883 The axis along which to integrate.
4885 Returns
4886 -------
4887 trapz : float or ndarray
4888 Definite integral of `y` = n-dimensional array as approximated along
4889 a single axis by the trapezoidal rule. If `y` is a 1-dimensional array,
4890 then the result is a float. If `n` is greater than 1, then the result
4891 is an `n`-1 dimensional array.
4893 See Also
4894 --------
4895 sum, cumsum
4897 Notes
4898 -----
4899 Image [2]_ illustrates trapezoidal rule -- y-axis locations of points
4900 will be taken from `y` array, by default x-axis distances between
4901 points will be 1.0, alternatively they can be provided with `x` array
4902 or with `dx` scalar. Return value will be equal to combined area under
4903 the red lines.
4906 References
4907 ----------
4908 .. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule
4910 .. [2] Illustration image:
4911 https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
4913 Examples
4914 --------
4915 Use the trapezoidal rule on evenly spaced points:
4917 >>> np.trapz([1, 2, 3])
4918 4.0
4920 The spacing between sample points can be selected by either the
4921 ``x`` or ``dx`` arguments:
4923 >>> np.trapz([1, 2, 3], x=[4, 6, 8])
4924 8.0
4925 >>> np.trapz([1, 2, 3], dx=2)
4926 8.0
4928 Using a decreasing ``x`` corresponds to integrating in reverse:
4930 >>> np.trapz([1, 2, 3], x=[8, 6, 4])
4931 -8.0
4933 More generally ``x`` is used to integrate along a parametric curve. We can
4934 estimate the integral :math:`\int_0^1 x^2 = 1/3` using:
4936 >>> x = np.linspace(0, 1, num=50)
4937 >>> y = x**2
4938 >>> np.trapz(y, x)
4939 0.33340274885464394
4941 Or estimate the area of a circle, noting we repeat the sample which closes
4942 the curve:
4944 >>> theta = np.linspace(0, 2 * np.pi, num=1000, endpoint=True)
4945 >>> np.trapz(np.cos(theta), x=np.sin(theta))
4946 3.141571941375841
4948 ``np.trapz`` can be applied along a specified axis to do multiple
4949 computations in one call:
4951 >>> a = np.arange(6).reshape(2, 3)
4952 >>> a
4953 array([[0, 1, 2],
4954 [3, 4, 5]])
4955 >>> np.trapz(a, axis=0)
4956 array([1.5, 2.5, 3.5])
4957 >>> np.trapz(a, axis=1)
4958 array([2., 8.])
4959 """
4960 y = asanyarray(y)
4961 if x is None:
4962 d = dx
4963 else:
4964 x = asanyarray(x)
4965 if x.ndim == 1:
4966 d = diff(x)
4967 # reshape to correct shape
4968 shape = [1]*y.ndim
4969 shape[axis] = d.shape[0]
4970 d = d.reshape(shape)
4971 else:
4972 d = diff(x, axis=axis)
4973 nd = y.ndim
4974 slice1 = [slice(None)]*nd
4975 slice2 = [slice(None)]*nd
4976 slice1[axis] = slice(1, None)
4977 slice2[axis] = slice(None, -1)
4978 try:
4979 ret = (d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0).sum(axis)
4980 except ValueError:
4981 # Operations didn't work, cast to ndarray
4982 d = np.asarray(d)
4983 y = np.asarray(y)
4984 ret = add.reduce(d * (y[tuple(slice1)]+y[tuple(slice2)])/2.0, axis)
4985 return ret
4988# __array_function__ has no __code__ or other attributes normal Python funcs we
4989# wrap everything into a C callable. SciPy however, tries to "clone" `trapz`
4990# into a new Python function which requires `__code__` and a few other
4991# attributes. So we create a dummy clone and copy over its attributes allowing
4992# SciPy <= 1.10 to work: https://github.com/scipy/scipy/issues/17811
4993assert not hasattr(trapz, "__code__")
4995def _fake_trapz(y, x=None, dx=1.0, axis=-1):
4996 return trapz(y, x=x, dx=dx, axis=axis)
4999trapz.__code__ = _fake_trapz.__code__
5000trapz.__globals__ = _fake_trapz.__globals__
5001trapz.__defaults__ = _fake_trapz.__defaults__
5002trapz.__closure__ = _fake_trapz.__closure__
5003trapz.__kwdefaults__ = _fake_trapz.__kwdefaults__
5006def _meshgrid_dispatcher(*xi, copy=None, sparse=None, indexing=None):
5007 return xi
5010# Based on scitools meshgrid
5011@array_function_dispatch(_meshgrid_dispatcher)
5012def meshgrid(*xi, copy=True, sparse=False, indexing='xy'):
5013 """
5014 Return a list of coordinate matrices from coordinate vectors.
5016 Make N-D coordinate arrays for vectorized evaluations of
5017 N-D scalar/vector fields over N-D grids, given
5018 one-dimensional coordinate arrays x1, x2,..., xn.
5020 .. versionchanged:: 1.9
5021 1-D and 0-D cases are allowed.
5023 Parameters
5024 ----------
5025 x1, x2,..., xn : array_like
5026 1-D arrays representing the coordinates of a grid.
5027 indexing : {'xy', 'ij'}, optional
5028 Cartesian ('xy', default) or matrix ('ij') indexing of output.
5029 See Notes for more details.
5031 .. versionadded:: 1.7.0
5032 sparse : bool, optional
5033 If True the shape of the returned coordinate array for dimension *i*
5034 is reduced from ``(N1, ..., Ni, ... Nn)`` to
5035 ``(1, ..., 1, Ni, 1, ..., 1)``. These sparse coordinate grids are
5036 intended to be use with :ref:`basics.broadcasting`. When all
5037 coordinates are used in an expression, broadcasting still leads to a
5038 fully-dimensonal result array.
5040 Default is False.
5042 .. versionadded:: 1.7.0
5043 copy : bool, optional
5044 If False, a view into the original arrays are returned in order to
5045 conserve memory. Default is True. Please note that
5046 ``sparse=False, copy=False`` will likely return non-contiguous
5047 arrays. Furthermore, more than one element of a broadcast array
5048 may refer to a single memory location. If you need to write to the
5049 arrays, make copies first.
5051 .. versionadded:: 1.7.0
5053 Returns
5054 -------
5055 X1, X2,..., XN : list of ndarrays
5056 For vectors `x1`, `x2`,..., `xn` with lengths ``Ni=len(xi)``,
5057 returns ``(N1, N2, N3,..., Nn)`` shaped arrays if indexing='ij'
5058 or ``(N2, N1, N3,..., Nn)`` shaped arrays if indexing='xy'
5059 with the elements of `xi` repeated to fill the matrix along
5060 the first dimension for `x1`, the second for `x2` and so on.
5062 Notes
5063 -----
5064 This function supports both indexing conventions through the indexing
5065 keyword argument. Giving the string 'ij' returns a meshgrid with
5066 matrix indexing, while 'xy' returns a meshgrid with Cartesian indexing.
5067 In the 2-D case with inputs of length M and N, the outputs are of shape
5068 (N, M) for 'xy' indexing and (M, N) for 'ij' indexing. In the 3-D case
5069 with inputs of length M, N and P, outputs are of shape (N, M, P) for
5070 'xy' indexing and (M, N, P) for 'ij' indexing. The difference is
5071 illustrated by the following code snippet::
5073 xv, yv = np.meshgrid(x, y, indexing='ij')
5074 for i in range(nx):
5075 for j in range(ny):
5076 # treat xv[i,j], yv[i,j]
5078 xv, yv = np.meshgrid(x, y, indexing='xy')
5079 for i in range(nx):
5080 for j in range(ny):
5081 # treat xv[j,i], yv[j,i]
5083 In the 1-D and 0-D case, the indexing and sparse keywords have no effect.
5085 See Also
5086 --------
5087 mgrid : Construct a multi-dimensional "meshgrid" using indexing notation.
5088 ogrid : Construct an open multi-dimensional "meshgrid" using indexing
5089 notation.
5090 how-to-index
5092 Examples
5093 --------
5094 >>> nx, ny = (3, 2)
5095 >>> x = np.linspace(0, 1, nx)
5096 >>> y = np.linspace(0, 1, ny)
5097 >>> xv, yv = np.meshgrid(x, y)
5098 >>> xv
5099 array([[0. , 0.5, 1. ],
5100 [0. , 0.5, 1. ]])
5101 >>> yv
5102 array([[0., 0., 0.],
5103 [1., 1., 1.]])
5105 The result of `meshgrid` is a coordinate grid:
5107 >>> import matplotlib.pyplot as plt
5108 >>> plt.plot(xv, yv, marker='o', color='k', linestyle='none')
5109 >>> plt.show()
5111 You can create sparse output arrays to save memory and computation time.
5113 >>> xv, yv = np.meshgrid(x, y, sparse=True)
5114 >>> xv
5115 array([[0. , 0.5, 1. ]])
5116 >>> yv
5117 array([[0.],
5118 [1.]])
5120 `meshgrid` is very useful to evaluate functions on a grid. If the
5121 function depends on all coordinates, both dense and sparse outputs can be
5122 used.
5124 >>> x = np.linspace(-5, 5, 101)
5125 >>> y = np.linspace(-5, 5, 101)
5126 >>> # full coordinate arrays
5127 >>> xx, yy = np.meshgrid(x, y)
5128 >>> zz = np.sqrt(xx**2 + yy**2)
5129 >>> xx.shape, yy.shape, zz.shape
5130 ((101, 101), (101, 101), (101, 101))
5131 >>> # sparse coordinate arrays
5132 >>> xs, ys = np.meshgrid(x, y, sparse=True)
5133 >>> zs = np.sqrt(xs**2 + ys**2)
5134 >>> xs.shape, ys.shape, zs.shape
5135 ((1, 101), (101, 1), (101, 101))
5136 >>> np.array_equal(zz, zs)
5137 True
5139 >>> h = plt.contourf(x, y, zs)
5140 >>> plt.axis('scaled')
5141 >>> plt.colorbar()
5142 >>> plt.show()
5143 """
5144 ndim = len(xi)
5146 if indexing not in ['xy', 'ij']:
5147 raise ValueError(
5148 "Valid values for `indexing` are 'xy' and 'ij'.")
5150 s0 = (1,) * ndim
5151 output = [np.asanyarray(x).reshape(s0[:i] + (-1,) + s0[i + 1:])
5152 for i, x in enumerate(xi)]
5154 if indexing == 'xy' and ndim > 1:
5155 # switch first and second axis
5156 output[0].shape = (1, -1) + s0[2:]
5157 output[1].shape = (-1, 1) + s0[2:]
5159 if not sparse:
5160 # Return the full N-D matrix (not only the 1-D vector)
5161 output = np.broadcast_arrays(*output, subok=True)
5163 if copy:
5164 output = [x.copy() for x in output]
5166 return output
5169def _delete_dispatcher(arr, obj, axis=None):
5170 return (arr, obj)
5173@array_function_dispatch(_delete_dispatcher)
5174def delete(arr, obj, axis=None):
5175 """
5176 Return a new array with sub-arrays along an axis deleted. For a one
5177 dimensional array, this returns those entries not returned by
5178 `arr[obj]`.
5180 Parameters
5181 ----------
5182 arr : array_like
5183 Input array.
5184 obj : slice, int or array of ints
5185 Indicate indices of sub-arrays to remove along the specified axis.
5187 .. versionchanged:: 1.19.0
5188 Boolean indices are now treated as a mask of elements to remove,
5189 rather than being cast to the integers 0 and 1.
5191 axis : int, optional
5192 The axis along which to delete the subarray defined by `obj`.
5193 If `axis` is None, `obj` is applied to the flattened array.
5195 Returns
5196 -------
5197 out : ndarray
5198 A copy of `arr` with the elements specified by `obj` removed. Note
5199 that `delete` does not occur in-place. If `axis` is None, `out` is
5200 a flattened array.
5202 See Also
5203 --------
5204 insert : Insert elements into an array.
5205 append : Append elements at the end of an array.
5207 Notes
5208 -----
5209 Often it is preferable to use a boolean mask. For example:
5211 >>> arr = np.arange(12) + 1
5212 >>> mask = np.ones(len(arr), dtype=bool)
5213 >>> mask[[0,2,4]] = False
5214 >>> result = arr[mask,...]
5216 Is equivalent to ``np.delete(arr, [0,2,4], axis=0)``, but allows further
5217 use of `mask`.
5219 Examples
5220 --------
5221 >>> arr = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]])
5222 >>> arr
5223 array([[ 1, 2, 3, 4],
5224 [ 5, 6, 7, 8],
5225 [ 9, 10, 11, 12]])
5226 >>> np.delete(arr, 1, 0)
5227 array([[ 1, 2, 3, 4],
5228 [ 9, 10, 11, 12]])
5230 >>> np.delete(arr, np.s_[::2], 1)
5231 array([[ 2, 4],
5232 [ 6, 8],
5233 [10, 12]])
5234 >>> np.delete(arr, [1,3,5], None)
5235 array([ 1, 3, 5, 7, 8, 9, 10, 11, 12])
5237 """
5238 wrap = None
5239 if type(arr) is not ndarray:
5240 try:
5241 wrap = arr.__array_wrap__
5242 except AttributeError:
5243 pass
5245 arr = asarray(arr)
5246 ndim = arr.ndim
5247 arrorder = 'F' if arr.flags.fnc else 'C'
5248 if axis is None:
5249 if ndim != 1:
5250 arr = arr.ravel()
5251 # needed for np.matrix, which is still not 1d after being ravelled
5252 ndim = arr.ndim
5253 axis = ndim - 1
5254 else:
5255 axis = normalize_axis_index(axis, ndim)
5257 slobj = [slice(None)]*ndim
5258 N = arr.shape[axis]
5259 newshape = list(arr.shape)
5261 if isinstance(obj, slice):
5262 start, stop, step = obj.indices(N)
5263 xr = range(start, stop, step)
5264 numtodel = len(xr)
5266 if numtodel <= 0:
5267 if wrap:
5268 return wrap(arr.copy(order=arrorder))
5269 else:
5270 return arr.copy(order=arrorder)
5272 # Invert if step is negative:
5273 if step < 0:
5274 step = -step
5275 start = xr[-1]
5276 stop = xr[0] + 1
5278 newshape[axis] -= numtodel
5279 new = empty(newshape, arr.dtype, arrorder)
5280 # copy initial chunk
5281 if start == 0:
5282 pass
5283 else:
5284 slobj[axis] = slice(None, start)
5285 new[tuple(slobj)] = arr[tuple(slobj)]
5286 # copy end chunk
5287 if stop == N:
5288 pass
5289 else:
5290 slobj[axis] = slice(stop-numtodel, None)
5291 slobj2 = [slice(None)]*ndim
5292 slobj2[axis] = slice(stop, None)
5293 new[tuple(slobj)] = arr[tuple(slobj2)]
5294 # copy middle pieces
5295 if step == 1:
5296 pass
5297 else: # use array indexing.
5298 keep = ones(stop-start, dtype=bool)
5299 keep[:stop-start:step] = False
5300 slobj[axis] = slice(start, stop-numtodel)
5301 slobj2 = [slice(None)]*ndim
5302 slobj2[axis] = slice(start, stop)
5303 arr = arr[tuple(slobj2)]
5304 slobj2[axis] = keep
5305 new[tuple(slobj)] = arr[tuple(slobj2)]
5306 if wrap:
5307 return wrap(new)
5308 else:
5309 return new
5311 if isinstance(obj, (int, integer)) and not isinstance(obj, bool):
5312 single_value = True
5313 else:
5314 single_value = False
5315 _obj = obj
5316 obj = np.asarray(obj)
5317 # `size == 0` to allow empty lists similar to indexing, but (as there)
5318 # is really too generic:
5319 if obj.size == 0 and not isinstance(_obj, np.ndarray):
5320 obj = obj.astype(intp)
5321 elif obj.size == 1 and obj.dtype.kind in "ui":
5322 # For a size 1 integer array we can use the single-value path
5323 # (most dtypes, except boolean, should just fail later).
5324 obj = obj.item()
5325 single_value = True
5327 if single_value:
5328 # optimization for a single value
5329 if (obj < -N or obj >= N):
5330 raise IndexError(
5331 "index %i is out of bounds for axis %i with "
5332 "size %i" % (obj, axis, N))
5333 if (obj < 0):
5334 obj += N
5335 newshape[axis] -= 1
5336 new = empty(newshape, arr.dtype, arrorder)
5337 slobj[axis] = slice(None, obj)
5338 new[tuple(slobj)] = arr[tuple(slobj)]
5339 slobj[axis] = slice(obj, None)
5340 slobj2 = [slice(None)]*ndim
5341 slobj2[axis] = slice(obj+1, None)
5342 new[tuple(slobj)] = arr[tuple(slobj2)]
5343 else:
5344 if obj.dtype == bool:
5345 if obj.shape != (N,):
5346 raise ValueError('boolean array argument obj to delete '
5347 'must be one dimensional and match the axis '
5348 'length of {}'.format(N))
5350 # optimization, the other branch is slower
5351 keep = ~obj
5352 else:
5353 keep = ones(N, dtype=bool)
5354 keep[obj,] = False
5356 slobj[axis] = keep
5357 new = arr[tuple(slobj)]
5359 if wrap:
5360 return wrap(new)
5361 else:
5362 return new
5365def _insert_dispatcher(arr, obj, values, axis=None):
5366 return (arr, obj, values)
5369@array_function_dispatch(_insert_dispatcher)
5370def insert(arr, obj, values, axis=None):
5371 """
5372 Insert values along the given axis before the given indices.
5374 Parameters
5375 ----------
5376 arr : array_like
5377 Input array.
5378 obj : int, slice or sequence of ints
5379 Object that defines the index or indices before which `values` is
5380 inserted.
5382 .. versionadded:: 1.8.0
5384 Support for multiple insertions when `obj` is a single scalar or a
5385 sequence with one element (similar to calling insert multiple
5386 times).
5387 values : array_like
5388 Values to insert into `arr`. If the type of `values` is different
5389 from that of `arr`, `values` is converted to the type of `arr`.
5390 `values` should be shaped so that ``arr[...,obj,...] = values``
5391 is legal.
5392 axis : int, optional
5393 Axis along which to insert `values`. If `axis` is None then `arr`
5394 is flattened first.
5396 Returns
5397 -------
5398 out : ndarray
5399 A copy of `arr` with `values` inserted. Note that `insert`
5400 does not occur in-place: a new array is returned. If
5401 `axis` is None, `out` is a flattened array.
5403 See Also
5404 --------
5405 append : Append elements at the end of an array.
5406 concatenate : Join a sequence of arrays along an existing axis.
5407 delete : Delete elements from an array.
5409 Notes
5410 -----
5411 Note that for higher dimensional inserts ``obj=0`` behaves very different
5412 from ``obj=[0]`` just like ``arr[:,0,:] = values`` is different from
5413 ``arr[:,[0],:] = values``.
5415 Examples
5416 --------
5417 >>> a = np.array([[1, 1], [2, 2], [3, 3]])
5418 >>> a
5419 array([[1, 1],
5420 [2, 2],
5421 [3, 3]])
5422 >>> np.insert(a, 1, 5)
5423 array([1, 5, 1, ..., 2, 3, 3])
5424 >>> np.insert(a, 1, 5, axis=1)
5425 array([[1, 5, 1],
5426 [2, 5, 2],
5427 [3, 5, 3]])
5429 Difference between sequence and scalars:
5431 >>> np.insert(a, [1], [[1],[2],[3]], axis=1)
5432 array([[1, 1, 1],
5433 [2, 2, 2],
5434 [3, 3, 3]])
5435 >>> np.array_equal(np.insert(a, 1, [1, 2, 3], axis=1),
5436 ... np.insert(a, [1], [[1],[2],[3]], axis=1))
5437 True
5439 >>> b = a.flatten()
5440 >>> b
5441 array([1, 1, 2, 2, 3, 3])
5442 >>> np.insert(b, [2, 2], [5, 6])
5443 array([1, 1, 5, ..., 2, 3, 3])
5445 >>> np.insert(b, slice(2, 4), [5, 6])
5446 array([1, 1, 5, ..., 2, 3, 3])
5448 >>> np.insert(b, [2, 2], [7.13, False]) # type casting
5449 array([1, 1, 7, ..., 2, 3, 3])
5451 >>> x = np.arange(8).reshape(2, 4)
5452 >>> idx = (1, 3)
5453 >>> np.insert(x, idx, 999, axis=1)
5454 array([[ 0, 999, 1, 2, 999, 3],
5455 [ 4, 999, 5, 6, 999, 7]])
5457 """
5458 wrap = None
5459 if type(arr) is not ndarray:
5460 try:
5461 wrap = arr.__array_wrap__
5462 except AttributeError:
5463 pass
5465 arr = asarray(arr)
5466 ndim = arr.ndim
5467 arrorder = 'F' if arr.flags.fnc else 'C'
5468 if axis is None:
5469 if ndim != 1:
5470 arr = arr.ravel()
5471 # needed for np.matrix, which is still not 1d after being ravelled
5472 ndim = arr.ndim
5473 axis = ndim - 1
5474 else:
5475 axis = normalize_axis_index(axis, ndim)
5476 slobj = [slice(None)]*ndim
5477 N = arr.shape[axis]
5478 newshape = list(arr.shape)
5480 if isinstance(obj, slice):
5481 # turn it into a range object
5482 indices = arange(*obj.indices(N), dtype=intp)
5483 else:
5484 # need to copy obj, because indices will be changed in-place
5485 indices = np.array(obj)
5486 if indices.dtype == bool:
5487 # See also delete
5488 # 2012-10-11, NumPy 1.8
5489 warnings.warn(
5490 "in the future insert will treat boolean arrays and "
5491 "array-likes as a boolean index instead of casting it to "
5492 "integer", FutureWarning, stacklevel=2)
5493 indices = indices.astype(intp)
5494 # Code after warning period:
5495 #if obj.ndim != 1:
5496 # raise ValueError('boolean array argument obj to insert '
5497 # 'must be one dimensional')
5498 #indices = np.flatnonzero(obj)
5499 elif indices.ndim > 1:
5500 raise ValueError(
5501 "index array argument obj to insert must be one dimensional "
5502 "or scalar")
5503 if indices.size == 1:
5504 index = indices.item()
5505 if index < -N or index > N:
5506 raise IndexError(f"index {obj} is out of bounds for axis {axis} "
5507 f"with size {N}")
5508 if (index < 0):
5509 index += N
5511 # There are some object array corner cases here, but we cannot avoid
5512 # that:
5513 values = array(values, copy=False, ndmin=arr.ndim, dtype=arr.dtype)
5514 if indices.ndim == 0:
5515 # broadcasting is very different here, since a[:,0,:] = ... behaves
5516 # very different from a[:,[0],:] = ...! This changes values so that
5517 # it works likes the second case. (here a[:,0:1,:])
5518 values = np.moveaxis(values, 0, axis)
5519 numnew = values.shape[axis]
5520 newshape[axis] += numnew
5521 new = empty(newshape, arr.dtype, arrorder)
5522 slobj[axis] = slice(None, index)
5523 new[tuple(slobj)] = arr[tuple(slobj)]
5524 slobj[axis] = slice(index, index+numnew)
5525 new[tuple(slobj)] = values
5526 slobj[axis] = slice(index+numnew, None)
5527 slobj2 = [slice(None)] * ndim
5528 slobj2[axis] = slice(index, None)
5529 new[tuple(slobj)] = arr[tuple(slobj2)]
5530 if wrap:
5531 return wrap(new)
5532 return new
5533 elif indices.size == 0 and not isinstance(obj, np.ndarray):
5534 # Can safely cast the empty list to intp
5535 indices = indices.astype(intp)
5537 indices[indices < 0] += N
5539 numnew = len(indices)
5540 order = indices.argsort(kind='mergesort') # stable sort
5541 indices[order] += np.arange(numnew)
5543 newshape[axis] += numnew
5544 old_mask = ones(newshape[axis], dtype=bool)
5545 old_mask[indices] = False
5547 new = empty(newshape, arr.dtype, arrorder)
5548 slobj2 = [slice(None)]*ndim
5549 slobj[axis] = indices
5550 slobj2[axis] = old_mask
5551 new[tuple(slobj)] = values
5552 new[tuple(slobj2)] = arr
5554 if wrap:
5555 return wrap(new)
5556 return new
5559def _append_dispatcher(arr, values, axis=None):
5560 return (arr, values)
5563@array_function_dispatch(_append_dispatcher)
5564def append(arr, values, axis=None):
5565 """
5566 Append values to the end of an array.
5568 Parameters
5569 ----------
5570 arr : array_like
5571 Values are appended to a copy of this array.
5572 values : array_like
5573 These values are appended to a copy of `arr`. It must be of the
5574 correct shape (the same shape as `arr`, excluding `axis`). If
5575 `axis` is not specified, `values` can be any shape and will be
5576 flattened before use.
5577 axis : int, optional
5578 The axis along which `values` are appended. If `axis` is not
5579 given, both `arr` and `values` are flattened before use.
5581 Returns
5582 -------
5583 append : ndarray
5584 A copy of `arr` with `values` appended to `axis`. Note that
5585 `append` does not occur in-place: a new array is allocated and
5586 filled. If `axis` is None, `out` is a flattened array.
5588 See Also
5589 --------
5590 insert : Insert elements into an array.
5591 delete : Delete elements from an array.
5593 Examples
5594 --------
5595 >>> np.append([1, 2, 3], [[4, 5, 6], [7, 8, 9]])
5596 array([1, 2, 3, ..., 7, 8, 9])
5598 When `axis` is specified, `values` must have the correct shape.
5600 >>> np.append([[1, 2, 3], [4, 5, 6]], [[7, 8, 9]], axis=0)
5601 array([[1, 2, 3],
5602 [4, 5, 6],
5603 [7, 8, 9]])
5604 >>> np.append([[1, 2, 3], [4, 5, 6]], [7, 8, 9], axis=0)
5605 Traceback (most recent call last):
5606 ...
5607 ValueError: all the input arrays must have same number of dimensions, but
5608 the array at index 0 has 2 dimension(s) and the array at index 1 has 1
5609 dimension(s)
5611 """
5612 arr = asanyarray(arr)
5613 if axis is None:
5614 if arr.ndim != 1:
5615 arr = arr.ravel()
5616 values = ravel(values)
5617 axis = arr.ndim-1
5618 return concatenate((arr, values), axis=axis)
5621def _digitize_dispatcher(x, bins, right=None):
5622 return (x, bins)
5625@array_function_dispatch(_digitize_dispatcher)
5626def digitize(x, bins, right=False):
5627 """
5628 Return the indices of the bins to which each value in input array belongs.
5630 ========= ============= ============================
5631 `right` order of bins returned index `i` satisfies
5632 ========= ============= ============================
5633 ``False`` increasing ``bins[i-1] <= x < bins[i]``
5634 ``True`` increasing ``bins[i-1] < x <= bins[i]``
5635 ``False`` decreasing ``bins[i-1] > x >= bins[i]``
5636 ``True`` decreasing ``bins[i-1] >= x > bins[i]``
5637 ========= ============= ============================
5639 If values in `x` are beyond the bounds of `bins`, 0 or ``len(bins)`` is
5640 returned as appropriate.
5642 Parameters
5643 ----------
5644 x : array_like
5645 Input array to be binned. Prior to NumPy 1.10.0, this array had to
5646 be 1-dimensional, but can now have any shape.
5647 bins : array_like
5648 Array of bins. It has to be 1-dimensional and monotonic.
5649 right : bool, optional
5650 Indicating whether the intervals include the right or the left bin
5651 edge. Default behavior is (right==False) indicating that the interval
5652 does not include the right edge. The left bin end is open in this
5653 case, i.e., bins[i-1] <= x < bins[i] is the default behavior for
5654 monotonically increasing bins.
5656 Returns
5657 -------
5658 indices : ndarray of ints
5659 Output array of indices, of same shape as `x`.
5661 Raises
5662 ------
5663 ValueError
5664 If `bins` is not monotonic.
5665 TypeError
5666 If the type of the input is complex.
5668 See Also
5669 --------
5670 bincount, histogram, unique, searchsorted
5672 Notes
5673 -----
5674 If values in `x` are such that they fall outside the bin range,
5675 attempting to index `bins` with the indices that `digitize` returns
5676 will result in an IndexError.
5678 .. versionadded:: 1.10.0
5680 `np.digitize` is implemented in terms of `np.searchsorted`. This means
5681 that a binary search is used to bin the values, which scales much better
5682 for larger number of bins than the previous linear search. It also removes
5683 the requirement for the input array to be 1-dimensional.
5685 For monotonically _increasing_ `bins`, the following are equivalent::
5687 np.digitize(x, bins, right=True)
5688 np.searchsorted(bins, x, side='left')
5690 Note that as the order of the arguments are reversed, the side must be too.
5691 The `searchsorted` call is marginally faster, as it does not do any
5692 monotonicity checks. Perhaps more importantly, it supports all dtypes.
5694 Examples
5695 --------
5696 >>> x = np.array([0.2, 6.4, 3.0, 1.6])
5697 >>> bins = np.array([0.0, 1.0, 2.5, 4.0, 10.0])
5698 >>> inds = np.digitize(x, bins)
5699 >>> inds
5700 array([1, 4, 3, 2])
5701 >>> for n in range(x.size):
5702 ... print(bins[inds[n]-1], "<=", x[n], "<", bins[inds[n]])
5703 ...
5704 0.0 <= 0.2 < 1.0
5705 4.0 <= 6.4 < 10.0
5706 2.5 <= 3.0 < 4.0
5707 1.0 <= 1.6 < 2.5
5709 >>> x = np.array([1.2, 10.0, 12.4, 15.5, 20.])
5710 >>> bins = np.array([0, 5, 10, 15, 20])
5711 >>> np.digitize(x,bins,right=True)
5712 array([1, 2, 3, 4, 4])
5713 >>> np.digitize(x,bins,right=False)
5714 array([1, 3, 3, 4, 5])
5715 """
5716 x = _nx.asarray(x)
5717 bins = _nx.asarray(bins)
5719 # here for compatibility, searchsorted below is happy to take this
5720 if np.issubdtype(x.dtype, _nx.complexfloating):
5721 raise TypeError("x may not be complex")
5723 mono = _monotonicity(bins)
5724 if mono == 0:
5725 raise ValueError("bins must be monotonically increasing or decreasing")
5727 # this is backwards because the arguments below are swapped
5728 side = 'left' if right else 'right'
5729 if mono == -1:
5730 # reverse the bins, and invert the results
5731 return len(bins) - _nx.searchsorted(bins[::-1], x, side=side)
5732 else:
5733 return _nx.searchsorted(bins, x, side=side)