Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/numpy/lib/function_base.py: 13%
1227 statements
« prev ^ index » next coverage.py v7.0.5, created at 2023-01-17 06:27 +0000
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 _insert, add_docstring, bincount, normalize_axis_index, _monotonicity,
28 interp as compiled_interp, interp_complex as compiled_interp_complex
29 )
30from numpy.core.umath import _add_newdoc_ufunc as add_newdoc_ufunc
32import builtins
34# needed in this module for compatibility
35from numpy.lib.histograms import histogram, histogramdd # noqa: F401
38array_function_dispatch = functools.partial(
39 overrides.array_function_dispatch, module='numpy')
42__all__ = [
43 'select', 'piecewise', 'trim_zeros', 'copy', 'iterable', 'percentile',
44 'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'disp', 'flip',
45 'rot90', 'extract', 'place', 'vectorize', 'asarray_chkfinite', 'average',
46 'bincount', 'digitize', 'cov', 'corrcoef',
47 'msort', 'median', 'sinc', 'hamming', 'hanning', 'bartlett',
48 'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc', 'add_docstring',
49 'meshgrid', 'delete', 'insert', 'append', 'interp', 'add_newdoc_ufunc',
50 'quantile'
51 ]
53# _QuantileMethods is a dictionary listing all the supported methods to
54# compute quantile/percentile.
55#
56# Below virtual_index refer to the index of the element where the percentile
57# would be found in the sorted sample.
58# When the sample contains exactly the percentile wanted, the virtual_index is
59# an integer to the index of this element.
60# When the percentile wanted is in between two elements, the virtual_index
61# is made of a integer part (a.k.a 'i' or 'left') and a fractional part
62# (a.k.a 'g' or 'gamma')
63#
64# Each method in _QuantileMethods has two properties
65# get_virtual_index : Callable
66# The function used to compute the virtual_index.
67# fix_gamma : Callable
68# A function used for discret methods to force the index to a specific value.
69_QuantileMethods = dict(
70 # --- HYNDMAN and FAN METHODS
71 # Discrete methods
72 inverted_cdf=dict(
73 get_virtual_index=lambda n, quantiles: _inverted_cdf(n, quantiles),
74 fix_gamma=lambda gamma, _: gamma, # should never be called
75 ),
76 averaged_inverted_cdf=dict(
77 get_virtual_index=lambda n, quantiles: (n * quantiles) - 1,
78 fix_gamma=lambda gamma, _: _get_gamma_mask(
79 shape=gamma.shape,
80 default_value=1.,
81 conditioned_value=0.5,
82 where=gamma == 0),
83 ),
84 closest_observation=dict(
85 get_virtual_index=lambda n, quantiles: _closest_observation(n,
86 quantiles),
87 fix_gamma=lambda gamma, _: gamma, # should never be called
88 ),
89 # Continuous methods
90 interpolated_inverted_cdf=dict(
91 get_virtual_index=lambda n, quantiles:
92 _compute_virtual_index(n, quantiles, 0, 1),
93 fix_gamma=lambda gamma, _: gamma,
94 ),
95 hazen=dict(
96 get_virtual_index=lambda n, quantiles:
97 _compute_virtual_index(n, quantiles, 0.5, 0.5),
98 fix_gamma=lambda gamma, _: gamma,
99 ),
100 weibull=dict(
101 get_virtual_index=lambda n, quantiles:
102 _compute_virtual_index(n, quantiles, 0, 0),
103 fix_gamma=lambda gamma, _: gamma,
104 ),
105 # Default method.
106 # To avoid some rounding issues, `(n-1) * quantiles` is preferred to
107 # `_compute_virtual_index(n, quantiles, 1, 1)`.
108 # They are mathematically equivalent.
109 linear=dict(
110 get_virtual_index=lambda n, quantiles: (n - 1) * quantiles,
111 fix_gamma=lambda gamma, _: gamma,
112 ),
113 median_unbiased=dict(
114 get_virtual_index=lambda n, quantiles:
115 _compute_virtual_index(n, quantiles, 1 / 3.0, 1 / 3.0),
116 fix_gamma=lambda gamma, _: gamma,
117 ),
118 normal_unbiased=dict(
119 get_virtual_index=lambda n, quantiles:
120 _compute_virtual_index(n, quantiles, 3 / 8.0, 3 / 8.0),
121 fix_gamma=lambda gamma, _: gamma,
122 ),
123 # --- OTHER METHODS
124 lower=dict(
125 get_virtual_index=lambda n, quantiles: np.floor(
126 (n - 1) * quantiles).astype(np.intp),
127 fix_gamma=lambda gamma, _: gamma,
128 # should never be called, index dtype is int
129 ),
130 higher=dict(
131 get_virtual_index=lambda n, quantiles: np.ceil(
132 (n - 1) * quantiles).astype(np.intp),
133 fix_gamma=lambda gamma, _: gamma,
134 # should never be called, index dtype is int
135 ),
136 midpoint=dict(
137 get_virtual_index=lambda n, quantiles: 0.5 * (
138 np.floor((n - 1) * quantiles)
139 + np.ceil((n - 1) * quantiles)),
140 fix_gamma=lambda gamma, index: _get_gamma_mask(
141 shape=gamma.shape,
142 default_value=0.5,
143 conditioned_value=0.,
144 where=index % 1 == 0),
145 ),
146 nearest=dict(
147 get_virtual_index=lambda n, quantiles: np.around(
148 (n - 1) * quantiles).astype(np.intp),
149 fix_gamma=lambda gamma, _: gamma,
150 # should never be called, index dtype is int
151 ))
154def _rot90_dispatcher(m, k=None, axes=None):
155 return (m,)
158@array_function_dispatch(_rot90_dispatcher)
159def rot90(m, k=1, axes=(0, 1)):
160 """
161 Rotate an array by 90 degrees in the plane specified by axes.
163 Rotation direction is from the first towards the second axis.
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 else:
1315 return outvals
1318def _diff_dispatcher(a, n=None, axis=None, prepend=None, append=None):
1319 return (a, prepend, append)
1322@array_function_dispatch(_diff_dispatcher)
1323def diff(a, n=1, axis=-1, prepend=np._NoValue, append=np._NoValue):
1324 """
1325 Calculate the n-th discrete difference along the given axis.
1327 The first difference is given by ``out[i] = a[i+1] - a[i]`` along
1328 the given axis, higher differences are calculated by using `diff`
1329 recursively.
1331 Parameters
1332 ----------
1333 a : array_like
1334 Input array
1335 n : int, optional
1336 The number of times values are differenced. If zero, the input
1337 is returned as-is.
1338 axis : int, optional
1339 The axis along which the difference is taken, default is the
1340 last axis.
1341 prepend, append : array_like, optional
1342 Values to prepend or append to `a` along axis prior to
1343 performing the difference. Scalar values are expanded to
1344 arrays with length 1 in the direction of axis and the shape
1345 of the input array in along all other axes. Otherwise the
1346 dimension and shape must match `a` except along axis.
1348 .. versionadded:: 1.16.0
1350 Returns
1351 -------
1352 diff : ndarray
1353 The n-th differences. The shape of the output is the same as `a`
1354 except along `axis` where the dimension is smaller by `n`. The
1355 type of the output is the same as the type of the difference
1356 between any two elements of `a`. This is the same as the type of
1357 `a` in most cases. A notable exception is `datetime64`, which
1358 results in a `timedelta64` output array.
1360 See Also
1361 --------
1362 gradient, ediff1d, cumsum
1364 Notes
1365 -----
1366 Type is preserved for boolean arrays, so the result will contain
1367 `False` when consecutive elements are the same and `True` when they
1368 differ.
1370 For unsigned integer arrays, the results will also be unsigned. This
1371 should not be surprising, as the result is consistent with
1372 calculating the difference directly:
1374 >>> u8_arr = np.array([1, 0], dtype=np.uint8)
1375 >>> np.diff(u8_arr)
1376 array([255], dtype=uint8)
1377 >>> u8_arr[1,...] - u8_arr[0,...]
1378 255
1380 If this is not desirable, then the array should be cast to a larger
1381 integer type first:
1383 >>> i16_arr = u8_arr.astype(np.int16)
1384 >>> np.diff(i16_arr)
1385 array([-1], dtype=int16)
1387 Examples
1388 --------
1389 >>> x = np.array([1, 2, 4, 7, 0])
1390 >>> np.diff(x)
1391 array([ 1, 2, 3, -7])
1392 >>> np.diff(x, n=2)
1393 array([ 1, 1, -10])
1395 >>> x = np.array([[1, 3, 6, 10], [0, 5, 6, 8]])
1396 >>> np.diff(x)
1397 array([[2, 3, 4],
1398 [5, 1, 2]])
1399 >>> np.diff(x, axis=0)
1400 array([[-1, 2, 0, -2]])
1402 >>> x = np.arange('1066-10-13', '1066-10-16', dtype=np.datetime64)
1403 >>> np.diff(x)
1404 array([1, 1], dtype='timedelta64[D]')
1406 """
1407 if n == 0:
1408 return a
1409 if n < 0:
1410 raise ValueError(
1411 "order must be non-negative but got " + repr(n))
1413 a = asanyarray(a)
1414 nd = a.ndim
1415 if nd == 0:
1416 raise ValueError("diff requires input that is at least one dimensional")
1417 axis = normalize_axis_index(axis, nd)
1419 combined = []
1420 if prepend is not np._NoValue:
1421 prepend = np.asanyarray(prepend)
1422 if prepend.ndim == 0:
1423 shape = list(a.shape)
1424 shape[axis] = 1
1425 prepend = np.broadcast_to(prepend, tuple(shape))
1426 combined.append(prepend)
1428 combined.append(a)
1430 if append is not np._NoValue:
1431 append = np.asanyarray(append)
1432 if append.ndim == 0:
1433 shape = list(a.shape)
1434 shape[axis] = 1
1435 append = np.broadcast_to(append, tuple(shape))
1436 combined.append(append)
1438 if len(combined) > 1:
1439 a = np.concatenate(combined, axis)
1441 slice1 = [slice(None)] * nd
1442 slice2 = [slice(None)] * nd
1443 slice1[axis] = slice(1, None)
1444 slice2[axis] = slice(None, -1)
1445 slice1 = tuple(slice1)
1446 slice2 = tuple(slice2)
1448 op = not_equal if a.dtype == np.bool_ else subtract
1449 for _ in range(n):
1450 a = op(a[slice1], a[slice2])
1452 return a
1455def _interp_dispatcher(x, xp, fp, left=None, right=None, period=None):
1456 return (x, xp, fp)
1459@array_function_dispatch(_interp_dispatcher)
1460def interp(x, xp, fp, left=None, right=None, period=None):
1461 """
1462 One-dimensional linear interpolation for monotonically increasing sample points.
1464 Returns the one-dimensional piecewise linear interpolant to a function
1465 with given discrete data points (`xp`, `fp`), evaluated at `x`.
1467 Parameters
1468 ----------
1469 x : array_like
1470 The x-coordinates at which to evaluate the interpolated values.
1472 xp : 1-D sequence of floats
1473 The x-coordinates of the data points, must be increasing if argument
1474 `period` is not specified. Otherwise, `xp` is internally sorted after
1475 normalizing the periodic boundaries with ``xp = xp % period``.
1477 fp : 1-D sequence of float or complex
1478 The y-coordinates of the data points, same length as `xp`.
1480 left : optional float or complex corresponding to fp
1481 Value to return for `x < xp[0]`, default is `fp[0]`.
1483 right : optional float or complex corresponding to fp
1484 Value to return for `x > xp[-1]`, default is `fp[-1]`.
1486 period : None or float, optional
1487 A period for the x-coordinates. This parameter allows the proper
1488 interpolation of angular x-coordinates. Parameters `left` and `right`
1489 are ignored if `period` is specified.
1491 .. versionadded:: 1.10.0
1493 Returns
1494 -------
1495 y : float or complex (corresponding to fp) or ndarray
1496 The interpolated values, same shape as `x`.
1498 Raises
1499 ------
1500 ValueError
1501 If `xp` and `fp` have different length
1502 If `xp` or `fp` are not 1-D sequences
1503 If `period == 0`
1505 See Also
1506 --------
1507 scipy.interpolate
1509 Warnings
1510 --------
1511 The x-coordinate sequence is expected to be increasing, but this is not
1512 explicitly enforced. However, if the sequence `xp` is non-increasing,
1513 interpolation results are meaningless.
1515 Note that, since NaN is unsortable, `xp` also cannot contain NaNs.
1517 A simple check for `xp` being strictly increasing is::
1519 np.all(np.diff(xp) > 0)
1521 Examples
1522 --------
1523 >>> xp = [1, 2, 3]
1524 >>> fp = [3, 2, 0]
1525 >>> np.interp(2.5, xp, fp)
1526 1.0
1527 >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
1528 array([3. , 3. , 2.5 , 0.56, 0. ])
1529 >>> UNDEF = -99.0
1530 >>> np.interp(3.14, xp, fp, right=UNDEF)
1531 -99.0
1533 Plot an interpolant to the sine function:
1535 >>> x = np.linspace(0, 2*np.pi, 10)
1536 >>> y = np.sin(x)
1537 >>> xvals = np.linspace(0, 2*np.pi, 50)
1538 >>> yinterp = np.interp(xvals, x, y)
1539 >>> import matplotlib.pyplot as plt
1540 >>> plt.plot(x, y, 'o')
1541 [<matplotlib.lines.Line2D object at 0x...>]
1542 >>> plt.plot(xvals, yinterp, '-x')
1543 [<matplotlib.lines.Line2D object at 0x...>]
1544 >>> plt.show()
1546 Interpolation with periodic x-coordinates:
1548 >>> x = [-180, -170, -185, 185, -10, -5, 0, 365]
1549 >>> xp = [190, -190, 350, -350]
1550 >>> fp = [5, 10, 3, 4]
1551 >>> np.interp(x, xp, fp, period=360)
1552 array([7.5 , 5. , 8.75, 6.25, 3. , 3.25, 3.5 , 3.75])
1554 Complex interpolation:
1556 >>> x = [1.5, 4.0]
1557 >>> xp = [2,3,5]
1558 >>> fp = [1.0j, 0, 2+3j]
1559 >>> np.interp(x, xp, fp)
1560 array([0.+1.j , 1.+1.5j])
1562 """
1564 fp = np.asarray(fp)
1566 if np.iscomplexobj(fp):
1567 interp_func = compiled_interp_complex
1568 input_dtype = np.complex128
1569 else:
1570 interp_func = compiled_interp
1571 input_dtype = np.float64
1573 if period is not None:
1574 if period == 0:
1575 raise ValueError("period must be a non-zero value")
1576 period = abs(period)
1577 left = None
1578 right = None
1580 x = np.asarray(x, dtype=np.float64)
1581 xp = np.asarray(xp, dtype=np.float64)
1582 fp = np.asarray(fp, dtype=input_dtype)
1584 if xp.ndim != 1 or fp.ndim != 1:
1585 raise ValueError("Data points must be 1-D sequences")
1586 if xp.shape[0] != fp.shape[0]:
1587 raise ValueError("fp and xp are not of the same length")
1588 # normalizing periodic boundaries
1589 x = x % period
1590 xp = xp % period
1591 asort_xp = np.argsort(xp)
1592 xp = xp[asort_xp]
1593 fp = fp[asort_xp]
1594 xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period))
1595 fp = np.concatenate((fp[-1:], fp, fp[0:1]))
1597 return interp_func(x, xp, fp, left, right)
1600def _angle_dispatcher(z, deg=None):
1601 return (z,)
1604@array_function_dispatch(_angle_dispatcher)
1605def angle(z, deg=False):
1606 """
1607 Return the angle of the complex argument.
1609 Parameters
1610 ----------
1611 z : array_like
1612 A complex number or sequence of complex numbers.
1613 deg : bool, optional
1614 Return angle in degrees if True, radians if False (default).
1616 Returns
1617 -------
1618 angle : ndarray or scalar
1619 The counterclockwise angle from the positive real axis on the complex
1620 plane in the range ``(-pi, pi]``, with dtype as numpy.float64.
1622 .. versionchanged:: 1.16.0
1623 This function works on subclasses of ndarray like `ma.array`.
1625 See Also
1626 --------
1627 arctan2
1628 absolute
1630 Notes
1631 -----
1632 Although the angle of the complex number 0 is undefined, ``numpy.angle(0)``
1633 returns the value 0.
1635 Examples
1636 --------
1637 >>> np.angle([1.0, 1.0j, 1+1j]) # in radians
1638 array([ 0. , 1.57079633, 0.78539816]) # may vary
1639 >>> np.angle(1+1j, deg=True) # in degrees
1640 45.0
1642 """
1643 z = asanyarray(z)
1644 if issubclass(z.dtype.type, _nx.complexfloating):
1645 zimag = z.imag
1646 zreal = z.real
1647 else:
1648 zimag = 0
1649 zreal = z
1651 a = arctan2(zimag, zreal)
1652 if deg:
1653 a *= 180/pi
1654 return a
1657def _unwrap_dispatcher(p, discont=None, axis=None, *, period=None):
1658 return (p,)
1661@array_function_dispatch(_unwrap_dispatcher)
1662def unwrap(p, discont=None, axis=-1, *, period=2*pi):
1663 r"""
1664 Unwrap by taking the complement of large deltas with respect to the period.
1666 This unwraps a signal `p` by changing elements which have an absolute
1667 difference from their predecessor of more than ``max(discont, period/2)``
1668 to their `period`-complementary values.
1670 For the default case where `period` is :math:`2\pi` and `discont` is
1671 :math:`\pi`, this unwraps a radian phase `p` such that adjacent differences
1672 are never greater than :math:`\pi` by adding :math:`2k\pi` for some
1673 integer :math:`k`.
1675 Parameters
1676 ----------
1677 p : array_like
1678 Input array.
1679 discont : float, optional
1680 Maximum discontinuity between values, default is ``period/2``.
1681 Values below ``period/2`` are treated as if they were ``period/2``.
1682 To have an effect different from the default, `discont` should be
1683 larger than ``period/2``.
1684 axis : int, optional
1685 Axis along which unwrap will operate, default is the last axis.
1686 period : float, optional
1687 Size of the range over which the input wraps. By default, it is
1688 ``2 pi``.
1690 .. versionadded:: 1.21.0
1692 Returns
1693 -------
1694 out : ndarray
1695 Output array.
1697 See Also
1698 --------
1699 rad2deg, deg2rad
1701 Notes
1702 -----
1703 If the discontinuity in `p` is smaller than ``period/2``,
1704 but larger than `discont`, no unwrapping is done because taking
1705 the complement would only make the discontinuity larger.
1707 Examples
1708 --------
1709 >>> phase = np.linspace(0, np.pi, num=5)
1710 >>> phase[3:] += np.pi
1711 >>> phase
1712 array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531]) # may vary
1713 >>> np.unwrap(phase)
1714 array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ]) # may vary
1715 >>> np.unwrap([0, 1, 2, -1, 0], period=4)
1716 array([0, 1, 2, 3, 4])
1717 >>> np.unwrap([ 1, 2, 3, 4, 5, 6, 1, 2, 3], period=6)
1718 array([1, 2, 3, 4, 5, 6, 7, 8, 9])
1719 >>> np.unwrap([2, 3, 4, 5, 2, 3, 4, 5], period=4)
1720 array([2, 3, 4, 5, 6, 7, 8, 9])
1721 >>> phase_deg = np.mod(np.linspace(0 ,720, 19), 360) - 180
1722 >>> np.unwrap(phase_deg, period=360)
1723 array([-180., -140., -100., -60., -20., 20., 60., 100., 140.,
1724 180., 220., 260., 300., 340., 380., 420., 460., 500.,
1725 540.])
1726 """
1727 p = asarray(p)
1728 nd = p.ndim
1729 dd = diff(p, axis=axis)
1730 if discont is None:
1731 discont = period/2
1732 slice1 = [slice(None, None)]*nd # full slices
1733 slice1[axis] = slice(1, None)
1734 slice1 = tuple(slice1)
1735 dtype = np.result_type(dd, period)
1736 if _nx.issubdtype(dtype, _nx.integer):
1737 interval_high, rem = divmod(period, 2)
1738 boundary_ambiguous = rem == 0
1739 else:
1740 interval_high = period / 2
1741 boundary_ambiguous = True
1742 interval_low = -interval_high
1743 ddmod = mod(dd - interval_low, period) + interval_low
1744 if boundary_ambiguous:
1745 # for `mask = (abs(dd) == period/2)`, the above line made
1746 # `ddmod[mask] == -period/2`. correct these such that
1747 # `ddmod[mask] == sign(dd[mask])*period/2`.
1748 _nx.copyto(ddmod, interval_high,
1749 where=(ddmod == interval_low) & (dd > 0))
1750 ph_correct = ddmod - dd
1751 _nx.copyto(ph_correct, 0, where=abs(dd) < discont)
1752 up = array(p, copy=True, dtype=dtype)
1753 up[slice1] = p[slice1] + ph_correct.cumsum(axis)
1754 return up
1757def _sort_complex(a):
1758 return (a,)
1761@array_function_dispatch(_sort_complex)
1762def sort_complex(a):
1763 """
1764 Sort a complex array using the real part first, then the imaginary part.
1766 Parameters
1767 ----------
1768 a : array_like
1769 Input array
1771 Returns
1772 -------
1773 out : complex ndarray
1774 Always returns a sorted complex array.
1776 Examples
1777 --------
1778 >>> np.sort_complex([5, 3, 6, 2, 1])
1779 array([1.+0.j, 2.+0.j, 3.+0.j, 5.+0.j, 6.+0.j])
1781 >>> np.sort_complex([1 + 2j, 2 - 1j, 3 - 2j, 3 - 3j, 3 + 5j])
1782 array([1.+2.j, 2.-1.j, 3.-3.j, 3.-2.j, 3.+5.j])
1784 """
1785 b = array(a, copy=True)
1786 b.sort()
1787 if not issubclass(b.dtype.type, _nx.complexfloating):
1788 if b.dtype.char in 'bhBH':
1789 return b.astype('F')
1790 elif b.dtype.char == 'g':
1791 return b.astype('G')
1792 else:
1793 return b.astype('D')
1794 else:
1795 return b
1798def _trim_zeros(filt, trim=None):
1799 return (filt,)
1802@array_function_dispatch(_trim_zeros)
1803def trim_zeros(filt, trim='fb'):
1804 """
1805 Trim the leading and/or trailing zeros from a 1-D array or sequence.
1807 Parameters
1808 ----------
1809 filt : 1-D array or sequence
1810 Input array.
1811 trim : str, optional
1812 A string with 'f' representing trim from front and 'b' to trim from
1813 back. Default is 'fb', trim zeros from both front and back of the
1814 array.
1816 Returns
1817 -------
1818 trimmed : 1-D array or sequence
1819 The result of trimming the input. The input data type is preserved.
1821 Examples
1822 --------
1823 >>> a = np.array((0, 0, 0, 1, 2, 3, 0, 2, 1, 0))
1824 >>> np.trim_zeros(a)
1825 array([1, 2, 3, 0, 2, 1])
1827 >>> np.trim_zeros(a, 'b')
1828 array([0, 0, 0, ..., 0, 2, 1])
1830 The input data type is preserved, list/tuple in means list/tuple out.
1832 >>> np.trim_zeros([0, 1, 2, 0])
1833 [1, 2]
1835 """
1837 first = 0
1838 trim = trim.upper()
1839 if 'F' in trim:
1840 for i in filt:
1841 if i != 0.:
1842 break
1843 else:
1844 first = first + 1
1845 last = len(filt)
1846 if 'B' in trim:
1847 for i in filt[::-1]:
1848 if i != 0.:
1849 break
1850 else:
1851 last = last - 1
1852 return filt[first:last]
1855def _extract_dispatcher(condition, arr):
1856 return (condition, arr)
1859@array_function_dispatch(_extract_dispatcher)
1860def extract(condition, arr):
1861 """
1862 Return the elements of an array that satisfy some condition.
1864 This is equivalent to ``np.compress(ravel(condition), ravel(arr))``. If
1865 `condition` is boolean ``np.extract`` is equivalent to ``arr[condition]``.
1867 Note that `place` does the exact opposite of `extract`.
1869 Parameters
1870 ----------
1871 condition : array_like
1872 An array whose nonzero or True entries indicate the elements of `arr`
1873 to extract.
1874 arr : array_like
1875 Input array of the same size as `condition`.
1877 Returns
1878 -------
1879 extract : ndarray
1880 Rank 1 array of values from `arr` where `condition` is True.
1882 See Also
1883 --------
1884 take, put, copyto, compress, place
1886 Examples
1887 --------
1888 >>> arr = np.arange(12).reshape((3, 4))
1889 >>> arr
1890 array([[ 0, 1, 2, 3],
1891 [ 4, 5, 6, 7],
1892 [ 8, 9, 10, 11]])
1893 >>> condition = np.mod(arr, 3)==0
1894 >>> condition
1895 array([[ True, False, False, True],
1896 [False, False, True, False],
1897 [False, True, False, False]])
1898 >>> np.extract(condition, arr)
1899 array([0, 3, 6, 9])
1902 If `condition` is boolean:
1904 >>> arr[condition]
1905 array([0, 3, 6, 9])
1907 """
1908 return _nx.take(ravel(arr), nonzero(ravel(condition))[0])
1911def _place_dispatcher(arr, mask, vals):
1912 return (arr, mask, vals)
1915@array_function_dispatch(_place_dispatcher)
1916def place(arr, mask, vals):
1917 """
1918 Change elements of an array based on conditional and input values.
1920 Similar to ``np.copyto(arr, vals, where=mask)``, the difference is that
1921 `place` uses the first N elements of `vals`, where N is the number of
1922 True values in `mask`, while `copyto` uses the elements where `mask`
1923 is True.
1925 Note that `extract` does the exact opposite of `place`.
1927 Parameters
1928 ----------
1929 arr : ndarray
1930 Array to put data into.
1931 mask : array_like
1932 Boolean mask array. Must have the same size as `a`.
1933 vals : 1-D sequence
1934 Values to put into `a`. Only the first N elements are used, where
1935 N is the number of True values in `mask`. If `vals` is smaller
1936 than N, it will be repeated, and if elements of `a` are to be masked,
1937 this sequence must be non-empty.
1939 See Also
1940 --------
1941 copyto, put, take, extract
1943 Examples
1944 --------
1945 >>> arr = np.arange(6).reshape(2, 3)
1946 >>> np.place(arr, arr>2, [44, 55])
1947 >>> arr
1948 array([[ 0, 1, 2],
1949 [44, 55, 44]])
1951 """
1952 if not isinstance(arr, np.ndarray):
1953 raise TypeError("argument 1 must be numpy.ndarray, "
1954 "not {name}".format(name=type(arr).__name__))
1956 return _insert(arr, mask, vals)
1959def disp(mesg, device=None, linefeed=True):
1960 """
1961 Display a message on a device.
1963 Parameters
1964 ----------
1965 mesg : str
1966 Message to display.
1967 device : object
1968 Device to write message. If None, defaults to ``sys.stdout`` which is
1969 very similar to ``print``. `device` needs to have ``write()`` and
1970 ``flush()`` methods.
1971 linefeed : bool, optional
1972 Option whether to print a line feed or not. Defaults to True.
1974 Raises
1975 ------
1976 AttributeError
1977 If `device` does not have a ``write()`` or ``flush()`` method.
1979 Examples
1980 --------
1981 Besides ``sys.stdout``, a file-like object can also be used as it has
1982 both required methods:
1984 >>> from io import StringIO
1985 >>> buf = StringIO()
1986 >>> np.disp(u'"Display" in a file', device=buf)
1987 >>> buf.getvalue()
1988 '"Display" in a file\\n'
1990 """
1991 if device is None:
1992 device = sys.stdout
1993 if linefeed:
1994 device.write('%s\n' % mesg)
1995 else:
1996 device.write('%s' % mesg)
1997 device.flush()
1998 return
2001# See https://docs.scipy.org/doc/numpy/reference/c-api.generalized-ufuncs.html
2002_DIMENSION_NAME = r'\w+'
2003_CORE_DIMENSION_LIST = '(?:{0:}(?:,{0:})*)?'.format(_DIMENSION_NAME)
2004_ARGUMENT = r'\({}\)'.format(_CORE_DIMENSION_LIST)
2005_ARGUMENT_LIST = '{0:}(?:,{0:})*'.format(_ARGUMENT)
2006_SIGNATURE = '^{0:}->{0:}$'.format(_ARGUMENT_LIST)
2009def _parse_gufunc_signature(signature):
2010 """
2011 Parse string signatures for a generalized universal function.
2013 Arguments
2014 ---------
2015 signature : string
2016 Generalized universal function signature, e.g., ``(m,n),(n,p)->(m,p)``
2017 for ``np.matmul``.
2019 Returns
2020 -------
2021 Tuple of input and output core dimensions parsed from the signature, each
2022 of the form List[Tuple[str, ...]].
2023 """
2024 signature = re.sub(r'\s+', '', signature)
2026 if not re.match(_SIGNATURE, signature):
2027 raise ValueError(
2028 'not a valid gufunc signature: {}'.format(signature))
2029 return tuple([tuple(re.findall(_DIMENSION_NAME, arg))
2030 for arg in re.findall(_ARGUMENT, arg_list)]
2031 for arg_list in signature.split('->'))
2034def _update_dim_sizes(dim_sizes, arg, core_dims):
2035 """
2036 Incrementally check and update core dimension sizes for a single argument.
2038 Arguments
2039 ---------
2040 dim_sizes : Dict[str, int]
2041 Sizes of existing core dimensions. Will be updated in-place.
2042 arg : ndarray
2043 Argument to examine.
2044 core_dims : Tuple[str, ...]
2045 Core dimensions for this argument.
2046 """
2047 if not core_dims:
2048 return
2050 num_core_dims = len(core_dims)
2051 if arg.ndim < num_core_dims:
2052 raise ValueError(
2053 '%d-dimensional argument does not have enough '
2054 'dimensions for all core dimensions %r'
2055 % (arg.ndim, core_dims))
2057 core_shape = arg.shape[-num_core_dims:]
2058 for dim, size in zip(core_dims, core_shape):
2059 if dim in dim_sizes:
2060 if size != dim_sizes[dim]:
2061 raise ValueError(
2062 'inconsistent size for core dimension %r: %r vs %r'
2063 % (dim, size, dim_sizes[dim]))
2064 else:
2065 dim_sizes[dim] = size
2068def _parse_input_dimensions(args, input_core_dims):
2069 """
2070 Parse broadcast and core dimensions for vectorize with a signature.
2072 Arguments
2073 ---------
2074 args : Tuple[ndarray, ...]
2075 Tuple of input arguments to examine.
2076 input_core_dims : List[Tuple[str, ...]]
2077 List of core dimensions corresponding to each input.
2079 Returns
2080 -------
2081 broadcast_shape : Tuple[int, ...]
2082 Common shape to broadcast all non-core dimensions to.
2083 dim_sizes : Dict[str, int]
2084 Common sizes for named core dimensions.
2085 """
2086 broadcast_args = []
2087 dim_sizes = {}
2088 for arg, core_dims in zip(args, input_core_dims):
2089 _update_dim_sizes(dim_sizes, arg, core_dims)
2090 ndim = arg.ndim - len(core_dims)
2091 dummy_array = np.lib.stride_tricks.as_strided(0, arg.shape[:ndim])
2092 broadcast_args.append(dummy_array)
2093 broadcast_shape = np.lib.stride_tricks._broadcast_shape(*broadcast_args)
2094 return broadcast_shape, dim_sizes
2097def _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims):
2098 """Helper for calculating broadcast shapes with core dimensions."""
2099 return [broadcast_shape + tuple(dim_sizes[dim] for dim in core_dims)
2100 for core_dims in list_of_core_dims]
2103def _create_arrays(broadcast_shape, dim_sizes, list_of_core_dims, dtypes,
2104 results=None):
2105 """Helper for creating output arrays in vectorize."""
2106 shapes = _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims)
2107 if dtypes is None:
2108 dtypes = [None] * len(shapes)
2109 if results is None:
2110 arrays = tuple(np.empty(shape=shape, dtype=dtype)
2111 for shape, dtype in zip(shapes, dtypes))
2112 else:
2113 arrays = tuple(np.empty_like(result, shape=shape, dtype=dtype)
2114 for result, shape, dtype
2115 in zip(results, shapes, dtypes))
2116 return arrays
2119@set_module('numpy')
2120class vectorize:
2121 """
2122 vectorize(pyfunc, otypes=None, doc=None, excluded=None, cache=False,
2123 signature=None)
2125 Generalized function class.
2127 Define a vectorized function which takes a nested sequence of objects or
2128 numpy arrays as inputs and returns a single numpy array or a tuple of numpy
2129 arrays. The vectorized function evaluates `pyfunc` over successive tuples
2130 of the input arrays like the python map function, except it uses the
2131 broadcasting rules of numpy.
2133 The data type of the output of `vectorized` is determined by calling
2134 the function with the first element of the input. This can be avoided
2135 by specifying the `otypes` argument.
2137 Parameters
2138 ----------
2139 pyfunc : callable
2140 A python function or method.
2141 otypes : str or list of dtypes, optional
2142 The output data type. It must be specified as either a string of
2143 typecode characters or a list of data type specifiers. There should
2144 be one data type specifier for each output.
2145 doc : str, optional
2146 The docstring for the function. If None, the docstring will be the
2147 ``pyfunc.__doc__``.
2148 excluded : set, optional
2149 Set of strings or integers representing the positional or keyword
2150 arguments for which the function will not be vectorized. These will be
2151 passed directly to `pyfunc` unmodified.
2153 .. versionadded:: 1.7.0
2155 cache : bool, optional
2156 If `True`, then cache the first function call that determines the number
2157 of outputs if `otypes` is not provided.
2159 .. versionadded:: 1.7.0
2161 signature : string, optional
2162 Generalized universal function signature, e.g., ``(m,n),(n)->(m)`` for
2163 vectorized matrix-vector multiplication. If provided, ``pyfunc`` will
2164 be called with (and expected to return) arrays with shapes given by the
2165 size of corresponding core dimensions. By default, ``pyfunc`` is
2166 assumed to take scalars as input and output.
2168 .. versionadded:: 1.12.0
2170 Returns
2171 -------
2172 vectorized : callable
2173 Vectorized function.
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 """
2271 def __init__(self, pyfunc, otypes=None, doc=None, excluded=None,
2272 cache=False, signature=None):
2273 self.pyfunc = pyfunc
2274 self.cache = cache
2275 self.signature = signature
2276 self._ufunc = {} # Caching to improve default performance
2278 if doc is None:
2279 self.__doc__ = pyfunc.__doc__
2280 else:
2281 self.__doc__ = doc
2283 if isinstance(otypes, str):
2284 for char in otypes:
2285 if char not in typecodes['All']:
2286 raise ValueError("Invalid otype specified: %s" % (char,))
2287 elif iterable(otypes):
2288 otypes = ''.join([_nx.dtype(x).char for x in otypes])
2289 elif otypes is not None:
2290 raise ValueError("Invalid otype specification")
2291 self.otypes = otypes
2293 # Excluded variable support
2294 if excluded is None:
2295 excluded = set()
2296 self.excluded = set(excluded)
2298 if signature is not None:
2299 self._in_and_out_core_dims = _parse_gufunc_signature(signature)
2300 else:
2301 self._in_and_out_core_dims = None
2303 def __call__(self, *args, **kwargs):
2304 """
2305 Return arrays with the results of `pyfunc` broadcast (vectorized) over
2306 `args` and `kwargs` not in `excluded`.
2307 """
2308 excluded = self.excluded
2309 if not kwargs and not excluded:
2310 func = self.pyfunc
2311 vargs = args
2312 else:
2313 # The wrapper accepts only positional arguments: we use `names` and
2314 # `inds` to mutate `the_args` and `kwargs` to pass to the original
2315 # function.
2316 nargs = len(args)
2318 names = [_n for _n in kwargs if _n not in excluded]
2319 inds = [_i for _i in range(nargs) if _i not in excluded]
2320 the_args = list(args)
2322 def func(*vargs):
2323 for _n, _i in enumerate(inds):
2324 the_args[_i] = vargs[_n]
2325 kwargs.update(zip(names, vargs[len(inds):]))
2326 return self.pyfunc(*the_args, **kwargs)
2328 vargs = [args[_i] for _i in inds]
2329 vargs.extend([kwargs[_n] for _n in names])
2331 return self._vectorize_call(func=func, args=vargs)
2333 def _get_ufunc_and_otypes(self, func, args):
2334 """Return (ufunc, otypes)."""
2335 # frompyfunc will fail if args is empty
2336 if not args:
2337 raise ValueError('args can not be empty')
2339 if self.otypes is not None:
2340 otypes = self.otypes
2342 # self._ufunc is a dictionary whose keys are the number of
2343 # arguments (i.e. len(args)) and whose values are ufuncs created
2344 # by frompyfunc. len(args) can be different for different calls if
2345 # self.pyfunc has parameters with default values. We only use the
2346 # cache when func is self.pyfunc, which occurs when the call uses
2347 # only positional arguments and no arguments are excluded.
2349 nin = len(args)
2350 nout = len(self.otypes)
2351 if func is not self.pyfunc or nin not in self._ufunc:
2352 ufunc = frompyfunc(func, nin, nout)
2353 else:
2354 ufunc = None # We'll get it from self._ufunc
2355 if func is self.pyfunc:
2356 ufunc = self._ufunc.setdefault(nin, ufunc)
2357 else:
2358 # Get number of outputs and output types by calling the function on
2359 # the first entries of args. We also cache the result to prevent
2360 # the subsequent call when the ufunc is evaluated.
2361 # Assumes that ufunc first evaluates the 0th elements in the input
2362 # arrays (the input values are not checked to ensure this)
2363 args = [asarray(arg) for arg in args]
2364 if builtins.any(arg.size == 0 for arg in args):
2365 raise ValueError('cannot call `vectorize` on size 0 inputs '
2366 'unless `otypes` is set')
2368 inputs = [arg.flat[0] for arg in args]
2369 outputs = func(*inputs)
2371 # Performance note: profiling indicates that -- for simple
2372 # functions at least -- this wrapping can almost double the
2373 # execution time.
2374 # Hence we make it optional.
2375 if self.cache:
2376 _cache = [outputs]
2378 def _func(*vargs):
2379 if _cache:
2380 return _cache.pop()
2381 else:
2382 return func(*vargs)
2383 else:
2384 _func = func
2386 if isinstance(outputs, tuple):
2387 nout = len(outputs)
2388 else:
2389 nout = 1
2390 outputs = (outputs,)
2392 otypes = ''.join([asarray(outputs[_k]).dtype.char
2393 for _k in range(nout)])
2395 # Performance note: profiling indicates that creating the ufunc is
2396 # not a significant cost compared with wrapping so it seems not
2397 # worth trying to cache this.
2398 ufunc = frompyfunc(_func, len(args), nout)
2400 return ufunc, otypes
2402 def _vectorize_call(self, func, args):
2403 """Vectorized call to `func` over positional `args`."""
2404 if self.signature is not None:
2405 res = self._vectorize_call_with_signature(func, args)
2406 elif not args:
2407 res = func()
2408 else:
2409 ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args)
2411 # Convert args to object arrays first
2412 inputs = [asanyarray(a, dtype=object) for a in args]
2414 outputs = ufunc(*inputs)
2416 if ufunc.nout == 1:
2417 res = asanyarray(outputs, dtype=otypes[0])
2418 else:
2419 res = tuple([asanyarray(x, dtype=t)
2420 for x, t in zip(outputs, otypes)])
2421 return res
2423 def _vectorize_call_with_signature(self, func, args):
2424 """Vectorized call over positional arguments with a signature."""
2425 input_core_dims, output_core_dims = self._in_and_out_core_dims
2427 if len(args) != len(input_core_dims):
2428 raise TypeError('wrong number of positional arguments: '
2429 'expected %r, got %r'
2430 % (len(input_core_dims), len(args)))
2431 args = tuple(asanyarray(arg) for arg in args)
2433 broadcast_shape, dim_sizes = _parse_input_dimensions(
2434 args, input_core_dims)
2435 input_shapes = _calculate_shapes(broadcast_shape, dim_sizes,
2436 input_core_dims)
2437 args = [np.broadcast_to(arg, shape, subok=True)
2438 for arg, shape in zip(args, input_shapes)]
2440 outputs = None
2441 otypes = self.otypes
2442 nout = len(output_core_dims)
2444 for index in np.ndindex(*broadcast_shape):
2445 results = func(*(arg[index] for arg in args))
2447 n_results = len(results) if isinstance(results, tuple) else 1
2449 if nout != n_results:
2450 raise ValueError(
2451 'wrong number of outputs from pyfunc: expected %r, got %r'
2452 % (nout, n_results))
2454 if nout == 1:
2455 results = (results,)
2457 if outputs is None:
2458 for result, core_dims in zip(results, output_core_dims):
2459 _update_dim_sizes(dim_sizes, result, core_dims)
2461 outputs = _create_arrays(broadcast_shape, dim_sizes,
2462 output_core_dims, otypes, results)
2464 for output, result in zip(outputs, results):
2465 output[index] = result
2467 if outputs is None:
2468 # did not call the function even once
2469 if otypes is None:
2470 raise ValueError('cannot call `vectorize` on size 0 inputs '
2471 'unless `otypes` is set')
2472 if builtins.any(dim not in dim_sizes
2473 for dims in output_core_dims
2474 for dim in dims):
2475 raise ValueError('cannot call `vectorize` with a signature '
2476 'including new output dimensions on size 0 '
2477 'inputs')
2478 outputs = _create_arrays(broadcast_shape, dim_sizes,
2479 output_core_dims, otypes)
2481 return outputs[0] if nout == 1 else outputs
2484def _cov_dispatcher(m, y=None, rowvar=None, bias=None, ddof=None,
2485 fweights=None, aweights=None, *, dtype=None):
2486 return (m, y, fweights, aweights)
2489@array_function_dispatch(_cov_dispatcher)
2490def cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None,
2491 aweights=None, *, dtype=None):
2492 """
2493 Estimate a covariance matrix, given data and weights.
2495 Covariance indicates the level to which two variables vary together.
2496 If we examine N-dimensional samples, :math:`X = [x_1, x_2, ... x_N]^T`,
2497 then the covariance matrix element :math:`C_{ij}` is the covariance of
2498 :math:`x_i` and :math:`x_j`. The element :math:`C_{ii}` is the variance
2499 of :math:`x_i`.
2501 See the notes for an outline of the algorithm.
2503 Parameters
2504 ----------
2505 m : array_like
2506 A 1-D or 2-D array containing multiple variables and observations.
2507 Each row of `m` represents a variable, and each column a single
2508 observation of all those variables. Also see `rowvar` below.
2509 y : array_like, optional
2510 An additional set of variables and observations. `y` has the same form
2511 as that of `m`.
2512 rowvar : bool, optional
2513 If `rowvar` is True (default), then each row represents a
2514 variable, with observations in the columns. Otherwise, the relationship
2515 is transposed: each column represents a variable, while the rows
2516 contain observations.
2517 bias : bool, optional
2518 Default normalization (False) is by ``(N - 1)``, where ``N`` is the
2519 number of observations given (unbiased estimate). If `bias` is True,
2520 then normalization is by ``N``. These values can be overridden by using
2521 the keyword ``ddof`` in numpy versions >= 1.5.
2522 ddof : int, optional
2523 If not ``None`` the default value implied by `bias` is overridden.
2524 Note that ``ddof=1`` will return the unbiased estimate, even if both
2525 `fweights` and `aweights` are specified, and ``ddof=0`` will return
2526 the simple average. See the notes for the details. The default value
2527 is ``None``.
2529 .. versionadded:: 1.5
2530 fweights : array_like, int, optional
2531 1-D array of integer frequency weights; the number of times each
2532 observation vector should be repeated.
2534 .. versionadded:: 1.10
2535 aweights : array_like, optional
2536 1-D array of observation vector weights. These relative weights are
2537 typically large for observations considered "important" and smaller for
2538 observations considered less "important". If ``ddof=0`` the array of
2539 weights can be used to assign probabilities to observation vectors.
2541 .. versionadded:: 1.10
2542 dtype : data-type, optional
2543 Data-type of the result. By default, the return data-type will have
2544 at least `numpy.float64` precision.
2546 .. versionadded:: 1.20
2548 Returns
2549 -------
2550 out : ndarray
2551 The covariance matrix of the variables.
2553 See Also
2554 --------
2555 corrcoef : Normalized covariance matrix
2557 Notes
2558 -----
2559 Assume that the observations are in the columns of the observation
2560 array `m` and let ``f = fweights`` and ``a = aweights`` for brevity. The
2561 steps to compute the weighted covariance are as follows::
2563 >>> m = np.arange(10, dtype=np.float64)
2564 >>> f = np.arange(10) * 2
2565 >>> a = np.arange(10) ** 2.
2566 >>> ddof = 1
2567 >>> w = f * a
2568 >>> v1 = np.sum(w)
2569 >>> v2 = np.sum(w * a)
2570 >>> m -= np.sum(m * w, axis=None, keepdims=True) / v1
2571 >>> cov = np.dot(m * w, m.T) * v1 / (v1**2 - ddof * v2)
2573 Note that when ``a == 1``, the normalization factor
2574 ``v1 / (v1**2 - ddof * v2)`` goes over to ``1 / (np.sum(f) - ddof)``
2575 as it should.
2577 Examples
2578 --------
2579 Consider two variables, :math:`x_0` and :math:`x_1`, which
2580 correlate perfectly, but in opposite directions:
2582 >>> x = np.array([[0, 2], [1, 1], [2, 0]]).T
2583 >>> x
2584 array([[0, 1, 2],
2585 [2, 1, 0]])
2587 Note how :math:`x_0` increases while :math:`x_1` decreases. The covariance
2588 matrix shows this clearly:
2590 >>> np.cov(x)
2591 array([[ 1., -1.],
2592 [-1., 1.]])
2594 Note that element :math:`C_{0,1}`, which shows the correlation between
2595 :math:`x_0` and :math:`x_1`, is negative.
2597 Further, note how `x` and `y` are combined:
2599 >>> x = [-2.1, -1, 4.3]
2600 >>> y = [3, 1.1, 0.12]
2601 >>> X = np.stack((x, y), axis=0)
2602 >>> np.cov(X)
2603 array([[11.71 , -4.286 ], # may vary
2604 [-4.286 , 2.144133]])
2605 >>> np.cov(x, y)
2606 array([[11.71 , -4.286 ], # may vary
2607 [-4.286 , 2.144133]])
2608 >>> np.cov(x)
2609 array(11.71)
2611 """
2612 # Check inputs
2613 if ddof is not None and ddof != int(ddof):
2614 raise ValueError(
2615 "ddof must be integer")
2617 # Handles complex arrays too
2618 m = np.asarray(m)
2619 if m.ndim > 2:
2620 raise ValueError("m has more than 2 dimensions")
2622 if y is not None:
2623 y = np.asarray(y)
2624 if y.ndim > 2:
2625 raise ValueError("y has more than 2 dimensions")
2627 if dtype is None:
2628 if y is None:
2629 dtype = np.result_type(m, np.float64)
2630 else:
2631 dtype = np.result_type(m, y, np.float64)
2633 X = array(m, ndmin=2, dtype=dtype)
2634 if not rowvar and X.shape[0] != 1:
2635 X = X.T
2636 if X.shape[0] == 0:
2637 return np.array([]).reshape(0, 0)
2638 if y is not None:
2639 y = array(y, copy=False, ndmin=2, dtype=dtype)
2640 if not rowvar and y.shape[0] != 1:
2641 y = y.T
2642 X = np.concatenate((X, y), axis=0)
2644 if ddof is None:
2645 if bias == 0:
2646 ddof = 1
2647 else:
2648 ddof = 0
2650 # Get the product of frequencies and weights
2651 w = None
2652 if fweights is not None:
2653 fweights = np.asarray(fweights, dtype=float)
2654 if not np.all(fweights == np.around(fweights)):
2655 raise TypeError(
2656 "fweights must be integer")
2657 if fweights.ndim > 1:
2658 raise RuntimeError(
2659 "cannot handle multidimensional fweights")
2660 if fweights.shape[0] != X.shape[1]:
2661 raise RuntimeError(
2662 "incompatible numbers of samples and fweights")
2663 if any(fweights < 0):
2664 raise ValueError(
2665 "fweights cannot be negative")
2666 w = fweights
2667 if aweights is not None:
2668 aweights = np.asarray(aweights, dtype=float)
2669 if aweights.ndim > 1:
2670 raise RuntimeError(
2671 "cannot handle multidimensional aweights")
2672 if aweights.shape[0] != X.shape[1]:
2673 raise RuntimeError(
2674 "incompatible numbers of samples and aweights")
2675 if any(aweights < 0):
2676 raise ValueError(
2677 "aweights cannot be negative")
2678 if w is None:
2679 w = aweights
2680 else:
2681 w *= aweights
2683 avg, w_sum = average(X, axis=1, weights=w, returned=True)
2684 w_sum = w_sum[0]
2686 # Determine the normalization
2687 if w is None:
2688 fact = X.shape[1] - ddof
2689 elif ddof == 0:
2690 fact = w_sum
2691 elif aweights is None:
2692 fact = w_sum - ddof
2693 else:
2694 fact = w_sum - ddof*sum(w*aweights)/w_sum
2696 if fact <= 0:
2697 warnings.warn("Degrees of freedom <= 0 for slice",
2698 RuntimeWarning, stacklevel=3)
2699 fact = 0.0
2701 X -= avg[:, None]
2702 if w is None:
2703 X_T = X.T
2704 else:
2705 X_T = (X*w).T
2706 c = dot(X, X_T.conj())
2707 c *= np.true_divide(1, fact)
2708 return c.squeeze()
2711def _corrcoef_dispatcher(x, y=None, rowvar=None, bias=None, ddof=None, *,
2712 dtype=None):
2713 return (x, y)
2716@array_function_dispatch(_corrcoef_dispatcher)
2717def corrcoef(x, y=None, rowvar=True, bias=np._NoValue, ddof=np._NoValue, *,
2718 dtype=None):
2719 """
2720 Return Pearson product-moment correlation coefficients.
2722 Please refer to the documentation for `cov` for more detail. The
2723 relationship between the correlation coefficient matrix, `R`, and the
2724 covariance matrix, `C`, is
2726 .. math:: R_{ij} = \\frac{ C_{ij} } { \\sqrt{ C_{ii} C_{jj} } }
2728 The values of `R` are between -1 and 1, inclusive.
2730 Parameters
2731 ----------
2732 x : array_like
2733 A 1-D or 2-D array containing multiple variables and observations.
2734 Each row of `x` represents a variable, and each column a single
2735 observation of all those variables. Also see `rowvar` below.
2736 y : array_like, optional
2737 An additional set of variables and observations. `y` has the same
2738 shape as `x`.
2739 rowvar : bool, optional
2740 If `rowvar` is True (default), then each row represents a
2741 variable, with observations in the columns. Otherwise, the relationship
2742 is transposed: each column represents a variable, while the rows
2743 contain observations.
2744 bias : _NoValue, optional
2745 Has no effect, do not use.
2747 .. deprecated:: 1.10.0
2748 ddof : _NoValue, optional
2749 Has no effect, do not use.
2751 .. deprecated:: 1.10.0
2752 dtype : data-type, optional
2753 Data-type of the result. By default, the return data-type will have
2754 at least `numpy.float64` precision.
2756 .. versionadded:: 1.20
2758 Returns
2759 -------
2760 R : ndarray
2761 The correlation coefficient matrix of the variables.
2763 See Also
2764 --------
2765 cov : Covariance matrix
2767 Notes
2768 -----
2769 Due to floating point rounding the resulting array may not be Hermitian,
2770 the diagonal elements may not be 1, and the elements may not satisfy the
2771 inequality abs(a) <= 1. The real and imaginary parts are clipped to the
2772 interval [-1, 1] in an attempt to improve on that situation but is not
2773 much help in the complex case.
2775 This function accepts but discards arguments `bias` and `ddof`. This is
2776 for backwards compatibility with previous versions of this function. These
2777 arguments had no effect on the return values of the function and can be
2778 safely ignored in this and previous versions of numpy.
2780 Examples
2781 --------
2782 In this example we generate two random arrays, ``xarr`` and ``yarr``, and
2783 compute the row-wise and column-wise Pearson correlation coefficients,
2784 ``R``. Since ``rowvar`` is true by default, we first find the row-wise
2785 Pearson correlation coefficients between the variables of ``xarr``.
2787 >>> import numpy as np
2788 >>> rng = np.random.default_rng(seed=42)
2789 >>> xarr = rng.random((3, 3))
2790 >>> xarr
2791 array([[0.77395605, 0.43887844, 0.85859792],
2792 [0.69736803, 0.09417735, 0.97562235],
2793 [0.7611397 , 0.78606431, 0.12811363]])
2794 >>> R1 = np.corrcoef(xarr)
2795 >>> R1
2796 array([[ 1. , 0.99256089, -0.68080986],
2797 [ 0.99256089, 1. , -0.76492172],
2798 [-0.68080986, -0.76492172, 1. ]])
2800 If we add another set of variables and observations ``yarr``, we can
2801 compute the row-wise Pearson correlation coefficients between the
2802 variables in ``xarr`` and ``yarr``.
2804 >>> yarr = rng.random((3, 3))
2805 >>> yarr
2806 array([[0.45038594, 0.37079802, 0.92676499],
2807 [0.64386512, 0.82276161, 0.4434142 ],
2808 [0.22723872, 0.55458479, 0.06381726]])
2809 >>> R2 = np.corrcoef(xarr, yarr)
2810 >>> R2
2811 array([[ 1. , 0.99256089, -0.68080986, 0.75008178, -0.934284 ,
2812 -0.99004057],
2813 [ 0.99256089, 1. , -0.76492172, 0.82502011, -0.97074098,
2814 -0.99981569],
2815 [-0.68080986, -0.76492172, 1. , -0.99507202, 0.89721355,
2816 0.77714685],
2817 [ 0.75008178, 0.82502011, -0.99507202, 1. , -0.93657855,
2818 -0.83571711],
2819 [-0.934284 , -0.97074098, 0.89721355, -0.93657855, 1. ,
2820 0.97517215],
2821 [-0.99004057, -0.99981569, 0.77714685, -0.83571711, 0.97517215,
2822 1. ]])
2824 Finally if we use the option ``rowvar=False``, the columns are now
2825 being treated as the variables and we will find the column-wise Pearson
2826 correlation coefficients between variables in ``xarr`` and ``yarr``.
2828 >>> R3 = np.corrcoef(xarr, yarr, rowvar=False)
2829 >>> R3
2830 array([[ 1. , 0.77598074, -0.47458546, -0.75078643, -0.9665554 ,
2831 0.22423734],
2832 [ 0.77598074, 1. , -0.92346708, -0.99923895, -0.58826587,
2833 -0.44069024],
2834 [-0.47458546, -0.92346708, 1. , 0.93773029, 0.23297648,
2835 0.75137473],
2836 [-0.75078643, -0.99923895, 0.93773029, 1. , 0.55627469,
2837 0.47536961],
2838 [-0.9665554 , -0.58826587, 0.23297648, 0.55627469, 1. ,
2839 -0.46666491],
2840 [ 0.22423734, -0.44069024, 0.75137473, 0.47536961, -0.46666491,
2841 1. ]])
2843 """
2844 if bias is not np._NoValue or ddof is not np._NoValue:
2845 # 2015-03-15, 1.10
2846 warnings.warn('bias and ddof have no effect and are deprecated',
2847 DeprecationWarning, stacklevel=3)
2848 c = cov(x, y, rowvar, dtype=dtype)
2849 try:
2850 d = diag(c)
2851 except ValueError:
2852 # scalar covariance
2853 # nan if incorrect value (nan, inf, 0), 1 otherwise
2854 return c / c
2855 stddev = sqrt(d.real)
2856 c /= stddev[:, None]
2857 c /= stddev[None, :]
2859 # Clip real and imaginary parts to [-1, 1]. This does not guarantee
2860 # abs(a[i,j]) <= 1 for complex arrays, but is the best we can do without
2861 # excessive work.
2862 np.clip(c.real, -1, 1, out=c.real)
2863 if np.iscomplexobj(c):
2864 np.clip(c.imag, -1, 1, out=c.imag)
2866 return c
2869@set_module('numpy')
2870def blackman(M):
2871 """
2872 Return the Blackman window.
2874 The Blackman window is a taper formed by using the first three
2875 terms of a summation of cosines. It was designed to have close to the
2876 minimal leakage possible. It is close to optimal, only slightly worse
2877 than a Kaiser window.
2879 Parameters
2880 ----------
2881 M : int
2882 Number of points in the output window. If zero or less, an empty
2883 array is returned.
2885 Returns
2886 -------
2887 out : ndarray
2888 The window, with the maximum value normalized to one (the value one
2889 appears only if the number of samples is odd).
2891 See Also
2892 --------
2893 bartlett, hamming, hanning, kaiser
2895 Notes
2896 -----
2897 The Blackman window is defined as
2899 .. math:: w(n) = 0.42 - 0.5 \\cos(2\\pi n/M) + 0.08 \\cos(4\\pi n/M)
2901 Most references to the Blackman window come from the signal processing
2902 literature, where it is used as one of many windowing functions for
2903 smoothing values. It is also known as an apodization (which means
2904 "removing the foot", i.e. smoothing discontinuities at the beginning
2905 and end of the sampled signal) or tapering function. It is known as a
2906 "near optimal" tapering function, almost as good (by some measures)
2907 as the kaiser window.
2909 References
2910 ----------
2911 Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra,
2912 Dover Publications, New York.
2914 Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing.
2915 Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471.
2917 Examples
2918 --------
2919 >>> import matplotlib.pyplot as plt
2920 >>> np.blackman(12)
2921 array([-1.38777878e-17, 3.26064346e-02, 1.59903635e-01, # may vary
2922 4.14397981e-01, 7.36045180e-01, 9.67046769e-01,
2923 9.67046769e-01, 7.36045180e-01, 4.14397981e-01,
2924 1.59903635e-01, 3.26064346e-02, -1.38777878e-17])
2926 Plot the window and the frequency response:
2928 >>> from numpy.fft import fft, fftshift
2929 >>> window = np.blackman(51)
2930 >>> plt.plot(window)
2931 [<matplotlib.lines.Line2D object at 0x...>]
2932 >>> plt.title("Blackman window")
2933 Text(0.5, 1.0, 'Blackman window')
2934 >>> plt.ylabel("Amplitude")
2935 Text(0, 0.5, 'Amplitude')
2936 >>> plt.xlabel("Sample")
2937 Text(0.5, 0, 'Sample')
2938 >>> plt.show()
2940 >>> plt.figure()
2941 <Figure size 640x480 with 0 Axes>
2942 >>> A = fft(window, 2048) / 25.5
2943 >>> mag = np.abs(fftshift(A))
2944 >>> freq = np.linspace(-0.5, 0.5, len(A))
2945 >>> with np.errstate(divide='ignore', invalid='ignore'):
2946 ... response = 20 * np.log10(mag)
2947 ...
2948 >>> response = np.clip(response, -100, 100)
2949 >>> plt.plot(freq, response)
2950 [<matplotlib.lines.Line2D object at 0x...>]
2951 >>> plt.title("Frequency response of Blackman window")
2952 Text(0.5, 1.0, 'Frequency response of Blackman window')
2953 >>> plt.ylabel("Magnitude [dB]")
2954 Text(0, 0.5, 'Magnitude [dB]')
2955 >>> plt.xlabel("Normalized frequency [cycles per sample]")
2956 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
2957 >>> _ = plt.axis('tight')
2958 >>> plt.show()
2960 """
2961 if M < 1:
2962 return array([], dtype=np.result_type(M, 0.0))
2963 if M == 1:
2964 return ones(1, dtype=np.result_type(M, 0.0))
2965 n = arange(1-M, M, 2)
2966 return 0.42 + 0.5*cos(pi*n/(M-1)) + 0.08*cos(2.0*pi*n/(M-1))
2969@set_module('numpy')
2970def bartlett(M):
2971 """
2972 Return the Bartlett window.
2974 The Bartlett window is very similar to a triangular window, except
2975 that the end points are at zero. It is often used in signal
2976 processing for tapering a signal, without generating too much
2977 ripple in the frequency domain.
2979 Parameters
2980 ----------
2981 M : int
2982 Number of points in the output window. If zero or less, an
2983 empty array is returned.
2985 Returns
2986 -------
2987 out : array
2988 The triangular window, with the maximum value normalized to one
2989 (the value one appears only if the number of samples is odd), with
2990 the first and last samples equal to zero.
2992 See Also
2993 --------
2994 blackman, hamming, hanning, kaiser
2996 Notes
2997 -----
2998 The Bartlett window is defined as
3000 .. math:: w(n) = \\frac{2}{M-1} \\left(
3001 \\frac{M-1}{2} - \\left|n - \\frac{M-1}{2}\\right|
3002 \\right)
3004 Most references to the Bartlett window come from the signal processing
3005 literature, where it is used as one of many windowing functions for
3006 smoothing values. Note that convolution with this window produces linear
3007 interpolation. It is also known as an apodization (which means "removing
3008 the foot", i.e. smoothing discontinuities at the beginning and end of the
3009 sampled signal) or tapering function. The Fourier transform of the
3010 Bartlett window is the product of two sinc functions. Note the excellent
3011 discussion in Kanasewich [2]_.
3013 References
3014 ----------
3015 .. [1] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra",
3016 Biometrika 37, 1-16, 1950.
3017 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
3018 The University of Alberta Press, 1975, pp. 109-110.
3019 .. [3] A.V. Oppenheim and R.W. Schafer, "Discrete-Time Signal
3020 Processing", Prentice-Hall, 1999, pp. 468-471.
3021 .. [4] Wikipedia, "Window function",
3022 https://en.wikipedia.org/wiki/Window_function
3023 .. [5] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3024 "Numerical Recipes", Cambridge University Press, 1986, page 429.
3026 Examples
3027 --------
3028 >>> import matplotlib.pyplot as plt
3029 >>> np.bartlett(12)
3030 array([ 0. , 0.18181818, 0.36363636, 0.54545455, 0.72727273, # may vary
3031 0.90909091, 0.90909091, 0.72727273, 0.54545455, 0.36363636,
3032 0.18181818, 0. ])
3034 Plot the window and its frequency response (requires SciPy and matplotlib):
3036 >>> from numpy.fft import fft, fftshift
3037 >>> window = np.bartlett(51)
3038 >>> plt.plot(window)
3039 [<matplotlib.lines.Line2D object at 0x...>]
3040 >>> plt.title("Bartlett window")
3041 Text(0.5, 1.0, 'Bartlett window')
3042 >>> plt.ylabel("Amplitude")
3043 Text(0, 0.5, 'Amplitude')
3044 >>> plt.xlabel("Sample")
3045 Text(0.5, 0, 'Sample')
3046 >>> plt.show()
3048 >>> plt.figure()
3049 <Figure size 640x480 with 0 Axes>
3050 >>> A = fft(window, 2048) / 25.5
3051 >>> mag = np.abs(fftshift(A))
3052 >>> freq = np.linspace(-0.5, 0.5, len(A))
3053 >>> with np.errstate(divide='ignore', invalid='ignore'):
3054 ... response = 20 * np.log10(mag)
3055 ...
3056 >>> response = np.clip(response, -100, 100)
3057 >>> plt.plot(freq, response)
3058 [<matplotlib.lines.Line2D object at 0x...>]
3059 >>> plt.title("Frequency response of Bartlett window")
3060 Text(0.5, 1.0, 'Frequency response of Bartlett window')
3061 >>> plt.ylabel("Magnitude [dB]")
3062 Text(0, 0.5, 'Magnitude [dB]')
3063 >>> plt.xlabel("Normalized frequency [cycles per sample]")
3064 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
3065 >>> _ = plt.axis('tight')
3066 >>> plt.show()
3068 """
3069 if M < 1:
3070 return array([], dtype=np.result_type(M, 0.0))
3071 if M == 1:
3072 return ones(1, dtype=np.result_type(M, 0.0))
3073 n = arange(1-M, M, 2)
3074 return where(less_equal(n, 0), 1 + n/(M-1), 1 - n/(M-1))
3077@set_module('numpy')
3078def hanning(M):
3079 """
3080 Return the Hanning window.
3082 The Hanning window is a taper formed by using a weighted cosine.
3084 Parameters
3085 ----------
3086 M : int
3087 Number of points in the output window. If zero or less, an
3088 empty array is returned.
3090 Returns
3091 -------
3092 out : ndarray, shape(M,)
3093 The window, with the maximum value normalized to one (the value
3094 one appears only if `M` is odd).
3096 See Also
3097 --------
3098 bartlett, blackman, hamming, kaiser
3100 Notes
3101 -----
3102 The Hanning window is defined as
3104 .. math:: w(n) = 0.5 - 0.5\\cos\\left(\\frac{2\\pi{n}}{M-1}\\right)
3105 \\qquad 0 \\leq n \\leq M-1
3107 The Hanning was named for Julius von Hann, an Austrian meteorologist.
3108 It is also known as the Cosine Bell. Some authors prefer that it be
3109 called a Hann window, to help avoid confusion with the very similar
3110 Hamming window.
3112 Most references to the Hanning window come from the signal processing
3113 literature, where it is used as one of many windowing functions for
3114 smoothing values. It is also known as an apodization (which means
3115 "removing the foot", i.e. smoothing discontinuities at the beginning
3116 and end of the sampled signal) or tapering function.
3118 References
3119 ----------
3120 .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
3121 spectra, Dover Publications, New York.
3122 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
3123 The University of Alberta Press, 1975, pp. 106-108.
3124 .. [3] Wikipedia, "Window function",
3125 https://en.wikipedia.org/wiki/Window_function
3126 .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3127 "Numerical Recipes", Cambridge University Press, 1986, page 425.
3129 Examples
3130 --------
3131 >>> np.hanning(12)
3132 array([0. , 0.07937323, 0.29229249, 0.57115742, 0.82743037,
3133 0.97974649, 0.97974649, 0.82743037, 0.57115742, 0.29229249,
3134 0.07937323, 0. ])
3136 Plot the window and its frequency response:
3138 >>> import matplotlib.pyplot as plt
3139 >>> from numpy.fft import fft, fftshift
3140 >>> window = np.hanning(51)
3141 >>> plt.plot(window)
3142 [<matplotlib.lines.Line2D object at 0x...>]
3143 >>> plt.title("Hann window")
3144 Text(0.5, 1.0, 'Hann window')
3145 >>> plt.ylabel("Amplitude")
3146 Text(0, 0.5, 'Amplitude')
3147 >>> plt.xlabel("Sample")
3148 Text(0.5, 0, 'Sample')
3149 >>> plt.show()
3151 >>> plt.figure()
3152 <Figure size 640x480 with 0 Axes>
3153 >>> A = fft(window, 2048) / 25.5
3154 >>> mag = np.abs(fftshift(A))
3155 >>> freq = np.linspace(-0.5, 0.5, len(A))
3156 >>> with np.errstate(divide='ignore', invalid='ignore'):
3157 ... response = 20 * np.log10(mag)
3158 ...
3159 >>> response = np.clip(response, -100, 100)
3160 >>> plt.plot(freq, response)
3161 [<matplotlib.lines.Line2D object at 0x...>]
3162 >>> plt.title("Frequency response of the Hann window")
3163 Text(0.5, 1.0, 'Frequency response of the Hann window')
3164 >>> plt.ylabel("Magnitude [dB]")
3165 Text(0, 0.5, 'Magnitude [dB]')
3166 >>> plt.xlabel("Normalized frequency [cycles per sample]")
3167 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
3168 >>> plt.axis('tight')
3169 ...
3170 >>> plt.show()
3172 """
3173 if M < 1:
3174 return array([], dtype=np.result_type(M, 0.0))
3175 if M == 1:
3176 return ones(1, dtype=np.result_type(M, 0.0))
3177 n = arange(1-M, M, 2)
3178 return 0.5 + 0.5*cos(pi*n/(M-1))
3181@set_module('numpy')
3182def hamming(M):
3183 """
3184 Return the Hamming window.
3186 The Hamming window is a taper formed by using a weighted cosine.
3188 Parameters
3189 ----------
3190 M : int
3191 Number of points in the output window. If zero or less, an
3192 empty array is returned.
3194 Returns
3195 -------
3196 out : ndarray
3197 The window, with the maximum value normalized to one (the value
3198 one appears only if the number of samples is odd).
3200 See Also
3201 --------
3202 bartlett, blackman, hanning, kaiser
3204 Notes
3205 -----
3206 The Hamming window is defined as
3208 .. math:: w(n) = 0.54 - 0.46\\cos\\left(\\frac{2\\pi{n}}{M-1}\\right)
3209 \\qquad 0 \\leq n \\leq M-1
3211 The Hamming was named for R. W. Hamming, an associate of J. W. Tukey
3212 and is described in Blackman and Tukey. It was recommended for
3213 smoothing the truncated autocovariance function in the time domain.
3214 Most references to the Hamming window come from the signal processing
3215 literature, where it is used as one of many windowing functions for
3216 smoothing values. It is also known as an apodization (which means
3217 "removing the foot", i.e. smoothing discontinuities at the beginning
3218 and end of the sampled signal) or tapering function.
3220 References
3221 ----------
3222 .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
3223 spectra, Dover Publications, New York.
3224 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
3225 University of Alberta Press, 1975, pp. 109-110.
3226 .. [3] Wikipedia, "Window function",
3227 https://en.wikipedia.org/wiki/Window_function
3228 .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3229 "Numerical Recipes", Cambridge University Press, 1986, page 425.
3231 Examples
3232 --------
3233 >>> np.hamming(12)
3234 array([ 0.08 , 0.15302337, 0.34890909, 0.60546483, 0.84123594, # may vary
3235 0.98136677, 0.98136677, 0.84123594, 0.60546483, 0.34890909,
3236 0.15302337, 0.08 ])
3238 Plot the window and the frequency response:
3240 >>> import matplotlib.pyplot as plt
3241 >>> from numpy.fft import fft, fftshift
3242 >>> window = np.hamming(51)
3243 >>> plt.plot(window)
3244 [<matplotlib.lines.Line2D object at 0x...>]
3245 >>> plt.title("Hamming window")
3246 Text(0.5, 1.0, 'Hamming window')
3247 >>> plt.ylabel("Amplitude")
3248 Text(0, 0.5, 'Amplitude')
3249 >>> plt.xlabel("Sample")
3250 Text(0.5, 0, 'Sample')
3251 >>> plt.show()
3253 >>> plt.figure()
3254 <Figure size 640x480 with 0 Axes>
3255 >>> A = fft(window, 2048) / 25.5
3256 >>> mag = np.abs(fftshift(A))
3257 >>> freq = np.linspace(-0.5, 0.5, len(A))
3258 >>> response = 20 * np.log10(mag)
3259 >>> response = np.clip(response, -100, 100)
3260 >>> plt.plot(freq, response)
3261 [<matplotlib.lines.Line2D object at 0x...>]
3262 >>> plt.title("Frequency response of Hamming window")
3263 Text(0.5, 1.0, 'Frequency response of Hamming window')
3264 >>> plt.ylabel("Magnitude [dB]")
3265 Text(0, 0.5, 'Magnitude [dB]')
3266 >>> plt.xlabel("Normalized frequency [cycles per sample]")
3267 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
3268 >>> plt.axis('tight')
3269 ...
3270 >>> plt.show()
3272 """
3273 if M < 1:
3274 return array([], dtype=np.result_type(M, 0.0))
3275 if M == 1:
3276 return ones(1, dtype=np.result_type(M, 0.0))
3277 n = arange(1-M, M, 2)
3278 return 0.54 + 0.46*cos(pi*n/(M-1))
3281## Code from cephes for i0
3283_i0A = [
3284 -4.41534164647933937950E-18,
3285 3.33079451882223809783E-17,
3286 -2.43127984654795469359E-16,
3287 1.71539128555513303061E-15,
3288 -1.16853328779934516808E-14,
3289 7.67618549860493561688E-14,
3290 -4.85644678311192946090E-13,
3291 2.95505266312963983461E-12,
3292 -1.72682629144155570723E-11,
3293 9.67580903537323691224E-11,
3294 -5.18979560163526290666E-10,
3295 2.65982372468238665035E-9,
3296 -1.30002500998624804212E-8,
3297 6.04699502254191894932E-8,
3298 -2.67079385394061173391E-7,
3299 1.11738753912010371815E-6,
3300 -4.41673835845875056359E-6,
3301 1.64484480707288970893E-5,
3302 -5.75419501008210370398E-5,
3303 1.88502885095841655729E-4,
3304 -5.76375574538582365885E-4,
3305 1.63947561694133579842E-3,
3306 -4.32430999505057594430E-3,
3307 1.05464603945949983183E-2,
3308 -2.37374148058994688156E-2,
3309 4.93052842396707084878E-2,
3310 -9.49010970480476444210E-2,
3311 1.71620901522208775349E-1,
3312 -3.04682672343198398683E-1,
3313 6.76795274409476084995E-1
3314 ]
3316_i0B = [
3317 -7.23318048787475395456E-18,
3318 -4.83050448594418207126E-18,
3319 4.46562142029675999901E-17,
3320 3.46122286769746109310E-17,
3321 -2.82762398051658348494E-16,
3322 -3.42548561967721913462E-16,
3323 1.77256013305652638360E-15,
3324 3.81168066935262242075E-15,
3325 -9.55484669882830764870E-15,
3326 -4.15056934728722208663E-14,
3327 1.54008621752140982691E-14,
3328 3.85277838274214270114E-13,
3329 7.18012445138366623367E-13,
3330 -1.79417853150680611778E-12,
3331 -1.32158118404477131188E-11,
3332 -3.14991652796324136454E-11,
3333 1.18891471078464383424E-11,
3334 4.94060238822496958910E-10,
3335 3.39623202570838634515E-9,
3336 2.26666899049817806459E-8,
3337 2.04891858946906374183E-7,
3338 2.89137052083475648297E-6,
3339 6.88975834691682398426E-5,
3340 3.36911647825569408990E-3,
3341 8.04490411014108831608E-1
3342 ]
3345def _chbevl(x, vals):
3346 b0 = vals[0]
3347 b1 = 0.0
3349 for i in range(1, len(vals)):
3350 b2 = b1
3351 b1 = b0
3352 b0 = x*b1 - b2 + vals[i]
3354 return 0.5*(b0 - b2)
3357def _i0_1(x):
3358 return exp(x) * _chbevl(x/2.0-2, _i0A)
3361def _i0_2(x):
3362 return exp(x) * _chbevl(32.0/x - 2.0, _i0B) / sqrt(x)
3365def _i0_dispatcher(x):
3366 return (x,)
3369@array_function_dispatch(_i0_dispatcher)
3370def i0(x):
3371 """
3372 Modified Bessel function of the first kind, order 0.
3374 Usually denoted :math:`I_0`.
3376 Parameters
3377 ----------
3378 x : array_like of float
3379 Argument of the Bessel function.
3381 Returns
3382 -------
3383 out : ndarray, shape = x.shape, dtype = float
3384 The modified Bessel function evaluated at each of the elements of `x`.
3386 See Also
3387 --------
3388 scipy.special.i0, scipy.special.iv, scipy.special.ive
3390 Notes
3391 -----
3392 The scipy implementation is recommended over this function: it is a
3393 proper ufunc written in C, and more than an order of magnitude faster.
3395 We use the algorithm published by Clenshaw [1]_ and referenced by
3396 Abramowitz and Stegun [2]_, for which the function domain is
3397 partitioned into the two intervals [0,8] and (8,inf), and Chebyshev
3398 polynomial expansions are employed in each interval. Relative error on
3399 the domain [0,30] using IEEE arithmetic is documented [3]_ as having a
3400 peak of 5.8e-16 with an rms of 1.4e-16 (n = 30000).
3402 References
3403 ----------
3404 .. [1] C. W. Clenshaw, "Chebyshev series for mathematical functions", in
3405 *National Physical Laboratory Mathematical Tables*, vol. 5, London:
3406 Her Majesty's Stationery Office, 1962.
3407 .. [2] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical
3408 Functions*, 10th printing, New York: Dover, 1964, pp. 379.
3409 https://personal.math.ubc.ca/~cbm/aands/page_379.htm
3410 .. [3] https://metacpan.org/pod/distribution/Math-Cephes/lib/Math/Cephes.pod#i0:-Modified-Bessel-function-of-order-zero
3412 Examples
3413 --------
3414 >>> np.i0(0.)
3415 array(1.0)
3416 >>> np.i0([0, 1, 2, 3])
3417 array([1. , 1.26606588, 2.2795853 , 4.88079259])
3419 """
3420 x = np.asanyarray(x)
3421 if x.dtype.kind == 'c':
3422 raise TypeError("i0 not supported for complex values")
3423 if x.dtype.kind != 'f':
3424 x = x.astype(float)
3425 x = np.abs(x)
3426 return piecewise(x, [x <= 8.0], [_i0_1, _i0_2])
3428## End of cephes code for i0
3431@set_module('numpy')
3432def kaiser(M, beta):
3433 """
3434 Return the Kaiser window.
3436 The Kaiser window is a taper formed by using a Bessel function.
3438 Parameters
3439 ----------
3440 M : int
3441 Number of points in the output window. If zero or less, an
3442 empty array is returned.
3443 beta : float
3444 Shape parameter for window.
3446 Returns
3447 -------
3448 out : array
3449 The window, with the maximum value normalized to one (the value
3450 one appears only if the number of samples is odd).
3452 See Also
3453 --------
3454 bartlett, blackman, hamming, hanning
3456 Notes
3457 -----
3458 The Kaiser window is defined as
3460 .. math:: w(n) = I_0\\left( \\beta \\sqrt{1-\\frac{4n^2}{(M-1)^2}}
3461 \\right)/I_0(\\beta)
3463 with
3465 .. math:: \\quad -\\frac{M-1}{2} \\leq n \\leq \\frac{M-1}{2},
3467 where :math:`I_0` is the modified zeroth-order Bessel function.
3469 The Kaiser was named for Jim Kaiser, who discovered a simple
3470 approximation to the DPSS window based on Bessel functions. The Kaiser
3471 window is a very good approximation to the Digital Prolate Spheroidal
3472 Sequence, or Slepian window, which is the transform which maximizes the
3473 energy in the main lobe of the window relative to total energy.
3475 The Kaiser can approximate many other windows by varying the beta
3476 parameter.
3478 ==== =======================
3479 beta Window shape
3480 ==== =======================
3481 0 Rectangular
3482 5 Similar to a Hamming
3483 6 Similar to a Hanning
3484 8.6 Similar to a Blackman
3485 ==== =======================
3487 A beta value of 14 is probably a good starting point. Note that as beta
3488 gets large, the window narrows, and so the number of samples needs to be
3489 large enough to sample the increasingly narrow spike, otherwise NaNs will
3490 get returned.
3492 Most references to the Kaiser window come from the signal processing
3493 literature, where it is used as one of many windowing functions for
3494 smoothing values. It is also known as an apodization (which means
3495 "removing the foot", i.e. smoothing discontinuities at the beginning
3496 and end of the sampled signal) or tapering function.
3498 References
3499 ----------
3500 .. [1] J. F. Kaiser, "Digital Filters" - Ch 7 in "Systems analysis by
3501 digital computer", Editors: F.F. Kuo and J.F. Kaiser, p 218-285.
3502 John Wiley and Sons, New York, (1966).
3503 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
3504 University of Alberta Press, 1975, pp. 177-178.
3505 .. [3] Wikipedia, "Window function",
3506 https://en.wikipedia.org/wiki/Window_function
3508 Examples
3509 --------
3510 >>> import matplotlib.pyplot as plt
3511 >>> np.kaiser(12, 14)
3512 array([7.72686684e-06, 3.46009194e-03, 4.65200189e-02, # may vary
3513 2.29737120e-01, 5.99885316e-01, 9.45674898e-01,
3514 9.45674898e-01, 5.99885316e-01, 2.29737120e-01,
3515 4.65200189e-02, 3.46009194e-03, 7.72686684e-06])
3518 Plot the window and the frequency response:
3520 >>> from numpy.fft import fft, fftshift
3521 >>> window = np.kaiser(51, 14)
3522 >>> plt.plot(window)
3523 [<matplotlib.lines.Line2D object at 0x...>]
3524 >>> plt.title("Kaiser window")
3525 Text(0.5, 1.0, 'Kaiser window')
3526 >>> plt.ylabel("Amplitude")
3527 Text(0, 0.5, 'Amplitude')
3528 >>> plt.xlabel("Sample")
3529 Text(0.5, 0, 'Sample')
3530 >>> plt.show()
3532 >>> plt.figure()
3533 <Figure size 640x480 with 0 Axes>
3534 >>> A = fft(window, 2048) / 25.5
3535 >>> mag = np.abs(fftshift(A))
3536 >>> freq = np.linspace(-0.5, 0.5, len(A))
3537 >>> response = 20 * np.log10(mag)
3538 >>> response = np.clip(response, -100, 100)
3539 >>> plt.plot(freq, response)
3540 [<matplotlib.lines.Line2D object at 0x...>]
3541 >>> plt.title("Frequency response of Kaiser window")
3542 Text(0.5, 1.0, 'Frequency response of Kaiser window')
3543 >>> plt.ylabel("Magnitude [dB]")
3544 Text(0, 0.5, 'Magnitude [dB]')
3545 >>> plt.xlabel("Normalized frequency [cycles per sample]")
3546 Text(0.5, 0, 'Normalized frequency [cycles per sample]')
3547 >>> plt.axis('tight')
3548 (-0.5, 0.5, -100.0, ...) # may vary
3549 >>> plt.show()
3551 """
3552 if M == 1:
3553 return np.ones(1, dtype=np.result_type(M, 0.0))
3554 n = arange(0, M)
3555 alpha = (M-1)/2.0
3556 return i0(beta * sqrt(1-((n-alpha)/alpha)**2.0))/i0(float(beta))
3559def _sinc_dispatcher(x):
3560 return (x,)
3563@array_function_dispatch(_sinc_dispatcher)
3564def sinc(x):
3565 r"""
3566 Return the normalized sinc function.
3568 The sinc function is equal to :math:`\sin(\pi x)/(\pi x)` for any argument
3569 :math:`x\ne 0`. ``sinc(0)`` takes the limit value 1, making ``sinc`` not
3570 only everywhere continuous but also infinitely differentiable.
3572 .. note::
3574 Note the normalization factor of ``pi`` used in the definition.
3575 This is the most commonly used definition in signal processing.
3576 Use ``sinc(x / np.pi)`` to obtain the unnormalized sinc function
3577 :math:`\sin(x)/x` that is more common in mathematics.
3579 Parameters
3580 ----------
3581 x : ndarray
3582 Array (possibly multi-dimensional) of values for which to calculate
3583 ``sinc(x)``.
3585 Returns
3586 -------
3587 out : ndarray
3588 ``sinc(x)``, which has the same shape as the input.
3590 Notes
3591 -----
3592 The name sinc is short for "sine cardinal" or "sinus cardinalis".
3594 The sinc function is used in various signal processing applications,
3595 including in anti-aliasing, in the construction of a Lanczos resampling
3596 filter, and in interpolation.
3598 For bandlimited interpolation of discrete-time signals, the ideal
3599 interpolation kernel is proportional to the sinc function.
3601 References
3602 ----------
3603 .. [1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram Web
3604 Resource. http://mathworld.wolfram.com/SincFunction.html
3605 .. [2] Wikipedia, "Sinc function",
3606 https://en.wikipedia.org/wiki/Sinc_function
3608 Examples
3609 --------
3610 >>> import matplotlib.pyplot as plt
3611 >>> x = np.linspace(-4, 4, 41)
3612 >>> np.sinc(x)
3613 array([-3.89804309e-17, -4.92362781e-02, -8.40918587e-02, # may vary
3614 -8.90384387e-02, -5.84680802e-02, 3.89804309e-17,
3615 6.68206631e-02, 1.16434881e-01, 1.26137788e-01,
3616 8.50444803e-02, -3.89804309e-17, -1.03943254e-01,
3617 -1.89206682e-01, -2.16236208e-01, -1.55914881e-01,
3618 3.89804309e-17, 2.33872321e-01, 5.04551152e-01,
3619 7.56826729e-01, 9.35489284e-01, 1.00000000e+00,
3620 9.35489284e-01, 7.56826729e-01, 5.04551152e-01,
3621 2.33872321e-01, 3.89804309e-17, -1.55914881e-01,
3622 -2.16236208e-01, -1.89206682e-01, -1.03943254e-01,
3623 -3.89804309e-17, 8.50444803e-02, 1.26137788e-01,
3624 1.16434881e-01, 6.68206631e-02, 3.89804309e-17,
3625 -5.84680802e-02, -8.90384387e-02, -8.40918587e-02,
3626 -4.92362781e-02, -3.89804309e-17])
3628 >>> plt.plot(x, np.sinc(x))
3629 [<matplotlib.lines.Line2D object at 0x...>]
3630 >>> plt.title("Sinc Function")
3631 Text(0.5, 1.0, 'Sinc Function')
3632 >>> plt.ylabel("Amplitude")
3633 Text(0, 0.5, 'Amplitude')
3634 >>> plt.xlabel("X")
3635 Text(0.5, 0, 'X')
3636 >>> plt.show()
3638 """
3639 x = np.asanyarray(x)
3640 y = pi * where(x == 0, 1.0e-20, x)
3641 return sin(y)/y
3644def _msort_dispatcher(a):
3645 return (a,)
3648@array_function_dispatch(_msort_dispatcher)
3649def msort(a):
3650 """
3651 Return a copy of an array sorted along the first axis.
3653 .. deprecated:: 1.24
3655 msort is deprecated, use ``np.sort(a, axis=0)`` instead.
3657 Parameters
3658 ----------
3659 a : array_like
3660 Array to be sorted.
3662 Returns
3663 -------
3664 sorted_array : ndarray
3665 Array of the same type and shape as `a`.
3667 See Also
3668 --------
3669 sort
3671 Notes
3672 -----
3673 ``np.msort(a)`` is equivalent to ``np.sort(a, axis=0)``.
3675 Examples
3676 --------
3677 >>> a = np.array([[1, 4], [3, 1]])
3678 >>> np.msort(a) # sort along the first axis
3679 array([[1, 1],
3680 [3, 4]])
3682 """
3683 # 2022-10-20 1.24
3684 warnings.warn(
3685 "msort is deprecated, use np.sort(a, axis=0) instead",
3686 DeprecationWarning,
3687 stacklevel=3,
3688 )
3689 b = array(a, subok=True, copy=True)
3690 b.sort(0)
3691 return b
3694def _ureduce(a, func, keepdims=False, **kwargs):
3695 """
3696 Internal Function.
3697 Call `func` with `a` as first argument swapping the axes to use extended
3698 axis on functions that don't support it natively.
3700 Returns result and a.shape with axis dims set to 1.
3702 Parameters
3703 ----------
3704 a : array_like
3705 Input array or object that can be converted to an array.
3706 func : callable
3707 Reduction function capable of receiving a single axis argument.
3708 It is called with `a` as first argument followed by `kwargs`.
3709 kwargs : keyword arguments
3710 additional keyword arguments to pass to `func`.
3712 Returns
3713 -------
3714 result : tuple
3715 Result of func(a, **kwargs) and a.shape with axis dims set to 1
3716 which can be used to reshape the result to the same shape a ufunc with
3717 keepdims=True would produce.
3719 """
3720 a = np.asanyarray(a)
3721 axis = kwargs.get('axis', None)
3722 out = kwargs.get('out', None)
3724 if keepdims is np._NoValue:
3725 keepdims = False
3727 nd = a.ndim
3728 if axis is not None:
3729 axis = _nx.normalize_axis_tuple(axis, nd)
3731 if keepdims:
3732 if out is not None:
3733 index_out = tuple(
3734 0 if i in axis else slice(None) for i in range(nd))
3735 kwargs['out'] = out[(Ellipsis, ) + index_out]
3737 if len(axis) == 1:
3738 kwargs['axis'] = axis[0]
3739 else:
3740 keep = set(range(nd)) - set(axis)
3741 nkeep = len(keep)
3742 # swap axis that should not be reduced to front
3743 for i, s in enumerate(sorted(keep)):
3744 a = a.swapaxes(i, s)
3745 # merge reduced axis
3746 a = a.reshape(a.shape[:nkeep] + (-1,))
3747 kwargs['axis'] = -1
3748 else:
3749 if keepdims:
3750 if out is not None:
3751 index_out = (0, ) * nd
3752 kwargs['out'] = out[(Ellipsis, ) + index_out]
3754 r = func(a, **kwargs)
3756 if out is not None:
3757 return out
3759 if keepdims:
3760 if axis is None:
3761 index_r = (np.newaxis, ) * nd
3762 else:
3763 index_r = tuple(
3764 np.newaxis if i in axis else slice(None)
3765 for i in range(nd))
3766 r = r[(Ellipsis, ) + index_r]
3768 return r
3771def _median_dispatcher(
3772 a, axis=None, out=None, overwrite_input=None, keepdims=None):
3773 return (a, out)
3776@array_function_dispatch(_median_dispatcher)
3777def median(a, axis=None, out=None, overwrite_input=False, keepdims=False):
3778 """
3779 Compute the median along the specified axis.
3781 Returns the median of the array elements.
3783 Parameters
3784 ----------
3785 a : array_like
3786 Input array or object that can be converted to an array.
3787 axis : {int, sequence of int, None}, optional
3788 Axis or axes along which the medians are computed. The default
3789 is to compute the median along a flattened version of the array.
3790 A sequence of axes is supported since version 1.9.0.
3791 out : ndarray, optional
3792 Alternative output array in which to place the result. It must
3793 have the same shape and buffer length as the expected output,
3794 but the type (of the output) will be cast if necessary.
3795 overwrite_input : bool, optional
3796 If True, then allow use of memory of input array `a` for
3797 calculations. The input array will be modified by the call to
3798 `median`. This will save memory when you do not need to preserve
3799 the contents of the input array. Treat the input as undefined,
3800 but it will probably be fully or partially sorted. Default is
3801 False. If `overwrite_input` is ``True`` and `a` is not already an
3802 `ndarray`, an error will be raised.
3803 keepdims : bool, optional
3804 If this is set to True, the axes which are reduced are left
3805 in the result as dimensions with size one. With this option,
3806 the result will broadcast correctly against the original `arr`.
3808 .. versionadded:: 1.9.0
3810 Returns
3811 -------
3812 median : ndarray
3813 A new array holding the result. If the input contains integers
3814 or floats smaller than ``float64``, then the output data-type is
3815 ``np.float64``. Otherwise, the data-type of the output is the
3816 same as that of the input. If `out` is specified, that array is
3817 returned instead.
3819 See Also
3820 --------
3821 mean, percentile
3823 Notes
3824 -----
3825 Given a vector ``V`` of length ``N``, the median of ``V`` is the
3826 middle value of a sorted copy of ``V``, ``V_sorted`` - i
3827 e., ``V_sorted[(N-1)/2]``, when ``N`` is odd, and the average of the
3828 two middle values of ``V_sorted`` when ``N`` is even.
3830 Examples
3831 --------
3832 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
3833 >>> a
3834 array([[10, 7, 4],
3835 [ 3, 2, 1]])
3836 >>> np.median(a)
3837 3.5
3838 >>> np.median(a, axis=0)
3839 array([6.5, 4.5, 2.5])
3840 >>> np.median(a, axis=1)
3841 array([7., 2.])
3842 >>> m = np.median(a, axis=0)
3843 >>> out = np.zeros_like(m)
3844 >>> np.median(a, axis=0, out=m)
3845 array([6.5, 4.5, 2.5])
3846 >>> m
3847 array([6.5, 4.5, 2.5])
3848 >>> b = a.copy()
3849 >>> np.median(b, axis=1, overwrite_input=True)
3850 array([7., 2.])
3851 >>> assert not np.all(a==b)
3852 >>> b = a.copy()
3853 >>> np.median(b, axis=None, overwrite_input=True)
3854 3.5
3855 >>> assert not np.all(a==b)
3857 """
3858 return _ureduce(a, func=_median, keepdims=keepdims, axis=axis, out=out,
3859 overwrite_input=overwrite_input)
3862def _median(a, axis=None, out=None, overwrite_input=False):
3863 # can't be reasonably be implemented in terms of percentile as we have to
3864 # call mean to not break astropy
3865 a = np.asanyarray(a)
3867 # Set the partition indexes
3868 if axis is None:
3869 sz = a.size
3870 else:
3871 sz = a.shape[axis]
3872 if sz % 2 == 0:
3873 szh = sz // 2
3874 kth = [szh - 1, szh]
3875 else:
3876 kth = [(sz - 1) // 2]
3877 # Check if the array contains any nan's
3878 if np.issubdtype(a.dtype, np.inexact):
3879 kth.append(-1)
3881 if overwrite_input:
3882 if axis is None:
3883 part = a.ravel()
3884 part.partition(kth)
3885 else:
3886 a.partition(kth, axis=axis)
3887 part = a
3888 else:
3889 part = partition(a, kth, axis=axis)
3891 if part.shape == ():
3892 # make 0-D arrays work
3893 return part.item()
3894 if axis is None:
3895 axis = 0
3897 indexer = [slice(None)] * part.ndim
3898 index = part.shape[axis] // 2
3899 if part.shape[axis] % 2 == 1:
3900 # index with slice to allow mean (below) to work
3901 indexer[axis] = slice(index, index+1)
3902 else:
3903 indexer[axis] = slice(index-1, index+1)
3904 indexer = tuple(indexer)
3906 # Use mean in both odd and even case to coerce data type,
3907 # using out array if needed.
3908 rout = mean(part[indexer], axis=axis, out=out)
3909 # Check if the array contains any nan's
3910 if np.issubdtype(a.dtype, np.inexact) and sz > 0:
3911 # If nans are possible, warn and replace by nans like mean would.
3912 rout = np.lib.utils._median_nancheck(part, rout, axis)
3914 return rout
3917def _percentile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
3918 method=None, keepdims=None, *, interpolation=None):
3919 return (a, q, out)
3922@array_function_dispatch(_percentile_dispatcher)
3923def percentile(a,
3924 q,
3925 axis=None,
3926 out=None,
3927 overwrite_input=False,
3928 method="linear",
3929 keepdims=False,
3930 *,
3931 interpolation=None):
3932 """
3933 Compute the q-th percentile of the data along the specified axis.
3935 Returns the q-th percentile(s) of the array elements.
3937 Parameters
3938 ----------
3939 a : array_like of real numbers
3940 Input array or object that can be converted to an array.
3941 q : array_like of float
3942 Percentile or sequence of percentiles to compute, which must be between
3943 0 and 100 inclusive.
3944 axis : {int, tuple of int, None}, optional
3945 Axis or axes along which the percentiles are computed. The
3946 default is to compute the percentile(s) along a flattened
3947 version of the array.
3949 .. versionchanged:: 1.9.0
3950 A tuple of axes is supported
3951 out : ndarray, optional
3952 Alternative output array in which to place the result. It must
3953 have the same shape and buffer length as the expected output,
3954 but the type (of the output) will be cast if necessary.
3955 overwrite_input : bool, optional
3956 If True, then allow the input array `a` to be modified by intermediate
3957 calculations, to save memory. In this case, the contents of the input
3958 `a` after this function completes is undefined.
3959 method : str, optional
3960 This parameter specifies the method to use for estimating the
3961 percentile. There are many different methods, some unique to NumPy.
3962 See the notes for explanation. The options sorted by their R type
3963 as summarized in the H&F paper [1]_ are:
3965 1. 'inverted_cdf'
3966 2. 'averaged_inverted_cdf'
3967 3. 'closest_observation'
3968 4. 'interpolated_inverted_cdf'
3969 5. 'hazen'
3970 6. 'weibull'
3971 7. 'linear' (default)
3972 8. 'median_unbiased'
3973 9. 'normal_unbiased'
3975 The first three methods are discontinuous. NumPy further defines the
3976 following discontinuous variations of the default 'linear' (7.) option:
3978 * 'lower'
3979 * 'higher',
3980 * 'midpoint'
3981 * 'nearest'
3983 .. versionchanged:: 1.22.0
3984 This argument was previously called "interpolation" and only
3985 offered the "linear" default and last four options.
3987 keepdims : bool, optional
3988 If this is set to True, the axes which are reduced are left in
3989 the result as dimensions with size one. With this option, the
3990 result will broadcast correctly against the original array `a`.
3992 .. versionadded:: 1.9.0
3994 interpolation : str, optional
3995 Deprecated name for the method keyword argument.
3997 .. deprecated:: 1.22.0
3999 Returns
4000 -------
4001 percentile : scalar or ndarray
4002 If `q` is a single percentile and `axis=None`, then the result
4003 is a scalar. If multiple percentiles are given, first axis of
4004 the result corresponds to the percentiles. The other axes are
4005 the axes that remain after the reduction of `a`. If the input
4006 contains integers or floats smaller than ``float64``, the output
4007 data-type is ``float64``. Otherwise, the output data-type is the
4008 same as that of the input. If `out` is specified, that array is
4009 returned instead.
4011 See Also
4012 --------
4013 mean
4014 median : equivalent to ``percentile(..., 50)``
4015 nanpercentile
4016 quantile : equivalent to percentile, except q in the range [0, 1].
4018 Notes
4019 -----
4020 Given a vector ``V`` of length ``n``, the q-th percentile of ``V`` is
4021 the value ``q/100`` of the way from the minimum to the maximum in a
4022 sorted copy of ``V``. The values and distances of the two nearest
4023 neighbors as well as the `method` parameter will determine the
4024 percentile if the normalized ranking does not match the location of
4025 ``q`` exactly. This function is the same as the median if ``q=50``, the
4026 same as the minimum if ``q=0`` and the same as the maximum if
4027 ``q=100``.
4029 The optional `method` parameter specifies the method to use when the
4030 desired percentile lies between two indexes ``i`` and ``j = i + 1``.
4031 In that case, we first determine ``i + g``, a virtual index that lies
4032 between ``i`` and ``j``, where ``i`` is the floor and ``g`` is the
4033 fractional part of the index. The final result is, then, an interpolation
4034 of ``a[i]`` and ``a[j]`` based on ``g``. During the computation of ``g``,
4035 ``i`` and ``j`` are modified using correction constants ``alpha`` and
4036 ``beta`` whose choices depend on the ``method`` used. Finally, note that
4037 since Python uses 0-based indexing, the code subtracts another 1 from the
4038 index internally.
4040 The following formula determines the virtual index ``i + g``, the location
4041 of the percentile in the sorted sample:
4043 .. math::
4044 i + g = (q / 100) * ( n - alpha - beta + 1 ) + alpha
4046 The different methods then work as follows
4048 inverted_cdf:
4049 method 1 of H&F [1]_.
4050 This method gives discontinuous results:
4052 * if g > 0 ; then take j
4053 * if g = 0 ; then take i
4055 averaged_inverted_cdf:
4056 method 2 of H&F [1]_.
4057 This method give discontinuous results:
4059 * if g > 0 ; then take j
4060 * if g = 0 ; then average between bounds
4062 closest_observation:
4063 method 3 of H&F [1]_.
4064 This method give discontinuous results:
4066 * if g > 0 ; then take j
4067 * if g = 0 and index is odd ; then take j
4068 * if g = 0 and index is even ; then take i
4070 interpolated_inverted_cdf:
4071 method 4 of H&F [1]_.
4072 This method give continuous results using:
4074 * alpha = 0
4075 * beta = 1
4077 hazen:
4078 method 5 of H&F [1]_.
4079 This method give continuous results using:
4081 * alpha = 1/2
4082 * beta = 1/2
4084 weibull:
4085 method 6 of H&F [1]_.
4086 This method give continuous results using:
4088 * alpha = 0
4089 * beta = 0
4091 linear:
4092 method 7 of H&F [1]_.
4093 This method give continuous results using:
4095 * alpha = 1
4096 * beta = 1
4098 median_unbiased:
4099 method 8 of H&F [1]_.
4100 This method is probably the best method if the sample
4101 distribution function is unknown (see reference).
4102 This method give continuous results using:
4104 * alpha = 1/3
4105 * beta = 1/3
4107 normal_unbiased:
4108 method 9 of H&F [1]_.
4109 This method is probably the best method if the sample
4110 distribution function is known to be normal.
4111 This method give continuous results using:
4113 * alpha = 3/8
4114 * beta = 3/8
4116 lower:
4117 NumPy method kept for backwards compatibility.
4118 Takes ``i`` as the interpolation point.
4120 higher:
4121 NumPy method kept for backwards compatibility.
4122 Takes ``j`` as the interpolation point.
4124 nearest:
4125 NumPy method kept for backwards compatibility.
4126 Takes ``i`` or ``j``, whichever is nearest.
4128 midpoint:
4129 NumPy method kept for backwards compatibility.
4130 Uses ``(i + j) / 2``.
4132 Examples
4133 --------
4134 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
4135 >>> a
4136 array([[10, 7, 4],
4137 [ 3, 2, 1]])
4138 >>> np.percentile(a, 50)
4139 3.5
4140 >>> np.percentile(a, 50, axis=0)
4141 array([6.5, 4.5, 2.5])
4142 >>> np.percentile(a, 50, axis=1)
4143 array([7., 2.])
4144 >>> np.percentile(a, 50, axis=1, keepdims=True)
4145 array([[7.],
4146 [2.]])
4148 >>> m = np.percentile(a, 50, axis=0)
4149 >>> out = np.zeros_like(m)
4150 >>> np.percentile(a, 50, axis=0, out=out)
4151 array([6.5, 4.5, 2.5])
4152 >>> m
4153 array([6.5, 4.5, 2.5])
4155 >>> b = a.copy()
4156 >>> np.percentile(b, 50, axis=1, overwrite_input=True)
4157 array([7., 2.])
4158 >>> assert not np.all(a == b)
4160 The different methods can be visualized graphically:
4162 .. plot::
4164 import matplotlib.pyplot as plt
4166 a = np.arange(4)
4167 p = np.linspace(0, 100, 6001)
4168 ax = plt.gca()
4169 lines = [
4170 ('linear', '-', 'C0'),
4171 ('inverted_cdf', ':', 'C1'),
4172 # Almost the same as `inverted_cdf`:
4173 ('averaged_inverted_cdf', '-.', 'C1'),
4174 ('closest_observation', ':', 'C2'),
4175 ('interpolated_inverted_cdf', '--', 'C1'),
4176 ('hazen', '--', 'C3'),
4177 ('weibull', '-.', 'C4'),
4178 ('median_unbiased', '--', 'C5'),
4179 ('normal_unbiased', '-.', 'C6'),
4180 ]
4181 for method, style, color in lines:
4182 ax.plot(
4183 p, np.percentile(a, p, method=method),
4184 label=method, linestyle=style, color=color)
4185 ax.set(
4186 title='Percentiles for different methods and data: ' + str(a),
4187 xlabel='Percentile',
4188 ylabel='Estimated percentile value',
4189 yticks=a)
4190 ax.legend(bbox_to_anchor=(1.03, 1))
4191 plt.tight_layout()
4192 plt.show()
4194 References
4195 ----------
4196 .. [1] R. J. Hyndman and Y. Fan,
4197 "Sample quantiles in statistical packages,"
4198 The American Statistician, 50(4), pp. 361-365, 1996
4200 """
4201 if interpolation is not None:
4202 method = _check_interpolation_as_method(
4203 method, interpolation, "percentile")
4205 a = np.asanyarray(a)
4206 if a.dtype.kind == "c":
4207 raise TypeError("a must be an array of real numbers")
4209 q = np.true_divide(q, 100)
4210 q = asanyarray(q) # undo any decay that the ufunc performed (see gh-13105)
4211 if not _quantile_is_valid(q):
4212 raise ValueError("Percentiles must be in the range [0, 100]")
4213 return _quantile_unchecked(
4214 a, q, axis, out, overwrite_input, method, keepdims)
4217def _quantile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
4218 method=None, keepdims=None, *, interpolation=None):
4219 return (a, q, out)
4222@array_function_dispatch(_quantile_dispatcher)
4223def quantile(a,
4224 q,
4225 axis=None,
4226 out=None,
4227 overwrite_input=False,
4228 method="linear",
4229 keepdims=False,
4230 *,
4231 interpolation=None):
4232 """
4233 Compute the q-th quantile of the data along the specified axis.
4235 .. versionadded:: 1.15.0
4237 Parameters
4238 ----------
4239 a : array_like of real numbers
4240 Input array or object that can be converted to an array.
4241 q : array_like of float
4242 Quantile or sequence of quantiles to compute, which must be between
4243 0 and 1 inclusive.
4244 axis : {int, tuple of int, None}, optional
4245 Axis or axes along which the quantiles are computed. The default is
4246 to compute the quantile(s) along a flattened version of the array.
4247 out : ndarray, optional
4248 Alternative output array in which to place the result. It must have
4249 the same shape and buffer length as the expected output, but the
4250 type (of the output) will be cast if necessary.
4251 overwrite_input : bool, optional
4252 If True, then allow the input array `a` to be modified by
4253 intermediate calculations, to save memory. In this case, the
4254 contents of the input `a` after this function completes is
4255 undefined.
4256 method : str, optional
4257 This parameter specifies the method to use for estimating the
4258 quantile. There are many different methods, some unique to NumPy.
4259 See the notes for explanation. The options sorted by their R type
4260 as summarized in the H&F paper [1]_ are:
4262 1. 'inverted_cdf'
4263 2. 'averaged_inverted_cdf'
4264 3. 'closest_observation'
4265 4. 'interpolated_inverted_cdf'
4266 5. 'hazen'
4267 6. 'weibull'
4268 7. 'linear' (default)
4269 8. 'median_unbiased'
4270 9. 'normal_unbiased'
4272 The first three methods are discontinuous. NumPy further defines the
4273 following discontinuous variations of the default 'linear' (7.) option:
4275 * 'lower'
4276 * 'higher',
4277 * 'midpoint'
4278 * 'nearest'
4280 .. versionchanged:: 1.22.0
4281 This argument was previously called "interpolation" and only
4282 offered the "linear" default and last four options.
4284 keepdims : bool, optional
4285 If this is set to True, the axes which are reduced are left in
4286 the result as dimensions with size one. With this option, the
4287 result will broadcast correctly against the original array `a`.
4289 interpolation : str, optional
4290 Deprecated name for the method keyword argument.
4292 .. deprecated:: 1.22.0
4294 Returns
4295 -------
4296 quantile : scalar or ndarray
4297 If `q` is a single quantile and `axis=None`, then the result
4298 is a scalar. If multiple quantiles are given, first axis of
4299 the result corresponds to the quantiles. The other axes are
4300 the axes that remain after the reduction of `a`. If the input
4301 contains integers or floats smaller than ``float64``, the output
4302 data-type is ``float64``. Otherwise, the output data-type is the
4303 same as that of the input. If `out` is specified, that array is
4304 returned instead.
4306 See Also
4307 --------
4308 mean
4309 percentile : equivalent to quantile, but with q in the range [0, 100].
4310 median : equivalent to ``quantile(..., 0.5)``
4311 nanquantile
4313 Notes
4314 -----
4315 Given a vector ``V`` of length ``n``, the q-th quantile of ``V`` is
4316 the value ``q`` of the way from the minimum to the maximum in a
4317 sorted copy of ``V``. The values and distances of the two nearest
4318 neighbors as well as the `method` parameter will determine the
4319 quantile if the normalized ranking does not match the location of
4320 ``q`` exactly. This function is the same as the median if ``q=0.5``, the
4321 same as the minimum if ``q=0.0`` and the same as the maximum if
4322 ``q=1.0``.
4324 The optional `method` parameter specifies the method to use when the
4325 desired quantile lies between two indexes ``i`` and ``j = i + 1``.
4326 In that case, we first determine ``i + g``, a virtual index that lies
4327 between ``i`` and ``j``, where ``i`` is the floor and ``g`` is the
4328 fractional part of the index. The final result is, then, an interpolation
4329 of ``a[i]`` and ``a[j]`` based on ``g``. During the computation of ``g``,
4330 ``i`` and ``j`` are modified using correction constants ``alpha`` and
4331 ``beta`` whose choices depend on the ``method`` used. Finally, note that
4332 since Python uses 0-based indexing, the code subtracts another 1 from the
4333 index internally.
4335 The following formula determines the virtual index ``i + g``, the location
4336 of the quantile in the sorted sample:
4338 .. math::
4339 i + g = q * ( n - alpha - beta + 1 ) + alpha
4341 The different methods then work as follows
4343 inverted_cdf:
4344 method 1 of H&F [1]_.
4345 This method gives discontinuous results:
4347 * if g > 0 ; then take j
4348 * if g = 0 ; then take i
4350 averaged_inverted_cdf:
4351 method 2 of H&F [1]_.
4352 This method gives discontinuous results:
4354 * if g > 0 ; then take j
4355 * if g = 0 ; then average between bounds
4357 closest_observation:
4358 method 3 of H&F [1]_.
4359 This method gives discontinuous results:
4361 * if g > 0 ; then take j
4362 * if g = 0 and index is odd ; then take j
4363 * if g = 0 and index is even ; then take i
4365 interpolated_inverted_cdf:
4366 method 4 of H&F [1]_.
4367 This method gives continuous results using:
4369 * alpha = 0
4370 * beta = 1
4372 hazen:
4373 method 5 of H&F [1]_.
4374 This method gives continuous results using:
4376 * alpha = 1/2
4377 * beta = 1/2
4379 weibull:
4380 method 6 of H&F [1]_.
4381 This method gives continuous results using:
4383 * alpha = 0
4384 * beta = 0
4386 linear:
4387 method 7 of H&F [1]_.
4388 This method gives continuous results using:
4390 * alpha = 1
4391 * beta = 1
4393 median_unbiased:
4394 method 8 of H&F [1]_.
4395 This method is probably the best method if the sample
4396 distribution function is unknown (see reference).
4397 This method gives continuous results using:
4399 * alpha = 1/3
4400 * beta = 1/3
4402 normal_unbiased:
4403 method 9 of H&F [1]_.
4404 This method is probably the best method if the sample
4405 distribution function is known to be normal.
4406 This method gives continuous results using:
4408 * alpha = 3/8
4409 * beta = 3/8
4411 lower:
4412 NumPy method kept for backwards compatibility.
4413 Takes ``i`` as the interpolation point.
4415 higher:
4416 NumPy method kept for backwards compatibility.
4417 Takes ``j`` as the interpolation point.
4419 nearest:
4420 NumPy method kept for backwards compatibility.
4421 Takes ``i`` or ``j``, whichever is nearest.
4423 midpoint:
4424 NumPy method kept for backwards compatibility.
4425 Uses ``(i + j) / 2``.
4427 Examples
4428 --------
4429 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
4430 >>> a
4431 array([[10, 7, 4],
4432 [ 3, 2, 1]])
4433 >>> np.quantile(a, 0.5)
4434 3.5
4435 >>> np.quantile(a, 0.5, axis=0)
4436 array([6.5, 4.5, 2.5])
4437 >>> np.quantile(a, 0.5, axis=1)
4438 array([7., 2.])
4439 >>> np.quantile(a, 0.5, axis=1, keepdims=True)
4440 array([[7.],
4441 [2.]])
4442 >>> m = np.quantile(a, 0.5, axis=0)
4443 >>> out = np.zeros_like(m)
4444 >>> np.quantile(a, 0.5, axis=0, out=out)
4445 array([6.5, 4.5, 2.5])
4446 >>> m
4447 array([6.5, 4.5, 2.5])
4448 >>> b = a.copy()
4449 >>> np.quantile(b, 0.5, axis=1, overwrite_input=True)
4450 array([7., 2.])
4451 >>> assert not np.all(a == b)
4453 See also `numpy.percentile` for a visualization of most methods.
4455 References
4456 ----------
4457 .. [1] R. J. Hyndman and Y. Fan,
4458 "Sample quantiles in statistical packages,"
4459 The American Statistician, 50(4), pp. 361-365, 1996
4461 """
4462 if interpolation is not None:
4463 method = _check_interpolation_as_method(
4464 method, interpolation, "quantile")
4466 a = np.asanyarray(a)
4467 if a.dtype.kind == "c":
4468 raise TypeError("a must be an array of real numbers")
4470 q = np.asanyarray(q)
4471 if not _quantile_is_valid(q):
4472 raise ValueError("Quantiles must be in the range [0, 1]")
4473 return _quantile_unchecked(
4474 a, q, axis, out, overwrite_input, method, keepdims)
4477def _quantile_unchecked(a,
4478 q,
4479 axis=None,
4480 out=None,
4481 overwrite_input=False,
4482 method="linear",
4483 keepdims=False):
4484 """Assumes that q is in [0, 1], and is an ndarray"""
4485 return _ureduce(a,
4486 func=_quantile_ureduce_func,
4487 q=q,
4488 keepdims=keepdims,
4489 axis=axis,
4490 out=out,
4491 overwrite_input=overwrite_input,
4492 method=method)
4495def _quantile_is_valid(q):
4496 # avoid expensive reductions, relevant for arrays with < O(1000) elements
4497 if q.ndim == 1 and q.size < 10:
4498 for i in range(q.size):
4499 if not (0.0 <= q[i] <= 1.0):
4500 return False
4501 else:
4502 if not (np.all(0 <= q) and np.all(q <= 1)):
4503 return False
4504 return True
4507def _check_interpolation_as_method(method, interpolation, fname):
4508 # Deprecated NumPy 1.22, 2021-11-08
4509 warnings.warn(
4510 f"the `interpolation=` argument to {fname} was renamed to "
4511 "`method=`, which has additional options.\n"
4512 "Users of the modes 'nearest', 'lower', 'higher', or "
4513 "'midpoint' are encouraged to review the method they used. "
4514 "(Deprecated NumPy 1.22)",
4515 DeprecationWarning, stacklevel=4)
4516 if method != "linear":
4517 # sanity check, we assume this basically never happens
4518 raise TypeError(
4519 "You shall not pass both `method` and `interpolation`!\n"
4520 "(`interpolation` is Deprecated in favor of `method`)")
4521 return interpolation
4524def _compute_virtual_index(n, quantiles, alpha: float, beta: float):
4525 """
4526 Compute the floating point indexes of an array for the linear
4527 interpolation of quantiles.
4528 n : array_like
4529 The sample sizes.
4530 quantiles : array_like
4531 The quantiles values.
4532 alpha : float
4533 A constant used to correct the index computed.
4534 beta : float
4535 A constant used to correct the index computed.
4537 alpha and beta values depend on the chosen method
4538 (see quantile documentation)
4540 Reference:
4541 Hyndman&Fan paper "Sample Quantiles in Statistical Packages",
4542 DOI: 10.1080/00031305.1996.10473566
4543 """
4544 return n * quantiles + (
4545 alpha + quantiles * (1 - alpha - beta)
4546 ) - 1
4549def _get_gamma(virtual_indexes, previous_indexes, method):
4550 """
4551 Compute gamma (a.k.a 'm' or 'weight') for the linear interpolation
4552 of quantiles.
4554 virtual_indexes : array_like
4555 The indexes where the percentile is supposed to be found in the sorted
4556 sample.
4557 previous_indexes : array_like
4558 The floor values of virtual_indexes.
4559 interpolation : dict
4560 The interpolation method chosen, which may have a specific rule
4561 modifying gamma.
4563 gamma is usually the fractional part of virtual_indexes but can be modified
4564 by the interpolation method.
4565 """
4566 gamma = np.asanyarray(virtual_indexes - previous_indexes)
4567 gamma = method["fix_gamma"](gamma, virtual_indexes)
4568 return np.asanyarray(gamma)
4571def _lerp(a, b, t, out=None):
4572 """
4573 Compute the linear interpolation weighted by gamma on each point of
4574 two same shape array.
4576 a : array_like
4577 Left bound.
4578 b : array_like
4579 Right bound.
4580 t : array_like
4581 The interpolation weight.
4582 out : array_like
4583 Output array.
4584 """
4585 diff_b_a = subtract(b, a)
4586 # asanyarray is a stop-gap until gh-13105
4587 lerp_interpolation = asanyarray(add(a, diff_b_a * t, out=out))
4588 subtract(b, diff_b_a * (1 - t), out=lerp_interpolation, where=t >= 0.5)
4589 if lerp_interpolation.ndim == 0 and out is None:
4590 lerp_interpolation = lerp_interpolation[()] # unpack 0d arrays
4591 return lerp_interpolation
4594def _get_gamma_mask(shape, default_value, conditioned_value, where):
4595 out = np.full(shape, default_value)
4596 np.copyto(out, conditioned_value, where=where, casting="unsafe")
4597 return out
4600def _discret_interpolation_to_boundaries(index, gamma_condition_fun):
4601 previous = np.floor(index)
4602 next = previous + 1
4603 gamma = index - previous
4604 res = _get_gamma_mask(shape=index.shape,
4605 default_value=next,
4606 conditioned_value=previous,
4607 where=gamma_condition_fun(gamma, index)
4608 ).astype(np.intp)
4609 # Some methods can lead to out-of-bound integers, clip them:
4610 res[res < 0] = 0
4611 return res
4614def _closest_observation(n, quantiles):
4615 gamma_fun = lambda gamma, index: (gamma == 0) & (np.floor(index) % 2 == 0)
4616 return _discret_interpolation_to_boundaries((n * quantiles) - 1 - 0.5,
4617 gamma_fun)
4620def _inverted_cdf(n, quantiles):
4621 gamma_fun = lambda gamma, _: (gamma == 0)
4622 return _discret_interpolation_to_boundaries((n * quantiles) - 1,
4623 gamma_fun)
4626def _quantile_ureduce_func(
4627 a: np.array,
4628 q: np.array,
4629 axis: int = None,
4630 out=None,
4631 overwrite_input: bool = False,
4632 method="linear",
4633) -> np.array:
4634 if q.ndim > 2:
4635 # The code below works fine for nd, but it might not have useful
4636 # semantics. For now, keep the supported dimensions the same as it was
4637 # before.
4638 raise ValueError("q must be a scalar or 1d")
4639 if overwrite_input:
4640 if axis is None:
4641 axis = 0
4642 arr = a.ravel()
4643 else:
4644 arr = a
4645 else:
4646 if axis is None:
4647 axis = 0
4648 arr = a.flatten()
4649 else:
4650 arr = a.copy()
4651 result = _quantile(arr,
4652 quantiles=q,
4653 axis=axis,
4654 method=method,
4655 out=out)
4656 return result
4659def _get_indexes(arr, virtual_indexes, valid_values_count):
4660 """
4661 Get the valid indexes of arr neighbouring virtual_indexes.
4662 Note
4663 This is a companion function to linear interpolation of
4664 Quantiles
4666 Returns
4667 -------
4668 (previous_indexes, next_indexes): Tuple
4669 A Tuple of virtual_indexes neighbouring indexes
4670 """
4671 previous_indexes = np.asanyarray(np.floor(virtual_indexes))
4672 next_indexes = np.asanyarray(previous_indexes + 1)
4673 indexes_above_bounds = virtual_indexes >= valid_values_count - 1
4674 # When indexes is above max index, take the max value of the array
4675 if indexes_above_bounds.any():
4676 previous_indexes[indexes_above_bounds] = -1
4677 next_indexes[indexes_above_bounds] = -1
4678 # When indexes is below min index, take the min value of the array
4679 indexes_below_bounds = virtual_indexes < 0
4680 if indexes_below_bounds.any():
4681 previous_indexes[indexes_below_bounds] = 0
4682 next_indexes[indexes_below_bounds] = 0
4683 if np.issubdtype(arr.dtype, np.inexact):
4684 # After the sort, slices having NaNs will have for last element a NaN
4685 virtual_indexes_nans = np.isnan(virtual_indexes)
4686 if virtual_indexes_nans.any():
4687 previous_indexes[virtual_indexes_nans] = -1
4688 next_indexes[virtual_indexes_nans] = -1
4689 previous_indexes = previous_indexes.astype(np.intp)
4690 next_indexes = next_indexes.astype(np.intp)
4691 return previous_indexes, next_indexes
4694def _quantile(
4695 arr: np.array,
4696 quantiles: np.array,
4697 axis: int = -1,
4698 method="linear",
4699 out=None,
4700):
4701 """
4702 Private function that doesn't support extended axis or keepdims.
4703 These methods are extended to this function using _ureduce
4704 See nanpercentile for parameter usage
4705 It computes the quantiles of the array for the given axis.
4706 A linear interpolation is performed based on the `interpolation`.
4708 By default, the method is "linear" where alpha == beta == 1 which
4709 performs the 7th method of Hyndman&Fan.
4710 With "median_unbiased" we get alpha == beta == 1/3
4711 thus the 8th method of Hyndman&Fan.
4712 """
4713 # --- Setup
4714 arr = np.asanyarray(arr)
4715 values_count = arr.shape[axis]
4716 # The dimensions of `q` are prepended to the output shape, so we need the
4717 # axis being sampled from `arr` to be last.
4718 DATA_AXIS = 0
4719 if axis != DATA_AXIS: # But moveaxis is slow, so only call it if axis!=0.
4720 arr = np.moveaxis(arr, axis, destination=DATA_AXIS)
4721 # --- Computation of indexes
4722 # Index where to find the value in the sorted array.
4723 # Virtual because it is a floating point value, not an valid index.
4724 # The nearest neighbours are used for interpolation
4725 try:
4726 method = _QuantileMethods[method]
4727 except KeyError:
4728 raise ValueError(
4729 f"{method!r} is not a valid method. Use one of: "
4730 f"{_QuantileMethods.keys()}") from None
4731 virtual_indexes = method["get_virtual_index"](values_count, quantiles)
4732 virtual_indexes = np.asanyarray(virtual_indexes)
4733 if np.issubdtype(virtual_indexes.dtype, np.integer):
4734 # No interpolation needed, take the points along axis
4735 if np.issubdtype(arr.dtype, np.inexact):
4736 # may contain nan, which would sort to the end
4737 arr.partition(concatenate((virtual_indexes.ravel(), [-1])), axis=0)
4738 slices_having_nans = np.isnan(arr[-1])
4739 else:
4740 # cannot contain nan
4741 arr.partition(virtual_indexes.ravel(), axis=0)
4742 slices_having_nans = np.array(False, dtype=bool)
4743 result = take(arr, virtual_indexes, axis=0, out=out)
4744 else:
4745 previous_indexes, next_indexes = _get_indexes(arr,
4746 virtual_indexes,
4747 values_count)
4748 # --- Sorting
4749 arr.partition(
4750 np.unique(np.concatenate(([0, -1],
4751 previous_indexes.ravel(),
4752 next_indexes.ravel(),
4753 ))),
4754 axis=DATA_AXIS)
4755 if np.issubdtype(arr.dtype, np.inexact):
4756 slices_having_nans = np.isnan(
4757 take(arr, indices=-1, axis=DATA_AXIS)
4758 )
4759 else:
4760 slices_having_nans = None
4761 # --- Get values from indexes
4762 previous = np.take(arr, previous_indexes, axis=DATA_AXIS)
4763 next = np.take(arr, next_indexes, axis=DATA_AXIS)
4764 # --- Linear interpolation
4765 gamma = _get_gamma(virtual_indexes, previous_indexes, method)
4766 result_shape = virtual_indexes.shape + (1,) * (arr.ndim - 1)
4767 gamma = gamma.reshape(result_shape)
4768 result = _lerp(previous,
4769 next,
4770 gamma,
4771 out=out)
4772 if np.any(slices_having_nans):
4773 if result.ndim == 0 and out is None:
4774 # can't write to a scalar
4775 result = arr.dtype.type(np.nan)
4776 else:
4777 result[..., slices_having_nans] = np.nan
4778 return result
4781def _trapz_dispatcher(y, x=None, dx=None, axis=None):
4782 return (y, x)
4785@array_function_dispatch(_trapz_dispatcher)
4786def trapz(y, x=None, dx=1.0, axis=-1):
4787 r"""
4788 Integrate along the given axis using the composite trapezoidal rule.
4790 If `x` is provided, the integration happens in sequence along its
4791 elements - they are not sorted.
4793 Integrate `y` (`x`) along each 1d slice on the given axis, compute
4794 :math:`\int y(x) dx`.
4795 When `x` is specified, this integrates along the parametric curve,
4796 computing :math:`\int_t y(t) dt =
4797 \int_t y(t) \left.\frac{dx}{dt}\right|_{x=x(t)} dt`.
4799 Parameters
4800 ----------
4801 y : array_like
4802 Input array to integrate.
4803 x : array_like, optional
4804 The sample points corresponding to the `y` values. If `x` is None,
4805 the sample points are assumed to be evenly spaced `dx` apart. The
4806 default is None.
4807 dx : scalar, optional
4808 The spacing between sample points when `x` is None. The default is 1.
4809 axis : int, optional
4810 The axis along which to integrate.
4812 Returns
4813 -------
4814 trapz : float or ndarray
4815 Definite integral of `y` = n-dimensional array as approximated along
4816 a single axis by the trapezoidal rule. If `y` is a 1-dimensional array,
4817 then the result is a float. If `n` is greater than 1, then the result
4818 is an `n`-1 dimensional array.
4820 See Also
4821 --------
4822 sum, cumsum
4824 Notes
4825 -----
4826 Image [2]_ illustrates trapezoidal rule -- y-axis locations of points
4827 will be taken from `y` array, by default x-axis distances between
4828 points will be 1.0, alternatively they can be provided with `x` array
4829 or with `dx` scalar. Return value will be equal to combined area under
4830 the red lines.
4833 References
4834 ----------
4835 .. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule
4837 .. [2] Illustration image:
4838 https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
4840 Examples
4841 --------
4842 Use the trapezoidal rule on evenly spaced points:
4844 >>> np.trapz([1, 2, 3])
4845 4.0
4847 The spacing between sample points can be selected by either the
4848 ``x`` or ``dx`` arguments:
4850 >>> np.trapz([1, 2, 3], x=[4, 6, 8])
4851 8.0
4852 >>> np.trapz([1, 2, 3], dx=2)
4853 8.0
4855 Using a decreasing ``x`` corresponds to integrating in reverse:
4857 >>> np.trapz([1, 2, 3], x=[8, 6, 4])
4858 -8.0
4860 More generally ``x`` is used to integrate along a parametric curve. We can
4861 estimate the integral :math:`\int_0^1 x^2 = 1/3` using:
4863 >>> x = np.linspace(0, 1, num=50)
4864 >>> y = x**2
4865 >>> np.trapz(y, x)
4866 0.33340274885464394
4868 Or estimate the area of a circle, noting we repeat the sample which closes
4869 the curve:
4871 >>> theta = np.linspace(0, 2 * np.pi, num=1000, endpoint=True)
4872 >>> np.trapz(np.cos(theta), x=np.sin(theta))
4873 3.141571941375841
4875 ``np.trapz`` can be applied along a specified axis to do multiple
4876 computations in one call:
4878 >>> a = np.arange(6).reshape(2, 3)
4879 >>> a
4880 array([[0, 1, 2],
4881 [3, 4, 5]])
4882 >>> np.trapz(a, axis=0)
4883 array([1.5, 2.5, 3.5])
4884 >>> np.trapz(a, axis=1)
4885 array([2., 8.])
4886 """
4887 y = asanyarray(y)
4888 if x is None:
4889 d = dx
4890 else:
4891 x = asanyarray(x)
4892 if x.ndim == 1:
4893 d = diff(x)
4894 # reshape to correct shape
4895 shape = [1]*y.ndim
4896 shape[axis] = d.shape[0]
4897 d = d.reshape(shape)
4898 else:
4899 d = diff(x, axis=axis)
4900 nd = y.ndim
4901 slice1 = [slice(None)]*nd
4902 slice2 = [slice(None)]*nd
4903 slice1[axis] = slice(1, None)
4904 slice2[axis] = slice(None, -1)
4905 try:
4906 ret = (d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0).sum(axis)
4907 except ValueError:
4908 # Operations didn't work, cast to ndarray
4909 d = np.asarray(d)
4910 y = np.asarray(y)
4911 ret = add.reduce(d * (y[tuple(slice1)]+y[tuple(slice2)])/2.0, axis)
4912 return ret
4915def _meshgrid_dispatcher(*xi, copy=None, sparse=None, indexing=None):
4916 return xi
4919# Based on scitools meshgrid
4920@array_function_dispatch(_meshgrid_dispatcher)
4921def meshgrid(*xi, copy=True, sparse=False, indexing='xy'):
4922 """
4923 Return coordinate matrices from coordinate vectors.
4925 Make N-D coordinate arrays for vectorized evaluations of
4926 N-D scalar/vector fields over N-D grids, given
4927 one-dimensional coordinate arrays x1, x2,..., xn.
4929 .. versionchanged:: 1.9
4930 1-D and 0-D cases are allowed.
4932 Parameters
4933 ----------
4934 x1, x2,..., xn : array_like
4935 1-D arrays representing the coordinates of a grid.
4936 indexing : {'xy', 'ij'}, optional
4937 Cartesian ('xy', default) or matrix ('ij') indexing of output.
4938 See Notes for more details.
4940 .. versionadded:: 1.7.0
4941 sparse : bool, optional
4942 If True the shape of the returned coordinate array for dimension *i*
4943 is reduced from ``(N1, ..., Ni, ... Nn)`` to
4944 ``(1, ..., 1, Ni, 1, ..., 1)``. These sparse coordinate grids are
4945 intended to be use with :ref:`basics.broadcasting`. When all
4946 coordinates are used in an expression, broadcasting still leads to a
4947 fully-dimensonal result array.
4949 Default is False.
4951 .. versionadded:: 1.7.0
4952 copy : bool, optional
4953 If False, a view into the original arrays are returned in order to
4954 conserve memory. Default is True. Please note that
4955 ``sparse=False, copy=False`` will likely return non-contiguous
4956 arrays. Furthermore, more than one element of a broadcast array
4957 may refer to a single memory location. If you need to write to the
4958 arrays, make copies first.
4960 .. versionadded:: 1.7.0
4962 Returns
4963 -------
4964 X1, X2,..., XN : ndarray
4965 For vectors `x1`, `x2`,..., `xn` with lengths ``Ni=len(xi)``,
4966 returns ``(N1, N2, N3,..., Nn)`` shaped arrays if indexing='ij'
4967 or ``(N2, N1, N3,..., Nn)`` shaped arrays if indexing='xy'
4968 with the elements of `xi` repeated to fill the matrix along
4969 the first dimension for `x1`, the second for `x2` and so on.
4971 Notes
4972 -----
4973 This function supports both indexing conventions through the indexing
4974 keyword argument. Giving the string 'ij' returns a meshgrid with
4975 matrix indexing, while 'xy' returns a meshgrid with Cartesian indexing.
4976 In the 2-D case with inputs of length M and N, the outputs are of shape
4977 (N, M) for 'xy' indexing and (M, N) for 'ij' indexing. In the 3-D case
4978 with inputs of length M, N and P, outputs are of shape (N, M, P) for
4979 'xy' indexing and (M, N, P) for 'ij' indexing. The difference is
4980 illustrated by the following code snippet::
4982 xv, yv = np.meshgrid(x, y, indexing='ij')
4983 for i in range(nx):
4984 for j in range(ny):
4985 # treat xv[i,j], yv[i,j]
4987 xv, yv = np.meshgrid(x, y, indexing='xy')
4988 for i in range(nx):
4989 for j in range(ny):
4990 # treat xv[j,i], yv[j,i]
4992 In the 1-D and 0-D case, the indexing and sparse keywords have no effect.
4994 See Also
4995 --------
4996 mgrid : Construct a multi-dimensional "meshgrid" using indexing notation.
4997 ogrid : Construct an open multi-dimensional "meshgrid" using indexing
4998 notation.
4999 how-to-index
5001 Examples
5002 --------
5003 >>> nx, ny = (3, 2)
5004 >>> x = np.linspace(0, 1, nx)
5005 >>> y = np.linspace(0, 1, ny)
5006 >>> xv, yv = np.meshgrid(x, y)
5007 >>> xv
5008 array([[0. , 0.5, 1. ],
5009 [0. , 0.5, 1. ]])
5010 >>> yv
5011 array([[0., 0., 0.],
5012 [1., 1., 1.]])
5014 The result of `meshgrid` is a coordinate grid:
5016 >>> import matplotlib.pyplot as plt
5017 >>> plt.plot(xv, yv, marker='o', color='k', linestyle='none')
5018 >>> plt.show()
5020 You can create sparse output arrays to save memory and computation time.
5022 >>> xv, yv = np.meshgrid(x, y, sparse=True)
5023 >>> xv
5024 array([[0. , 0.5, 1. ]])
5025 >>> yv
5026 array([[0.],
5027 [1.]])
5029 `meshgrid` is very useful to evaluate functions on a grid. If the
5030 function depends on all coordinates, both dense and sparse outputs can be
5031 used.
5033 >>> x = np.linspace(-5, 5, 101)
5034 >>> y = np.linspace(-5, 5, 101)
5035 >>> # full coordinate arrays
5036 >>> xx, yy = np.meshgrid(x, y)
5037 >>> zz = np.sqrt(xx**2 + yy**2)
5038 >>> xx.shape, yy.shape, zz.shape
5039 ((101, 101), (101, 101), (101, 101))
5040 >>> # sparse coordinate arrays
5041 >>> xs, ys = np.meshgrid(x, y, sparse=True)
5042 >>> zs = np.sqrt(xs**2 + ys**2)
5043 >>> xs.shape, ys.shape, zs.shape
5044 ((1, 101), (101, 1), (101, 101))
5045 >>> np.array_equal(zz, zs)
5046 True
5048 >>> h = plt.contourf(x, y, zs)
5049 >>> plt.axis('scaled')
5050 >>> plt.colorbar()
5051 >>> plt.show()
5052 """
5053 ndim = len(xi)
5055 if indexing not in ['xy', 'ij']:
5056 raise ValueError(
5057 "Valid values for `indexing` are 'xy' and 'ij'.")
5059 s0 = (1,) * ndim
5060 output = [np.asanyarray(x).reshape(s0[:i] + (-1,) + s0[i + 1:])
5061 for i, x in enumerate(xi)]
5063 if indexing == 'xy' and ndim > 1:
5064 # switch first and second axis
5065 output[0].shape = (1, -1) + s0[2:]
5066 output[1].shape = (-1, 1) + s0[2:]
5068 if not sparse:
5069 # Return the full N-D matrix (not only the 1-D vector)
5070 output = np.broadcast_arrays(*output, subok=True)
5072 if copy:
5073 output = [x.copy() for x in output]
5075 return output
5078def _delete_dispatcher(arr, obj, axis=None):
5079 return (arr, obj)
5082@array_function_dispatch(_delete_dispatcher)
5083def delete(arr, obj, axis=None):
5084 """
5085 Return a new array with sub-arrays along an axis deleted. For a one
5086 dimensional array, this returns those entries not returned by
5087 `arr[obj]`.
5089 Parameters
5090 ----------
5091 arr : array_like
5092 Input array.
5093 obj : slice, int or array of ints
5094 Indicate indices of sub-arrays to remove along the specified axis.
5096 .. versionchanged:: 1.19.0
5097 Boolean indices are now treated as a mask of elements to remove,
5098 rather than being cast to the integers 0 and 1.
5100 axis : int, optional
5101 The axis along which to delete the subarray defined by `obj`.
5102 If `axis` is None, `obj` is applied to the flattened array.
5104 Returns
5105 -------
5106 out : ndarray
5107 A copy of `arr` with the elements specified by `obj` removed. Note
5108 that `delete` does not occur in-place. If `axis` is None, `out` is
5109 a flattened array.
5111 See Also
5112 --------
5113 insert : Insert elements into an array.
5114 append : Append elements at the end of an array.
5116 Notes
5117 -----
5118 Often it is preferable to use a boolean mask. For example:
5120 >>> arr = np.arange(12) + 1
5121 >>> mask = np.ones(len(arr), dtype=bool)
5122 >>> mask[[0,2,4]] = False
5123 >>> result = arr[mask,...]
5125 Is equivalent to ``np.delete(arr, [0,2,4], axis=0)``, but allows further
5126 use of `mask`.
5128 Examples
5129 --------
5130 >>> arr = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]])
5131 >>> arr
5132 array([[ 1, 2, 3, 4],
5133 [ 5, 6, 7, 8],
5134 [ 9, 10, 11, 12]])
5135 >>> np.delete(arr, 1, 0)
5136 array([[ 1, 2, 3, 4],
5137 [ 9, 10, 11, 12]])
5139 >>> np.delete(arr, np.s_[::2], 1)
5140 array([[ 2, 4],
5141 [ 6, 8],
5142 [10, 12]])
5143 >>> np.delete(arr, [1,3,5], None)
5144 array([ 1, 3, 5, 7, 8, 9, 10, 11, 12])
5146 """
5147 wrap = None
5148 if type(arr) is not ndarray:
5149 try:
5150 wrap = arr.__array_wrap__
5151 except AttributeError:
5152 pass
5154 arr = asarray(arr)
5155 ndim = arr.ndim
5156 arrorder = 'F' if arr.flags.fnc else 'C'
5157 if axis is None:
5158 if ndim != 1:
5159 arr = arr.ravel()
5160 # needed for np.matrix, which is still not 1d after being ravelled
5161 ndim = arr.ndim
5162 axis = ndim - 1
5163 else:
5164 axis = normalize_axis_index(axis, ndim)
5166 slobj = [slice(None)]*ndim
5167 N = arr.shape[axis]
5168 newshape = list(arr.shape)
5170 if isinstance(obj, slice):
5171 start, stop, step = obj.indices(N)
5172 xr = range(start, stop, step)
5173 numtodel = len(xr)
5175 if numtodel <= 0:
5176 if wrap:
5177 return wrap(arr.copy(order=arrorder))
5178 else:
5179 return arr.copy(order=arrorder)
5181 # Invert if step is negative:
5182 if step < 0:
5183 step = -step
5184 start = xr[-1]
5185 stop = xr[0] + 1
5187 newshape[axis] -= numtodel
5188 new = empty(newshape, arr.dtype, arrorder)
5189 # copy initial chunk
5190 if start == 0:
5191 pass
5192 else:
5193 slobj[axis] = slice(None, start)
5194 new[tuple(slobj)] = arr[tuple(slobj)]
5195 # copy end chunk
5196 if stop == N:
5197 pass
5198 else:
5199 slobj[axis] = slice(stop-numtodel, None)
5200 slobj2 = [slice(None)]*ndim
5201 slobj2[axis] = slice(stop, None)
5202 new[tuple(slobj)] = arr[tuple(slobj2)]
5203 # copy middle pieces
5204 if step == 1:
5205 pass
5206 else: # use array indexing.
5207 keep = ones(stop-start, dtype=bool)
5208 keep[:stop-start:step] = False
5209 slobj[axis] = slice(start, stop-numtodel)
5210 slobj2 = [slice(None)]*ndim
5211 slobj2[axis] = slice(start, stop)
5212 arr = arr[tuple(slobj2)]
5213 slobj2[axis] = keep
5214 new[tuple(slobj)] = arr[tuple(slobj2)]
5215 if wrap:
5216 return wrap(new)
5217 else:
5218 return new
5220 if isinstance(obj, (int, integer)) and not isinstance(obj, bool):
5221 single_value = True
5222 else:
5223 single_value = False
5224 _obj = obj
5225 obj = np.asarray(obj)
5226 # `size == 0` to allow empty lists similar to indexing, but (as there)
5227 # is really too generic:
5228 if obj.size == 0 and not isinstance(_obj, np.ndarray):
5229 obj = obj.astype(intp)
5230 elif obj.size == 1 and obj.dtype.kind in "ui":
5231 # For a size 1 integer array we can use the single-value path
5232 # (most dtypes, except boolean, should just fail later).
5233 obj = obj.item()
5234 single_value = True
5236 if single_value:
5237 # optimization for a single value
5238 if (obj < -N or obj >= N):
5239 raise IndexError(
5240 "index %i is out of bounds for axis %i with "
5241 "size %i" % (obj, axis, N))
5242 if (obj < 0):
5243 obj += N
5244 newshape[axis] -= 1
5245 new = empty(newshape, arr.dtype, arrorder)
5246 slobj[axis] = slice(None, obj)
5247 new[tuple(slobj)] = arr[tuple(slobj)]
5248 slobj[axis] = slice(obj, None)
5249 slobj2 = [slice(None)]*ndim
5250 slobj2[axis] = slice(obj+1, None)
5251 new[tuple(slobj)] = arr[tuple(slobj2)]
5252 else:
5253 if obj.dtype == bool:
5254 if obj.shape != (N,):
5255 raise ValueError('boolean array argument obj to delete '
5256 'must be one dimensional and match the axis '
5257 'length of {}'.format(N))
5259 # optimization, the other branch is slower
5260 keep = ~obj
5261 else:
5262 keep = ones(N, dtype=bool)
5263 keep[obj,] = False
5265 slobj[axis] = keep
5266 new = arr[tuple(slobj)]
5268 if wrap:
5269 return wrap(new)
5270 else:
5271 return new
5274def _insert_dispatcher(arr, obj, values, axis=None):
5275 return (arr, obj, values)
5278@array_function_dispatch(_insert_dispatcher)
5279def insert(arr, obj, values, axis=None):
5280 """
5281 Insert values along the given axis before the given indices.
5283 Parameters
5284 ----------
5285 arr : array_like
5286 Input array.
5287 obj : int, slice or sequence of ints
5288 Object that defines the index or indices before which `values` is
5289 inserted.
5291 .. versionadded:: 1.8.0
5293 Support for multiple insertions when `obj` is a single scalar or a
5294 sequence with one element (similar to calling insert multiple
5295 times).
5296 values : array_like
5297 Values to insert into `arr`. If the type of `values` is different
5298 from that of `arr`, `values` is converted to the type of `arr`.
5299 `values` should be shaped so that ``arr[...,obj,...] = values``
5300 is legal.
5301 axis : int, optional
5302 Axis along which to insert `values`. If `axis` is None then `arr`
5303 is flattened first.
5305 Returns
5306 -------
5307 out : ndarray
5308 A copy of `arr` with `values` inserted. Note that `insert`
5309 does not occur in-place: a new array is returned. If
5310 `axis` is None, `out` is a flattened array.
5312 See Also
5313 --------
5314 append : Append elements at the end of an array.
5315 concatenate : Join a sequence of arrays along an existing axis.
5316 delete : Delete elements from an array.
5318 Notes
5319 -----
5320 Note that for higher dimensional inserts ``obj=0`` behaves very different
5321 from ``obj=[0]`` just like ``arr[:,0,:] = values`` is different from
5322 ``arr[:,[0],:] = values``.
5324 Examples
5325 --------
5326 >>> a = np.array([[1, 1], [2, 2], [3, 3]])
5327 >>> a
5328 array([[1, 1],
5329 [2, 2],
5330 [3, 3]])
5331 >>> np.insert(a, 1, 5)
5332 array([1, 5, 1, ..., 2, 3, 3])
5333 >>> np.insert(a, 1, 5, axis=1)
5334 array([[1, 5, 1],
5335 [2, 5, 2],
5336 [3, 5, 3]])
5338 Difference between sequence and scalars:
5340 >>> np.insert(a, [1], [[1],[2],[3]], axis=1)
5341 array([[1, 1, 1],
5342 [2, 2, 2],
5343 [3, 3, 3]])
5344 >>> np.array_equal(np.insert(a, 1, [1, 2, 3], axis=1),
5345 ... np.insert(a, [1], [[1],[2],[3]], axis=1))
5346 True
5348 >>> b = a.flatten()
5349 >>> b
5350 array([1, 1, 2, 2, 3, 3])
5351 >>> np.insert(b, [2, 2], [5, 6])
5352 array([1, 1, 5, ..., 2, 3, 3])
5354 >>> np.insert(b, slice(2, 4), [5, 6])
5355 array([1, 1, 5, ..., 2, 3, 3])
5357 >>> np.insert(b, [2, 2], [7.13, False]) # type casting
5358 array([1, 1, 7, ..., 2, 3, 3])
5360 >>> x = np.arange(8).reshape(2, 4)
5361 >>> idx = (1, 3)
5362 >>> np.insert(x, idx, 999, axis=1)
5363 array([[ 0, 999, 1, 2, 999, 3],
5364 [ 4, 999, 5, 6, 999, 7]])
5366 """
5367 wrap = None
5368 if type(arr) is not ndarray:
5369 try:
5370 wrap = arr.__array_wrap__
5371 except AttributeError:
5372 pass
5374 arr = asarray(arr)
5375 ndim = arr.ndim
5376 arrorder = 'F' if arr.flags.fnc else 'C'
5377 if axis is None:
5378 if ndim != 1:
5379 arr = arr.ravel()
5380 # needed for np.matrix, which is still not 1d after being ravelled
5381 ndim = arr.ndim
5382 axis = ndim - 1
5383 else:
5384 axis = normalize_axis_index(axis, ndim)
5385 slobj = [slice(None)]*ndim
5386 N = arr.shape[axis]
5387 newshape = list(arr.shape)
5389 if isinstance(obj, slice):
5390 # turn it into a range object
5391 indices = arange(*obj.indices(N), dtype=intp)
5392 else:
5393 # need to copy obj, because indices will be changed in-place
5394 indices = np.array(obj)
5395 if indices.dtype == bool:
5396 # See also delete
5397 # 2012-10-11, NumPy 1.8
5398 warnings.warn(
5399 "in the future insert will treat boolean arrays and "
5400 "array-likes as a boolean index instead of casting it to "
5401 "integer", FutureWarning, stacklevel=3)
5402 indices = indices.astype(intp)
5403 # Code after warning period:
5404 #if obj.ndim != 1:
5405 # raise ValueError('boolean array argument obj to insert '
5406 # 'must be one dimensional')
5407 #indices = np.flatnonzero(obj)
5408 elif indices.ndim > 1:
5409 raise ValueError(
5410 "index array argument obj to insert must be one dimensional "
5411 "or scalar")
5412 if indices.size == 1:
5413 index = indices.item()
5414 if index < -N or index > N:
5415 raise IndexError(f"index {obj} is out of bounds for axis {axis} "
5416 f"with size {N}")
5417 if (index < 0):
5418 index += N
5420 # There are some object array corner cases here, but we cannot avoid
5421 # that:
5422 values = array(values, copy=False, ndmin=arr.ndim, dtype=arr.dtype)
5423 if indices.ndim == 0:
5424 # broadcasting is very different here, since a[:,0,:] = ... behaves
5425 # very different from a[:,[0],:] = ...! This changes values so that
5426 # it works likes the second case. (here a[:,0:1,:])
5427 values = np.moveaxis(values, 0, axis)
5428 numnew = values.shape[axis]
5429 newshape[axis] += numnew
5430 new = empty(newshape, arr.dtype, arrorder)
5431 slobj[axis] = slice(None, index)
5432 new[tuple(slobj)] = arr[tuple(slobj)]
5433 slobj[axis] = slice(index, index+numnew)
5434 new[tuple(slobj)] = values
5435 slobj[axis] = slice(index+numnew, None)
5436 slobj2 = [slice(None)] * ndim
5437 slobj2[axis] = slice(index, None)
5438 new[tuple(slobj)] = arr[tuple(slobj2)]
5439 if wrap:
5440 return wrap(new)
5441 return new
5442 elif indices.size == 0 and not isinstance(obj, np.ndarray):
5443 # Can safely cast the empty list to intp
5444 indices = indices.astype(intp)
5446 indices[indices < 0] += N
5448 numnew = len(indices)
5449 order = indices.argsort(kind='mergesort') # stable sort
5450 indices[order] += np.arange(numnew)
5452 newshape[axis] += numnew
5453 old_mask = ones(newshape[axis], dtype=bool)
5454 old_mask[indices] = False
5456 new = empty(newshape, arr.dtype, arrorder)
5457 slobj2 = [slice(None)]*ndim
5458 slobj[axis] = indices
5459 slobj2[axis] = old_mask
5460 new[tuple(slobj)] = values
5461 new[tuple(slobj2)] = arr
5463 if wrap:
5464 return wrap(new)
5465 return new
5468def _append_dispatcher(arr, values, axis=None):
5469 return (arr, values)
5472@array_function_dispatch(_append_dispatcher)
5473def append(arr, values, axis=None):
5474 """
5475 Append values to the end of an array.
5477 Parameters
5478 ----------
5479 arr : array_like
5480 Values are appended to a copy of this array.
5481 values : array_like
5482 These values are appended to a copy of `arr`. It must be of the
5483 correct shape (the same shape as `arr`, excluding `axis`). If
5484 `axis` is not specified, `values` can be any shape and will be
5485 flattened before use.
5486 axis : int, optional
5487 The axis along which `values` are appended. If `axis` is not
5488 given, both `arr` and `values` are flattened before use.
5490 Returns
5491 -------
5492 append : ndarray
5493 A copy of `arr` with `values` appended to `axis`. Note that
5494 `append` does not occur in-place: a new array is allocated and
5495 filled. If `axis` is None, `out` is a flattened array.
5497 See Also
5498 --------
5499 insert : Insert elements into an array.
5500 delete : Delete elements from an array.
5502 Examples
5503 --------
5504 >>> np.append([1, 2, 3], [[4, 5, 6], [7, 8, 9]])
5505 array([1, 2, 3, ..., 7, 8, 9])
5507 When `axis` is specified, `values` must have the correct shape.
5509 >>> np.append([[1, 2, 3], [4, 5, 6]], [[7, 8, 9]], axis=0)
5510 array([[1, 2, 3],
5511 [4, 5, 6],
5512 [7, 8, 9]])
5513 >>> np.append([[1, 2, 3], [4, 5, 6]], [7, 8, 9], axis=0)
5514 Traceback (most recent call last):
5515 ...
5516 ValueError: all the input arrays must have same number of dimensions, but
5517 the array at index 0 has 2 dimension(s) and the array at index 1 has 1
5518 dimension(s)
5520 """
5521 arr = asanyarray(arr)
5522 if axis is None:
5523 if arr.ndim != 1:
5524 arr = arr.ravel()
5525 values = ravel(values)
5526 axis = arr.ndim-1
5527 return concatenate((arr, values), axis=axis)
5530def _digitize_dispatcher(x, bins, right=None):
5531 return (x, bins)
5534@array_function_dispatch(_digitize_dispatcher)
5535def digitize(x, bins, right=False):
5536 """
5537 Return the indices of the bins to which each value in input array belongs.
5539 ========= ============= ============================
5540 `right` order of bins returned index `i` satisfies
5541 ========= ============= ============================
5542 ``False`` increasing ``bins[i-1] <= x < bins[i]``
5543 ``True`` increasing ``bins[i-1] < x <= bins[i]``
5544 ``False`` decreasing ``bins[i-1] > x >= bins[i]``
5545 ``True`` decreasing ``bins[i-1] >= x > bins[i]``
5546 ========= ============= ============================
5548 If values in `x` are beyond the bounds of `bins`, 0 or ``len(bins)`` is
5549 returned as appropriate.
5551 Parameters
5552 ----------
5553 x : array_like
5554 Input array to be binned. Prior to NumPy 1.10.0, this array had to
5555 be 1-dimensional, but can now have any shape.
5556 bins : array_like
5557 Array of bins. It has to be 1-dimensional and monotonic.
5558 right : bool, optional
5559 Indicating whether the intervals include the right or the left bin
5560 edge. Default behavior is (right==False) indicating that the interval
5561 does not include the right edge. The left bin end is open in this
5562 case, i.e., bins[i-1] <= x < bins[i] is the default behavior for
5563 monotonically increasing bins.
5565 Returns
5566 -------
5567 indices : ndarray of ints
5568 Output array of indices, of same shape as `x`.
5570 Raises
5571 ------
5572 ValueError
5573 If `bins` is not monotonic.
5574 TypeError
5575 If the type of the input is complex.
5577 See Also
5578 --------
5579 bincount, histogram, unique, searchsorted
5581 Notes
5582 -----
5583 If values in `x` are such that they fall outside the bin range,
5584 attempting to index `bins` with the indices that `digitize` returns
5585 will result in an IndexError.
5587 .. versionadded:: 1.10.0
5589 `np.digitize` is implemented in terms of `np.searchsorted`. This means
5590 that a binary search is used to bin the values, which scales much better
5591 for larger number of bins than the previous linear search. It also removes
5592 the requirement for the input array to be 1-dimensional.
5594 For monotonically _increasing_ `bins`, the following are equivalent::
5596 np.digitize(x, bins, right=True)
5597 np.searchsorted(bins, x, side='left')
5599 Note that as the order of the arguments are reversed, the side must be too.
5600 The `searchsorted` call is marginally faster, as it does not do any
5601 monotonicity checks. Perhaps more importantly, it supports all dtypes.
5603 Examples
5604 --------
5605 >>> x = np.array([0.2, 6.4, 3.0, 1.6])
5606 >>> bins = np.array([0.0, 1.0, 2.5, 4.0, 10.0])
5607 >>> inds = np.digitize(x, bins)
5608 >>> inds
5609 array([1, 4, 3, 2])
5610 >>> for n in range(x.size):
5611 ... print(bins[inds[n]-1], "<=", x[n], "<", bins[inds[n]])
5612 ...
5613 0.0 <= 0.2 < 1.0
5614 4.0 <= 6.4 < 10.0
5615 2.5 <= 3.0 < 4.0
5616 1.0 <= 1.6 < 2.5
5618 >>> x = np.array([1.2, 10.0, 12.4, 15.5, 20.])
5619 >>> bins = np.array([0, 5, 10, 15, 20])
5620 >>> np.digitize(x,bins,right=True)
5621 array([1, 2, 3, 4, 4])
5622 >>> np.digitize(x,bins,right=False)
5623 array([1, 3, 3, 4, 5])
5624 """
5625 x = _nx.asarray(x)
5626 bins = _nx.asarray(bins)
5628 # here for compatibility, searchsorted below is happy to take this
5629 if np.issubdtype(x.dtype, _nx.complexfloating):
5630 raise TypeError("x may not be complex")
5632 mono = _monotonicity(bins)
5633 if mono == 0:
5634 raise ValueError("bins must be monotonically increasing or decreasing")
5636 # this is backwards because the arguments below are swapped
5637 side = 'left' if right else 'right'
5638 if mono == -1:
5639 # reverse the bins, and invert the results
5640 return len(bins) - _nx.searchsorted(bins[::-1], x, side=side)
5641 else:
5642 return _nx.searchsorted(bins, x, side=side)