Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1import builtins
2import collections.abc
3import functools
4import re
5import warnings
7import numpy as np
8import numpy._core.numeric as _nx
9from numpy._core import overrides, transpose
10from numpy._core._multiarray_umath import _array_converter
11from numpy._core.fromnumeric import any, mean, nonzero, partition, ravel, sum
12from numpy._core.multiarray import (
13 _monotonicity,
14 _place,
15 bincount,
16 interp as compiled_interp,
17 interp_complex as compiled_interp_complex,
18 normalize_axis_index,
19)
20from numpy._core.numeric import (
21 absolute,
22 arange,
23 array,
24 asanyarray,
25 asarray,
26 concatenate,
27 dot,
28 empty,
29 integer,
30 intp,
31 isscalar,
32 ndarray,
33 ones,
34 take,
35 where,
36 zeros_like,
37)
38from numpy._core.numerictypes import typecodes
39from numpy._core.umath import (
40 add,
41 arctan2,
42 cos,
43 exp,
44 floor,
45 frompyfunc,
46 less_equal,
47 minimum,
48 mod,
49 not_equal,
50 pi,
51 sin,
52 sqrt,
53 subtract,
54)
55from numpy._utils import set_module
57# needed in this module for compatibility
58from numpy.lib._histograms_impl import histogram, histogramdd # noqa: F401
59from numpy.lib._twodim_base_impl import diag
61array_function_dispatch = functools.partial(
62 overrides.array_function_dispatch, module='numpy')
65__all__ = [
66 'select', 'piecewise', 'trim_zeros', 'copy', 'iterable', 'percentile',
67 'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'flip',
68 'rot90', 'extract', 'place', 'vectorize', 'asarray_chkfinite', 'average',
69 'bincount', 'digitize', 'cov', 'corrcoef',
70 'median', 'sinc', 'hamming', 'hanning', 'bartlett',
71 'blackman', 'kaiser', 'trapezoid', 'i0',
72 'meshgrid', 'delete', 'insert', 'append', 'interp',
73 'quantile'
74 ]
76# _QuantileMethods is a dictionary listing all the supported methods to
77# compute quantile/percentile.
78#
79# Below virtual_index refers to the index of the element where the percentile
80# would be found in the sorted sample.
81# When the sample contains exactly the percentile wanted, the virtual_index is
82# an integer to the index of this element.
83# When the percentile wanted is in between two elements, the virtual_index
84# is made of a integer part (a.k.a 'i' or 'left') and a fractional part
85# (a.k.a 'g' or 'gamma')
86#
87# Each method in _QuantileMethods has two properties
88# get_virtual_index : Callable
89# The function used to compute the virtual_index.
90# fix_gamma : Callable
91# A function used for discrete methods to force the index to a specific value.
92_QuantileMethods = {
93 # --- HYNDMAN and FAN METHODS
94 # Discrete methods
95 'inverted_cdf': {
96 'get_virtual_index': lambda n, quantiles: _inverted_cdf(n, quantiles), # noqa: PLW0108
97 'fix_gamma': None, # should never be called
98 },
99 'averaged_inverted_cdf': {
100 'get_virtual_index': lambda n, quantiles: (n * quantiles) - 1,
101 'fix_gamma': lambda gamma, _: _get_gamma_mask(
102 shape=gamma.shape,
103 default_value=1.,
104 conditioned_value=0.5,
105 where=gamma == 0),
106 },
107 'closest_observation': {
108 'get_virtual_index': lambda n, quantiles: _closest_observation(n, quantiles), # noqa: PLW0108
109 'fix_gamma': None, # should never be called
110 },
111 # Continuous methods
112 'interpolated_inverted_cdf': {
113 'get_virtual_index': lambda n, quantiles:
114 _compute_virtual_index(n, quantiles, 0, 1),
115 'fix_gamma': lambda gamma, _: gamma,
116 },
117 'hazen': {
118 'get_virtual_index': lambda n, quantiles:
119 _compute_virtual_index(n, quantiles, 0.5, 0.5),
120 'fix_gamma': lambda gamma, _: gamma,
121 },
122 'weibull': {
123 'get_virtual_index': lambda n, quantiles:
124 _compute_virtual_index(n, quantiles, 0, 0),
125 'fix_gamma': lambda gamma, _: gamma,
126 },
127 # Default method.
128 # To avoid some rounding issues, `(n-1) * quantiles` is preferred to
129 # `_compute_virtual_index(n, quantiles, 1, 1)`.
130 # They are mathematically equivalent.
131 'linear': {
132 'get_virtual_index': lambda n, quantiles: (n - 1) * quantiles,
133 'fix_gamma': lambda gamma, _: gamma,
134 },
135 'median_unbiased': {
136 'get_virtual_index': lambda n, quantiles:
137 _compute_virtual_index(n, quantiles, 1 / 3.0, 1 / 3.0),
138 'fix_gamma': lambda gamma, _: gamma,
139 },
140 'normal_unbiased': {
141 'get_virtual_index': lambda n, quantiles:
142 _compute_virtual_index(n, quantiles, 3 / 8.0, 3 / 8.0),
143 'fix_gamma': lambda gamma, _: gamma,
144 },
145 # --- OTHER METHODS
146 'lower': {
147 'get_virtual_index': lambda n, quantiles: np.floor(
148 (n - 1) * quantiles).astype(np.intp),
149 'fix_gamma': None, # should never be called, index dtype is int
150 },
151 'higher': {
152 'get_virtual_index': lambda n, quantiles: np.ceil(
153 (n - 1) * quantiles).astype(np.intp),
154 'fix_gamma': None, # should never be called, index dtype is int
155 },
156 'midpoint': {
157 'get_virtual_index': lambda n, quantiles: 0.5 * (
158 np.floor((n - 1) * quantiles)
159 + np.ceil((n - 1) * quantiles)),
160 'fix_gamma': lambda gamma, index: _get_gamma_mask(
161 shape=gamma.shape,
162 default_value=0.5,
163 conditioned_value=0.,
164 where=index % 1 == 0),
165 },
166 'nearest': {
167 'get_virtual_index': lambda n, quantiles: np.around(
168 (n - 1) * quantiles).astype(np.intp),
169 'fix_gamma': None,
170 # should never be called, index dtype is int
171 }}
174def _rot90_dispatcher(m, k=None, axes=None):
175 return (m,)
178@array_function_dispatch(_rot90_dispatcher)
179def rot90(m, k=1, axes=(0, 1)):
180 """
181 Rotate an array by 90 degrees in the plane specified by axes.
183 Rotation direction is from the first towards the second axis.
184 This means for a 2D array with the default `k` and `axes`, the
185 rotation will be counterclockwise.
187 Parameters
188 ----------
189 m : array_like
190 Array of two or more dimensions.
191 k : integer
192 Number of times the array is rotated by 90 degrees.
193 axes : (2,) array_like
194 The array is rotated in the plane defined by the axes.
195 Axes must be different.
197 Returns
198 -------
199 y : ndarray
200 A rotated view of `m`.
202 See Also
203 --------
204 flip : Reverse the order of elements in an array along the given axis.
205 fliplr : Flip an array horizontally.
206 flipud : Flip an array vertically.
208 Notes
209 -----
210 ``rot90(m, k=1, axes=(1,0))`` is the reverse of
211 ``rot90(m, k=1, axes=(0,1))``
213 ``rot90(m, k=1, axes=(1,0))`` is equivalent to
214 ``rot90(m, k=-1, axes=(0,1))``
216 Examples
217 --------
218 >>> import numpy as np
219 >>> m = np.array([[1,2],[3,4]], int)
220 >>> m
221 array([[1, 2],
222 [3, 4]])
223 >>> np.rot90(m)
224 array([[2, 4],
225 [1, 3]])
226 >>> np.rot90(m, 2)
227 array([[4, 3],
228 [2, 1]])
229 >>> m = np.arange(8).reshape((2,2,2))
230 >>> np.rot90(m, 1, (1,2))
231 array([[[1, 3],
232 [0, 2]],
233 [[5, 7],
234 [4, 6]]])
236 """
237 axes = tuple(axes)
238 if len(axes) != 2:
239 raise ValueError("len(axes) must be 2.")
241 m = asanyarray(m)
243 if axes[0] == axes[1] or absolute(axes[0] - axes[1]) == m.ndim:
244 raise ValueError("Axes must be different.")
246 if (axes[0] >= m.ndim or axes[0] < -m.ndim
247 or axes[1] >= m.ndim or axes[1] < -m.ndim):
248 raise ValueError(f"Axes={axes} out of range for array of ndim={m.ndim}.")
250 k %= 4
252 if k == 0:
253 return m[:]
254 if k == 2:
255 return flip(flip(m, axes[0]), axes[1])
257 axes_list = arange(0, m.ndim)
258 (axes_list[axes[0]], axes_list[axes[1]]) = (axes_list[axes[1]],
259 axes_list[axes[0]])
261 if k == 1:
262 return transpose(flip(m, axes[1]), axes_list)
263 else:
264 # k == 3
265 return flip(transpose(m, axes_list), axes[1])
268def _flip_dispatcher(m, axis=None):
269 return (m,)
272@array_function_dispatch(_flip_dispatcher)
273def flip(m, axis=None):
274 """
275 Reverse the order of elements in an array along the given axis.
277 The shape of the array is preserved, but the elements are reordered.
279 Parameters
280 ----------
281 m : array_like
282 Input array.
283 axis : None or int or tuple of ints, optional
284 Axis or axes along which to flip over. The default,
285 axis=None, will flip over all of the axes of the input array.
286 If axis is negative it counts from the last to the first axis.
288 If axis is a tuple of ints, flipping is performed on all of the axes
289 specified in the tuple.
291 Returns
292 -------
293 out : array_like
294 A view of `m` with the entries of axis reversed. Since a view is
295 returned, this operation is done in constant time.
297 See Also
298 --------
299 flipud : Flip an array vertically (axis=0).
300 fliplr : Flip an array horizontally (axis=1).
302 Notes
303 -----
304 flip(m, 0) is equivalent to flipud(m).
306 flip(m, 1) is equivalent to fliplr(m).
308 flip(m, n) corresponds to ``m[...,::-1,...]`` with ``::-1`` at position n.
310 flip(m) corresponds to ``m[::-1,::-1,...,::-1]`` with ``::-1`` at all
311 positions.
313 flip(m, (0, 1)) corresponds to ``m[::-1,::-1,...]`` with ``::-1`` at
314 position 0 and position 1.
316 Examples
317 --------
318 >>> import numpy as np
319 >>> A = np.arange(8).reshape((2,2,2))
320 >>> A
321 array([[[0, 1],
322 [2, 3]],
323 [[4, 5],
324 [6, 7]]])
325 >>> np.flip(A, 0)
326 array([[[4, 5],
327 [6, 7]],
328 [[0, 1],
329 [2, 3]]])
330 >>> np.flip(A, 1)
331 array([[[2, 3],
332 [0, 1]],
333 [[6, 7],
334 [4, 5]]])
335 >>> np.flip(A)
336 array([[[7, 6],
337 [5, 4]],
338 [[3, 2],
339 [1, 0]]])
340 >>> np.flip(A, (0, 2))
341 array([[[5, 4],
342 [7, 6]],
343 [[1, 0],
344 [3, 2]]])
345 >>> rng = np.random.default_rng()
346 >>> A = rng.normal(size=(3,4,5))
347 >>> np.all(np.flip(A,2) == A[:,:,::-1,...])
348 True
349 """
350 if not hasattr(m, 'ndim'):
351 m = asarray(m)
352 if axis is None:
353 indexer = (np.s_[::-1],) * m.ndim
354 else:
355 axis = _nx.normalize_axis_tuple(axis, m.ndim)
356 indexer = [np.s_[:]] * m.ndim
357 for ax in axis:
358 indexer[ax] = np.s_[::-1]
359 indexer = tuple(indexer)
360 return m[indexer]
363@set_module('numpy')
364def iterable(y):
365 """
366 Check whether or not an object can be iterated over.
368 Parameters
369 ----------
370 y : object
371 Input object.
373 Returns
374 -------
375 b : bool
376 Return ``True`` if the object has an iterator method or is a
377 sequence and ``False`` otherwise.
380 Examples
381 --------
382 >>> import numpy as np
383 >>> np.iterable([1, 2, 3])
384 True
385 >>> np.iterable(2)
386 False
388 Notes
389 -----
390 In most cases, the results of ``np.iterable(obj)`` are consistent with
391 ``isinstance(obj, collections.abc.Iterable)``. One notable exception is
392 the treatment of 0-dimensional arrays::
394 >>> from collections.abc import Iterable
395 >>> a = np.array(1.0) # 0-dimensional numpy array
396 >>> isinstance(a, Iterable)
397 True
398 >>> np.iterable(a)
399 False
401 """
402 try:
403 iter(y)
404 except TypeError:
405 return False
406 return True
409def _weights_are_valid(weights, a, axis):
410 """Validate weights array.
412 We assume, weights is not None.
413 """
414 wgt = np.asanyarray(weights)
416 # Sanity checks
417 if a.shape != wgt.shape:
418 if axis is None:
419 raise TypeError(
420 "Axis must be specified when shapes of a and weights "
421 "differ.")
422 if wgt.shape != tuple(a.shape[ax] for ax in axis):
423 raise ValueError(
424 "Shape of weights must be consistent with "
425 "shape of a along specified axis.")
427 # setup wgt to broadcast along axis
428 wgt = wgt.transpose(np.argsort(axis))
429 wgt = wgt.reshape(tuple((s if ax in axis else 1)
430 for ax, s in enumerate(a.shape)))
431 return wgt
434def _average_dispatcher(a, axis=None, weights=None, returned=None, *,
435 keepdims=None):
436 return (a, weights)
439@array_function_dispatch(_average_dispatcher)
440def average(a, axis=None, weights=None, returned=False, *,
441 keepdims=np._NoValue):
442 """
443 Compute the weighted average along the specified axis.
445 Parameters
446 ----------
447 a : array_like
448 Array containing data to be averaged. If `a` is not an array, a
449 conversion is attempted.
450 axis : None or int or tuple of ints, optional
451 Axis or axes along which to average `a`. The default,
452 `axis=None`, will average over all of the elements of the input array.
453 If axis is negative it counts from the last to the first axis.
454 If axis is a tuple of ints, averaging is performed on all of the axes
455 specified in the tuple instead of a single axis or all the axes as
456 before.
457 weights : array_like, optional
458 An array of weights associated with the values in `a`. Each value in
459 `a` contributes to the average according to its associated weight.
460 The array of weights must be the same shape as `a` if no axis is
461 specified, otherwise the weights must have dimensions and shape
462 consistent with `a` along the specified axis.
463 If `weights=None`, then all data in `a` are assumed to have a
464 weight equal to one.
465 The calculation is::
467 avg = sum(a * weights) / sum(weights)
469 where the sum is over all included elements.
470 The only constraint on the values of `weights` is that `sum(weights)`
471 must not be 0.
472 returned : bool, optional
473 Default is `False`. If `True`, the tuple (`average`, `sum_of_weights`)
474 is returned, otherwise only the average is returned.
475 If `weights=None`, `sum_of_weights` is equivalent to the number of
476 elements over which the average is taken.
477 keepdims : bool, optional
478 If this is set to True, the axes which are reduced are left
479 in the result as dimensions with size one. With this option,
480 the result will broadcast correctly against the original `a`.
481 *Note:* `keepdims` will not work with instances of `numpy.matrix`
482 or other classes whose methods do not support `keepdims`.
484 .. versionadded:: 1.23.0
486 Returns
487 -------
488 retval, [sum_of_weights] : array_type or double
489 Return the average along the specified axis. When `returned` is `True`,
490 return a tuple with the average as the first element and the sum
491 of the weights as the second element. `sum_of_weights` is of the
492 same type as `retval`. The result dtype follows a general pattern.
493 If `weights` is None, the result dtype will be that of `a` , or ``float64``
494 if `a` is integral. Otherwise, if `weights` is not None and `a` is non-
495 integral, the result type will be the type of lowest precision capable of
496 representing values of both `a` and `weights`. If `a` happens to be
497 integral, the previous rules still applies but the result dtype will
498 at least be ``float64``.
500 Raises
501 ------
502 ZeroDivisionError
503 When all weights along axis are zero. See `numpy.ma.average` for a
504 version robust to this type of error.
505 TypeError
506 When `weights` does not have the same shape as `a`, and `axis=None`.
507 ValueError
508 When `weights` does not have dimensions and shape consistent with `a`
509 along specified `axis`.
511 See Also
512 --------
513 mean
515 ma.average : average for masked arrays -- useful if your data contains
516 "missing" values
517 numpy.result_type : Returns the type that results from applying the
518 numpy type promotion rules to the arguments.
520 Examples
521 --------
522 >>> import numpy as np
523 >>> data = np.arange(1, 5)
524 >>> data
525 array([1, 2, 3, 4])
526 >>> np.average(data)
527 2.5
528 >>> np.average(np.arange(1, 11), weights=np.arange(10, 0, -1))
529 4.0
531 >>> data = np.arange(6).reshape((3, 2))
532 >>> data
533 array([[0, 1],
534 [2, 3],
535 [4, 5]])
536 >>> np.average(data, axis=1, weights=[1./4, 3./4])
537 array([0.75, 2.75, 4.75])
538 >>> np.average(data, weights=[1./4, 3./4])
539 Traceback (most recent call last):
540 ...
541 TypeError: Axis must be specified when shapes of a and weights differ.
543 With ``keepdims=True``, the following result has shape (3, 1).
545 >>> np.average(data, axis=1, keepdims=True)
546 array([[0.5],
547 [2.5],
548 [4.5]])
550 >>> data = np.arange(8).reshape((2, 2, 2))
551 >>> data
552 array([[[0, 1],
553 [2, 3]],
554 [[4, 5],
555 [6, 7]]])
556 >>> np.average(data, axis=(0, 1), weights=[[1./4, 3./4], [1., 1./2]])
557 array([3.4, 4.4])
558 >>> np.average(data, axis=0, weights=[[1./4, 3./4], [1., 1./2]])
559 Traceback (most recent call last):
560 ...
561 ValueError: Shape of weights must be consistent
562 with shape of a along specified axis.
563 """
564 a = np.asanyarray(a)
566 if axis is not None:
567 axis = _nx.normalize_axis_tuple(axis, a.ndim, argname="axis")
569 if keepdims is np._NoValue:
570 # Don't pass on the keepdims argument if one wasn't given.
571 keepdims_kw = {}
572 else:
573 keepdims_kw = {'keepdims': keepdims}
575 if weights is None:
576 avg = a.mean(axis, **keepdims_kw)
577 avg_as_array = np.asanyarray(avg)
578 scl = avg_as_array.dtype.type(a.size / avg_as_array.size)
579 else:
580 wgt = _weights_are_valid(weights=weights, a=a, axis=axis)
582 if issubclass(a.dtype.type, (np.integer, np.bool)):
583 result_dtype = np.result_type(a.dtype, wgt.dtype, 'f8')
584 else:
585 result_dtype = np.result_type(a.dtype, wgt.dtype)
587 scl = wgt.sum(axis=axis, dtype=result_dtype, **keepdims_kw)
588 if np.any(scl == 0.0):
589 raise ZeroDivisionError(
590 "Weights sum to zero, can't be normalized")
592 avg = avg_as_array = np.multiply(a, wgt,
593 dtype=result_dtype).sum(axis, **keepdims_kw) / scl
595 if returned:
596 if scl.shape != avg_as_array.shape:
597 scl = np.broadcast_to(scl, avg_as_array.shape, subok=True).copy()
598 return avg, scl
599 else:
600 return avg
603@set_module('numpy')
604def asarray_chkfinite(a, dtype=None, order=None):
605 """Convert the input to an array, checking for NaNs or Infs.
607 Parameters
608 ----------
609 a : array_like
610 Input data, in any form that can be converted to an array. This
611 includes lists, lists of tuples, tuples, tuples of tuples, tuples
612 of lists and ndarrays. Success requires no NaNs or Infs.
613 dtype : data-type, optional
614 By default, the data-type is inferred from the input data.
615 order : {'C', 'F', 'A', 'K'}, optional
616 The memory layout of the output.
617 'C' gives a row-major layout (C-style),
618 'F' gives a column-major layout (Fortran-style).
619 'C' and 'F' will copy if needed to ensure the output format.
620 'A' (any) is equivalent to 'F' if input a is non-contiguous or
621 Fortran-contiguous, otherwise, it is equivalent to 'C'.
622 Unlike 'C' or 'F', 'A' does not ensure that the result is contiguous.
623 'K' (keep) preserves the input order for the output.
624 'C' is the default.
626 Returns
627 -------
628 out : ndarray
629 Array interpretation of `a`. No copy is performed if the input
630 is already an ndarray. If `a` is a subclass of ndarray, a base
631 class ndarray is returned.
633 Raises
634 ------
635 ValueError
636 Raises ValueError if `a` contains NaN (Not a Number) or Inf (Infinity).
638 See Also
639 --------
640 asarray : Create and array.
641 asanyarray : Similar function which passes through subclasses.
642 ascontiguousarray : Convert input to a contiguous array.
643 asfortranarray : Convert input to an ndarray with column-major
644 memory order.
645 fromiter : Create an array from an iterator.
646 fromfunction : Construct an array by executing a function on grid
647 positions.
649 Examples
650 --------
651 >>> import numpy as np
653 Convert a list into an array. If all elements are finite, then
654 ``asarray_chkfinite`` is identical to ``asarray``.
656 >>> a = [1, 2]
657 >>> np.asarray_chkfinite(a, dtype=float)
658 array([1., 2.])
660 Raises ValueError if array_like contains Nans or Infs.
662 >>> a = [1, 2, np.inf]
663 >>> try:
664 ... np.asarray_chkfinite(a)
665 ... except ValueError:
666 ... print('ValueError')
667 ...
668 ValueError
670 """
671 a = asarray(a, dtype=dtype, order=order)
672 if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
673 raise ValueError(
674 "array must not contain infs or NaNs")
675 return a
678def _piecewise_dispatcher(x, condlist, funclist, *args, **kw):
679 yield x
680 # support the undocumented behavior of allowing scalars
681 if np.iterable(condlist):
682 yield from condlist
685@array_function_dispatch(_piecewise_dispatcher)
686def piecewise(x, condlist, funclist, *args, **kw):
687 """
688 Evaluate a piecewise-defined function.
690 Given a set of conditions and corresponding functions, evaluate each
691 function on the input data wherever its condition is true.
693 Parameters
694 ----------
695 x : ndarray or scalar
696 The input domain.
697 condlist : list of bool arrays or bool scalars
698 Each boolean array corresponds to a function in `funclist`. Wherever
699 `condlist[i]` is True, `funclist[i](x)` is used as the output value.
701 Each boolean array in `condlist` selects a piece of `x`,
702 and should therefore be of the same shape as `x`.
704 The length of `condlist` must correspond to that of `funclist`.
705 If one extra function is given, i.e. if
706 ``len(funclist) == len(condlist) + 1``, then that extra function
707 is the default value, used wherever all conditions are false.
708 funclist : list of callables, f(x,*args,**kw), or scalars
709 Each function is evaluated over `x` wherever its corresponding
710 condition is True. It should take a 1d array as input and give a 1d
711 array or a scalar value as output. If, instead of a callable,
712 a scalar is provided then a constant function (``lambda x: scalar``) is
713 assumed.
714 args : tuple, optional
715 Any further arguments given to `piecewise` are passed to the functions
716 upon execution, i.e., if called ``piecewise(..., ..., 1, 'a')``, then
717 each function is called as ``f(x, 1, 'a')``.
718 kw : dict, optional
719 Keyword arguments used in calling `piecewise` are passed to the
720 functions upon execution, i.e., if called
721 ``piecewise(..., ..., alpha=1)``, then each function is called as
722 ``f(x, alpha=1)``.
724 Returns
725 -------
726 out : ndarray
727 The output is the same shape and type as x and is found by
728 calling the functions in `funclist` on the appropriate portions of `x`,
729 as defined by the boolean arrays in `condlist`. Portions not covered
730 by any condition have a default value of 0.
733 See Also
734 --------
735 choose, select, where
737 Notes
738 -----
739 This is similar to choose or select, except that functions are
740 evaluated on elements of `x` that satisfy the corresponding condition from
741 `condlist`.
743 The result is::
745 |--
746 |funclist[0](x[condlist[0]])
747 out = |funclist[1](x[condlist[1]])
748 |...
749 |funclist[n2](x[condlist[n2]])
750 |--
752 Examples
753 --------
754 >>> import numpy as np
756 Define the signum function, which is -1 for ``x < 0`` and +1 for ``x >= 0``.
758 >>> x = np.linspace(-2.5, 2.5, 6)
759 >>> np.piecewise(x, [x < 0, x >= 0], [-1, 1])
760 array([-1., -1., -1., 1., 1., 1.])
762 Define the absolute value, which is ``-x`` for ``x <0`` and ``x`` for
763 ``x >= 0``.
765 >>> np.piecewise(x, [x < 0, x >= 0], [lambda x: -x, lambda x: x])
766 array([2.5, 1.5, 0.5, 0.5, 1.5, 2.5])
768 Apply the same function to a scalar value.
770 >>> y = -2
771 >>> np.piecewise(y, [y < 0, y >= 0], [lambda x: -x, lambda x: x])
772 array(2)
774 """
775 x = asanyarray(x)
776 n2 = len(funclist)
778 # undocumented: single condition is promoted to a list of one condition
779 if isscalar(condlist) or (
780 not isinstance(condlist[0], (list, ndarray)) and x.ndim != 0):
781 condlist = [condlist]
783 condlist = asarray(condlist, dtype=bool)
784 n = len(condlist)
786 if n == n2 - 1: # compute the "otherwise" condition.
787 condelse = ~np.any(condlist, axis=0, keepdims=True)
788 condlist = np.concatenate([condlist, condelse], axis=0)
789 n += 1
790 elif n != n2:
791 raise ValueError(
792 f"with {n} condition(s), either {n} or {n + 1} functions are expected"
793 )
795 y = zeros_like(x)
796 for cond, func in zip(condlist, funclist):
797 if not isinstance(func, collections.abc.Callable):
798 y[cond] = func
799 else:
800 vals = x[cond]
801 if vals.size > 0:
802 y[cond] = func(vals, *args, **kw)
804 return y
807def _select_dispatcher(condlist, choicelist, default=None):
808 yield from condlist
809 yield from choicelist
812@array_function_dispatch(_select_dispatcher)
813def select(condlist, choicelist, default=0):
814 """
815 Return an array drawn from elements in choicelist, depending on conditions.
817 Parameters
818 ----------
819 condlist : list of bool ndarrays
820 The list of conditions which determine from which array in `choicelist`
821 the output elements are taken. When multiple conditions are satisfied,
822 the first one encountered in `condlist` is used.
823 choicelist : list of ndarrays
824 The list of arrays from which the output elements are taken. It has
825 to be of the same length as `condlist`.
826 default : array_like, optional
827 The element inserted in `output` when all conditions evaluate to False.
829 Returns
830 -------
831 output : ndarray
832 The output at position m is the m-th element of the array in
833 `choicelist` where the m-th element of the corresponding array in
834 `condlist` is True.
836 See Also
837 --------
838 where : Return elements from one of two arrays depending on condition.
839 take, choose, compress, diag, diagonal
841 Examples
842 --------
843 >>> import numpy as np
845 Beginning with an array of integers from 0 to 5 (inclusive),
846 elements less than ``3`` are negated, elements greater than ``3``
847 are squared, and elements not meeting either of these conditions
848 (exactly ``3``) are replaced with a `default` value of ``42``.
850 >>> x = np.arange(6)
851 >>> condlist = [x<3, x>3]
852 >>> choicelist = [-x, x**2]
853 >>> np.select(condlist, choicelist, 42)
854 array([ 0, -1, -2, 42, 16, 25])
856 When multiple conditions are satisfied, the first one encountered in
857 `condlist` is used.
859 >>> condlist = [x<=4, x>3]
860 >>> choicelist = [x, x**2]
861 >>> np.select(condlist, choicelist, 55)
862 array([ 0, 1, 2, 3, 4, 25])
864 """
865 # Check the size of condlist and choicelist are the same, or abort.
866 if len(condlist) != len(choicelist):
867 raise ValueError(
868 'list of cases must be same length as list of conditions')
870 # Now that the dtype is known, handle the deprecated select([], []) case
871 if len(condlist) == 0:
872 raise ValueError("select with an empty condition list is not possible")
874 # TODO: This preserves the Python int, float, complex manually to get the
875 # right `result_type` with NEP 50. Most likely we will grow a better
876 # way to spell this (and this can be replaced).
877 choicelist = [
878 choice if type(choice) in (int, float, complex) else np.asarray(choice)
879 for choice in choicelist]
880 choicelist.append(default if type(default) in (int, float, complex)
881 else np.asarray(default))
883 try:
884 dtype = np.result_type(*choicelist)
885 except TypeError as e:
886 msg = f'Choicelist and default value do not have a common dtype: {e}'
887 raise TypeError(msg) from None
889 # Convert conditions to arrays and broadcast conditions and choices
890 # as the shape is needed for the result. Doing it separately optimizes
891 # for example when all choices are scalars.
892 condlist = np.broadcast_arrays(*condlist)
893 choicelist = np.broadcast_arrays(*choicelist)
895 # If cond array is not an ndarray in boolean format or scalar bool, abort.
896 for i, cond in enumerate(condlist):
897 if cond.dtype.type is not np.bool:
898 raise TypeError(
899 f'invalid entry {i} in condlist: should be boolean ndarray')
901 if choicelist[0].ndim == 0:
902 # This may be common, so avoid the call.
903 result_shape = condlist[0].shape
904 else:
905 result_shape = np.broadcast_arrays(condlist[0], choicelist[0])[0].shape
907 result = np.full(result_shape, choicelist[-1], dtype)
909 # Use np.copyto to burn each choicelist array onto result, using the
910 # corresponding condlist as a boolean mask. This is done in reverse
911 # order since the first choice should take precedence.
912 choicelist = choicelist[-2::-1]
913 condlist = condlist[::-1]
914 for choice, cond in zip(choicelist, condlist):
915 np.copyto(result, choice, where=cond)
917 return result
920def _copy_dispatcher(a, order=None, subok=None):
921 return (a,)
924@array_function_dispatch(_copy_dispatcher)
925def copy(a, order='K', subok=False):
926 """
927 Return an array copy of the given object.
929 Parameters
930 ----------
931 a : array_like
932 Input data.
933 order : {'C', 'F', 'A', 'K'}, optional
934 Controls the memory layout of the copy. 'C' means C-order,
935 'F' means F-order, 'A' means 'F' if `a` is Fortran contiguous,
936 'C' otherwise. 'K' means match the layout of `a` as closely
937 as possible. (Note that this function and :meth:`ndarray.copy` are very
938 similar, but have different default values for their order=
939 arguments.)
940 subok : bool, optional
941 If True, then sub-classes will be passed-through, otherwise the
942 returned array will be forced to be a base-class array (defaults to False).
944 Returns
945 -------
946 arr : ndarray
947 Array interpretation of `a`.
949 See Also
950 --------
951 ndarray.copy : Preferred method for creating an array copy
953 Notes
954 -----
955 This is equivalent to:
957 >>> np.array(a, copy=True) #doctest: +SKIP
959 The copy made of the data is shallow, i.e., for arrays with object dtype,
960 the new array will point to the same objects.
961 See Examples from `ndarray.copy`.
963 Examples
964 --------
965 >>> import numpy as np
967 Create an array x, with a reference y and a copy z:
969 >>> x = np.array([1, 2, 3])
970 >>> y = x
971 >>> z = np.copy(x)
973 Note that, when we modify x, y changes, but not z:
975 >>> x[0] = 10
976 >>> x[0] == y[0]
977 True
978 >>> x[0] == z[0]
979 False
981 Note that, np.copy clears previously set WRITEABLE=False flag.
983 >>> a = np.array([1, 2, 3])
984 >>> a.flags["WRITEABLE"] = False
985 >>> b = np.copy(a)
986 >>> b.flags["WRITEABLE"]
987 True
988 >>> b[0] = 3
989 >>> b
990 array([3, 2, 3])
991 """
992 return array(a, order=order, subok=subok, copy=True)
994# Basic operations
997def _gradient_dispatcher(f, *varargs, axis=None, edge_order=None):
998 yield f
999 yield from varargs
1002@array_function_dispatch(_gradient_dispatcher)
1003def gradient(f, *varargs, axis=None, edge_order=1):
1004 """
1005 Return the gradient of an N-dimensional array.
1007 The gradient is computed using second order accurate central differences
1008 in the interior points and either first or second order accurate one-sides
1009 (forward or backwards) differences at the boundaries.
1010 The returned gradient hence has the same shape as the input array.
1012 Parameters
1013 ----------
1014 f : array_like
1015 An N-dimensional array containing samples of a scalar function.
1016 varargs : list of scalar or array, optional
1017 Spacing between f values. Default unitary spacing for all dimensions.
1018 Spacing can be specified using:
1020 1. single scalar to specify a sample distance for all dimensions.
1021 2. N scalars to specify a constant sample distance for each dimension.
1022 i.e. `dx`, `dy`, `dz`, ...
1023 3. N arrays to specify the coordinates of the values along each
1024 dimension of F. The length of the array must match the size of
1025 the corresponding dimension
1026 4. Any combination of N scalars/arrays with the meaning of 2. and 3.
1028 If `axis` is given, the number of varargs must equal the number of axes
1029 specified in the axis parameter.
1030 Default: 1. (see Examples below).
1032 edge_order : {1, 2}, optional
1033 Gradient is calculated using N-th order accurate differences
1034 at the boundaries. Default: 1.
1035 axis : None or int or tuple of ints, optional
1036 Gradient is calculated only along the given axis or axes
1037 The default (axis = None) is to calculate the gradient for all the axes
1038 of the input array. axis may be negative, in which case it counts from
1039 the last to the first axis.
1041 Returns
1042 -------
1043 gradient : ndarray or tuple of ndarray
1044 A tuple of ndarrays (or a single ndarray if there is only one
1045 dimension) corresponding to the derivatives of f with respect
1046 to each dimension. Each derivative has the same shape as f.
1048 Examples
1049 --------
1050 >>> import numpy as np
1051 >>> f = np.array([1, 2, 4, 7, 11, 16])
1052 >>> np.gradient(f)
1053 array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
1054 >>> np.gradient(f, 2)
1055 array([0.5 , 0.75, 1.25, 1.75, 2.25, 2.5 ])
1057 Spacing can be also specified with an array that represents the coordinates
1058 of the values F along the dimensions.
1059 For instance a uniform spacing:
1061 >>> x = np.arange(f.size)
1062 >>> np.gradient(f, x)
1063 array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
1065 Or a non uniform one:
1067 >>> x = np.array([0., 1., 1.5, 3.5, 4., 6.])
1068 >>> np.gradient(f, x)
1069 array([1. , 3. , 3.5, 6.7, 6.9, 2.5])
1071 For two dimensional arrays, the return will be two arrays ordered by
1072 axis. In this example the first array stands for the gradient in
1073 rows and the second one in columns direction:
1075 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]]))
1076 (array([[ 2., 2., -1.],
1077 [ 2., 2., -1.]]),
1078 array([[1. , 2.5, 4. ],
1079 [1. , 1. , 1. ]]))
1081 In this example the spacing is also specified:
1082 uniform for axis=0 and non uniform for axis=1
1084 >>> dx = 2.
1085 >>> y = [1., 1.5, 3.5]
1086 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]]), dx, y)
1087 (array([[ 1. , 1. , -0.5],
1088 [ 1. , 1. , -0.5]]),
1089 array([[2. , 2. , 2. ],
1090 [2. , 1.7, 0.5]]))
1092 It is possible to specify how boundaries are treated using `edge_order`
1094 >>> x = np.array([0, 1, 2, 3, 4])
1095 >>> f = x**2
1096 >>> np.gradient(f, edge_order=1)
1097 array([1., 2., 4., 6., 7.])
1098 >>> np.gradient(f, edge_order=2)
1099 array([0., 2., 4., 6., 8.])
1101 The `axis` keyword can be used to specify a subset of axes of which the
1102 gradient is calculated
1104 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]]), axis=0)
1105 array([[ 2., 2., -1.],
1106 [ 2., 2., -1.]])
1108 The `varargs` argument defines the spacing between sample points in the
1109 input array. It can take two forms:
1111 1. An array, specifying coordinates, which may be unevenly spaced:
1113 >>> x = np.array([0., 2., 3., 6., 8.])
1114 >>> y = x ** 2
1115 >>> np.gradient(y, x, edge_order=2)
1116 array([ 0., 4., 6., 12., 16.])
1118 2. A scalar, representing the fixed sample distance:
1120 >>> dx = 2
1121 >>> x = np.array([0., 2., 4., 6., 8.])
1122 >>> y = x ** 2
1123 >>> np.gradient(y, dx, edge_order=2)
1124 array([ 0., 4., 8., 12., 16.])
1126 It's possible to provide different data for spacing along each dimension.
1127 The number of arguments must match the number of dimensions in the input
1128 data.
1130 >>> dx = 2
1131 >>> dy = 3
1132 >>> x = np.arange(0, 6, dx)
1133 >>> y = np.arange(0, 9, dy)
1134 >>> xs, ys = np.meshgrid(x, y)
1135 >>> zs = xs + 2 * ys
1136 >>> np.gradient(zs, dy, dx) # Passing two scalars
1137 (array([[2., 2., 2.],
1138 [2., 2., 2.],
1139 [2., 2., 2.]]),
1140 array([[1., 1., 1.],
1141 [1., 1., 1.],
1142 [1., 1., 1.]]))
1144 Mixing scalars and arrays is also allowed:
1146 >>> np.gradient(zs, y, dx) # Passing one array and one scalar
1147 (array([[2., 2., 2.],
1148 [2., 2., 2.],
1149 [2., 2., 2.]]),
1150 array([[1., 1., 1.],
1151 [1., 1., 1.],
1152 [1., 1., 1.]]))
1154 Notes
1155 -----
1156 Assuming that :math:`f\\in C^{3}` (i.e., :math:`f` has at least 3 continuous
1157 derivatives) and let :math:`h_{*}` be a non-homogeneous stepsize, we
1158 minimize the "consistency error" :math:`\\eta_{i}` between the true gradient
1159 and its estimate from a linear combination of the neighboring grid-points:
1161 .. math::
1163 \\eta_{i} = f_{i}^{\\left(1\\right)} -
1164 \\left[ \\alpha f\\left(x_{i}\\right) +
1165 \\beta f\\left(x_{i} + h_{d}\\right) +
1166 \\gamma f\\left(x_{i}-h_{s}\\right)
1167 \\right]
1169 By substituting :math:`f(x_{i} + h_{d})` and :math:`f(x_{i} - h_{s})`
1170 with their Taylor series expansion, this translates into solving
1171 the following the linear system:
1173 .. math::
1175 \\left\\{
1176 \\begin{array}{r}
1177 \\alpha+\\beta+\\gamma=0 \\\\
1178 \\beta h_{d}-\\gamma h_{s}=1 \\\\
1179 \\beta h_{d}^{2}+\\gamma h_{s}^{2}=0
1180 \\end{array}
1181 \\right.
1183 The resulting approximation of :math:`f_{i}^{(1)}` is the following:
1185 .. math::
1187 \\hat f_{i}^{(1)} =
1188 \\frac{
1189 h_{s}^{2}f\\left(x_{i} + h_{d}\\right)
1190 + \\left(h_{d}^{2} - h_{s}^{2}\\right)f\\left(x_{i}\\right)
1191 - h_{d}^{2}f\\left(x_{i}-h_{s}\\right)}
1192 { h_{s}h_{d}\\left(h_{d} + h_{s}\\right)}
1193 + \\mathcal{O}\\left(\\frac{h_{d}h_{s}^{2}
1194 + h_{s}h_{d}^{2}}{h_{d}
1195 + h_{s}}\\right)
1197 It is worth noting that if :math:`h_{s}=h_{d}`
1198 (i.e., data are evenly spaced)
1199 we find the standard second order approximation:
1201 .. math::
1203 \\hat f_{i}^{(1)}=
1204 \\frac{f\\left(x_{i+1}\\right) - f\\left(x_{i-1}\\right)}{2h}
1205 + \\mathcal{O}\\left(h^{2}\\right)
1207 With a similar procedure the forward/backward approximations used for
1208 boundaries can be derived.
1210 References
1211 ----------
1212 .. [1] Quarteroni A., Sacco R., Saleri F. (2007) Numerical Mathematics
1213 (Texts in Applied Mathematics). New York: Springer.
1214 .. [2] Durran D. R. (1999) Numerical Methods for Wave Equations
1215 in Geophysical Fluid Dynamics. New York: Springer.
1216 .. [3] Fornberg B. (1988) Generation of Finite Difference Formulas on
1217 Arbitrarily Spaced Grids,
1218 Mathematics of Computation 51, no. 184 : 699-706.
1219 `PDF <https://www.ams.org/journals/mcom/1988-51-184/
1220 S0025-5718-1988-0935077-0/S0025-5718-1988-0935077-0.pdf>`_.
1221 """
1222 f = np.asanyarray(f)
1223 N = f.ndim # number of dimensions
1225 if axis is None:
1226 axes = tuple(range(N))
1227 else:
1228 axes = _nx.normalize_axis_tuple(axis, N)
1230 len_axes = len(axes)
1231 n = len(varargs)
1232 if n == 0:
1233 # no spacing argument - use 1 in all axes
1234 dx = [1.0] * len_axes
1235 elif n == 1 and np.ndim(varargs[0]) == 0:
1236 # single scalar for all axes
1237 dx = varargs * len_axes
1238 elif n == len_axes:
1239 # scalar or 1d array for each axis
1240 dx = list(varargs)
1241 for i, distances in enumerate(dx):
1242 distances = np.asanyarray(distances)
1243 if distances.ndim == 0:
1244 continue
1245 elif distances.ndim != 1:
1246 raise ValueError("distances must be either scalars or 1d")
1247 if len(distances) != f.shape[axes[i]]:
1248 raise ValueError("when 1d, distances must match "
1249 "the length of the corresponding dimension")
1250 if np.issubdtype(distances.dtype, np.integer):
1251 # Convert numpy integer types to float64 to avoid modular
1252 # arithmetic in np.diff(distances).
1253 distances = distances.astype(np.float64)
1254 diffx = np.diff(distances)
1255 # if distances are constant reduce to the scalar case
1256 # since it brings a consistent speedup
1257 if (diffx == diffx[0]).all():
1258 diffx = diffx[0]
1259 dx[i] = diffx
1260 else:
1261 raise TypeError("invalid number of arguments")
1263 if edge_order > 2:
1264 raise ValueError("'edge_order' greater than 2 not supported")
1266 # use central differences on interior and one-sided differences on the
1267 # endpoints. This preserves second order-accuracy over the full domain.
1269 outvals = []
1271 # create slice objects --- initially all are [:, :, ..., :]
1272 slice1 = [slice(None)] * N
1273 slice2 = [slice(None)] * N
1274 slice3 = [slice(None)] * N
1275 slice4 = [slice(None)] * N
1277 otype = f.dtype
1278 if otype.type is np.datetime64:
1279 # the timedelta dtype with the same unit information
1280 otype = np.dtype(otype.name.replace('datetime', 'timedelta'))
1281 # view as timedelta to allow addition
1282 f = f.view(otype)
1283 elif otype.type is np.timedelta64:
1284 pass
1285 elif np.issubdtype(otype, np.inexact):
1286 pass
1287 else:
1288 # All other types convert to floating point.
1289 # First check if f is a numpy integer type; if so, convert f to float64
1290 # to avoid modular arithmetic when computing the changes in f.
1291 if np.issubdtype(otype, np.integer):
1292 f = f.astype(np.float64)
1293 otype = np.float64
1295 for axis, ax_dx in zip(axes, dx):
1296 if f.shape[axis] < edge_order + 1:
1297 raise ValueError(
1298 "Shape of array too small to calculate a numerical gradient, "
1299 "at least (edge_order + 1) elements are required.")
1300 # result allocation
1301 out = np.empty_like(f, dtype=otype)
1303 # spacing for the current axis
1304 uniform_spacing = np.ndim(ax_dx) == 0
1306 # Numerical differentiation: 2nd order interior
1307 slice1[axis] = slice(1, -1)
1308 slice2[axis] = slice(None, -2)
1309 slice3[axis] = slice(1, -1)
1310 slice4[axis] = slice(2, None)
1312 if uniform_spacing:
1313 out[tuple(slice1)] = (f[tuple(slice4)] - f[tuple(slice2)]) / (2. * ax_dx)
1314 else:
1315 dx1 = ax_dx[0:-1]
1316 dx2 = ax_dx[1:]
1317 a = -(dx2) / (dx1 * (dx1 + dx2))
1318 b = (dx2 - dx1) / (dx1 * dx2)
1319 c = dx1 / (dx2 * (dx1 + dx2))
1320 # fix the shape for broadcasting
1321 shape = np.ones(N, dtype=int)
1322 shape[axis] = -1
1324 a = a.reshape(shape)
1325 b = b.reshape(shape)
1326 c = c.reshape(shape)
1327 # 1D equivalent -- out[1:-1] = a * f[:-2] + b * f[1:-1] + c * f[2:]
1328 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] \
1329 + c * f[tuple(slice4)]
1331 # Numerical differentiation: 1st order edges
1332 if edge_order == 1:
1333 slice1[axis] = 0
1334 slice2[axis] = 1
1335 slice3[axis] = 0
1336 dx_0 = ax_dx if uniform_spacing else ax_dx[0]
1337 # 1D equivalent -- out[0] = (f[1] - f[0]) / (x[1] - x[0])
1338 out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0
1340 slice1[axis] = -1
1341 slice2[axis] = -1
1342 slice3[axis] = -2
1343 dx_n = ax_dx if uniform_spacing else ax_dx[-1]
1344 # 1D equivalent -- out[-1] = (f[-1] - f[-2]) / (x[-1] - x[-2])
1345 out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_n
1347 # Numerical differentiation: 2nd order edges
1348 else:
1349 slice1[axis] = 0
1350 slice2[axis] = 0
1351 slice3[axis] = 1
1352 slice4[axis] = 2
1353 if uniform_spacing:
1354 a = -1.5 / ax_dx
1355 b = 2. / ax_dx
1356 c = -0.5 / ax_dx
1357 else:
1358 dx1 = ax_dx[0]
1359 dx2 = ax_dx[1]
1360 a = -(2. * dx1 + dx2) / (dx1 * (dx1 + dx2))
1361 b = (dx1 + dx2) / (dx1 * dx2)
1362 c = - dx1 / (dx2 * (dx1 + dx2))
1363 # 1D equivalent -- out[0] = a * f[0] + b * f[1] + c * f[2]
1364 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] \
1365 + c * f[tuple(slice4)]
1367 slice1[axis] = -1
1368 slice2[axis] = -3
1369 slice3[axis] = -2
1370 slice4[axis] = -1
1371 if uniform_spacing:
1372 a = 0.5 / ax_dx
1373 b = -2. / ax_dx
1374 c = 1.5 / ax_dx
1375 else:
1376 dx1 = ax_dx[-2]
1377 dx2 = ax_dx[-1]
1378 a = (dx2) / (dx1 * (dx1 + dx2))
1379 b = - (dx2 + dx1) / (dx1 * dx2)
1380 c = (2. * dx2 + dx1) / (dx2 * (dx1 + dx2))
1381 # 1D equivalent -- out[-1] = a * f[-3] + b * f[-2] + c * f[-1]
1382 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] \
1383 + c * f[tuple(slice4)]
1385 outvals.append(out)
1387 # reset the slice object in this dimension to ":"
1388 slice1[axis] = slice(None)
1389 slice2[axis] = slice(None)
1390 slice3[axis] = slice(None)
1391 slice4[axis] = slice(None)
1393 if len_axes == 1:
1394 return outvals[0]
1395 return tuple(outvals)
1398def _diff_dispatcher(a, n=None, axis=None, prepend=None, append=None):
1399 return (a, prepend, append)
1402@array_function_dispatch(_diff_dispatcher)
1403def diff(a, n=1, axis=-1, prepend=np._NoValue, append=np._NoValue):
1404 """
1405 Calculate the n-th discrete difference along the given axis.
1407 The first difference is given by ``out[i] = a[i+1] - a[i]`` along
1408 the given axis, higher differences are calculated by using `diff`
1409 recursively.
1411 Parameters
1412 ----------
1413 a : array_like
1414 Input array
1415 n : int, optional
1416 The number of times values are differenced. If zero, the input
1417 is returned as-is.
1418 axis : int, optional
1419 The axis along which the difference is taken, default is the
1420 last axis.
1421 prepend, append : array_like, optional
1422 Values to prepend or append to `a` along axis prior to
1423 performing the difference. Scalar values are expanded to
1424 arrays with length 1 in the direction of axis and the shape
1425 of the input array in along all other axes. Otherwise the
1426 dimension and shape must match `a` except along axis.
1428 Returns
1429 -------
1430 diff : ndarray
1431 The n-th differences. The shape of the output is the same as `a`
1432 except along `axis` where the dimension is smaller by `n`. The
1433 type of the output is the same as the type of the difference
1434 between any two elements of `a`. This is the same as the type of
1435 `a` in most cases. A notable exception is `datetime64`, which
1436 results in a `timedelta64` output array.
1438 See Also
1439 --------
1440 gradient, ediff1d, cumsum
1442 Notes
1443 -----
1444 Type is preserved for boolean arrays, so the result will contain
1445 `False` when consecutive elements are the same and `True` when they
1446 differ.
1448 For unsigned integer arrays, the results will also be unsigned. This
1449 should not be surprising, as the result is consistent with
1450 calculating the difference directly:
1452 >>> u8_arr = np.array([1, 0], dtype=np.uint8)
1453 >>> np.diff(u8_arr)
1454 array([255], dtype=uint8)
1455 >>> u8_arr[1,...] - u8_arr[0,...]
1456 np.uint8(255)
1458 If this is not desirable, then the array should be cast to a larger
1459 integer type first:
1461 >>> i16_arr = u8_arr.astype(np.int16)
1462 >>> np.diff(i16_arr)
1463 array([-1], dtype=int16)
1465 Examples
1466 --------
1467 >>> import numpy as np
1468 >>> x = np.array([1, 2, 4, 7, 0])
1469 >>> np.diff(x)
1470 array([ 1, 2, 3, -7])
1471 >>> np.diff(x, n=2)
1472 array([ 1, 1, -10])
1474 >>> x = np.array([[1, 3, 6, 10], [0, 5, 6, 8]])
1475 >>> np.diff(x)
1476 array([[2, 3, 4],
1477 [5, 1, 2]])
1478 >>> np.diff(x, axis=0)
1479 array([[-1, 2, 0, -2]])
1481 >>> x = np.arange('1066-10-13', '1066-10-16', dtype=np.datetime64)
1482 >>> np.diff(x)
1483 array([1, 1], dtype='timedelta64[D]')
1485 """
1486 if n == 0:
1487 return a
1488 if n < 0:
1489 raise ValueError(
1490 "order must be non-negative but got " + repr(n))
1492 a = asanyarray(a)
1493 nd = a.ndim
1494 if nd == 0:
1495 raise ValueError("diff requires input that is at least one dimensional")
1496 axis = normalize_axis_index(axis, nd)
1498 combined = []
1499 if prepend is not np._NoValue:
1500 prepend = np.asanyarray(prepend)
1501 if prepend.ndim == 0:
1502 shape = list(a.shape)
1503 shape[axis] = 1
1504 prepend = np.broadcast_to(prepend, tuple(shape))
1505 combined.append(prepend)
1507 combined.append(a)
1509 if append is not np._NoValue:
1510 append = np.asanyarray(append)
1511 if append.ndim == 0:
1512 shape = list(a.shape)
1513 shape[axis] = 1
1514 append = np.broadcast_to(append, tuple(shape))
1515 combined.append(append)
1517 if len(combined) > 1:
1518 a = np.concatenate(combined, axis)
1520 slice1 = [slice(None)] * nd
1521 slice2 = [slice(None)] * nd
1522 slice1[axis] = slice(1, None)
1523 slice2[axis] = slice(None, -1)
1524 slice1 = tuple(slice1)
1525 slice2 = tuple(slice2)
1527 op = not_equal if a.dtype == np.bool else subtract
1528 for _ in range(n):
1529 a = op(a[slice1], a[slice2])
1531 return a
1534def _interp_dispatcher(x, xp, fp, left=None, right=None, period=None):
1535 return (x, xp, fp)
1538@array_function_dispatch(_interp_dispatcher)
1539def interp(x, xp, fp, left=None, right=None, period=None):
1540 """
1541 One-dimensional linear interpolation for monotonically increasing sample points.
1543 Returns the one-dimensional piecewise linear interpolant to a function
1544 with given discrete data points (`xp`, `fp`), evaluated at `x`.
1546 Parameters
1547 ----------
1548 x : array_like
1549 The x-coordinates at which to evaluate the interpolated values.
1551 xp : 1-D sequence of floats
1552 The x-coordinates of the data points, must be increasing if argument
1553 `period` is not specified. Otherwise, `xp` is internally sorted after
1554 normalizing the periodic boundaries with ``xp = xp % period``.
1556 fp : 1-D sequence of float or complex
1557 The y-coordinates of the data points, same length as `xp`.
1559 left : optional float or complex corresponding to fp
1560 Value to return for `x < xp[0]`, default is `fp[0]`.
1562 right : optional float or complex corresponding to fp
1563 Value to return for `x > xp[-1]`, default is `fp[-1]`.
1565 period : None or float, optional
1566 A period for the x-coordinates. This parameter allows the proper
1567 interpolation of angular x-coordinates. Parameters `left` and `right`
1568 are ignored if `period` is specified.
1570 Returns
1571 -------
1572 y : float or complex (corresponding to fp) or ndarray
1573 The interpolated values, same shape as `x`.
1575 Raises
1576 ------
1577 ValueError
1578 If `xp` and `fp` have different length
1579 If `xp` or `fp` are not 1-D sequences
1580 If `period == 0`
1582 See Also
1583 --------
1584 scipy.interpolate
1586 Warnings
1587 --------
1588 The x-coordinate sequence is expected to be increasing, but this is not
1589 explicitly enforced. However, if the sequence `xp` is non-increasing,
1590 interpolation results are meaningless.
1592 Note that, since NaN is unsortable, `xp` also cannot contain NaNs.
1594 A simple check for `xp` being strictly increasing is::
1596 np.all(np.diff(xp) > 0)
1598 Examples
1599 --------
1600 >>> import numpy as np
1601 >>> xp = [1, 2, 3]
1602 >>> fp = [3, 2, 0]
1603 >>> np.interp(2.5, xp, fp)
1604 1.0
1605 >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
1606 array([3. , 3. , 2.5 , 0.56, 0. ])
1607 >>> UNDEF = -99.0
1608 >>> np.interp(3.14, xp, fp, right=UNDEF)
1609 -99.0
1611 Plot an interpolant to the sine function:
1613 >>> x = np.linspace(0, 2*np.pi, 10)
1614 >>> y = np.sin(x)
1615 >>> xvals = np.linspace(0, 2*np.pi, 50)
1616 >>> yinterp = np.interp(xvals, x, y)
1617 >>> import matplotlib.pyplot as plt
1618 >>> plt.plot(x, y, 'o')
1619 [<matplotlib.lines.Line2D object at 0x...>]
1620 >>> plt.plot(xvals, yinterp, '-x')
1621 [<matplotlib.lines.Line2D object at 0x...>]
1622 >>> plt.show()
1624 Interpolation with periodic x-coordinates:
1626 >>> x = [-180, -170, -185, 185, -10, -5, 0, 365]
1627 >>> xp = [190, -190, 350, -350]
1628 >>> fp = [5, 10, 3, 4]
1629 >>> np.interp(x, xp, fp, period=360)
1630 array([7.5 , 5. , 8.75, 6.25, 3. , 3.25, 3.5 , 3.75])
1632 Complex interpolation:
1634 >>> x = [1.5, 4.0]
1635 >>> xp = [2,3,5]
1636 >>> fp = [1.0j, 0, 2+3j]
1637 >>> np.interp(x, xp, fp)
1638 array([0.+1.j , 1.+1.5j])
1640 """
1642 fp = np.asarray(fp)
1644 if np.iscomplexobj(fp):
1645 interp_func = compiled_interp_complex
1646 input_dtype = np.complex128
1647 else:
1648 interp_func = compiled_interp
1649 input_dtype = np.float64
1651 if period is not None:
1652 if period == 0:
1653 raise ValueError("period must be a non-zero value")
1654 period = abs(period)
1655 left = None
1656 right = None
1658 x = np.asarray(x, dtype=np.float64)
1659 xp = np.asarray(xp, dtype=np.float64)
1660 fp = np.asarray(fp, dtype=input_dtype)
1662 if xp.ndim != 1 or fp.ndim != 1:
1663 raise ValueError("Data points must be 1-D sequences")
1664 if xp.shape[0] != fp.shape[0]:
1665 raise ValueError("fp and xp are not of the same length")
1666 # normalizing periodic boundaries
1667 x = x % period
1668 xp = xp % period
1669 asort_xp = np.argsort(xp)
1670 xp = xp[asort_xp]
1671 fp = fp[asort_xp]
1672 xp = np.concatenate((xp[-1:] - period, xp, xp[0:1] + period))
1673 fp = np.concatenate((fp[-1:], fp, fp[0:1]))
1675 return interp_func(x, xp, fp, left, right)
1678def _angle_dispatcher(z, deg=None):
1679 return (z,)
1682@array_function_dispatch(_angle_dispatcher)
1683def angle(z, deg=False):
1684 """
1685 Return the angle of the complex argument.
1687 Parameters
1688 ----------
1689 z : array_like
1690 A complex number or sequence of complex numbers.
1691 deg : bool, optional
1692 Return angle in degrees if True, radians if False (default).
1694 Returns
1695 -------
1696 angle : ndarray or scalar
1697 The counterclockwise angle from the positive real axis on the complex
1698 plane in the range ``(-pi, pi]``, with dtype as numpy.float64.
1700 See Also
1701 --------
1702 arctan2
1703 absolute
1705 Notes
1706 -----
1707 This function passes the imaginary and real parts of the argument to
1708 `arctan2` to compute the result; consequently, it follows the convention
1709 of `arctan2` when the magnitude of the argument is zero. See example.
1711 Examples
1712 --------
1713 >>> import numpy as np
1714 >>> np.angle([1.0, 1.0j, 1+1j]) # in radians
1715 array([ 0. , 1.57079633, 0.78539816]) # may vary
1716 >>> np.angle(1+1j, deg=True) # in degrees
1717 45.0
1718 >>> np.angle([0., -0., complex(0., -0.), complex(-0., -0.)]) # convention
1719 array([ 0. , 3.14159265, -0. , -3.14159265])
1721 """
1722 z = asanyarray(z)
1723 if issubclass(z.dtype.type, _nx.complexfloating):
1724 zimag = z.imag
1725 zreal = z.real
1726 else:
1727 zimag = 0
1728 zreal = z
1730 a = arctan2(zimag, zreal)
1731 if deg:
1732 a *= 180 / pi
1733 return a
1736def _unwrap_dispatcher(p, discont=None, axis=None, *, period=None):
1737 return (p,)
1740@array_function_dispatch(_unwrap_dispatcher)
1741def unwrap(p, discont=None, axis=-1, *, period=2 * pi):
1742 r"""
1743 Unwrap by taking the complement of large deltas with respect to the period.
1745 This unwraps a signal `p` by changing elements which have an absolute
1746 difference from their predecessor of more than ``max(discont, period/2)``
1747 to their `period`-complementary values.
1749 For the default case where `period` is :math:`2\pi` and `discont` is
1750 :math:`\pi`, this unwraps a radian phase `p` such that adjacent differences
1751 are never greater than :math:`\pi` by adding :math:`2k\pi` for some
1752 integer :math:`k`.
1754 Parameters
1755 ----------
1756 p : array_like
1757 Input array.
1758 discont : float, optional
1759 Maximum discontinuity between values, default is ``period/2``.
1760 Values below ``period/2`` are treated as if they were ``period/2``.
1761 To have an effect different from the default, `discont` should be
1762 larger than ``period/2``.
1763 axis : int, optional
1764 Axis along which unwrap will operate, default is the last axis.
1765 period : float, optional
1766 Size of the range over which the input wraps. By default, it is
1767 ``2 pi``.
1769 .. versionadded:: 1.21.0
1771 Returns
1772 -------
1773 out : ndarray
1774 Output array.
1776 See Also
1777 --------
1778 rad2deg, deg2rad
1780 Notes
1781 -----
1782 If the discontinuity in `p` is smaller than ``period/2``,
1783 but larger than `discont`, no unwrapping is done because taking
1784 the complement would only make the discontinuity larger.
1786 Examples
1787 --------
1788 >>> import numpy as np
1790 >>> phase = np.linspace(0, np.pi, num=5)
1791 >>> phase[3:] += np.pi
1792 >>> phase
1793 array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531]) # may vary
1794 >>> np.unwrap(phase)
1795 array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ]) # may vary
1796 >>> np.unwrap([0, 1, 2, -1, 0], period=4)
1797 array([0, 1, 2, 3, 4])
1798 >>> np.unwrap([ 1, 2, 3, 4, 5, 6, 1, 2, 3], period=6)
1799 array([1, 2, 3, 4, 5, 6, 7, 8, 9])
1800 >>> np.unwrap([2, 3, 4, 5, 2, 3, 4, 5], period=4)
1801 array([2, 3, 4, 5, 6, 7, 8, 9])
1802 >>> phase_deg = np.mod(np.linspace(0 ,720, 19), 360) - 180
1803 >>> np.unwrap(phase_deg, period=360)
1804 array([-180., -140., -100., -60., -20., 20., 60., 100., 140.,
1805 180., 220., 260., 300., 340., 380., 420., 460., 500.,
1806 540.])
1808 This example plots the unwrapping of the wrapped input signal `w`.
1809 First generate `w`, then apply `unwrap` to get `u`.
1811 >>> t = np.linspace(0, 25, 801)
1812 >>> w = np.mod(1.5 * np.sin(1.1 * t + 0.26) * (1 - t / 6 + (t / 23) ** 3), 2.0) - 1
1813 >>> u = np.unwrap(w, period=2.0)
1815 Plot `w` and `u`.
1817 >>> import matplotlib.pyplot as plt
1818 >>> plt.plot(t, w, label='w (a signal wrapped to [-1, 1])')
1819 >>> plt.plot(t, u, linewidth=2.5, alpha=0.5, label='unwrap(w, period=2)')
1820 >>> plt.xlabel('t')
1821 >>> plt.grid(alpha=0.6)
1822 >>> plt.legend(framealpha=1, shadow=True)
1823 >>> plt.show()
1824 """
1825 p = asarray(p)
1826 nd = p.ndim
1827 dd = diff(p, axis=axis)
1828 if discont is None:
1829 discont = period / 2
1830 slice1 = [slice(None, None)] * nd # full slices
1831 slice1[axis] = slice(1, None)
1832 slice1 = tuple(slice1)
1833 dtype = np.result_type(dd, period)
1834 if _nx.issubdtype(dtype, _nx.integer):
1835 interval_high, rem = divmod(period, 2)
1836 boundary_ambiguous = rem == 0
1837 else:
1838 interval_high = period / 2
1839 boundary_ambiguous = True
1840 interval_low = -interval_high
1841 ddmod = mod(dd - interval_low, period) + interval_low
1842 if boundary_ambiguous:
1843 # for `mask = (abs(dd) == period/2)`, the above line made
1844 # `ddmod[mask] == -period/2`. correct these such that
1845 # `ddmod[mask] == sign(dd[mask])*period/2`.
1846 _nx.copyto(ddmod, interval_high,
1847 where=(ddmod == interval_low) & (dd > 0))
1848 ph_correct = ddmod - dd
1849 _nx.copyto(ph_correct, 0, where=abs(dd) < discont)
1850 up = array(p, copy=True, dtype=dtype)
1851 up[slice1] = p[slice1] + ph_correct.cumsum(axis)
1852 return up
1855def _sort_complex(a):
1856 return (a,)
1859@array_function_dispatch(_sort_complex)
1860def sort_complex(a):
1861 """
1862 Sort a complex array using the real part first, then the imaginary part.
1864 Parameters
1865 ----------
1866 a : array_like
1867 Input array
1869 Returns
1870 -------
1871 out : complex ndarray
1872 Always returns a sorted complex array.
1874 Examples
1875 --------
1876 >>> import numpy as np
1877 >>> np.sort_complex([5, 3, 6, 2, 1])
1878 array([1.+0.j, 2.+0.j, 3.+0.j, 5.+0.j, 6.+0.j])
1880 >>> np.sort_complex([1 + 2j, 2 - 1j, 3 - 2j, 3 - 3j, 3 + 5j])
1881 array([1.+2.j, 2.-1.j, 3.-3.j, 3.-2.j, 3.+5.j])
1883 """
1884 b = array(a, copy=True)
1885 b.sort()
1886 if not issubclass(b.dtype.type, _nx.complexfloating):
1887 if b.dtype.char in 'bhBH':
1888 return b.astype('F')
1889 elif b.dtype.char == 'g':
1890 return b.astype('G')
1891 else:
1892 return b.astype('D')
1893 else:
1894 return b
1897def _arg_trim_zeros(filt):
1898 """Return indices of the first and last non-zero element.
1900 Parameters
1901 ----------
1902 filt : array_like
1903 Input array.
1905 Returns
1906 -------
1907 start, stop : ndarray
1908 Two arrays containing the indices of the first and last non-zero
1909 element in each dimension.
1911 See also
1912 --------
1913 trim_zeros
1915 Examples
1916 --------
1917 >>> import numpy as np
1918 >>> _arg_trim_zeros(np.array([0, 0, 1, 1, 0]))
1919 (array([2]), array([3]))
1920 """
1921 nonzero = (
1922 np.argwhere(filt)
1923 if filt.dtype != np.object_
1924 # Historically, `trim_zeros` treats `None` in an object array
1925 # as non-zero while argwhere doesn't, account for that
1926 else np.argwhere(filt != 0)
1927 )
1928 if nonzero.size == 0:
1929 start = stop = np.array([], dtype=np.intp)
1930 else:
1931 start = nonzero.min(axis=0)
1932 stop = nonzero.max(axis=0)
1933 return start, stop
1936def _trim_zeros(filt, trim=None, axis=None):
1937 return (filt,)
1940@array_function_dispatch(_trim_zeros)
1941def trim_zeros(filt, trim='fb', axis=None):
1942 """Remove values along a dimension which are zero along all other.
1944 Parameters
1945 ----------
1946 filt : array_like
1947 Input array.
1948 trim : {"fb", "f", "b"}, optional
1949 A string with 'f' representing trim from front and 'b' to trim from
1950 back. By default, zeros are trimmed on both sides.
1951 Front and back refer to the edges of a dimension, with "front" referring
1952 to the side with the lowest index 0, and "back" referring to the highest
1953 index (or index -1).
1954 axis : int or sequence, optional
1955 If None, `filt` is cropped such that the smallest bounding box is
1956 returned that still contains all values which are not zero.
1957 If an axis is specified, `filt` will be sliced in that dimension only
1958 on the sides specified by `trim`. The remaining area will be the
1959 smallest that still contains all values wich are not zero.
1961 .. versionadded:: 2.2.0
1963 Returns
1964 -------
1965 trimmed : ndarray or sequence
1966 The result of trimming the input. The number of dimensions and the
1967 input data type are preserved.
1969 Notes
1970 -----
1971 For all-zero arrays, the first axis is trimmed first.
1973 Examples
1974 --------
1975 >>> import numpy as np
1976 >>> a = np.array((0, 0, 0, 1, 2, 3, 0, 2, 1, 0))
1977 >>> np.trim_zeros(a)
1978 array([1, 2, 3, 0, 2, 1])
1980 >>> np.trim_zeros(a, trim='b')
1981 array([0, 0, 0, ..., 0, 2, 1])
1983 Multiple dimensions are supported.
1985 >>> b = np.array([[0, 0, 2, 3, 0, 0],
1986 ... [0, 1, 0, 3, 0, 0],
1987 ... [0, 0, 0, 0, 0, 0]])
1988 >>> np.trim_zeros(b)
1989 array([[0, 2, 3],
1990 [1, 0, 3]])
1992 >>> np.trim_zeros(b, axis=-1)
1993 array([[0, 2, 3],
1994 [1, 0, 3],
1995 [0, 0, 0]])
1997 The input data type is preserved, list/tuple in means list/tuple out.
1999 >>> np.trim_zeros([0, 1, 2, 0])
2000 [1, 2]
2002 """
2003 filt_ = np.asarray(filt)
2005 trim = trim.lower()
2006 if trim not in {"fb", "bf", "f", "b"}:
2007 raise ValueError(f"unexpected character(s) in `trim`: {trim!r}")
2008 if axis is None:
2009 axis_tuple = tuple(range(filt_.ndim))
2010 else:
2011 axis_tuple = _nx.normalize_axis_tuple(axis, filt_.ndim, argname="axis")
2013 if not axis_tuple:
2014 # No trimming requested -> return input unmodified.
2015 return filt
2017 start, stop = _arg_trim_zeros(filt_)
2018 stop += 1 # Adjust for slicing
2020 if start.size == 0:
2021 # filt is all-zero -> assign same values to start and stop so that
2022 # resulting slice will be empty
2023 start = stop = np.zeros(filt_.ndim, dtype=np.intp)
2024 else:
2025 if 'f' not in trim:
2026 start = (None,) * filt_.ndim
2027 if 'b' not in trim:
2028 stop = (None,) * filt_.ndim
2030 sl = tuple(slice(start[ax], stop[ax]) if ax in axis_tuple else slice(None)
2031 for ax in range(filt_.ndim))
2032 if len(sl) == 1:
2033 # filt is 1D -> avoid multi-dimensional slicing to preserve
2034 # non-array input types
2035 return filt[sl[0]]
2036 return filt[sl]
2039def _extract_dispatcher(condition, arr):
2040 return (condition, arr)
2043@array_function_dispatch(_extract_dispatcher)
2044def extract(condition, arr):
2045 """
2046 Return the elements of an array that satisfy some condition.
2048 This is equivalent to ``np.compress(ravel(condition), ravel(arr))``. If
2049 `condition` is boolean ``np.extract`` is equivalent to ``arr[condition]``.
2051 Note that `place` does the exact opposite of `extract`.
2053 Parameters
2054 ----------
2055 condition : array_like
2056 An array whose nonzero or True entries indicate the elements of `arr`
2057 to extract.
2058 arr : array_like
2059 Input array of the same size as `condition`.
2061 Returns
2062 -------
2063 extract : ndarray
2064 Rank 1 array of values from `arr` where `condition` is True.
2066 See Also
2067 --------
2068 take, put, copyto, compress, place
2070 Examples
2071 --------
2072 >>> import numpy as np
2073 >>> arr = np.arange(12).reshape((3, 4))
2074 >>> arr
2075 array([[ 0, 1, 2, 3],
2076 [ 4, 5, 6, 7],
2077 [ 8, 9, 10, 11]])
2078 >>> condition = np.mod(arr, 3)==0
2079 >>> condition
2080 array([[ True, False, False, True],
2081 [False, False, True, False],
2082 [False, True, False, False]])
2083 >>> np.extract(condition, arr)
2084 array([0, 3, 6, 9])
2087 If `condition` is boolean:
2089 >>> arr[condition]
2090 array([0, 3, 6, 9])
2092 """
2093 return _nx.take(ravel(arr), nonzero(ravel(condition))[0])
2096def _place_dispatcher(arr, mask, vals):
2097 return (arr, mask, vals)
2100@array_function_dispatch(_place_dispatcher)
2101def place(arr, mask, vals):
2102 """
2103 Change elements of an array based on conditional and input values.
2105 Similar to ``np.copyto(arr, vals, where=mask)``, the difference is that
2106 `place` uses the first N elements of `vals`, where N is the number of
2107 True values in `mask`, while `copyto` uses the elements where `mask`
2108 is True.
2110 Note that `extract` does the exact opposite of `place`.
2112 Parameters
2113 ----------
2114 arr : ndarray
2115 Array to put data into.
2116 mask : array_like
2117 Boolean mask array. Must have the same size as `a`.
2118 vals : 1-D sequence
2119 Values to put into `a`. Only the first N elements are used, where
2120 N is the number of True values in `mask`. If `vals` is smaller
2121 than N, it will be repeated, and if elements of `a` are to be masked,
2122 this sequence must be non-empty.
2124 See Also
2125 --------
2126 copyto, put, take, extract
2128 Examples
2129 --------
2130 >>> import numpy as np
2131 >>> arr = np.arange(6).reshape(2, 3)
2132 >>> np.place(arr, arr>2, [44, 55])
2133 >>> arr
2134 array([[ 0, 1, 2],
2135 [44, 55, 44]])
2137 """
2138 return _place(arr, mask, vals)
2141# See https://docs.scipy.org/doc/numpy/reference/c-api.generalized-ufuncs.html
2142_DIMENSION_NAME = r'\w+'
2143_CORE_DIMENSION_LIST = f'(?:{_DIMENSION_NAME}(?:,{_DIMENSION_NAME})*)?'
2144_ARGUMENT = fr'\({_CORE_DIMENSION_LIST}\)'
2145_ARGUMENT_LIST = f'{_ARGUMENT}(?:,{_ARGUMENT})*'
2146_SIGNATURE = f'^{_ARGUMENT_LIST}->{_ARGUMENT_LIST}$'
2149def _parse_gufunc_signature(signature):
2150 """
2151 Parse string signatures for a generalized universal function.
2153 Arguments
2154 ---------
2155 signature : string
2156 Generalized universal function signature, e.g., ``(m,n),(n,p)->(m,p)``
2157 for ``np.matmul``.
2159 Returns
2160 -------
2161 Tuple of input and output core dimensions parsed from the signature, each
2162 of the form List[Tuple[str, ...]].
2163 """
2164 signature = re.sub(r'\s+', '', signature)
2166 if not re.match(_SIGNATURE, signature):
2167 raise ValueError(
2168 f'not a valid gufunc signature: {signature}')
2169 return tuple([tuple(re.findall(_DIMENSION_NAME, arg))
2170 for arg in re.findall(_ARGUMENT, arg_list)]
2171 for arg_list in signature.split('->'))
2174def _update_dim_sizes(dim_sizes, arg, core_dims):
2175 """
2176 Incrementally check and update core dimension sizes for a single argument.
2178 Arguments
2179 ---------
2180 dim_sizes : Dict[str, int]
2181 Sizes of existing core dimensions. Will be updated in-place.
2182 arg : ndarray
2183 Argument to examine.
2184 core_dims : Tuple[str, ...]
2185 Core dimensions for this argument.
2186 """
2187 if not core_dims:
2188 return
2190 num_core_dims = len(core_dims)
2191 if arg.ndim < num_core_dims:
2192 raise ValueError(
2193 '%d-dimensional argument does not have enough '
2194 'dimensions for all core dimensions %r'
2195 % (arg.ndim, core_dims))
2197 core_shape = arg.shape[-num_core_dims:]
2198 for dim, size in zip(core_dims, core_shape):
2199 if dim in dim_sizes:
2200 if size != dim_sizes[dim]:
2201 raise ValueError(
2202 'inconsistent size for core dimension %r: %r vs %r'
2203 % (dim, size, dim_sizes[dim]))
2204 else:
2205 dim_sizes[dim] = size
2208def _parse_input_dimensions(args, input_core_dims):
2209 """
2210 Parse broadcast and core dimensions for vectorize with a signature.
2212 Arguments
2213 ---------
2214 args : Tuple[ndarray, ...]
2215 Tuple of input arguments to examine.
2216 input_core_dims : List[Tuple[str, ...]]
2217 List of core dimensions corresponding to each input.
2219 Returns
2220 -------
2221 broadcast_shape : Tuple[int, ...]
2222 Common shape to broadcast all non-core dimensions to.
2223 dim_sizes : Dict[str, int]
2224 Common sizes for named core dimensions.
2225 """
2226 broadcast_args = []
2227 dim_sizes = {}
2228 for arg, core_dims in zip(args, input_core_dims):
2229 _update_dim_sizes(dim_sizes, arg, core_dims)
2230 ndim = arg.ndim - len(core_dims)
2231 dummy_array = np.lib.stride_tricks.as_strided(0, arg.shape[:ndim])
2232 broadcast_args.append(dummy_array)
2233 broadcast_shape = np.lib._stride_tricks_impl._broadcast_shape(
2234 *broadcast_args
2235 )
2236 return broadcast_shape, dim_sizes
2239def _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims):
2240 """Helper for calculating broadcast shapes with core dimensions."""
2241 return [broadcast_shape + tuple(dim_sizes[dim] for dim in core_dims)
2242 for core_dims in list_of_core_dims]
2245def _create_arrays(broadcast_shape, dim_sizes, list_of_core_dims, dtypes,
2246 results=None):
2247 """Helper for creating output arrays in vectorize."""
2248 shapes = _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims)
2249 if dtypes is None:
2250 dtypes = [None] * len(shapes)
2251 if results is None:
2252 arrays = tuple(np.empty(shape=shape, dtype=dtype)
2253 for shape, dtype in zip(shapes, dtypes))
2254 else:
2255 arrays = tuple(np.empty_like(result, shape=shape, dtype=dtype)
2256 for result, shape, dtype
2257 in zip(results, shapes, dtypes))
2258 return arrays
2261def _get_vectorize_dtype(dtype):
2262 if dtype.char in "SU":
2263 return dtype.char
2264 return dtype
2267@set_module('numpy')
2268class vectorize:
2269 """
2270 vectorize(pyfunc=np._NoValue, otypes=None, doc=None, excluded=None,
2271 cache=False, signature=None)
2273 Returns an object that acts like pyfunc, but takes arrays as input.
2275 Define a vectorized function which takes a nested sequence of objects or
2276 numpy arrays as inputs and returns a single numpy array or a tuple of numpy
2277 arrays. The vectorized function evaluates `pyfunc` over successive tuples
2278 of the input arrays like the python map function, except it uses the
2279 broadcasting rules of numpy.
2281 The data type of the output of `vectorized` is determined by calling
2282 the function with the first element of the input. This can be avoided
2283 by specifying the `otypes` argument.
2285 Parameters
2286 ----------
2287 pyfunc : callable, optional
2288 A python function or method.
2289 Can be omitted to produce a decorator with keyword arguments.
2290 otypes : str or list of dtypes, optional
2291 The output data type. It must be specified as either a string of
2292 typecode characters or a list of data type specifiers. There should
2293 be one data type specifier for each output.
2294 doc : str, optional
2295 The docstring for the function. If None, the docstring will be the
2296 ``pyfunc.__doc__``.
2297 excluded : set, optional
2298 Set of strings or integers representing the positional or keyword
2299 arguments for which the function will not be vectorized. These will be
2300 passed directly to `pyfunc` unmodified.
2302 cache : bool, optional
2303 If neither `otypes` nor `signature` are provided, and `cache` is ``True``, then
2304 cache the number of outputs.
2306 signature : string, optional
2307 Generalized universal function signature, e.g., ``(m,n),(n)->(m)`` for
2308 vectorized matrix-vector multiplication. If provided, ``pyfunc`` will
2309 be called with (and expected to return) arrays with shapes given by the
2310 size of corresponding core dimensions. By default, ``pyfunc`` is
2311 assumed to take scalars as input and output.
2313 Returns
2314 -------
2315 out : callable
2316 A vectorized function if ``pyfunc`` was provided,
2317 a decorator otherwise.
2319 See Also
2320 --------
2321 frompyfunc : Takes an arbitrary Python function and returns a ufunc
2323 Notes
2324 -----
2325 The `vectorize` function is provided primarily for convenience, not for
2326 performance. The implementation is essentially a for loop.
2328 If neither `otypes` nor `signature` are specified, then a call to the function with
2329 the first argument will be used to determine the number of outputs. The results of
2330 this call will be cached if `cache` is `True` to prevent calling the function
2331 twice. However, to implement the cache, the original function must be wrapped
2332 which will slow down subsequent calls, so only do this if your function is
2333 expensive.
2335 The new keyword argument interface and `excluded` argument support
2336 further degrades performance.
2338 References
2339 ----------
2340 .. [1] :doc:`/reference/c-api/generalized-ufuncs`
2342 Examples
2343 --------
2344 >>> import numpy as np
2345 >>> def myfunc(a, b):
2346 ... "Return a-b if a>b, otherwise return a+b"
2347 ... if a > b:
2348 ... return a - b
2349 ... else:
2350 ... return a + b
2352 >>> vfunc = np.vectorize(myfunc)
2353 >>> vfunc([1, 2, 3, 4], 2)
2354 array([3, 4, 1, 2])
2356 The docstring is taken from the input function to `vectorize` unless it
2357 is specified:
2359 >>> vfunc.__doc__
2360 'Return a-b if a>b, otherwise return a+b'
2361 >>> vfunc = np.vectorize(myfunc, doc='Vectorized `myfunc`')
2362 >>> vfunc.__doc__
2363 'Vectorized `myfunc`'
2365 The output type is determined by evaluating the first element of the input,
2366 unless it is specified:
2368 >>> out = vfunc([1, 2, 3, 4], 2)
2369 >>> type(out[0])
2370 <class 'numpy.int64'>
2371 >>> vfunc = np.vectorize(myfunc, otypes=[float])
2372 >>> out = vfunc([1, 2, 3, 4], 2)
2373 >>> type(out[0])
2374 <class 'numpy.float64'>
2376 The `excluded` argument can be used to prevent vectorizing over certain
2377 arguments. This can be useful for array-like arguments of a fixed length
2378 such as the coefficients for a polynomial as in `polyval`:
2380 >>> def mypolyval(p, x):
2381 ... _p = list(p)
2382 ... res = _p.pop(0)
2383 ... while _p:
2384 ... res = res*x + _p.pop(0)
2385 ... return res
2387 Here, we exclude the zeroth argument from vectorization whether it is
2388 passed by position or keyword.
2390 >>> vpolyval = np.vectorize(mypolyval, excluded={0, 'p'})
2391 >>> vpolyval([1, 2, 3], x=[0, 1])
2392 array([3, 6])
2393 >>> vpolyval(p=[1, 2, 3], x=[0, 1])
2394 array([3, 6])
2396 The `signature` argument allows for vectorizing functions that act on
2397 non-scalar arrays of fixed length. For example, you can use it for a
2398 vectorized calculation of Pearson correlation coefficient and its p-value:
2400 >>> import scipy.stats
2401 >>> pearsonr = np.vectorize(scipy.stats.pearsonr,
2402 ... signature='(n),(n)->(),()')
2403 >>> pearsonr([[0, 1, 2, 3]], [[1, 2, 3, 4], [4, 3, 2, 1]])
2404 (array([ 1., -1.]), array([ 0., 0.]))
2406 Or for a vectorized convolution:
2408 >>> convolve = np.vectorize(np.convolve, signature='(n),(m)->(k)')
2409 >>> convolve(np.eye(4), [1, 2, 1])
2410 array([[1., 2., 1., 0., 0., 0.],
2411 [0., 1., 2., 1., 0., 0.],
2412 [0., 0., 1., 2., 1., 0.],
2413 [0., 0., 0., 1., 2., 1.]])
2415 Decorator syntax is supported. The decorator can be called as
2416 a function to provide keyword arguments:
2418 >>> @np.vectorize
2419 ... def identity(x):
2420 ... return x
2421 ...
2422 >>> identity([0, 1, 2])
2423 array([0, 1, 2])
2424 >>> @np.vectorize(otypes=[float])
2425 ... def as_float(x):
2426 ... return x
2427 ...
2428 >>> as_float([0, 1, 2])
2429 array([0., 1., 2.])
2430 """
2431 def __init__(self, pyfunc=np._NoValue, otypes=None, doc=None,
2432 excluded=None, cache=False, signature=None):
2434 if (pyfunc != np._NoValue) and (not callable(pyfunc)):
2435 # Splitting the error message to keep
2436 # the length below 79 characters.
2437 part1 = "When used as a decorator, "
2438 part2 = "only accepts keyword arguments."
2439 raise TypeError(part1 + part2)
2441 self.pyfunc = pyfunc
2442 self.cache = cache
2443 self.signature = signature
2444 if pyfunc != np._NoValue and hasattr(pyfunc, '__name__'):
2445 self.__name__ = pyfunc.__name__
2447 self._ufunc = {} # Caching to improve default performance
2448 self._doc = None
2449 self.__doc__ = doc
2450 if doc is None and hasattr(pyfunc, '__doc__'):
2451 self.__doc__ = pyfunc.__doc__
2452 else:
2453 self._doc = doc
2455 if isinstance(otypes, str):
2456 for char in otypes:
2457 if char not in typecodes['All']:
2458 raise ValueError(f"Invalid otype specified: {char}")
2459 elif iterable(otypes):
2460 otypes = [_get_vectorize_dtype(_nx.dtype(x)) for x in otypes]
2461 elif otypes is not None:
2462 raise ValueError("Invalid otype specification")
2463 self.otypes = otypes
2465 # Excluded variable support
2466 if excluded is None:
2467 excluded = set()
2468 self.excluded = set(excluded)
2470 if signature is not None:
2471 self._in_and_out_core_dims = _parse_gufunc_signature(signature)
2472 else:
2473 self._in_and_out_core_dims = None
2475 def _init_stage_2(self, pyfunc, *args, **kwargs):
2476 self.__name__ = pyfunc.__name__
2477 self.pyfunc = pyfunc
2478 if self._doc is None:
2479 self.__doc__ = pyfunc.__doc__
2480 else:
2481 self.__doc__ = self._doc
2483 def _call_as_normal(self, *args, **kwargs):
2484 """
2485 Return arrays with the results of `pyfunc` broadcast (vectorized) over
2486 `args` and `kwargs` not in `excluded`.
2487 """
2488 excluded = self.excluded
2489 if not kwargs and not excluded:
2490 func = self.pyfunc
2491 vargs = args
2492 else:
2493 # The wrapper accepts only positional arguments: we use `names` and
2494 # `inds` to mutate `the_args` and `kwargs` to pass to the original
2495 # function.
2496 nargs = len(args)
2498 names = [_n for _n in kwargs if _n not in excluded]
2499 inds = [_i for _i in range(nargs) if _i not in excluded]
2500 the_args = list(args)
2502 def func(*vargs):
2503 for _n, _i in enumerate(inds):
2504 the_args[_i] = vargs[_n]
2505 kwargs.update(zip(names, vargs[len(inds):]))
2506 return self.pyfunc(*the_args, **kwargs)
2508 vargs = [args[_i] for _i in inds]
2509 vargs.extend([kwargs[_n] for _n in names])
2511 return self._vectorize_call(func=func, args=vargs)
2513 def __call__(self, *args, **kwargs):
2514 if self.pyfunc is np._NoValue:
2515 self._init_stage_2(*args, **kwargs)
2516 return self
2518 return self._call_as_normal(*args, **kwargs)
2520 def _get_ufunc_and_otypes(self, func, args):
2521 """Return (ufunc, otypes)."""
2522 # frompyfunc will fail if args is empty
2523 if not args:
2524 raise ValueError('args can not be empty')
2526 if self.otypes is not None:
2527 otypes = self.otypes
2529 # self._ufunc is a dictionary whose keys are the number of
2530 # arguments (i.e. len(args)) and whose values are ufuncs created
2531 # by frompyfunc. len(args) can be different for different calls if
2532 # self.pyfunc has parameters with default values. We only use the
2533 # cache when func is self.pyfunc, which occurs when the call uses
2534 # only positional arguments and no arguments are excluded.
2536 nin = len(args)
2537 nout = len(self.otypes)
2538 if func is not self.pyfunc or nin not in self._ufunc:
2539 ufunc = frompyfunc(func, nin, nout)
2540 else:
2541 ufunc = None # We'll get it from self._ufunc
2542 if func is self.pyfunc:
2543 ufunc = self._ufunc.setdefault(nin, ufunc)
2544 else:
2545 # Get number of outputs and output types by calling the function on
2546 # the first entries of args. We also cache the result to prevent
2547 # the subsequent call when the ufunc is evaluated.
2548 # Assumes that ufunc first evaluates the 0th elements in the input
2549 # arrays (the input values are not checked to ensure this)
2550 args = [asarray(a) for a in args]
2551 if builtins.any(arg.size == 0 for arg in args):
2552 raise ValueError('cannot call `vectorize` on size 0 inputs '
2553 'unless `otypes` is set')
2555 inputs = [arg.flat[0] for arg in args]
2556 outputs = func(*inputs)
2558 # Performance note: profiling indicates that -- for simple
2559 # functions at least -- this wrapping can almost double the
2560 # execution time.
2561 # Hence we make it optional.
2562 if self.cache:
2563 _cache = [outputs]
2565 def _func(*vargs):
2566 if _cache:
2567 return _cache.pop()
2568 else:
2569 return func(*vargs)
2570 else:
2571 _func = func
2573 if isinstance(outputs, tuple):
2574 nout = len(outputs)
2575 else:
2576 nout = 1
2577 outputs = (outputs,)
2579 otypes = ''.join([asarray(outputs[_k]).dtype.char
2580 for _k in range(nout)])
2582 # Performance note: profiling indicates that creating the ufunc is
2583 # not a significant cost compared with wrapping so it seems not
2584 # worth trying to cache this.
2585 ufunc = frompyfunc(_func, len(args), nout)
2587 return ufunc, otypes
2589 def _vectorize_call(self, func, args):
2590 """Vectorized call to `func` over positional `args`."""
2591 if self.signature is not None:
2592 res = self._vectorize_call_with_signature(func, args)
2593 elif not args:
2594 res = func()
2595 else:
2596 ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args)
2597 # gh-29196: `dtype=object` should eventually be removed
2598 args = [asanyarray(a, dtype=object) for a in args]
2599 outputs = ufunc(*args, out=...)
2601 if ufunc.nout == 1:
2602 res = asanyarray(outputs, dtype=otypes[0])
2603 else:
2604 res = tuple(asanyarray(x, dtype=t)
2605 for x, t in zip(outputs, otypes))
2606 return res
2608 def _vectorize_call_with_signature(self, func, args):
2609 """Vectorized call over positional arguments with a signature."""
2610 input_core_dims, output_core_dims = self._in_and_out_core_dims
2612 if len(args) != len(input_core_dims):
2613 raise TypeError('wrong number of positional arguments: '
2614 'expected %r, got %r'
2615 % (len(input_core_dims), len(args)))
2616 args = tuple(asanyarray(arg) for arg in args)
2618 broadcast_shape, dim_sizes = _parse_input_dimensions(
2619 args, input_core_dims)
2620 input_shapes = _calculate_shapes(broadcast_shape, dim_sizes,
2621 input_core_dims)
2622 args = [np.broadcast_to(arg, shape, subok=True)
2623 for arg, shape in zip(args, input_shapes)]
2625 outputs = None
2626 otypes = self.otypes
2627 nout = len(output_core_dims)
2629 for index in np.ndindex(*broadcast_shape):
2630 results = func(*(arg[index] for arg in args))
2632 n_results = len(results) if isinstance(results, tuple) else 1
2634 if nout != n_results:
2635 raise ValueError(
2636 'wrong number of outputs from pyfunc: expected %r, got %r'
2637 % (nout, n_results))
2639 if nout == 1:
2640 results = (results,)
2642 if outputs is None:
2643 for result, core_dims in zip(results, output_core_dims):
2644 _update_dim_sizes(dim_sizes, result, core_dims)
2646 outputs = _create_arrays(broadcast_shape, dim_sizes,
2647 output_core_dims, otypes, results)
2649 for output, result in zip(outputs, results):
2650 output[index] = result
2652 if outputs is None:
2653 # did not call the function even once
2654 if otypes is None:
2655 raise ValueError('cannot call `vectorize` on size 0 inputs '
2656 'unless `otypes` is set')
2657 if builtins.any(dim not in dim_sizes
2658 for dims in output_core_dims
2659 for dim in dims):
2660 raise ValueError('cannot call `vectorize` with a signature '
2661 'including new output dimensions on size 0 '
2662 'inputs')
2663 outputs = _create_arrays(broadcast_shape, dim_sizes,
2664 output_core_dims, otypes)
2666 return outputs[0] if nout == 1 else outputs
2669def _cov_dispatcher(m, y=None, rowvar=None, bias=None, ddof=None,
2670 fweights=None, aweights=None, *, dtype=None):
2671 return (m, y, fweights, aweights)
2674@array_function_dispatch(_cov_dispatcher)
2675def cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None,
2676 aweights=None, *, dtype=None):
2677 """
2678 Estimate a covariance matrix, given data and weights.
2680 Covariance indicates the level to which two variables vary together.
2681 If we examine N-dimensional samples, :math:`X = [x_1, x_2, ..., x_N]^T`,
2682 then the covariance matrix element :math:`C_{ij}` is the covariance of
2683 :math:`x_i` and :math:`x_j`. The element :math:`C_{ii}` is the variance
2684 of :math:`x_i`.
2686 See the notes for an outline of the algorithm.
2688 Parameters
2689 ----------
2690 m : array_like
2691 A 1-D or 2-D array containing multiple variables and observations.
2692 Each row of `m` represents a variable, and each column a single
2693 observation of all those variables. Also see `rowvar` below.
2694 y : array_like, optional
2695 An additional set of variables and observations. `y` has the same form
2696 as that of `m`.
2697 rowvar : bool, optional
2698 If `rowvar` is True (default), then each row represents a
2699 variable, with observations in the columns. Otherwise, the relationship
2700 is transposed: each column represents a variable, while the rows
2701 contain observations.
2702 bias : bool, optional
2703 Default normalization (False) is by ``(N - 1)``, where ``N`` is the
2704 number of observations given (unbiased estimate). If `bias` is True,
2705 then normalization is by ``N``. These values can be overridden by using
2706 the keyword ``ddof`` in numpy versions >= 1.5.
2707 ddof : int, optional
2708 If not ``None`` the default value implied by `bias` is overridden.
2709 Note that ``ddof=1`` will return the unbiased estimate, even if both
2710 `fweights` and `aweights` are specified, and ``ddof=0`` will return
2711 the simple average. See the notes for the details. The default value
2712 is ``None``.
2713 fweights : array_like, int, optional
2714 1-D array of integer frequency weights; the number of times each
2715 observation vector should be repeated.
2716 aweights : array_like, optional
2717 1-D array of observation vector weights. These relative weights are
2718 typically large for observations considered "important" and smaller for
2719 observations considered less "important". If ``ddof=0`` the array of
2720 weights can be used to assign probabilities to observation vectors.
2721 dtype : data-type, optional
2722 Data-type of the result. By default, the return data-type will have
2723 at least `numpy.float64` precision.
2725 .. versionadded:: 1.20
2727 Returns
2728 -------
2729 out : ndarray
2730 The covariance matrix of the variables.
2732 See Also
2733 --------
2734 corrcoef : Normalized covariance matrix
2736 Notes
2737 -----
2738 Assume that the observations are in the columns of the observation
2739 array `m` and let ``f = fweights`` and ``a = aweights`` for brevity. The
2740 steps to compute the weighted covariance are as follows::
2742 >>> m = np.arange(10, dtype=np.float64)
2743 >>> f = np.arange(10) * 2
2744 >>> a = np.arange(10) ** 2.
2745 >>> ddof = 1
2746 >>> w = f * a
2747 >>> v1 = np.sum(w)
2748 >>> v2 = np.sum(w * a)
2749 >>> m -= np.sum(m * w, axis=None, keepdims=True) / v1
2750 >>> cov = np.dot(m * w, m.T) * v1 / (v1**2 - ddof * v2)
2752 Note that when ``a == 1``, the normalization factor
2753 ``v1 / (v1**2 - ddof * v2)`` goes over to ``1 / (np.sum(f) - ddof)``
2754 as it should.
2756 Examples
2757 --------
2758 >>> import numpy as np
2760 Consider two variables, :math:`x_0` and :math:`x_1`, which
2761 correlate perfectly, but in opposite directions:
2763 >>> x = np.array([[0, 2], [1, 1], [2, 0]]).T
2764 >>> x
2765 array([[0, 1, 2],
2766 [2, 1, 0]])
2768 Note how :math:`x_0` increases while :math:`x_1` decreases. The covariance
2769 matrix shows this clearly:
2771 >>> np.cov(x)
2772 array([[ 1., -1.],
2773 [-1., 1.]])
2775 Note that element :math:`C_{0,1}`, which shows the correlation between
2776 :math:`x_0` and :math:`x_1`, is negative.
2778 Further, note how `x` and `y` are combined:
2780 >>> x = [-2.1, -1, 4.3]
2781 >>> y = [3, 1.1, 0.12]
2782 >>> X = np.stack((x, y), axis=0)
2783 >>> np.cov(X)
2784 array([[11.71 , -4.286 ], # may vary
2785 [-4.286 , 2.144133]])
2786 >>> np.cov(x, y)
2787 array([[11.71 , -4.286 ], # may vary
2788 [-4.286 , 2.144133]])
2789 >>> np.cov(x)
2790 array(11.71)
2792 """
2793 # Check inputs
2794 if ddof is not None and ddof != int(ddof):
2795 raise ValueError(
2796 "ddof must be integer")
2798 # Handles complex arrays too
2799 m = np.asarray(m)
2800 if m.ndim > 2:
2801 raise ValueError("m has more than 2 dimensions")
2803 if y is not None:
2804 y = np.asarray(y)
2805 if y.ndim > 2:
2806 raise ValueError("y has more than 2 dimensions")
2808 if dtype is None:
2809 if y is None:
2810 dtype = np.result_type(m, np.float64)
2811 else:
2812 dtype = np.result_type(m, y, np.float64)
2814 X = array(m, ndmin=2, dtype=dtype)
2815 if not rowvar and m.ndim != 1:
2816 X = X.T
2817 if X.shape[0] == 0:
2818 return np.array([]).reshape(0, 0)
2819 if y is not None:
2820 y = array(y, copy=None, ndmin=2, dtype=dtype)
2821 if not rowvar and y.shape[0] != 1:
2822 y = y.T
2823 X = np.concatenate((X, y), axis=0)
2825 if ddof is None:
2826 if bias == 0:
2827 ddof = 1
2828 else:
2829 ddof = 0
2831 # Get the product of frequencies and weights
2832 w = None
2833 if fweights is not None:
2834 fweights = np.asarray(fweights, dtype=float)
2835 if not np.all(fweights == np.around(fweights)):
2836 raise TypeError(
2837 "fweights must be integer")
2838 if fweights.ndim > 1:
2839 raise RuntimeError(
2840 "cannot handle multidimensional fweights")
2841 if fweights.shape[0] != X.shape[1]:
2842 raise RuntimeError(
2843 "incompatible numbers of samples and fweights")
2844 if any(fweights < 0):
2845 raise ValueError(
2846 "fweights cannot be negative")
2847 w = fweights
2848 if aweights is not None:
2849 aweights = np.asarray(aweights, dtype=float)
2850 if aweights.ndim > 1:
2851 raise RuntimeError(
2852 "cannot handle multidimensional aweights")
2853 if aweights.shape[0] != X.shape[1]:
2854 raise RuntimeError(
2855 "incompatible numbers of samples and aweights")
2856 if any(aweights < 0):
2857 raise ValueError(
2858 "aweights cannot be negative")
2859 if w is None:
2860 w = aweights
2861 else:
2862 w *= aweights
2864 avg, w_sum = average(X, axis=1, weights=w, returned=True)
2865 w_sum = w_sum[0]
2867 # Determine the normalization
2868 if w is None:
2869 fact = X.shape[1] - ddof
2870 elif ddof == 0:
2871 fact = w_sum
2872 elif aweights is None:
2873 fact = w_sum - ddof
2874 else:
2875 fact = w_sum - ddof * sum(w * aweights) / w_sum
2877 if fact <= 0:
2878 warnings.warn("Degrees of freedom <= 0 for slice",
2879 RuntimeWarning, stacklevel=2)
2880 fact = 0.0
2882 X -= avg[:, None]
2883 if w is None:
2884 X_T = X.T
2885 else:
2886 X_T = (X * w).T
2887 c = dot(X, X_T.conj())
2888 c *= np.true_divide(1, fact)
2889 return c.squeeze()
2892def _corrcoef_dispatcher(x, y=None, rowvar=None, *,
2893 dtype=None):
2894 return (x, y)
2897@array_function_dispatch(_corrcoef_dispatcher)
2898def corrcoef(x, y=None, rowvar=True, *,
2899 dtype=None):
2900 """
2901 Return Pearson product-moment correlation coefficients.
2903 Please refer to the documentation for `cov` for more detail. The
2904 relationship between the correlation coefficient matrix, `R`, and the
2905 covariance matrix, `C`, is
2907 .. math:: R_{ij} = \\frac{ C_{ij} } { \\sqrt{ C_{ii} C_{jj} } }
2909 The values of `R` are between -1 and 1, inclusive.
2911 Parameters
2912 ----------
2913 x : array_like
2914 A 1-D or 2-D array containing multiple variables and observations.
2915 Each row of `x` represents a variable, and each column a single
2916 observation of all those variables. Also see `rowvar` below.
2917 y : array_like, optional
2918 An additional set of variables and observations. `y` has the same
2919 shape as `x`.
2920 rowvar : bool, optional
2921 If `rowvar` is True (default), then each row represents a
2922 variable, with observations in the columns. Otherwise, the relationship
2923 is transposed: each column represents a variable, while the rows
2924 contain observations.
2926 dtype : data-type, optional
2927 Data-type of the result. By default, the return data-type will have
2928 at least `numpy.float64` precision.
2930 .. versionadded:: 1.20
2932 Returns
2933 -------
2934 R : ndarray
2935 The correlation coefficient matrix of the variables.
2937 See Also
2938 --------
2939 cov : Covariance matrix
2941 Notes
2942 -----
2943 Due to floating point rounding the resulting array may not be Hermitian,
2944 the diagonal elements may not be 1, and the elements may not satisfy the
2945 inequality abs(a) <= 1. The real and imaginary parts are clipped to the
2946 interval [-1, 1] in an attempt to improve on that situation but is not
2947 much help in the complex case.
2949 Examples
2950 --------
2951 >>> import numpy as np
2953 In this example we generate two random arrays, ``xarr`` and ``yarr``, and
2954 compute the row-wise and column-wise Pearson correlation coefficients,
2955 ``R``. Since ``rowvar`` is true by default, we first find the row-wise
2956 Pearson correlation coefficients between the variables of ``xarr``.
2958 >>> import numpy as np
2959 >>> rng = np.random.default_rng(seed=42)
2960 >>> xarr = rng.random((3, 3))
2961 >>> xarr
2962 array([[0.77395605, 0.43887844, 0.85859792],
2963 [0.69736803, 0.09417735, 0.97562235],
2964 [0.7611397 , 0.78606431, 0.12811363]])
2965 >>> R1 = np.corrcoef(xarr)
2966 >>> R1
2967 array([[ 1. , 0.99256089, -0.68080986],
2968 [ 0.99256089, 1. , -0.76492172],
2969 [-0.68080986, -0.76492172, 1. ]])
2971 If we add another set of variables and observations ``yarr``, we can
2972 compute the row-wise Pearson correlation coefficients between the
2973 variables in ``xarr`` and ``yarr``.
2975 >>> yarr = rng.random((3, 3))
2976 >>> yarr
2977 array([[0.45038594, 0.37079802, 0.92676499],
2978 [0.64386512, 0.82276161, 0.4434142 ],
2979 [0.22723872, 0.55458479, 0.06381726]])
2980 >>> R2 = np.corrcoef(xarr, yarr)
2981 >>> R2
2982 array([[ 1. , 0.99256089, -0.68080986, 0.75008178, -0.934284 ,
2983 -0.99004057],
2984 [ 0.99256089, 1. , -0.76492172, 0.82502011, -0.97074098,
2985 -0.99981569],
2986 [-0.68080986, -0.76492172, 1. , -0.99507202, 0.89721355,
2987 0.77714685],
2988 [ 0.75008178, 0.82502011, -0.99507202, 1. , -0.93657855,
2989 -0.83571711],
2990 [-0.934284 , -0.97074098, 0.89721355, -0.93657855, 1. ,
2991 0.97517215],
2992 [-0.99004057, -0.99981569, 0.77714685, -0.83571711, 0.97517215,
2993 1. ]])
2995 Finally if we use the option ``rowvar=False``, the columns are now
2996 being treated as the variables and we will find the column-wise Pearson
2997 correlation coefficients between variables in ``xarr`` and ``yarr``.
2999 >>> R3 = np.corrcoef(xarr, yarr, rowvar=False)
3000 >>> R3
3001 array([[ 1. , 0.77598074, -0.47458546, -0.75078643, -0.9665554 ,
3002 0.22423734],
3003 [ 0.77598074, 1. , -0.92346708, -0.99923895, -0.58826587,
3004 -0.44069024],
3005 [-0.47458546, -0.92346708, 1. , 0.93773029, 0.23297648,
3006 0.75137473],
3007 [-0.75078643, -0.99923895, 0.93773029, 1. , 0.55627469,
3008 0.47536961],
3009 [-0.9665554 , -0.58826587, 0.23297648, 0.55627469, 1. ,
3010 -0.46666491],
3011 [ 0.22423734, -0.44069024, 0.75137473, 0.47536961, -0.46666491,
3012 1. ]])
3014 """
3015 c = cov(x, y, rowvar, dtype=dtype)
3016 try:
3017 d = diag(c)
3018 except ValueError:
3019 # scalar covariance
3020 # nan if incorrect value (nan, inf, 0), 1 otherwise
3021 return c / c
3022 stddev = sqrt(d.real)
3023 c /= stddev[:, None]
3024 c /= stddev[None, :]
3026 # Clip real and imaginary parts to [-1, 1]. This does not guarantee
3027 # abs(a[i,j]) <= 1 for complex arrays, but is the best we can do without
3028 # excessive work.
3029 np.clip(c.real, -1, 1, out=c.real)
3030 if np.iscomplexobj(c):
3031 np.clip(c.imag, -1, 1, out=c.imag)
3033 return c
3036@set_module('numpy')
3037def blackman(M):
3038 """
3039 Return the Blackman window.
3041 The Blackman window is a taper formed by using the first three
3042 terms of a summation of cosines. It was designed to have close to the
3043 minimal leakage possible. It is close to optimal, only slightly worse
3044 than a Kaiser window.
3046 Parameters
3047 ----------
3048 M : int
3049 Number of points in the output window. If zero or less, an empty
3050 array is returned.
3052 Returns
3053 -------
3054 out : ndarray
3055 The window, with the maximum value normalized to one (the value one
3056 appears only if the number of samples is odd).
3058 See Also
3059 --------
3060 bartlett, hamming, hanning, kaiser
3062 Notes
3063 -----
3064 The Blackman window is defined as
3066 .. math:: w(n) = 0.42 - 0.5 \\cos(2\\pi n/M) + 0.08 \\cos(4\\pi n/M)
3068 Most references to the Blackman window come from the signal processing
3069 literature, where it is used as one of many windowing functions for
3070 smoothing values. It is also known as an apodization (which means
3071 "removing the foot", i.e. smoothing discontinuities at the beginning
3072 and end of the sampled signal) or tapering function. It is known as a
3073 "near optimal" tapering function, almost as good (by some measures)
3074 as the kaiser window.
3076 References
3077 ----------
3078 Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra,
3079 Dover Publications, New York.
3081 Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing.
3082 Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471.
3084 Examples
3085 --------
3086 >>> import numpy as np
3087 >>> import matplotlib.pyplot as plt
3088 >>> np.blackman(12)
3089 array([-1.38777878e-17, 3.26064346e-02, 1.59903635e-01, # may vary
3090 4.14397981e-01, 7.36045180e-01, 9.67046769e-01,
3091 9.67046769e-01, 7.36045180e-01, 4.14397981e-01,
3092 1.59903635e-01, 3.26064346e-02, -1.38777878e-17])
3094 Plot the window and the frequency response.
3096 .. plot::
3097 :include-source:
3099 import matplotlib.pyplot as plt
3100 from numpy.fft import fft, fftshift
3101 window = np.blackman(51)
3102 plt.plot(window)
3103 plt.title("Blackman window")
3104 plt.ylabel("Amplitude")
3105 plt.xlabel("Sample")
3106 plt.show() # doctest: +SKIP
3108 plt.figure()
3109 A = fft(window, 2048) / 25.5
3110 mag = np.abs(fftshift(A))
3111 freq = np.linspace(-0.5, 0.5, len(A))
3112 with np.errstate(divide='ignore', invalid='ignore'):
3113 response = 20 * np.log10(mag)
3114 response = np.clip(response, -100, 100)
3115 plt.plot(freq, response)
3116 plt.title("Frequency response of Blackman window")
3117 plt.ylabel("Magnitude [dB]")
3118 plt.xlabel("Normalized frequency [cycles per sample]")
3119 plt.axis('tight')
3120 plt.show()
3122 """
3123 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3124 # to double is safe for a range.
3125 values = np.array([0.0, M])
3126 M = values[1]
3128 if M < 1:
3129 return array([], dtype=values.dtype)
3130 if M == 1:
3131 return ones(1, dtype=values.dtype)
3132 n = arange(1 - M, M, 2)
3133 return 0.42 + 0.5 * cos(pi * n / (M - 1)) + 0.08 * cos(2.0 * pi * n / (M - 1))
3136@set_module('numpy')
3137def bartlett(M):
3138 """
3139 Return the Bartlett window.
3141 The Bartlett window is very similar to a triangular window, except
3142 that the end points are at zero. It is often used in signal
3143 processing for tapering a signal, without generating too much
3144 ripple in the frequency domain.
3146 Parameters
3147 ----------
3148 M : int
3149 Number of points in the output window. If zero or less, an
3150 empty array is returned.
3152 Returns
3153 -------
3154 out : array
3155 The triangular window, with the maximum value normalized to one
3156 (the value one appears only if the number of samples is odd), with
3157 the first and last samples equal to zero.
3159 See Also
3160 --------
3161 blackman, hamming, hanning, kaiser
3163 Notes
3164 -----
3165 The Bartlett window is defined as
3167 .. math:: w(n) = \\frac{2}{M-1} \\left(
3168 \\frac{M-1}{2} - \\left|n - \\frac{M-1}{2}\\right|
3169 \\right)
3171 Most references to the Bartlett window come from the signal processing
3172 literature, where it is used as one of many windowing functions for
3173 smoothing values. Note that convolution with this window produces linear
3174 interpolation. It is also known as an apodization (which means "removing
3175 the foot", i.e. smoothing discontinuities at the beginning and end of the
3176 sampled signal) or tapering function. The Fourier transform of the
3177 Bartlett window is the product of two sinc functions. Note the excellent
3178 discussion in Kanasewich [2]_.
3180 References
3181 ----------
3182 .. [1] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra",
3183 Biometrika 37, 1-16, 1950.
3184 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
3185 The University of Alberta Press, 1975, pp. 109-110.
3186 .. [3] A.V. Oppenheim and R.W. Schafer, "Discrete-Time Signal
3187 Processing", Prentice-Hall, 1999, pp. 468-471.
3188 .. [4] Wikipedia, "Window function",
3189 https://en.wikipedia.org/wiki/Window_function
3190 .. [5] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3191 "Numerical Recipes", Cambridge University Press, 1986, page 429.
3193 Examples
3194 --------
3195 >>> import numpy as np
3196 >>> import matplotlib.pyplot as plt
3197 >>> np.bartlett(12)
3198 array([ 0. , 0.18181818, 0.36363636, 0.54545455, 0.72727273, # may vary
3199 0.90909091, 0.90909091, 0.72727273, 0.54545455, 0.36363636,
3200 0.18181818, 0. ])
3202 Plot the window and its frequency response (requires SciPy and matplotlib).
3204 .. plot::
3205 :include-source:
3207 import matplotlib.pyplot as plt
3208 from numpy.fft import fft, fftshift
3209 window = np.bartlett(51)
3210 plt.plot(window)
3211 plt.title("Bartlett window")
3212 plt.ylabel("Amplitude")
3213 plt.xlabel("Sample")
3214 plt.show()
3215 plt.figure()
3216 A = fft(window, 2048) / 25.5
3217 mag = np.abs(fftshift(A))
3218 freq = np.linspace(-0.5, 0.5, len(A))
3219 with np.errstate(divide='ignore', invalid='ignore'):
3220 response = 20 * np.log10(mag)
3221 response = np.clip(response, -100, 100)
3222 plt.plot(freq, response)
3223 plt.title("Frequency response of Bartlett window")
3224 plt.ylabel("Magnitude [dB]")
3225 plt.xlabel("Normalized frequency [cycles per sample]")
3226 plt.axis('tight')
3227 plt.show()
3229 """
3230 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3231 # to double is safe for a range.
3232 values = np.array([0.0, M])
3233 M = values[1]
3235 if M < 1:
3236 return array([], dtype=values.dtype)
3237 if M == 1:
3238 return ones(1, dtype=values.dtype)
3239 n = arange(1 - M, M, 2)
3240 return where(less_equal(n, 0), 1 + n / (M - 1), 1 - n / (M - 1))
3243@set_module('numpy')
3244def hanning(M):
3245 """
3246 Return the Hanning window.
3248 The Hanning window is a taper formed by using a weighted cosine.
3250 Parameters
3251 ----------
3252 M : int
3253 Number of points in the output window. If zero or less, an
3254 empty array is returned.
3256 Returns
3257 -------
3258 out : ndarray, shape(M,)
3259 The window, with the maximum value normalized to one (the value
3260 one appears only if `M` is odd).
3262 See Also
3263 --------
3264 bartlett, blackman, hamming, kaiser
3266 Notes
3267 -----
3268 The Hanning window is defined as
3270 .. math:: w(n) = 0.5 - 0.5\\cos\\left(\\frac{2\\pi{n}}{M-1}\\right)
3271 \\qquad 0 \\leq n \\leq M-1
3273 The Hanning was named for Julius von Hann, an Austrian meteorologist.
3274 It is also known as the Cosine Bell. Some authors prefer that it be
3275 called a Hann window, to help avoid confusion with the very similar
3276 Hamming window.
3278 Most references to the Hanning window come from the signal processing
3279 literature, where it is used as one of many windowing functions for
3280 smoothing values. It is also known as an apodization (which means
3281 "removing the foot", i.e. smoothing discontinuities at the beginning
3282 and end of the sampled signal) or tapering function.
3284 References
3285 ----------
3286 .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
3287 spectra, Dover Publications, New York.
3288 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
3289 The University of Alberta Press, 1975, pp. 106-108.
3290 .. [3] Wikipedia, "Window function",
3291 https://en.wikipedia.org/wiki/Window_function
3292 .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3293 "Numerical Recipes", Cambridge University Press, 1986, page 425.
3295 Examples
3296 --------
3297 >>> import numpy as np
3298 >>> np.hanning(12)
3299 array([0. , 0.07937323, 0.29229249, 0.57115742, 0.82743037,
3300 0.97974649, 0.97974649, 0.82743037, 0.57115742, 0.29229249,
3301 0.07937323, 0. ])
3303 Plot the window and its frequency response.
3305 .. plot::
3306 :include-source:
3308 import matplotlib.pyplot as plt
3309 from numpy.fft import fft, fftshift
3310 window = np.hanning(51)
3311 plt.plot(window)
3312 plt.title("Hann window")
3313 plt.ylabel("Amplitude")
3314 plt.xlabel("Sample")
3315 plt.show()
3317 plt.figure()
3318 A = fft(window, 2048) / 25.5
3319 mag = np.abs(fftshift(A))
3320 freq = np.linspace(-0.5, 0.5, len(A))
3321 with np.errstate(divide='ignore', invalid='ignore'):
3322 response = 20 * np.log10(mag)
3323 response = np.clip(response, -100, 100)
3324 plt.plot(freq, response)
3325 plt.title("Frequency response of the Hann window")
3326 plt.ylabel("Magnitude [dB]")
3327 plt.xlabel("Normalized frequency [cycles per sample]")
3328 plt.axis('tight')
3329 plt.show()
3331 """
3332 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3333 # to double is safe for a range.
3334 values = np.array([0.0, M])
3335 M = values[1]
3337 if M < 1:
3338 return array([], dtype=values.dtype)
3339 if M == 1:
3340 return ones(1, dtype=values.dtype)
3341 n = arange(1 - M, M, 2)
3342 return 0.5 + 0.5 * cos(pi * n / (M - 1))
3345@set_module('numpy')
3346def hamming(M):
3347 """
3348 Return the Hamming window.
3350 The Hamming window is a taper formed by using a weighted cosine.
3352 Parameters
3353 ----------
3354 M : int
3355 Number of points in the output window. If zero or less, an
3356 empty array is returned.
3358 Returns
3359 -------
3360 out : ndarray
3361 The window, with the maximum value normalized to one (the value
3362 one appears only if the number of samples is odd).
3364 See Also
3365 --------
3366 bartlett, blackman, hanning, kaiser
3368 Notes
3369 -----
3370 The Hamming window is defined as
3372 .. math:: w(n) = 0.54 - 0.46\\cos\\left(\\frac{2\\pi{n}}{M-1}\\right)
3373 \\qquad 0 \\leq n \\leq M-1
3375 The Hamming was named for R. W. Hamming, an associate of J. W. Tukey
3376 and is described in Blackman and Tukey. It was recommended for
3377 smoothing the truncated autocovariance function in the time domain.
3378 Most references to the Hamming window come from the signal processing
3379 literature, where it is used as one of many windowing functions for
3380 smoothing values. It is also known as an apodization (which means
3381 "removing the foot", i.e. smoothing discontinuities at the beginning
3382 and end of the sampled signal) or tapering function.
3384 References
3385 ----------
3386 .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
3387 spectra, Dover Publications, New York.
3388 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
3389 University of Alberta Press, 1975, pp. 109-110.
3390 .. [3] Wikipedia, "Window function",
3391 https://en.wikipedia.org/wiki/Window_function
3392 .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3393 "Numerical Recipes", Cambridge University Press, 1986, page 425.
3395 Examples
3396 --------
3397 >>> import numpy as np
3398 >>> np.hamming(12)
3399 array([ 0.08 , 0.15302337, 0.34890909, 0.60546483, 0.84123594, # may vary
3400 0.98136677, 0.98136677, 0.84123594, 0.60546483, 0.34890909,
3401 0.15302337, 0.08 ])
3403 Plot the window and the frequency response.
3405 .. plot::
3406 :include-source:
3408 import matplotlib.pyplot as plt
3409 from numpy.fft import fft, fftshift
3410 window = np.hamming(51)
3411 plt.plot(window)
3412 plt.title("Hamming window")
3413 plt.ylabel("Amplitude")
3414 plt.xlabel("Sample")
3415 plt.show()
3417 plt.figure()
3418 A = fft(window, 2048) / 25.5
3419 mag = np.abs(fftshift(A))
3420 freq = np.linspace(-0.5, 0.5, len(A))
3421 response = 20 * np.log10(mag)
3422 response = np.clip(response, -100, 100)
3423 plt.plot(freq, response)
3424 plt.title("Frequency response of Hamming window")
3425 plt.ylabel("Magnitude [dB]")
3426 plt.xlabel("Normalized frequency [cycles per sample]")
3427 plt.axis('tight')
3428 plt.show()
3430 """
3431 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3432 # to double is safe for a range.
3433 values = np.array([0.0, M])
3434 M = values[1]
3436 if M < 1:
3437 return array([], dtype=values.dtype)
3438 if M == 1:
3439 return ones(1, dtype=values.dtype)
3440 n = arange(1 - M, M, 2)
3441 return 0.54 + 0.46 * cos(pi * n / (M - 1))
3444## Code from cephes for i0
3446_i0A = [
3447 -4.41534164647933937950E-18,
3448 3.33079451882223809783E-17,
3449 -2.43127984654795469359E-16,
3450 1.71539128555513303061E-15,
3451 -1.16853328779934516808E-14,
3452 7.67618549860493561688E-14,
3453 -4.85644678311192946090E-13,
3454 2.95505266312963983461E-12,
3455 -1.72682629144155570723E-11,
3456 9.67580903537323691224E-11,
3457 -5.18979560163526290666E-10,
3458 2.65982372468238665035E-9,
3459 -1.30002500998624804212E-8,
3460 6.04699502254191894932E-8,
3461 -2.67079385394061173391E-7,
3462 1.11738753912010371815E-6,
3463 -4.41673835845875056359E-6,
3464 1.64484480707288970893E-5,
3465 -5.75419501008210370398E-5,
3466 1.88502885095841655729E-4,
3467 -5.76375574538582365885E-4,
3468 1.63947561694133579842E-3,
3469 -4.32430999505057594430E-3,
3470 1.05464603945949983183E-2,
3471 -2.37374148058994688156E-2,
3472 4.93052842396707084878E-2,
3473 -9.49010970480476444210E-2,
3474 1.71620901522208775349E-1,
3475 -3.04682672343198398683E-1,
3476 6.76795274409476084995E-1
3477 ]
3479_i0B = [
3480 -7.23318048787475395456E-18,
3481 -4.83050448594418207126E-18,
3482 4.46562142029675999901E-17,
3483 3.46122286769746109310E-17,
3484 -2.82762398051658348494E-16,
3485 -3.42548561967721913462E-16,
3486 1.77256013305652638360E-15,
3487 3.81168066935262242075E-15,
3488 -9.55484669882830764870E-15,
3489 -4.15056934728722208663E-14,
3490 1.54008621752140982691E-14,
3491 3.85277838274214270114E-13,
3492 7.18012445138366623367E-13,
3493 -1.79417853150680611778E-12,
3494 -1.32158118404477131188E-11,
3495 -3.14991652796324136454E-11,
3496 1.18891471078464383424E-11,
3497 4.94060238822496958910E-10,
3498 3.39623202570838634515E-9,
3499 2.26666899049817806459E-8,
3500 2.04891858946906374183E-7,
3501 2.89137052083475648297E-6,
3502 6.88975834691682398426E-5,
3503 3.36911647825569408990E-3,
3504 8.04490411014108831608E-1
3505 ]
3508def _chbevl(x, vals):
3509 b0 = vals[0]
3510 b1 = 0.0
3512 for i in range(1, len(vals)):
3513 b2 = b1
3514 b1 = b0
3515 b0 = x * b1 - b2 + vals[i]
3517 return 0.5 * (b0 - b2)
3520def _i0_1(x):
3521 return exp(x) * _chbevl(x / 2.0 - 2, _i0A)
3524def _i0_2(x):
3525 return exp(x) * _chbevl(32.0 / x - 2.0, _i0B) / sqrt(x)
3528def _i0_dispatcher(x):
3529 return (x,)
3532@array_function_dispatch(_i0_dispatcher)
3533def i0(x):
3534 """
3535 Modified Bessel function of the first kind, order 0.
3537 Usually denoted :math:`I_0`.
3539 Parameters
3540 ----------
3541 x : array_like of float
3542 Argument of the Bessel function.
3544 Returns
3545 -------
3546 out : ndarray, shape = x.shape, dtype = float
3547 The modified Bessel function evaluated at each of the elements of `x`.
3549 See Also
3550 --------
3551 scipy.special.i0, scipy.special.iv, scipy.special.ive
3553 Notes
3554 -----
3555 The scipy implementation is recommended over this function: it is a
3556 proper ufunc written in C, and more than an order of magnitude faster.
3558 We use the algorithm published by Clenshaw [1]_ and referenced by
3559 Abramowitz and Stegun [2]_, for which the function domain is
3560 partitioned into the two intervals [0,8] and (8,inf), and Chebyshev
3561 polynomial expansions are employed in each interval. Relative error on
3562 the domain [0,30] using IEEE arithmetic is documented [3]_ as having a
3563 peak of 5.8e-16 with an rms of 1.4e-16 (n = 30000).
3565 References
3566 ----------
3567 .. [1] C. W. Clenshaw, "Chebyshev series for mathematical functions", in
3568 *National Physical Laboratory Mathematical Tables*, vol. 5, London:
3569 Her Majesty's Stationery Office, 1962.
3570 .. [2] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical
3571 Functions*, 10th printing, New York: Dover, 1964, pp. 379.
3572 https://personal.math.ubc.ca/~cbm/aands/page_379.htm
3573 .. [3] https://metacpan.org/pod/distribution/Math-Cephes/lib/Math/Cephes.pod#i0:-Modified-Bessel-function-of-order-zero
3575 Examples
3576 --------
3577 >>> import numpy as np
3578 >>> np.i0(0.)
3579 array(1.0)
3580 >>> np.i0([0, 1, 2, 3])
3581 array([1. , 1.26606588, 2.2795853 , 4.88079259])
3583 """
3584 x = np.asanyarray(x)
3585 if x.dtype.kind == 'c':
3586 raise TypeError("i0 not supported for complex values")
3587 if x.dtype.kind != 'f':
3588 x = x.astype(float)
3589 x = np.abs(x)
3590 return piecewise(x, [x <= 8.0], [_i0_1, _i0_2])
3592## End of cephes code for i0
3595@set_module('numpy')
3596def kaiser(M, beta):
3597 """
3598 Return the Kaiser window.
3600 The Kaiser window is a taper formed by using a Bessel function.
3602 Parameters
3603 ----------
3604 M : int
3605 Number of points in the output window. If zero or less, an
3606 empty array is returned.
3607 beta : float
3608 Shape parameter for window.
3610 Returns
3611 -------
3612 out : array
3613 The window, with the maximum value normalized to one (the value
3614 one appears only if the number of samples is odd).
3616 See Also
3617 --------
3618 bartlett, blackman, hamming, hanning
3620 Notes
3621 -----
3622 The Kaiser window is defined as
3624 .. math:: w(n) = I_0\\left( \\beta \\sqrt{1-\\frac{4n^2}{(M-1)^2}}
3625 \\right)/I_0(\\beta)
3627 with
3629 .. math:: \\quad -\\frac{M-1}{2} \\leq n \\leq \\frac{M-1}{2},
3631 where :math:`I_0` is the modified zeroth-order Bessel function.
3633 The Kaiser was named for Jim Kaiser, who discovered a simple
3634 approximation to the DPSS window based on Bessel functions. The Kaiser
3635 window is a very good approximation to the Digital Prolate Spheroidal
3636 Sequence, or Slepian window, which is the transform which maximizes the
3637 energy in the main lobe of the window relative to total energy.
3639 The Kaiser can approximate many other windows by varying the beta
3640 parameter.
3642 ==== =======================
3643 beta Window shape
3644 ==== =======================
3645 0 Rectangular
3646 5 Similar to a Hamming
3647 6 Similar to a Hanning
3648 8.6 Similar to a Blackman
3649 ==== =======================
3651 A beta value of 14 is probably a good starting point. Note that as beta
3652 gets large, the window narrows, and so the number of samples needs to be
3653 large enough to sample the increasingly narrow spike, otherwise NaNs will
3654 get returned.
3656 Most references to the Kaiser window come from the signal processing
3657 literature, where it is used as one of many windowing functions for
3658 smoothing values. It is also known as an apodization (which means
3659 "removing the foot", i.e. smoothing discontinuities at the beginning
3660 and end of the sampled signal) or tapering function.
3662 References
3663 ----------
3664 .. [1] J. F. Kaiser, "Digital Filters" - Ch 7 in "Systems analysis by
3665 digital computer", Editors: F.F. Kuo and J.F. Kaiser, p 218-285.
3666 John Wiley and Sons, New York, (1966).
3667 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
3668 University of Alberta Press, 1975, pp. 177-178.
3669 .. [3] Wikipedia, "Window function",
3670 https://en.wikipedia.org/wiki/Window_function
3672 Examples
3673 --------
3674 >>> import numpy as np
3675 >>> import matplotlib.pyplot as plt
3676 >>> np.kaiser(12, 14)
3677 array([7.72686684e-06, 3.46009194e-03, 4.65200189e-02, # may vary
3678 2.29737120e-01, 5.99885316e-01, 9.45674898e-01,
3679 9.45674898e-01, 5.99885316e-01, 2.29737120e-01,
3680 4.65200189e-02, 3.46009194e-03, 7.72686684e-06])
3683 Plot the window and the frequency response.
3685 .. plot::
3686 :include-source:
3688 import matplotlib.pyplot as plt
3689 from numpy.fft import fft, fftshift
3690 window = np.kaiser(51, 14)
3691 plt.plot(window)
3692 plt.title("Kaiser window")
3693 plt.ylabel("Amplitude")
3694 plt.xlabel("Sample")
3695 plt.show()
3697 plt.figure()
3698 A = fft(window, 2048) / 25.5
3699 mag = np.abs(fftshift(A))
3700 freq = np.linspace(-0.5, 0.5, len(A))
3701 response = 20 * np.log10(mag)
3702 response = np.clip(response, -100, 100)
3703 plt.plot(freq, response)
3704 plt.title("Frequency response of Kaiser window")
3705 plt.ylabel("Magnitude [dB]")
3706 plt.xlabel("Normalized frequency [cycles per sample]")
3707 plt.axis('tight')
3708 plt.show()
3710 """
3711 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3712 # to double is safe for a range. (Simplified result_type with 0.0
3713 # strongly typed. result-type is not/less order sensitive, but that mainly
3714 # matters for integers anyway.)
3715 values = np.array([0.0, M, beta])
3716 M = values[1]
3717 beta = values[2]
3719 if M == 1:
3720 return np.ones(1, dtype=values.dtype)
3721 n = arange(0, M)
3722 alpha = (M - 1) / 2.0
3723 return i0(beta * sqrt(1 - ((n - alpha) / alpha)**2.0)) / i0(beta)
3726def _sinc_dispatcher(x):
3727 return (x,)
3730@array_function_dispatch(_sinc_dispatcher)
3731def sinc(x):
3732 r"""
3733 Return the normalized sinc function.
3735 The sinc function is equal to :math:`\sin(\pi x)/(\pi x)` for any argument
3736 :math:`x\ne 0`. ``sinc(0)`` takes the limit value 1, making ``sinc`` not
3737 only everywhere continuous but also infinitely differentiable.
3739 .. note::
3741 Note the normalization factor of ``pi`` used in the definition.
3742 This is the most commonly used definition in signal processing.
3743 Use ``sinc(x / np.pi)`` to obtain the unnormalized sinc function
3744 :math:`\sin(x)/x` that is more common in mathematics.
3746 Parameters
3747 ----------
3748 x : ndarray
3749 Array (possibly multi-dimensional) of values for which to calculate
3750 ``sinc(x)``.
3752 Returns
3753 -------
3754 out : ndarray
3755 ``sinc(x)``, which has the same shape as the input.
3757 Notes
3758 -----
3759 The name sinc is short for "sine cardinal" or "sinus cardinalis".
3761 The sinc function is used in various signal processing applications,
3762 including in anti-aliasing, in the construction of a Lanczos resampling
3763 filter, and in interpolation.
3765 For bandlimited interpolation of discrete-time signals, the ideal
3766 interpolation kernel is proportional to the sinc function.
3768 References
3769 ----------
3770 .. [1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram Web
3771 Resource. https://mathworld.wolfram.com/SincFunction.html
3772 .. [2] Wikipedia, "Sinc function",
3773 https://en.wikipedia.org/wiki/Sinc_function
3775 Examples
3776 --------
3777 >>> import numpy as np
3778 >>> import matplotlib.pyplot as plt
3779 >>> x = np.linspace(-4, 4, 41)
3780 >>> np.sinc(x)
3781 array([-3.89804309e-17, -4.92362781e-02, -8.40918587e-02, # may vary
3782 -8.90384387e-02, -5.84680802e-02, 3.89804309e-17,
3783 6.68206631e-02, 1.16434881e-01, 1.26137788e-01,
3784 8.50444803e-02, -3.89804309e-17, -1.03943254e-01,
3785 -1.89206682e-01, -2.16236208e-01, -1.55914881e-01,
3786 3.89804309e-17, 2.33872321e-01, 5.04551152e-01,
3787 7.56826729e-01, 9.35489284e-01, 1.00000000e+00,
3788 9.35489284e-01, 7.56826729e-01, 5.04551152e-01,
3789 2.33872321e-01, 3.89804309e-17, -1.55914881e-01,
3790 -2.16236208e-01, -1.89206682e-01, -1.03943254e-01,
3791 -3.89804309e-17, 8.50444803e-02, 1.26137788e-01,
3792 1.16434881e-01, 6.68206631e-02, 3.89804309e-17,
3793 -5.84680802e-02, -8.90384387e-02, -8.40918587e-02,
3794 -4.92362781e-02, -3.89804309e-17])
3796 >>> plt.plot(x, np.sinc(x))
3797 [<matplotlib.lines.Line2D object at 0x...>]
3798 >>> plt.title("Sinc Function")
3799 Text(0.5, 1.0, 'Sinc Function')
3800 >>> plt.ylabel("Amplitude")
3801 Text(0, 0.5, 'Amplitude')
3802 >>> plt.xlabel("X")
3803 Text(0.5, 0, 'X')
3804 >>> plt.show()
3806 """
3807 x = np.asanyarray(x)
3808 x = pi * x
3809 # Hope that 1e-20 is sufficient for objects...
3810 eps = np.finfo(x.dtype).eps if x.dtype.kind == "f" else 1e-20
3811 y = where(x, x, eps)
3812 return sin(y) / y
3815def _ureduce(a, func, keepdims=False, **kwargs):
3816 """
3817 Internal Function.
3818 Call `func` with `a` as first argument swapping the axes to use extended
3819 axis on functions that don't support it natively.
3821 Returns result and a.shape with axis dims set to 1.
3823 Parameters
3824 ----------
3825 a : array_like
3826 Input array or object that can be converted to an array.
3827 func : callable
3828 Reduction function capable of receiving a single axis argument.
3829 It is called with `a` as first argument followed by `kwargs`.
3830 kwargs : keyword arguments
3831 additional keyword arguments to pass to `func`.
3833 Returns
3834 -------
3835 result : tuple
3836 Result of func(a, **kwargs) and a.shape with axis dims set to 1
3837 which can be used to reshape the result to the same shape a ufunc with
3838 keepdims=True would produce.
3840 """
3841 a = np.asanyarray(a)
3842 axis = kwargs.get('axis')
3843 out = kwargs.get('out')
3845 if keepdims is np._NoValue:
3846 keepdims = False
3848 nd = a.ndim
3849 if axis is not None:
3850 axis = _nx.normalize_axis_tuple(axis, nd)
3852 if keepdims and out is not None:
3853 index_out = tuple(
3854 0 if i in axis else slice(None) for i in range(nd))
3855 kwargs['out'] = out[(Ellipsis, ) + index_out]
3857 if len(axis) == 1:
3858 kwargs['axis'] = axis[0]
3859 else:
3860 keep = sorted(set(range(nd)) - set(axis))
3861 nkeep = len(keep)
3863 def reshape_arr(a):
3864 # move axis that should not be reduced to front
3865 a = np.moveaxis(a, keep, range(nkeep))
3866 # merge reduced axis
3867 return a.reshape(a.shape[:nkeep] + (-1,))
3869 a = reshape_arr(a)
3871 weights = kwargs.get("weights")
3872 if weights is not None:
3873 kwargs["weights"] = reshape_arr(weights)
3875 kwargs['axis'] = -1
3876 elif keepdims and out is not None:
3877 index_out = (0, ) * nd
3878 kwargs['out'] = out[(Ellipsis, ) + index_out]
3880 r = func(a, **kwargs)
3882 if out is not None:
3883 return out
3885 if keepdims:
3886 if axis is None:
3887 index_r = (np.newaxis, ) * nd
3888 else:
3889 index_r = tuple(
3890 np.newaxis if i in axis else slice(None)
3891 for i in range(nd))
3892 r = r[(Ellipsis, ) + index_r]
3894 return r
3897def _median_dispatcher(
3898 a, axis=None, out=None, overwrite_input=None, keepdims=None):
3899 return (a, out)
3902@array_function_dispatch(_median_dispatcher)
3903def median(a, axis=None, out=None, overwrite_input=False, keepdims=False):
3904 """
3905 Compute the median along the specified axis.
3907 Returns the median of the array elements.
3909 Parameters
3910 ----------
3911 a : array_like
3912 Input array or object that can be converted to an array.
3913 axis : {int, sequence of int, None}, optional
3914 Axis or axes along which the medians are computed. The default,
3915 axis=None, will compute the median along a flattened version of
3916 the array. If a sequence of axes, the array is first flattened
3917 along the given axes, then the median is computed along the
3918 resulting flattened axis.
3919 out : ndarray, optional
3920 Alternative output array in which to place the result. It must
3921 have the same shape and buffer length as the expected output,
3922 but the type (of the output) will be cast if necessary.
3923 overwrite_input : bool, optional
3924 If True, then allow use of memory of input array `a` for
3925 calculations. The input array will be modified by the call to
3926 `median`. This will save memory when you do not need to preserve
3927 the contents of the input array. Treat the input as undefined,
3928 but it will probably be fully or partially sorted. Default is
3929 False. If `overwrite_input` is ``True`` and `a` is not already an
3930 `ndarray`, an error will be raised.
3931 keepdims : bool, optional
3932 If this is set to True, the axes which are reduced are left
3933 in the result as dimensions with size one. With this option,
3934 the result will broadcast correctly against the original `arr`.
3936 Returns
3937 -------
3938 median : ndarray
3939 A new array holding the result. If the input contains integers
3940 or floats smaller than ``float64``, then the output data-type is
3941 ``np.float64``. Otherwise, the data-type of the output is the
3942 same as that of the input. If `out` is specified, that array is
3943 returned instead.
3945 See Also
3946 --------
3947 mean, percentile
3949 Notes
3950 -----
3951 Given a vector ``V`` of length ``N``, the median of ``V`` is the
3952 middle value of a sorted copy of ``V``, ``V_sorted`` - i
3953 e., ``V_sorted[(N-1)/2]``, when ``N`` is odd, and the average of the
3954 two middle values of ``V_sorted`` when ``N`` is even.
3956 Examples
3957 --------
3958 >>> import numpy as np
3959 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
3960 >>> a
3961 array([[10, 7, 4],
3962 [ 3, 2, 1]])
3963 >>> np.median(a)
3964 np.float64(3.5)
3965 >>> np.median(a, axis=0)
3966 array([6.5, 4.5, 2.5])
3967 >>> np.median(a, axis=1)
3968 array([7., 2.])
3969 >>> np.median(a, axis=(0, 1))
3970 np.float64(3.5)
3971 >>> m = np.median(a, axis=0)
3972 >>> out = np.zeros_like(m)
3973 >>> np.median(a, axis=0, out=m)
3974 array([6.5, 4.5, 2.5])
3975 >>> m
3976 array([6.5, 4.5, 2.5])
3977 >>> b = a.copy()
3978 >>> np.median(b, axis=1, overwrite_input=True)
3979 array([7., 2.])
3980 >>> assert not np.all(a==b)
3981 >>> b = a.copy()
3982 >>> np.median(b, axis=None, overwrite_input=True)
3983 np.float64(3.5)
3984 >>> assert not np.all(a==b)
3986 """
3987 return _ureduce(a, func=_median, keepdims=keepdims, axis=axis, out=out,
3988 overwrite_input=overwrite_input)
3991def _median(a, axis=None, out=None, overwrite_input=False):
3992 # can't be reasonably be implemented in terms of percentile as we have to
3993 # call mean to not break astropy
3994 a = np.asanyarray(a)
3996 # Set the partition indexes
3997 if axis is None:
3998 sz = a.size
3999 else:
4000 sz = a.shape[axis]
4001 if sz % 2 == 0:
4002 szh = sz // 2
4003 kth = [szh - 1, szh]
4004 else:
4005 kth = [(sz - 1) // 2]
4007 # We have to check for NaNs (as of writing 'M' doesn't actually work).
4008 supports_nans = np.issubdtype(a.dtype, np.inexact) or a.dtype.kind in 'Mm'
4009 if supports_nans:
4010 kth.append(-1)
4012 if overwrite_input:
4013 if axis is None:
4014 part = a.ravel()
4015 part.partition(kth)
4016 else:
4017 a.partition(kth, axis=axis)
4018 part = a
4019 else:
4020 part = partition(a, kth, axis=axis)
4022 if part.shape == ():
4023 # make 0-D arrays work
4024 return part.item()
4025 if axis is None:
4026 axis = 0
4028 indexer = [slice(None)] * part.ndim
4029 index = part.shape[axis] // 2
4030 if part.shape[axis] % 2 == 1:
4031 # index with slice to allow mean (below) to work
4032 indexer[axis] = slice(index, index + 1)
4033 else:
4034 indexer[axis] = slice(index - 1, index + 1)
4035 indexer = tuple(indexer)
4037 # Use mean in both odd and even case to coerce data type,
4038 # using out array if needed.
4039 rout = mean(part[indexer], axis=axis, out=out)
4040 if supports_nans and sz > 0:
4041 # If nans are possible, warn and replace by nans like mean would.
4042 rout = np.lib._utils_impl._median_nancheck(part, rout, axis)
4044 return rout
4047def _percentile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
4048 method=None, keepdims=None, *, weights=None):
4049 return (a, q, out, weights)
4052@array_function_dispatch(_percentile_dispatcher)
4053def percentile(a,
4054 q,
4055 axis=None,
4056 out=None,
4057 overwrite_input=False,
4058 method="linear",
4059 keepdims=False,
4060 *,
4061 weights=None):
4062 """
4063 Compute the q-th percentile of the data along the specified axis.
4065 Returns the q-th percentile(s) of the array elements.
4067 Parameters
4068 ----------
4069 a : array_like of real numbers
4070 Input array or object that can be converted to an array.
4071 q : array_like of float
4072 Percentage or sequence of percentages for the percentiles to compute.
4073 Values must be between 0 and 100 inclusive.
4074 axis : {int, tuple of int, None}, optional
4075 Axis or axes along which the percentiles are computed. The
4076 default is to compute the percentile(s) along a flattened
4077 version of the array.
4078 out : ndarray, optional
4079 Alternative output array in which to place the result. It must
4080 have the same shape and buffer length as the expected output,
4081 but the type (of the output) will be cast if necessary.
4082 overwrite_input : bool, optional
4083 If True, then allow the input array `a` to be modified by intermediate
4084 calculations, to save memory. In this case, the contents of the input
4085 `a` after this function completes is undefined.
4086 method : str, optional
4087 This parameter specifies the method to use for estimating the
4088 percentile. There are many different methods, some unique to NumPy.
4089 See the notes for explanation. The options sorted by their R type
4090 as summarized in the H&F paper [1]_ are:
4092 1. 'inverted_cdf'
4093 2. 'averaged_inverted_cdf'
4094 3. 'closest_observation'
4095 4. 'interpolated_inverted_cdf'
4096 5. 'hazen'
4097 6. 'weibull'
4098 7. 'linear' (default)
4099 8. 'median_unbiased'
4100 9. 'normal_unbiased'
4102 The first three methods are discontinuous. NumPy further defines the
4103 following discontinuous variations of the default 'linear' (7.) option:
4105 * 'lower'
4106 * 'higher',
4107 * 'midpoint'
4108 * 'nearest'
4110 .. versionchanged:: 1.22.0
4111 This argument was previously called "interpolation" and only
4112 offered the "linear" default and last four options.
4114 keepdims : bool, optional
4115 If this is set to True, the axes which are reduced are left in
4116 the result as dimensions with size one. With this option, the
4117 result will broadcast correctly against the original array `a`.
4119 weights : array_like, optional
4120 An array of weights associated with the values in `a`. Each value in
4121 `a` contributes to the percentile according to its associated weight.
4122 The weights array can either be 1-D (in which case its length must be
4123 the size of `a` along the given axis) or of the same shape as `a`.
4124 If `weights=None`, then all data in `a` are assumed to have a
4125 weight equal to one.
4126 Only `method="inverted_cdf"` supports weights.
4127 See the notes for more details.
4129 .. versionadded:: 2.0.0
4131 Returns
4132 -------
4133 percentile : scalar or ndarray
4134 If `q` is a single percentile and `axis=None`, then the result
4135 is a scalar. If multiple percentiles are given, first axis of
4136 the result corresponds to the percentiles. The other axes are
4137 the axes that remain after the reduction of `a`. If the input
4138 contains integers or floats smaller than ``float64``, the output
4139 data-type is ``float64``. Otherwise, the output data-type is the
4140 same as that of the input. If `out` is specified, that array is
4141 returned instead.
4143 See Also
4144 --------
4145 mean
4146 median : equivalent to ``percentile(..., 50)``
4147 nanpercentile
4148 quantile : equivalent to percentile, except q in the range [0, 1].
4150 Notes
4151 -----
4152 The behavior of `numpy.percentile` with percentage `q` is
4153 that of `numpy.quantile` with argument ``q/100``.
4154 For more information, please see `numpy.quantile`.
4156 Examples
4157 --------
4158 >>> import numpy as np
4159 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
4160 >>> a
4161 array([[10, 7, 4],
4162 [ 3, 2, 1]])
4163 >>> np.percentile(a, 50)
4164 3.5
4165 >>> np.percentile(a, 50, axis=0)
4166 array([6.5, 4.5, 2.5])
4167 >>> np.percentile(a, 50, axis=1)
4168 array([7., 2.])
4169 >>> np.percentile(a, 50, axis=1, keepdims=True)
4170 array([[7.],
4171 [2.]])
4173 >>> m = np.percentile(a, 50, axis=0)
4174 >>> out = np.zeros_like(m)
4175 >>> np.percentile(a, 50, axis=0, out=out)
4176 array([6.5, 4.5, 2.5])
4177 >>> m
4178 array([6.5, 4.5, 2.5])
4180 >>> b = a.copy()
4181 >>> np.percentile(b, 50, axis=1, overwrite_input=True)
4182 array([7., 2.])
4183 >>> assert not np.all(a == b)
4185 The different methods can be visualized graphically:
4187 .. plot::
4189 import matplotlib.pyplot as plt
4191 a = np.arange(4)
4192 p = np.linspace(0, 100, 6001)
4193 ax = plt.gca()
4194 lines = [
4195 ('linear', '-', 'C0'),
4196 ('inverted_cdf', ':', 'C1'),
4197 # Almost the same as `inverted_cdf`:
4198 ('averaged_inverted_cdf', '-.', 'C1'),
4199 ('closest_observation', ':', 'C2'),
4200 ('interpolated_inverted_cdf', '--', 'C1'),
4201 ('hazen', '--', 'C3'),
4202 ('weibull', '-.', 'C4'),
4203 ('median_unbiased', '--', 'C5'),
4204 ('normal_unbiased', '-.', 'C6'),
4205 ]
4206 for method, style, color in lines:
4207 ax.plot(
4208 p, np.percentile(a, p, method=method),
4209 label=method, linestyle=style, color=color)
4210 ax.set(
4211 title='Percentiles for different methods and data: ' + str(a),
4212 xlabel='Percentile',
4213 ylabel='Estimated percentile value',
4214 yticks=a)
4215 ax.legend(bbox_to_anchor=(1.03, 1))
4216 plt.tight_layout()
4217 plt.show()
4219 References
4220 ----------
4221 .. [1] R. J. Hyndman and Y. Fan,
4222 "Sample quantiles in statistical packages,"
4223 The American Statistician, 50(4), pp. 361-365, 1996
4225 """
4226 a = np.asanyarray(a)
4227 if a.dtype.kind == "c":
4228 raise TypeError("a must be an array of real numbers")
4230 weak_q = type(q) in (int, float) # use weak promotion for final result type
4231 q = np.true_divide(q, 100, out=...)
4232 if not _quantile_is_valid(q):
4233 raise ValueError("Percentiles must be in the range [0, 100]")
4235 if weights is not None:
4236 if method != "inverted_cdf":
4237 msg = ("Only method 'inverted_cdf' supports weights. "
4238 f"Got: {method}.")
4239 raise ValueError(msg)
4240 if axis is not None:
4241 axis = _nx.normalize_axis_tuple(axis, a.ndim, argname="axis")
4242 weights = _weights_are_valid(weights=weights, a=a, axis=axis)
4243 if np.any(weights < 0):
4244 raise ValueError("Weights must be non-negative.")
4246 return _quantile_unchecked(
4247 a, q, axis, out, overwrite_input, method, keepdims, weights, weak_q)
4250def _quantile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
4251 method=None, keepdims=None, *, weights=None):
4252 return (a, q, out, weights)
4255@array_function_dispatch(_quantile_dispatcher)
4256def quantile(a,
4257 q,
4258 axis=None,
4259 out=None,
4260 overwrite_input=False,
4261 method="linear",
4262 keepdims=False,
4263 *,
4264 weights=None):
4265 """
4266 Compute the q-th quantile of the data along the specified axis.
4268 Parameters
4269 ----------
4270 a : array_like of real numbers
4271 Input array or object that can be converted to an array.
4272 q : array_like of float
4273 Probability or sequence of probabilities of the quantiles to compute.
4274 Values must be between 0 and 1 inclusive.
4275 axis : {int, tuple of int, None}, optional
4276 Axis or axes along which the quantiles are computed. The default is
4277 to compute the quantile(s) along a flattened version of the array.
4278 out : ndarray, optional
4279 Alternative output array in which to place the result. It must have
4280 the same shape and buffer length as the expected output, but the
4281 type (of the output) will be cast if necessary.
4282 overwrite_input : bool, optional
4283 If True, then allow the input array `a` to be modified by
4284 intermediate calculations, to save memory. In this case, the
4285 contents of the input `a` after this function completes is
4286 undefined.
4287 method : str, optional
4288 This parameter specifies the method to use for estimating the
4289 quantile. There are many different methods, some unique to NumPy.
4290 The recommended options, numbered as they appear in [1]_, are:
4292 1. 'inverted_cdf'
4293 2. 'averaged_inverted_cdf'
4294 3. 'closest_observation'
4295 4. 'interpolated_inverted_cdf'
4296 5. 'hazen'
4297 6. 'weibull'
4298 7. 'linear' (default)
4299 8. 'median_unbiased'
4300 9. 'normal_unbiased'
4302 The first three methods are discontinuous. For backward compatibility
4303 with previous versions of NumPy, the following discontinuous variations
4304 of the default 'linear' (7.) option are available:
4306 * 'lower'
4307 * 'higher',
4308 * 'midpoint'
4309 * 'nearest'
4311 See Notes for details.
4313 .. versionchanged:: 1.22.0
4314 This argument was previously called "interpolation" and only
4315 offered the "linear" default and last four options.
4317 keepdims : bool, optional
4318 If this is set to True, the axes which are reduced are left in
4319 the result as dimensions with size one. With this option, the
4320 result will broadcast correctly against the original array `a`.
4322 weights : array_like, optional
4323 An array of weights associated with the values in `a`. Each value in
4324 `a` contributes to the quantile according to its associated weight.
4325 The weights array can either be 1-D (in which case its length must be
4326 the size of `a` along the given axis) or of the same shape as `a`.
4327 If `weights=None`, then all data in `a` are assumed to have a
4328 weight equal to one.
4329 Only `method="inverted_cdf"` supports weights.
4330 See the notes for more details.
4332 .. versionadded:: 2.0.0
4334 Returns
4335 -------
4336 quantile : scalar or ndarray
4337 If `q` is a single probability and `axis=None`, then the result
4338 is a scalar. If multiple probability levels are given, first axis
4339 of the result corresponds to the quantiles. The other axes are
4340 the axes that remain after the reduction of `a`. If the input
4341 contains integers or floats smaller than ``float64``, the output
4342 data-type is ``float64``. Otherwise, the output data-type is the
4343 same as that of the input. If `out` is specified, that array is
4344 returned instead.
4346 See Also
4347 --------
4348 mean
4349 percentile : equivalent to quantile, but with q in the range [0, 100].
4350 median : equivalent to ``quantile(..., 0.5)``
4351 nanquantile
4353 Notes
4354 -----
4355 Given a sample `a` from an underlying distribution, `quantile` provides a
4356 nonparametric estimate of the inverse cumulative distribution function.
4358 By default, this is done by interpolating between adjacent elements in
4359 ``y``, a sorted copy of `a`::
4361 (1-g)*y[j] + g*y[j+1]
4363 where the index ``j`` and coefficient ``g`` are the integral and
4364 fractional components of ``q * (n-1)``, and ``n`` is the number of
4365 elements in the sample.
4367 This is a special case of Equation 1 of H&F [1]_. More generally,
4369 - ``j = (q*n + m - 1) // 1``, and
4370 - ``g = (q*n + m - 1) % 1``,
4372 where ``m`` may be defined according to several different conventions.
4373 The preferred convention may be selected using the ``method`` parameter:
4375 =============================== =============== ===============
4376 ``method`` number in H&F ``m``
4377 =============================== =============== ===============
4378 ``interpolated_inverted_cdf`` 4 ``0``
4379 ``hazen`` 5 ``1/2``
4380 ``weibull`` 6 ``q``
4381 ``linear`` (default) 7 ``1 - q``
4382 ``median_unbiased`` 8 ``q/3 + 1/3``
4383 ``normal_unbiased`` 9 ``q/4 + 3/8``
4384 =============================== =============== ===============
4386 Note that indices ``j`` and ``j + 1`` are clipped to the range ``0`` to
4387 ``n - 1`` when the results of the formula would be outside the allowed
4388 range of non-negative indices. The ``- 1`` in the formulas for ``j`` and
4389 ``g`` accounts for Python's 0-based indexing.
4391 The table above includes only the estimators from H&F that are continuous
4392 functions of probability `q` (estimators 4-9). NumPy also provides the
4393 three discontinuous estimators from H&F (estimators 1-3), where ``j`` is
4394 defined as above, ``m`` is defined as follows, and ``g`` is a function
4395 of the real-valued ``index = q*n + m - 1`` and ``j``.
4397 1. ``inverted_cdf``: ``m = 0`` and ``g = int(index - j > 0)``
4398 2. ``averaged_inverted_cdf``: ``m = 0`` and
4399 ``g = (1 + int(index - j > 0)) / 2``
4400 3. ``closest_observation``: ``m = -1/2`` and
4401 ``g = 1 - int((index == j) & (j%2 == 1))``
4403 For backward compatibility with previous versions of NumPy, `quantile`
4404 provides four additional discontinuous estimators. Like
4405 ``method='linear'``, all have ``m = 1 - q`` so that ``j = q*(n-1) // 1``,
4406 but ``g`` is defined as follows.
4408 - ``lower``: ``g = 0``
4409 - ``midpoint``: ``g = 0.5``
4410 - ``higher``: ``g = 1``
4411 - ``nearest``: ``g = (q*(n-1) % 1) > 0.5``
4413 **Weighted quantiles:**
4414 More formally, the quantile at probability level :math:`q` of a cumulative
4415 distribution function :math:`F(y)=P(Y \\leq y)` with probability measure
4416 :math:`P` is defined as any number :math:`x` that fulfills the
4417 *coverage conditions*
4419 .. math:: P(Y < x) \\leq q \\quad\\text{and}\\quad P(Y \\leq x) \\geq q
4421 with random variable :math:`Y\\sim P`.
4422 Sample quantiles, the result of `quantile`, provide nonparametric
4423 estimation of the underlying population counterparts, represented by the
4424 unknown :math:`F`, given a data vector `a` of length ``n``.
4426 Some of the estimators above arise when one considers :math:`F` as the
4427 empirical distribution function of the data, i.e.
4428 :math:`F(y) = \\frac{1}{n} \\sum_i 1_{a_i \\leq y}`.
4429 Then, different methods correspond to different choices of :math:`x` that
4430 fulfill the above coverage conditions. Methods that follow this approach
4431 are ``inverted_cdf`` and ``averaged_inverted_cdf``.
4433 For weighted quantiles, the coverage conditions still hold. The
4434 empirical cumulative distribution is simply replaced by its weighted
4435 version, i.e.
4436 :math:`P(Y \\leq t) = \\frac{1}{\\sum_i w_i} \\sum_i w_i 1_{x_i \\leq t}`.
4437 Only ``method="inverted_cdf"`` supports weights.
4439 Examples
4440 --------
4441 >>> import numpy as np
4442 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
4443 >>> a
4444 array([[10, 7, 4],
4445 [ 3, 2, 1]])
4446 >>> np.quantile(a, 0.5)
4447 3.5
4448 >>> np.quantile(a, 0.5, axis=0)
4449 array([6.5, 4.5, 2.5])
4450 >>> np.quantile(a, 0.5, axis=1)
4451 array([7., 2.])
4452 >>> np.quantile(a, 0.5, axis=1, keepdims=True)
4453 array([[7.],
4454 [2.]])
4455 >>> m = np.quantile(a, 0.5, axis=0)
4456 >>> out = np.zeros_like(m)
4457 >>> np.quantile(a, 0.5, axis=0, out=out)
4458 array([6.5, 4.5, 2.5])
4459 >>> m
4460 array([6.5, 4.5, 2.5])
4461 >>> b = a.copy()
4462 >>> np.quantile(b, 0.5, axis=1, overwrite_input=True)
4463 array([7., 2.])
4464 >>> assert not np.all(a == b)
4466 See also `numpy.percentile` for a visualization of most methods.
4468 References
4469 ----------
4470 .. [1] R. J. Hyndman and Y. Fan,
4471 "Sample quantiles in statistical packages,"
4472 The American Statistician, 50(4), pp. 361-365, 1996
4474 """
4475 a = np.asanyarray(a)
4476 if a.dtype.kind == "c":
4477 raise TypeError("a must be an array of real numbers")
4479 weak_q = type(q) in (int, float) # use weak promotion for final result type
4480 q = np.asanyarray(q)
4482 if not _quantile_is_valid(q):
4483 raise ValueError("Quantiles must be in the range [0, 1]")
4485 if weights is not None:
4486 if method != "inverted_cdf":
4487 msg = ("Only method 'inverted_cdf' supports weights. "
4488 f"Got: {method}.")
4489 raise ValueError(msg)
4490 if axis is not None:
4491 axis = _nx.normalize_axis_tuple(axis, a.ndim, argname="axis")
4492 weights = _weights_are_valid(weights=weights, a=a, axis=axis)
4493 if np.any(weights < 0):
4494 raise ValueError("Weights must be non-negative.")
4496 return _quantile_unchecked(
4497 a, q, axis, out, overwrite_input, method, keepdims, weights, weak_q)
4500def _quantile_unchecked(a,
4501 q,
4502 axis=None,
4503 out=None,
4504 overwrite_input=False,
4505 method="linear",
4506 keepdims=False,
4507 weights=None,
4508 weak_q=False):
4509 """Assumes that q is in [0, 1], and is an ndarray"""
4510 return _ureduce(a,
4511 func=_quantile_ureduce_func,
4512 q=q,
4513 weights=weights,
4514 keepdims=keepdims,
4515 axis=axis,
4516 out=out,
4517 overwrite_input=overwrite_input,
4518 method=method,
4519 weak_q=weak_q)
4522def _quantile_is_valid(q):
4523 # avoid expensive reductions, relevant for arrays with < O(1000) elements
4524 if q.ndim == 1 and q.size < 10:
4525 for i in range(q.size):
4526 if not (0.0 <= q[i] <= 1.0):
4527 return False
4528 elif not (q.min() >= 0 and q.max() <= 1):
4529 return False
4530 return True
4533def _compute_virtual_index(n, quantiles, alpha: float, beta: float):
4534 """
4535 Compute the floating point indexes of an array for the linear
4536 interpolation of quantiles.
4537 n : array_like
4538 The sample sizes.
4539 quantiles : array_like
4540 The quantiles values.
4541 alpha : float
4542 A constant used to correct the index computed.
4543 beta : float
4544 A constant used to correct the index computed.
4546 alpha and beta values depend on the chosen method
4547 (see quantile documentation)
4549 Reference:
4550 Hyndman&Fan paper "Sample Quantiles in Statistical Packages",
4551 DOI: 10.1080/00031305.1996.10473566
4552 """
4553 return n * quantiles + (
4554 alpha + quantiles * (1 - alpha - beta)
4555 ) - 1
4558def _get_gamma(virtual_indexes, previous_indexes, method):
4559 """
4560 Compute gamma (a.k.a 'm' or 'weight') for the linear interpolation
4561 of quantiles.
4563 virtual_indexes : array_like
4564 The indexes where the percentile is supposed to be found in the sorted
4565 sample.
4566 previous_indexes : array_like
4567 The floor values of virtual_indexes.
4568 method : dict
4569 The interpolation method chosen, which may have a specific rule
4570 modifying gamma.
4572 gamma is usually the fractional part of virtual_indexes but can be modified
4573 by the interpolation method.
4574 """
4575 gamma = np.asanyarray(virtual_indexes - previous_indexes)
4576 gamma = method["fix_gamma"](gamma, virtual_indexes)
4577 # Ensure both that we have an array, and that we keep the dtype
4578 # (which may have been matched to the input array).
4579 return np.asanyarray(gamma, dtype=virtual_indexes.dtype)
4582def _lerp(a, b, t, out=None):
4583 """
4584 Compute the linear interpolation weighted by gamma on each point of
4585 two same shape array.
4587 a : array_like
4588 Left bound.
4589 b : array_like
4590 Right bound.
4591 t : array_like
4592 The interpolation weight.
4593 out : array_like
4594 Output array.
4595 """
4596 diff_b_a = b - a
4597 lerp_interpolation = add(a, diff_b_a * t, out=... if out is None else out)
4598 subtract(b, diff_b_a * (1 - t), out=lerp_interpolation, where=t >= 0.5,
4599 casting='unsafe', dtype=type(lerp_interpolation.dtype))
4600 if lerp_interpolation.ndim == 0 and out is None:
4601 lerp_interpolation = lerp_interpolation[()] # unpack 0d arrays
4602 return lerp_interpolation
4605def _get_gamma_mask(shape, default_value, conditioned_value, where):
4606 out = np.full(shape, default_value)
4607 np.copyto(out, conditioned_value, where=where, casting="unsafe")
4608 return out
4611def _discrete_interpolation_to_boundaries(index, gamma_condition_fun):
4612 previous = np.floor(index)
4613 next = previous + 1
4614 gamma = index - previous
4615 res = _get_gamma_mask(shape=index.shape,
4616 default_value=next,
4617 conditioned_value=previous,
4618 where=gamma_condition_fun(gamma, index)
4619 ).astype(np.intp)
4620 # Some methods can lead to out-of-bound integers, clip them:
4621 res[res < 0] = 0
4622 return res
4625def _closest_observation(n, quantiles):
4626 # "choose the nearest even order statistic at g=0" (H&F (1996) pp. 362).
4627 # Order is 1-based so for zero-based indexing round to nearest odd index.
4628 gamma_fun = lambda gamma, index: (gamma == 0) & (np.floor(index) % 2 == 1)
4629 return _discrete_interpolation_to_boundaries((n * quantiles) - 1 - 0.5,
4630 gamma_fun)
4633def _inverted_cdf(n, quantiles):
4634 gamma_fun = lambda gamma, _: (gamma == 0)
4635 return _discrete_interpolation_to_boundaries((n * quantiles) - 1,
4636 gamma_fun)
4639def _quantile_ureduce_func(
4640 a: np.ndarray,
4641 q: np.ndarray,
4642 weights: np.ndarray | None,
4643 axis: int | None = None,
4644 out: np.ndarray | None = None,
4645 overwrite_input: bool = False,
4646 method: str = "linear",
4647 weak_q: bool = False,
4648) -> np.ndarray:
4649 if q.ndim > 2:
4650 # The code below works fine for nd, but it might not have useful
4651 # semantics. For now, keep the supported dimensions the same as it was
4652 # before.
4653 raise ValueError("q must be a scalar or 1d")
4654 if overwrite_input:
4655 if axis is None:
4656 axis = 0
4657 arr = a.ravel()
4658 wgt = None if weights is None else weights.ravel()
4659 else:
4660 arr = a
4661 wgt = weights
4662 elif axis is None:
4663 axis = 0
4664 arr = a.flatten()
4665 wgt = None if weights is None else weights.flatten()
4666 else:
4667 arr = a.copy()
4668 wgt = weights
4669 result = _quantile(arr,
4670 quantiles=q,
4671 axis=axis,
4672 method=method,
4673 out=out,
4674 weights=wgt,
4675 weak_q=weak_q)
4676 return result
4679def _get_indexes(arr, virtual_indexes, valid_values_count):
4680 """
4681 Get the valid indexes of arr neighbouring virtual_indexes.
4682 Note
4683 This is a companion function to linear interpolation of
4684 Quantiles
4686 Returns
4687 -------
4688 (previous_indexes, next_indexes): Tuple
4689 A Tuple of virtual_indexes neighbouring indexes
4690 """
4691 previous_indexes = floor(virtual_indexes, out=...)
4692 next_indexes = add(previous_indexes, 1, out=...)
4693 indexes_above_bounds = virtual_indexes >= valid_values_count - 1
4694 # When indexes is above max index, take the max value of the array
4695 if indexes_above_bounds.any():
4696 previous_indexes[indexes_above_bounds] = -1
4697 next_indexes[indexes_above_bounds] = -1
4698 # When indexes is below min index, take the min value of the array
4699 indexes_below_bounds = virtual_indexes < 0
4700 if indexes_below_bounds.any():
4701 previous_indexes[indexes_below_bounds] = 0
4702 next_indexes[indexes_below_bounds] = 0
4703 if np.issubdtype(arr.dtype, np.inexact):
4704 # After the sort, slices having NaNs will have for last element a NaN
4705 virtual_indexes_nans = np.isnan(virtual_indexes)
4706 if virtual_indexes_nans.any():
4707 previous_indexes[virtual_indexes_nans] = -1
4708 next_indexes[virtual_indexes_nans] = -1
4709 previous_indexes = previous_indexes.astype(np.intp)
4710 next_indexes = next_indexes.astype(np.intp)
4711 return previous_indexes, next_indexes
4714def _quantile(
4715 arr: "np.typing.ArrayLike",
4716 quantiles: np.ndarray,
4717 axis: int = -1,
4718 method: str = "linear",
4719 out: np.ndarray | None = None,
4720 weights: "np.typing.ArrayLike | None" = None,
4721 weak_q: bool = False,
4722) -> np.ndarray:
4723 """
4724 Private function that doesn't support extended axis or keepdims.
4725 These methods are extended to this function using _ureduce
4726 See nanpercentile for parameter usage
4727 It computes the quantiles of the array for the given axis.
4728 A linear interpolation is performed based on the `method`.
4730 By default, the method is "linear" where alpha == beta == 1 which
4731 performs the 7th method of Hyndman&Fan.
4732 With "median_unbiased" we get alpha == beta == 1/3
4733 thus the 8th method of Hyndman&Fan.
4734 """
4735 # --- Setup
4736 arr = np.asanyarray(arr)
4737 values_count = arr.shape[axis]
4738 # The dimensions of `q` are prepended to the output shape, so we need the
4739 # axis being sampled from `arr` to be last.
4740 if axis != 0: # But moveaxis is slow, so only call it if necessary.
4741 arr = np.moveaxis(arr, axis, destination=0)
4742 supports_nans = (
4743 np.issubdtype(arr.dtype, np.inexact) or arr.dtype.kind in 'Mm'
4744 )
4746 if weights is None:
4747 # --- Computation of indexes
4748 # Index where to find the value in the sorted array.
4749 # Virtual because it is a floating point value, not a valid index.
4750 # The nearest neighbours are used for interpolation
4751 try:
4752 method_props = _QuantileMethods[method]
4753 except KeyError:
4754 raise ValueError(
4755 f"{method!r} is not a valid method. Use one of: "
4756 f"{_QuantileMethods.keys()}") from None
4757 virtual_indexes = method_props["get_virtual_index"](values_count,
4758 quantiles)
4759 virtual_indexes = np.asanyarray(virtual_indexes)
4761 if method_props["fix_gamma"] is None:
4762 supports_integers = True
4763 else:
4764 int_virtual_indices = np.issubdtype(virtual_indexes.dtype,
4765 np.integer)
4766 supports_integers = method == 'linear' and int_virtual_indices
4768 if supports_integers:
4769 # No interpolation needed, take the points along axis
4770 if supports_nans:
4771 # may contain nan, which would sort to the end
4772 arr.partition(
4773 concatenate((virtual_indexes.ravel(), [-1])), axis=0,
4774 )
4775 slices_having_nans = np.isnan(arr[-1, ...])
4776 else:
4777 # cannot contain nan
4778 arr.partition(virtual_indexes.ravel(), axis=0)
4779 slices_having_nans = np.array(False, dtype=bool)
4780 result = take(arr, virtual_indexes, axis=0, out=out)
4781 else:
4782 previous_indexes, next_indexes = _get_indexes(arr,
4783 virtual_indexes,
4784 values_count)
4785 # --- Sorting
4786 arr.partition(
4787 np.unique(np.concatenate(([0, -1],
4788 previous_indexes.ravel(),
4789 next_indexes.ravel(),
4790 ))),
4791 axis=0)
4792 if supports_nans:
4793 slices_having_nans = np.isnan(arr[-1, ...])
4794 else:
4795 slices_having_nans = None
4796 # --- Get values from indexes
4797 previous = arr[previous_indexes]
4798 next = arr[next_indexes]
4799 # --- Linear interpolation
4800 gamma = _get_gamma(virtual_indexes, previous_indexes,
4801 method_props)
4802 if weak_q:
4803 gamma = float(gamma)
4804 else:
4805 result_shape = virtual_indexes.shape + (1,) * (arr.ndim - 1)
4806 gamma = gamma.reshape(result_shape)
4807 result = _lerp(previous,
4808 next,
4809 gamma,
4810 out=out)
4811 else:
4812 # Weighted case
4813 # This implements method="inverted_cdf", the only supported weighted
4814 # method, which needs to sort anyway.
4815 weights = np.asanyarray(weights)
4816 if axis != 0:
4817 weights = np.moveaxis(weights, axis, destination=0)
4818 index_array = np.argsort(arr, axis=0)
4820 # arr = arr[index_array, ...] # but this adds trailing dimensions of
4821 # 1.
4822 arr = np.take_along_axis(arr, index_array, axis=0)
4823 if weights.shape == arr.shape:
4824 weights = np.take_along_axis(weights, index_array, axis=0)
4825 else:
4826 # weights is 1d
4827 weights = weights.reshape(-1)[index_array, ...]
4829 if supports_nans:
4830 # may contain nan, which would sort to the end
4831 slices_having_nans = np.isnan(arr[-1, ...])
4832 else:
4833 # cannot contain nan
4834 slices_having_nans = np.array(False, dtype=bool)
4836 # We use the weights to calculate the empirical cumulative
4837 # distribution function cdf
4838 cdf = weights.cumsum(axis=0, dtype=np.float64)
4839 cdf /= cdf[-1, ...] # normalization to 1
4840 if np.isnan(cdf[-1]).any():
4841 # Above calculations should normally warn for the zero/inf case.
4842 raise ValueError("Weights included NaN, inf or were all zero.")
4843 # Search index i such that
4844 # sum(weights[j], j=0..i-1) < quantile <= sum(weights[j], j=0..i)
4845 # is then equivalent to
4846 # cdf[i-1] < quantile <= cdf[i]
4847 # Unfortunately, searchsorted only accepts 1-d arrays as first
4848 # argument, so we will need to iterate over dimensions.
4850 # Without the following cast, searchsorted can return surprising
4851 # results, e.g.
4852 # np.searchsorted(np.array([0.2, 0.4, 0.6, 0.8, 1.]),
4853 # np.array(0.4, dtype=np.float32), side="left")
4854 # returns 2 instead of 1 because 0.4 is not binary representable.
4855 if quantiles.dtype.kind == "f":
4856 cdf = cdf.astype(quantiles.dtype)
4857 # Weights must be non-negative, so we might have zero weights at the
4858 # beginning leading to some leading zeros in cdf. The call to
4859 # np.searchsorted for quantiles=0 will then pick the first element,
4860 # but should pick the first one larger than zero. We
4861 # therefore simply set 0 values in cdf to -1.
4862 if np.any(cdf[0, ...] == 0):
4863 cdf[cdf == 0] = -1
4865 def find_cdf_1d(arr, cdf):
4866 indices = np.searchsorted(cdf, quantiles, side="left")
4867 # We might have reached the maximum with i = len(arr), e.g. for
4868 # quantiles = 1, and need to cut it to len(arr) - 1.
4869 indices = minimum(indices, values_count - 1)
4870 result = take(arr, indices, axis=0)
4871 return result
4873 r_shape = arr.shape[1:]
4874 if quantiles.ndim > 0:
4875 r_shape = quantiles.shape + r_shape
4876 if out is None:
4877 result = np.empty_like(arr, shape=r_shape)
4878 else:
4879 if out.shape != r_shape:
4880 msg = (f"Wrong shape of argument 'out', shape={r_shape} is "
4881 f"required; got shape={out.shape}.")
4882 raise ValueError(msg)
4883 result = out
4885 # See apply_along_axis, which we do for axis=0. Note that Ni = (,)
4886 # always, so we remove it here.
4887 Nk = arr.shape[1:]
4888 for kk in np.ndindex(Nk):
4889 result[(...,) + kk] = find_cdf_1d(
4890 arr[np.s_[:, ] + kk], cdf[np.s_[:, ] + kk]
4891 )
4893 # Make result the same as in unweighted inverted_cdf.
4894 if result.shape == () and result.dtype == np.dtype("O"):
4895 result = result.item()
4897 if np.any(slices_having_nans):
4898 if result.ndim == 0 and out is None:
4899 # can't write to a scalar, but indexing will be correct
4900 result = arr[-1]
4901 else:
4902 np.copyto(result, arr[-1, ...], where=slices_having_nans)
4903 return result
4906def _trapezoid_dispatcher(y, x=None, dx=None, axis=None):
4907 return (y, x)
4910@array_function_dispatch(_trapezoid_dispatcher)
4911def trapezoid(y, x=None, dx=1.0, axis=-1):
4912 r"""
4913 Integrate along the given axis using the composite trapezoidal rule.
4915 If `x` is provided, the integration happens in sequence along its
4916 elements - they are not sorted.
4918 Integrate `y` (`x`) along each 1d slice on the given axis, compute
4919 :math:`\int y(x) dx`.
4920 When `x` is specified, this integrates along the parametric curve,
4921 computing :math:`\int_t y(t) dt =
4922 \int_t y(t) \left.\frac{dx}{dt}\right|_{x=x(t)} dt`.
4924 .. versionadded:: 2.0.0
4926 Parameters
4927 ----------
4928 y : array_like
4929 Input array to integrate.
4930 x : array_like, optional
4931 The sample points corresponding to the `y` values. If `x` is None,
4932 the sample points are assumed to be evenly spaced `dx` apart. The
4933 default is None.
4934 dx : scalar, optional
4935 The spacing between sample points when `x` is None. The default is 1.
4936 axis : int, optional
4937 The axis along which to integrate.
4939 Returns
4940 -------
4941 trapezoid : float or ndarray
4942 Definite integral of `y` = n-dimensional array as approximated along
4943 a single axis by the trapezoidal rule. If `y` is a 1-dimensional array,
4944 then the result is a float. If `n` is greater than 1, then the result
4945 is an `n`-1 dimensional array.
4947 See Also
4948 --------
4949 sum, cumsum
4951 Notes
4952 -----
4953 Image [2]_ illustrates trapezoidal rule -- y-axis locations of points
4954 will be taken from `y` array, by default x-axis distances between
4955 points will be 1.0, alternatively they can be provided with `x` array
4956 or with `dx` scalar. Return value will be equal to combined area under
4957 the red lines.
4960 References
4961 ----------
4962 .. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule
4964 .. [2] Illustration image:
4965 https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
4967 Examples
4968 --------
4969 >>> import numpy as np
4971 Use the trapezoidal rule on evenly spaced points:
4973 >>> np.trapezoid([1, 2, 3])
4974 4.0
4976 The spacing between sample points can be selected by either the
4977 ``x`` or ``dx`` arguments:
4979 >>> np.trapezoid([1, 2, 3], x=[4, 6, 8])
4980 8.0
4981 >>> np.trapezoid([1, 2, 3], dx=2)
4982 8.0
4984 Using a decreasing ``x`` corresponds to integrating in reverse:
4986 >>> np.trapezoid([1, 2, 3], x=[8, 6, 4])
4987 -8.0
4989 More generally ``x`` is used to integrate along a parametric curve. We can
4990 estimate the integral :math:`\int_0^1 x^2 = 1/3` using:
4992 >>> x = np.linspace(0, 1, num=50)
4993 >>> y = x**2
4994 >>> np.trapezoid(y, x)
4995 0.33340274885464394
4997 Or estimate the area of a circle, noting we repeat the sample which closes
4998 the curve:
5000 >>> theta = np.linspace(0, 2 * np.pi, num=1000, endpoint=True)
5001 >>> np.trapezoid(np.cos(theta), x=np.sin(theta))
5002 3.141571941375841
5004 ``np.trapezoid`` can be applied along a specified axis to do multiple
5005 computations in one call:
5007 >>> a = np.arange(6).reshape(2, 3)
5008 >>> a
5009 array([[0, 1, 2],
5010 [3, 4, 5]])
5011 >>> np.trapezoid(a, axis=0)
5012 array([1.5, 2.5, 3.5])
5013 >>> np.trapezoid(a, axis=1)
5014 array([2., 8.])
5015 """
5017 y = asanyarray(y)
5018 if x is None:
5019 d = dx
5020 else:
5021 x = asanyarray(x)
5022 if x.ndim == 1:
5023 d = diff(x)
5024 # reshape to correct shape
5025 shape = [1] * y.ndim
5026 shape[axis] = d.shape[0]
5027 d = d.reshape(shape)
5028 else:
5029 d = diff(x, axis=axis)
5030 nd = y.ndim
5031 slice1 = [slice(None)] * nd
5032 slice2 = [slice(None)] * nd
5033 slice1[axis] = slice(1, None)
5034 slice2[axis] = slice(None, -1)
5035 try:
5036 ret = (d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0).sum(axis)
5037 except ValueError:
5038 # Operations didn't work, cast to ndarray
5039 d = np.asarray(d)
5040 y = np.asarray(y)
5041 ret = add.reduce(d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0, axis)
5042 return ret
5045def _meshgrid_dispatcher(*xi, copy=None, sparse=None, indexing=None):
5046 return xi
5049# Based on scitools meshgrid
5050@array_function_dispatch(_meshgrid_dispatcher)
5051def meshgrid(*xi, copy=True, sparse=False, indexing='xy'):
5052 """
5053 Return a tuple of coordinate matrices from coordinate vectors.
5055 Make N-D coordinate arrays for vectorized evaluations of
5056 N-D scalar/vector fields over N-D grids, given
5057 one-dimensional coordinate arrays x1, x2,..., xn.
5059 Parameters
5060 ----------
5061 x1, x2,..., xn : array_like
5062 1-D arrays representing the coordinates of a grid.
5063 indexing : {'xy', 'ij'}, optional
5064 Cartesian ('xy', default) or matrix ('ij') indexing of output.
5065 See Notes for more details.
5066 sparse : bool, optional
5067 If True the shape of the returned coordinate array for dimension *i*
5068 is reduced from ``(N1, ..., Ni, ... Nn)`` to
5069 ``(1, ..., 1, Ni, 1, ..., 1)``. These sparse coordinate grids are
5070 intended to be used with :ref:`basics.broadcasting`. When all
5071 coordinates are used in an expression, broadcasting still leads to a
5072 fully-dimensonal result array.
5074 Default is False.
5076 copy : bool, optional
5077 If False, a view into the original arrays are returned in order to
5078 conserve memory. Default is True. Please note that
5079 ``sparse=False, copy=False`` will likely return non-contiguous
5080 arrays. Furthermore, more than one element of a broadcast array
5081 may refer to a single memory location. If you need to write to the
5082 arrays, make copies first.
5084 Returns
5085 -------
5086 X1, X2,..., XN : tuple of ndarrays
5087 For vectors `x1`, `x2`,..., `xn` with lengths ``Ni=len(xi)``,
5088 returns ``(N1, N2, N3,..., Nn)`` shaped arrays if indexing='ij'
5089 or ``(N2, N1, N3,..., Nn)`` shaped arrays if indexing='xy'
5090 with the elements of `xi` repeated to fill the matrix along
5091 the first dimension for `x1`, the second for `x2` and so on.
5093 Notes
5094 -----
5095 This function supports both indexing conventions through the indexing
5096 keyword argument. Giving the string 'ij' returns a meshgrid with
5097 matrix indexing, while 'xy' returns a meshgrid with Cartesian indexing.
5098 In the 2-D case with inputs of length M and N, the outputs are of shape
5099 (N, M) for 'xy' indexing and (M, N) for 'ij' indexing. In the 3-D case
5100 with inputs of length M, N and P, outputs are of shape (N, M, P) for
5101 'xy' indexing and (M, N, P) for 'ij' indexing. The difference is
5102 illustrated by the following code snippet::
5104 xv, yv = np.meshgrid(x, y, indexing='ij')
5105 for i in range(nx):
5106 for j in range(ny):
5107 # treat xv[i,j], yv[i,j]
5109 xv, yv = np.meshgrid(x, y, indexing='xy')
5110 for i in range(nx):
5111 for j in range(ny):
5112 # treat xv[j,i], yv[j,i]
5114 In the 1-D and 0-D case, the indexing and sparse keywords have no effect.
5116 See Also
5117 --------
5118 mgrid : Construct a multi-dimensional "meshgrid" using indexing notation.
5119 ogrid : Construct an open multi-dimensional "meshgrid" using indexing
5120 notation.
5121 :ref:`how-to-index`
5123 Examples
5124 --------
5125 >>> import numpy as np
5126 >>> nx, ny = (3, 2)
5127 >>> x = np.linspace(0, 1, nx)
5128 >>> y = np.linspace(0, 1, ny)
5129 >>> xv, yv = np.meshgrid(x, y)
5130 >>> xv
5131 array([[0. , 0.5, 1. ],
5132 [0. , 0.5, 1. ]])
5133 >>> yv
5134 array([[0., 0., 0.],
5135 [1., 1., 1.]])
5137 The result of `meshgrid` is a coordinate grid:
5139 >>> import matplotlib.pyplot as plt
5140 >>> plt.plot(xv, yv, marker='o', color='k', linestyle='none')
5141 >>> plt.show()
5143 You can create sparse output arrays to save memory and computation time.
5145 >>> xv, yv = np.meshgrid(x, y, sparse=True)
5146 >>> xv
5147 array([[0. , 0.5, 1. ]])
5148 >>> yv
5149 array([[0.],
5150 [1.]])
5152 `meshgrid` is very useful to evaluate functions on a grid. If the
5153 function depends on all coordinates, both dense and sparse outputs can be
5154 used.
5156 >>> x = np.linspace(-5, 5, 101)
5157 >>> y = np.linspace(-5, 5, 101)
5158 >>> # full coordinate arrays
5159 >>> xx, yy = np.meshgrid(x, y)
5160 >>> zz = np.sqrt(xx**2 + yy**2)
5161 >>> xx.shape, yy.shape, zz.shape
5162 ((101, 101), (101, 101), (101, 101))
5163 >>> # sparse coordinate arrays
5164 >>> xs, ys = np.meshgrid(x, y, sparse=True)
5165 >>> zs = np.sqrt(xs**2 + ys**2)
5166 >>> xs.shape, ys.shape, zs.shape
5167 ((1, 101), (101, 1), (101, 101))
5168 >>> np.array_equal(zz, zs)
5169 True
5171 >>> h = plt.contourf(x, y, zs)
5172 >>> plt.axis('scaled')
5173 >>> plt.colorbar()
5174 >>> plt.show()
5175 """
5176 ndim = len(xi)
5178 if indexing not in ['xy', 'ij']:
5179 raise ValueError(
5180 "Valid values for `indexing` are 'xy' and 'ij'.")
5182 s0 = (1,) * ndim
5183 output = [np.asanyarray(x).reshape(s0[:i] + (-1,) + s0[i + 1:])
5184 for i, x in enumerate(xi)]
5186 if indexing == 'xy' and ndim > 1:
5187 # switch first and second axis
5188 output[0].shape = (1, -1) + s0[2:]
5189 output[1].shape = (-1, 1) + s0[2:]
5191 if not sparse:
5192 # Return the full N-D matrix (not only the 1-D vector)
5193 output = np.broadcast_arrays(*output, subok=True)
5195 if copy:
5196 output = tuple(x.copy() for x in output)
5198 return output
5201def _delete_dispatcher(arr, obj, axis=None):
5202 return (arr, obj)
5205@array_function_dispatch(_delete_dispatcher)
5206def delete(arr, obj, axis=None):
5207 """
5208 Return a new array with sub-arrays along an axis deleted. For a one
5209 dimensional array, this returns those entries not returned by
5210 `arr[obj]`.
5212 Parameters
5213 ----------
5214 arr : array_like
5215 Input array.
5216 obj : slice, int, array-like of ints or bools
5217 Indicate indices of sub-arrays to remove along the specified axis.
5219 .. versionchanged:: 1.19.0
5220 Boolean indices are now treated as a mask of elements to remove,
5221 rather than being cast to the integers 0 and 1.
5223 axis : int, optional
5224 The axis along which to delete the subarray defined by `obj`.
5225 If `axis` is None, `obj` is applied to the flattened array.
5227 Returns
5228 -------
5229 out : ndarray
5230 A copy of `arr` with the elements specified by `obj` removed. Note
5231 that `delete` does not occur in-place. If `axis` is None, `out` is
5232 a flattened array.
5234 See Also
5235 --------
5236 insert : Insert elements into an array.
5237 append : Append elements at the end of an array.
5239 Notes
5240 -----
5241 Often it is preferable to use a boolean mask. For example:
5243 >>> arr = np.arange(12) + 1
5244 >>> mask = np.ones(len(arr), dtype=bool)
5245 >>> mask[[0,2,4]] = False
5246 >>> result = arr[mask,...]
5248 Is equivalent to ``np.delete(arr, [0,2,4], axis=0)``, but allows further
5249 use of `mask`.
5251 Examples
5252 --------
5253 >>> import numpy as np
5254 >>> arr = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]])
5255 >>> arr
5256 array([[ 1, 2, 3, 4],
5257 [ 5, 6, 7, 8],
5258 [ 9, 10, 11, 12]])
5259 >>> np.delete(arr, 1, 0)
5260 array([[ 1, 2, 3, 4],
5261 [ 9, 10, 11, 12]])
5263 >>> np.delete(arr, np.s_[::2], 1)
5264 array([[ 2, 4],
5265 [ 6, 8],
5266 [10, 12]])
5267 >>> np.delete(arr, [1,3,5], None)
5268 array([ 1, 3, 5, 7, 8, 9, 10, 11, 12])
5270 """
5271 conv = _array_converter(arr)
5272 arr, = conv.as_arrays(subok=False)
5274 ndim = arr.ndim
5275 arrorder = 'F' if arr.flags.fnc else 'C'
5276 if axis is None:
5277 if ndim != 1:
5278 arr = arr.ravel()
5279 # needed for np.matrix, which is still not 1d after being ravelled
5280 ndim = arr.ndim
5281 axis = ndim - 1
5282 else:
5283 axis = normalize_axis_index(axis, ndim)
5285 slobj = [slice(None)] * ndim
5286 N = arr.shape[axis]
5287 newshape = list(arr.shape)
5289 if isinstance(obj, slice):
5290 start, stop, step = obj.indices(N)
5291 xr = range(start, stop, step)
5292 numtodel = len(xr)
5294 if numtodel <= 0:
5295 return conv.wrap(arr.copy(order=arrorder), to_scalar=False)
5297 # Invert if step is negative:
5298 if step < 0:
5299 step = -step
5300 start = xr[-1]
5301 stop = xr[0] + 1
5303 newshape[axis] -= numtodel
5304 new = empty(newshape, arr.dtype, arrorder)
5305 # copy initial chunk
5306 if start == 0:
5307 pass
5308 else:
5309 slobj[axis] = slice(None, start)
5310 new[tuple(slobj)] = arr[tuple(slobj)]
5311 # copy end chunk
5312 if stop == N:
5313 pass
5314 else:
5315 slobj[axis] = slice(stop - numtodel, None)
5316 slobj2 = [slice(None)] * ndim
5317 slobj2[axis] = slice(stop, None)
5318 new[tuple(slobj)] = arr[tuple(slobj2)]
5319 # copy middle pieces
5320 if step == 1:
5321 pass
5322 else: # use array indexing.
5323 keep = ones(stop - start, dtype=bool)
5324 keep[:stop - start:step] = False
5325 slobj[axis] = slice(start, stop - numtodel)
5326 slobj2 = [slice(None)] * ndim
5327 slobj2[axis] = slice(start, stop)
5328 arr = arr[tuple(slobj2)]
5329 slobj2[axis] = keep
5330 new[tuple(slobj)] = arr[tuple(slobj2)]
5332 return conv.wrap(new, to_scalar=False)
5334 if isinstance(obj, (int, integer)) and not isinstance(obj, bool):
5335 single_value = True
5336 else:
5337 single_value = False
5338 _obj = obj
5339 obj = np.asarray(obj)
5340 # `size == 0` to allow empty lists similar to indexing, but (as there)
5341 # is really too generic:
5342 if obj.size == 0 and not isinstance(_obj, np.ndarray):
5343 obj = obj.astype(intp)
5344 elif obj.size == 1 and obj.dtype.kind in "ui":
5345 # For a size 1 integer array we can use the single-value path
5346 # (most dtypes, except boolean, should just fail later).
5347 obj = obj.item()
5348 single_value = True
5350 if single_value:
5351 # optimization for a single value
5352 if (obj < -N or obj >= N):
5353 raise IndexError(
5354 f"index {obj} is out of bounds for axis {axis} with "
5355 f"size {N}")
5356 if (obj < 0):
5357 obj += N
5358 newshape[axis] -= 1
5359 new = empty(newshape, arr.dtype, arrorder)
5360 slobj[axis] = slice(None, obj)
5361 new[tuple(slobj)] = arr[tuple(slobj)]
5362 slobj[axis] = slice(obj, None)
5363 slobj2 = [slice(None)] * ndim
5364 slobj2[axis] = slice(obj + 1, None)
5365 new[tuple(slobj)] = arr[tuple(slobj2)]
5366 else:
5367 if obj.dtype == bool:
5368 if obj.shape != (N,):
5369 raise ValueError('boolean array argument obj to delete '
5370 'must be one dimensional and match the axis '
5371 f'length of {N}')
5373 # optimization, the other branch is slower
5374 keep = ~obj
5375 else:
5376 keep = ones(N, dtype=bool)
5377 keep[obj,] = False
5379 slobj[axis] = keep
5380 new = arr[tuple(slobj)]
5382 return conv.wrap(new, to_scalar=False)
5385def _insert_dispatcher(arr, obj, values, axis=None):
5386 return (arr, obj, values)
5389@array_function_dispatch(_insert_dispatcher)
5390def insert(arr, obj, values, axis=None):
5391 """
5392 Insert values along the given axis before the given indices.
5394 Parameters
5395 ----------
5396 arr : array_like
5397 Input array.
5398 obj : slice, int, array-like of ints or bools
5399 Object that defines the index or indices before which `values` is
5400 inserted.
5402 .. versionchanged:: 2.1.2
5403 Boolean indices are now treated as a mask of elements to insert,
5404 rather than being cast to the integers 0 and 1.
5406 Support for multiple insertions when `obj` is a single scalar or a
5407 sequence with one element (similar to calling insert multiple
5408 times).
5409 values : array_like
5410 Values to insert into `arr`. If the type of `values` is different
5411 from that of `arr`, `values` is converted to the type of `arr`.
5412 `values` should be shaped so that ``arr[...,obj,...] = values``
5413 is legal.
5414 axis : int, optional
5415 Axis along which to insert `values`. If `axis` is None then `arr`
5416 is flattened first.
5418 Returns
5419 -------
5420 out : ndarray
5421 A copy of `arr` with `values` inserted. Note that `insert`
5422 does not occur in-place: a new array is returned. If
5423 `axis` is None, `out` is a flattened array.
5425 See Also
5426 --------
5427 append : Append elements at the end of an array.
5428 concatenate : Join a sequence of arrays along an existing axis.
5429 delete : Delete elements from an array.
5431 Notes
5432 -----
5433 Note that for higher dimensional inserts ``obj=0`` behaves very different
5434 from ``obj=[0]`` just like ``arr[:,0,:] = values`` is different from
5435 ``arr[:,[0],:] = values``. This is because of the difference between basic
5436 and advanced :ref:`indexing <basics.indexing>`.
5438 Examples
5439 --------
5440 >>> import numpy as np
5441 >>> a = np.arange(6).reshape(3, 2)
5442 >>> a
5443 array([[0, 1],
5444 [2, 3],
5445 [4, 5]])
5446 >>> np.insert(a, 1, 6)
5447 array([0, 6, 1, 2, 3, 4, 5])
5448 >>> np.insert(a, 1, 6, axis=1)
5449 array([[0, 6, 1],
5450 [2, 6, 3],
5451 [4, 6, 5]])
5453 Difference between sequence and scalars,
5454 showing how ``obj=[1]`` behaves different from ``obj=1``:
5456 >>> np.insert(a, [1], [[7],[8],[9]], axis=1)
5457 array([[0, 7, 1],
5458 [2, 8, 3],
5459 [4, 9, 5]])
5460 >>> np.insert(a, 1, [[7],[8],[9]], axis=1)
5461 array([[0, 7, 8, 9, 1],
5462 [2, 7, 8, 9, 3],
5463 [4, 7, 8, 9, 5]])
5464 >>> np.array_equal(np.insert(a, 1, [7, 8, 9], axis=1),
5465 ... np.insert(a, [1], [[7],[8],[9]], axis=1))
5466 True
5468 >>> b = a.flatten()
5469 >>> b
5470 array([0, 1, 2, 3, 4, 5])
5471 >>> np.insert(b, [2, 2], [6, 7])
5472 array([0, 1, 6, 7, 2, 3, 4, 5])
5474 >>> np.insert(b, slice(2, 4), [7, 8])
5475 array([0, 1, 7, 2, 8, 3, 4, 5])
5477 >>> np.insert(b, [2, 2], [7.13, False]) # type casting
5478 array([0, 1, 7, 0, 2, 3, 4, 5])
5480 >>> x = np.arange(8).reshape(2, 4)
5481 >>> idx = (1, 3)
5482 >>> np.insert(x, idx, 999, axis=1)
5483 array([[ 0, 999, 1, 2, 999, 3],
5484 [ 4, 999, 5, 6, 999, 7]])
5486 """
5487 conv = _array_converter(arr)
5488 arr, = conv.as_arrays(subok=False)
5490 ndim = arr.ndim
5491 arrorder = 'F' if arr.flags.fnc else 'C'
5492 if axis is None:
5493 if ndim != 1:
5494 arr = arr.ravel()
5495 # needed for np.matrix, which is still not 1d after being ravelled
5496 ndim = arr.ndim
5497 axis = ndim - 1
5498 else:
5499 axis = normalize_axis_index(axis, ndim)
5500 slobj = [slice(None)] * ndim
5501 N = arr.shape[axis]
5502 newshape = list(arr.shape)
5504 if isinstance(obj, slice):
5505 # turn it into a range object
5506 indices = arange(*obj.indices(N), dtype=intp)
5507 else:
5508 # need to copy obj, because indices will be changed in-place
5509 indices = np.array(obj)
5510 if indices.dtype == bool:
5511 if obj.ndim != 1:
5512 raise ValueError('boolean array argument obj to insert '
5513 'must be one dimensional')
5514 indices = np.flatnonzero(obj)
5515 elif indices.ndim > 1:
5516 raise ValueError(
5517 "index array argument obj to insert must be one dimensional "
5518 "or scalar")
5519 if indices.size == 1:
5520 index = indices.item()
5521 if index < -N or index > N:
5522 raise IndexError(f"index {obj} is out of bounds for axis {axis} "
5523 f"with size {N}")
5524 if (index < 0):
5525 index += N
5527 # There are some object array corner cases here, but we cannot avoid
5528 # that:
5529 values = array(values, copy=None, ndmin=arr.ndim, dtype=arr.dtype)
5530 if indices.ndim == 0:
5531 # broadcasting is very different here, since a[:,0,:] = ... behaves
5532 # very different from a[:,[0],:] = ...! This changes values so that
5533 # it works likes the second case. (here a[:,0:1,:])
5534 values = np.moveaxis(values, 0, axis)
5535 numnew = values.shape[axis]
5536 newshape[axis] += numnew
5537 new = empty(newshape, arr.dtype, arrorder)
5538 slobj[axis] = slice(None, index)
5539 new[tuple(slobj)] = arr[tuple(slobj)]
5540 slobj[axis] = slice(index, index + numnew)
5541 new[tuple(slobj)] = values
5542 slobj[axis] = slice(index + numnew, None)
5543 slobj2 = [slice(None)] * ndim
5544 slobj2[axis] = slice(index, None)
5545 new[tuple(slobj)] = arr[tuple(slobj2)]
5547 return conv.wrap(new, to_scalar=False)
5549 elif indices.size == 0 and not isinstance(obj, np.ndarray):
5550 # Can safely cast the empty list to intp
5551 indices = indices.astype(intp)
5553 indices[indices < 0] += N
5555 numnew = len(indices)
5556 order = indices.argsort(kind='mergesort') # stable sort
5557 indices[order] += np.arange(numnew)
5559 newshape[axis] += numnew
5560 old_mask = ones(newshape[axis], dtype=bool)
5561 old_mask[indices] = False
5563 new = empty(newshape, arr.dtype, arrorder)
5564 slobj2 = [slice(None)] * ndim
5565 slobj[axis] = indices
5566 slobj2[axis] = old_mask
5567 new[tuple(slobj)] = values
5568 new[tuple(slobj2)] = arr
5570 return conv.wrap(new, to_scalar=False)
5573def _append_dispatcher(arr, values, axis=None):
5574 return (arr, values)
5577@array_function_dispatch(_append_dispatcher)
5578def append(arr, values, axis=None):
5579 """
5580 Append values to the end of an array.
5582 Parameters
5583 ----------
5584 arr : array_like
5585 Values are appended to a copy of this array.
5586 values : array_like
5587 These values are appended to a copy of `arr`. It must be of the
5588 correct shape (the same shape as `arr`, excluding `axis`). If
5589 `axis` is not specified, `values` can be any shape and will be
5590 flattened before use.
5591 axis : int, optional
5592 The axis along which `values` are appended. If `axis` is not
5593 given, both `arr` and `values` are flattened before use.
5595 Returns
5596 -------
5597 append : ndarray
5598 A copy of `arr` with `values` appended to `axis`. Note that
5599 `append` does not occur in-place: a new array is allocated and
5600 filled. If `axis` is None, `out` is a flattened array.
5602 See Also
5603 --------
5604 insert : Insert elements into an array.
5605 delete : Delete elements from an array.
5607 Examples
5608 --------
5609 >>> import numpy as np
5610 >>> np.append([1, 2, 3], [[4, 5, 6], [7, 8, 9]])
5611 array([1, 2, 3, ..., 7, 8, 9])
5613 When `axis` is specified, `values` must have the correct shape.
5615 >>> np.append([[1, 2, 3], [4, 5, 6]], [[7, 8, 9]], axis=0)
5616 array([[1, 2, 3],
5617 [4, 5, 6],
5618 [7, 8, 9]])
5620 >>> np.append([[1, 2, 3], [4, 5, 6]], [7, 8, 9], axis=0)
5621 Traceback (most recent call last):
5622 ...
5623 ValueError: all the input arrays must have same number of dimensions, but
5624 the array at index 0 has 2 dimension(s) and the array at index 1 has 1
5625 dimension(s)
5627 >>> a = np.array([1, 2], dtype=int)
5628 >>> c = np.append(a, [])
5629 >>> c
5630 array([1., 2.])
5631 >>> c.dtype
5632 float64
5634 Default dtype for empty ndarrays is `float64` thus making the output of dtype
5635 `float64` when appended with dtype `int64`
5637 """
5638 arr = asanyarray(arr)
5639 if axis is None:
5640 if arr.ndim != 1:
5641 arr = arr.ravel()
5642 values = ravel(values)
5643 axis = arr.ndim - 1
5644 return concatenate((arr, values), axis=axis)
5647def _digitize_dispatcher(x, bins, right=None):
5648 return (x, bins)
5651@array_function_dispatch(_digitize_dispatcher)
5652def digitize(x, bins, right=False):
5653 """
5654 Return the indices of the bins to which each value in input array belongs.
5656 ========= ============= ============================
5657 `right` order of bins returned index `i` satisfies
5658 ========= ============= ============================
5659 ``False`` increasing ``bins[i-1] <= x < bins[i]``
5660 ``True`` increasing ``bins[i-1] < x <= bins[i]``
5661 ``False`` decreasing ``bins[i-1] > x >= bins[i]``
5662 ``True`` decreasing ``bins[i-1] >= x > bins[i]``
5663 ========= ============= ============================
5665 If values in `x` are beyond the bounds of `bins`, 0 or ``len(bins)`` is
5666 returned as appropriate.
5668 Parameters
5669 ----------
5670 x : array_like
5671 Input array to be binned. Prior to NumPy 1.10.0, this array had to
5672 be 1-dimensional, but can now have any shape.
5673 bins : array_like
5674 Array of bins. It has to be 1-dimensional and monotonic.
5675 right : bool, optional
5676 Indicating whether the intervals include the right or the left bin
5677 edge. Default behavior is (right==False) indicating that the interval
5678 does not include the right edge. The left bin end is open in this
5679 case, i.e., bins[i-1] <= x < bins[i] is the default behavior for
5680 monotonically increasing bins.
5682 Returns
5683 -------
5684 indices : ndarray of ints
5685 Output array of indices, of same shape as `x`.
5687 Raises
5688 ------
5689 ValueError
5690 If `bins` is not monotonic.
5691 TypeError
5692 If the type of the input is complex.
5694 See Also
5695 --------
5696 bincount, histogram, unique, searchsorted
5698 Notes
5699 -----
5700 If values in `x` are such that they fall outside the bin range,
5701 attempting to index `bins` with the indices that `digitize` returns
5702 will result in an IndexError.
5704 .. versionadded:: 1.10.0
5706 `numpy.digitize` is implemented in terms of `numpy.searchsorted`.
5707 This means that a binary search is used to bin the values, which scales
5708 much better for larger number of bins than the previous linear search.
5709 It also removes the requirement for the input array to be 1-dimensional.
5711 For monotonically *increasing* `bins`, the following are equivalent::
5713 np.digitize(x, bins, right=True)
5714 np.searchsorted(bins, x, side='left')
5716 Note that as the order of the arguments are reversed, the side must be too.
5717 The `searchsorted` call is marginally faster, as it does not do any
5718 monotonicity checks. Perhaps more importantly, it supports all dtypes.
5720 Examples
5721 --------
5722 >>> import numpy as np
5723 >>> x = np.array([0.2, 6.4, 3.0, 1.6])
5724 >>> bins = np.array([0.0, 1.0, 2.5, 4.0, 10.0])
5725 >>> inds = np.digitize(x, bins)
5726 >>> inds
5727 array([1, 4, 3, 2])
5728 >>> for n in range(x.size):
5729 ... print(bins[inds[n]-1], "<=", x[n], "<", bins[inds[n]])
5730 ...
5731 0.0 <= 0.2 < 1.0
5732 4.0 <= 6.4 < 10.0
5733 2.5 <= 3.0 < 4.0
5734 1.0 <= 1.6 < 2.5
5736 >>> x = np.array([1.2, 10.0, 12.4, 15.5, 20.])
5737 >>> bins = np.array([0, 5, 10, 15, 20])
5738 >>> np.digitize(x,bins,right=True)
5739 array([1, 2, 3, 4, 4])
5740 >>> np.digitize(x,bins,right=False)
5741 array([1, 3, 3, 4, 5])
5742 """
5743 x = _nx.asarray(x)
5744 bins = _nx.asarray(bins)
5746 # here for compatibility, searchsorted below is happy to take this
5747 if np.issubdtype(x.dtype, _nx.complexfloating):
5748 raise TypeError("x may not be complex")
5750 mono = _monotonicity(bins)
5751 if mono == 0:
5752 raise ValueError("bins must be monotonically increasing or decreasing")
5754 # this is backwards because the arguments below are swapped
5755 side = 'left' if right else 'right'
5756 if mono == -1:
5757 # reverse the bins, and invert the results
5758 return len(bins) - _nx.searchsorted(bins[::-1], x, side=side)
5759 else:
5760 return _nx.searchsorted(bins, x, side=side)