Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/numpy/lib/_function_base_impl.py: 13%
1291 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-04-09 06:12 +0000
1import builtins
2import collections.abc
3import functools
4import re
5import sys
6import warnings
8import numpy as np
9import numpy._core.numeric as _nx
10from numpy._core import transpose, overrides
11from numpy._core.numeric import (
12 ones, zeros_like, arange, concatenate, array, asarray, asanyarray, empty,
13 ndarray, take, dot, where, intp, integer, isscalar, absolute
14 )
15from numpy._core.umath import (
16 pi, add, arctan2, frompyfunc, cos, less_equal, sqrt, sin,
17 mod, exp, not_equal, subtract, minimum
18 )
19from numpy._core.fromnumeric import (
20 ravel, nonzero, partition, mean, any, sum
21 )
22from numpy._core.numerictypes import typecodes
23from numpy.lib._twodim_base_impl import diag
24from numpy._core.multiarray import (
25 _place, bincount, normalize_axis_index, _monotonicity,
26 interp as compiled_interp, interp_complex as compiled_interp_complex
27 )
28from numpy._core._multiarray_umath import _array_converter
29from numpy._utils import set_module
31# needed in this module for compatibility
32from numpy.lib._histograms_impl import histogram, histogramdd # noqa: F401
35array_function_dispatch = functools.partial(
36 overrides.array_function_dispatch, module='numpy')
39__all__ = [
40 'select', 'piecewise', 'trim_zeros', 'copy', 'iterable', 'percentile',
41 'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'flip',
42 'rot90', 'extract', 'place', 'vectorize', 'asarray_chkfinite', 'average',
43 'bincount', 'digitize', 'cov', 'corrcoef',
44 'median', 'sinc', 'hamming', 'hanning', 'bartlett',
45 'blackman', 'kaiser', 'trapezoid', 'trapz', 'i0',
46 'meshgrid', 'delete', 'insert', 'append', 'interp',
47 'quantile'
48 ]
50# _QuantileMethods is a dictionary listing all the supported methods to
51# compute quantile/percentile.
52#
53# Below virtual_index refers to the index of the element where the percentile
54# would be found in the sorted sample.
55# When the sample contains exactly the percentile wanted, the virtual_index is
56# an integer to the index of this element.
57# When the percentile wanted is in between two elements, the virtual_index
58# is made of a integer part (a.k.a 'i' or 'left') and a fractional part
59# (a.k.a 'g' or 'gamma')
60#
61# Each method in _QuantileMethods has two properties
62# get_virtual_index : Callable
63# The function used to compute the virtual_index.
64# fix_gamma : Callable
65# A function used for discret methods to force the index to a specific value.
66_QuantileMethods = dict(
67 # --- HYNDMAN and FAN METHODS
68 # Discrete methods
69 inverted_cdf=dict(
70 get_virtual_index=lambda n, quantiles: _inverted_cdf(n, quantiles),
71 fix_gamma=None, # should never be called
72 ),
73 averaged_inverted_cdf=dict(
74 get_virtual_index=lambda n, quantiles: (n * quantiles) - 1,
75 fix_gamma=lambda gamma, _: _get_gamma_mask(
76 shape=gamma.shape,
77 default_value=1.,
78 conditioned_value=0.5,
79 where=gamma == 0),
80 ),
81 closest_observation=dict(
82 get_virtual_index=lambda n, quantiles: _closest_observation(n,
83 quantiles),
84 fix_gamma=None, # should never be called
85 ),
86 # Continuous methods
87 interpolated_inverted_cdf=dict(
88 get_virtual_index=lambda n, quantiles:
89 _compute_virtual_index(n, quantiles, 0, 1),
90 fix_gamma=lambda gamma, _: gamma,
91 ),
92 hazen=dict(
93 get_virtual_index=lambda n, quantiles:
94 _compute_virtual_index(n, quantiles, 0.5, 0.5),
95 fix_gamma=lambda gamma, _: gamma,
96 ),
97 weibull=dict(
98 get_virtual_index=lambda n, quantiles:
99 _compute_virtual_index(n, quantiles, 0, 0),
100 fix_gamma=lambda gamma, _: gamma,
101 ),
102 # Default method.
103 # To avoid some rounding issues, `(n-1) * quantiles` is preferred to
104 # `_compute_virtual_index(n, quantiles, 1, 1)`.
105 # They are mathematically equivalent.
106 linear=dict(
107 get_virtual_index=lambda n, quantiles: (n - 1) * quantiles,
108 fix_gamma=lambda gamma, _: gamma,
109 ),
110 median_unbiased=dict(
111 get_virtual_index=lambda n, quantiles:
112 _compute_virtual_index(n, quantiles, 1 / 3.0, 1 / 3.0),
113 fix_gamma=lambda gamma, _: gamma,
114 ),
115 normal_unbiased=dict(
116 get_virtual_index=lambda n, quantiles:
117 _compute_virtual_index(n, quantiles, 3 / 8.0, 3 / 8.0),
118 fix_gamma=lambda gamma, _: gamma,
119 ),
120 # --- OTHER METHODS
121 lower=dict(
122 get_virtual_index=lambda n, quantiles: np.floor(
123 (n - 1) * quantiles).astype(np.intp),
124 fix_gamma=None, # should never be called, index dtype is int
125 ),
126 higher=dict(
127 get_virtual_index=lambda n, quantiles: np.ceil(
128 (n - 1) * quantiles).astype(np.intp),
129 fix_gamma=None, # should never be called, index dtype is int
130 ),
131 midpoint=dict(
132 get_virtual_index=lambda n, quantiles: 0.5 * (
133 np.floor((n - 1) * quantiles)
134 + np.ceil((n - 1) * quantiles)),
135 fix_gamma=lambda gamma, index: _get_gamma_mask(
136 shape=gamma.shape,
137 default_value=0.5,
138 conditioned_value=0.,
139 where=index % 1 == 0),
140 ),
141 nearest=dict(
142 get_virtual_index=lambda n, quantiles: np.around(
143 (n - 1) * quantiles).astype(np.intp),
144 fix_gamma=None,
145 # should never be called, index dtype is int
146 ))
149def _rot90_dispatcher(m, k=None, axes=None):
150 return (m,)
153@array_function_dispatch(_rot90_dispatcher)
154def rot90(m, k=1, axes=(0, 1)):
155 """
156 Rotate an array by 90 degrees in the plane specified by axes.
158 Rotation direction is from the first towards the second axis.
159 This means for a 2D array with the default `k` and `axes`, the
160 rotation will be counterclockwise.
162 Parameters
163 ----------
164 m : array_like
165 Array of two or more dimensions.
166 k : integer
167 Number of times the array is rotated by 90 degrees.
168 axes : (2,) array_like
169 The array is rotated in the plane defined by the axes.
170 Axes must be different.
172 .. versionadded:: 1.12.0
174 Returns
175 -------
176 y : ndarray
177 A rotated view of `m`.
179 See Also
180 --------
181 flip : Reverse the order of elements in an array along the given axis.
182 fliplr : Flip an array horizontally.
183 flipud : Flip an array vertically.
185 Notes
186 -----
187 ``rot90(m, k=1, axes=(1,0))`` is the reverse of
188 ``rot90(m, k=1, axes=(0,1))``
190 ``rot90(m, k=1, axes=(1,0))`` is equivalent to
191 ``rot90(m, k=-1, axes=(0,1))``
193 Examples
194 --------
195 >>> m = np.array([[1,2],[3,4]], int)
196 >>> m
197 array([[1, 2],
198 [3, 4]])
199 >>> np.rot90(m)
200 array([[2, 4],
201 [1, 3]])
202 >>> np.rot90(m, 2)
203 array([[4, 3],
204 [2, 1]])
205 >>> m = np.arange(8).reshape((2,2,2))
206 >>> np.rot90(m, 1, (1,2))
207 array([[[1, 3],
208 [0, 2]],
209 [[5, 7],
210 [4, 6]]])
212 """
213 axes = tuple(axes)
214 if len(axes) != 2:
215 raise ValueError("len(axes) must be 2.")
217 m = asanyarray(m)
219 if axes[0] == axes[1] or absolute(axes[0] - axes[1]) == m.ndim:
220 raise ValueError("Axes must be different.")
222 if (axes[0] >= m.ndim or axes[0] < -m.ndim
223 or axes[1] >= m.ndim or axes[1] < -m.ndim):
224 raise ValueError("Axes={} out of range for array of ndim={}."
225 .format(axes, m.ndim))
227 k %= 4
229 if k == 0:
230 return m[:]
231 if k == 2:
232 return flip(flip(m, axes[0]), axes[1])
234 axes_list = arange(0, m.ndim)
235 (axes_list[axes[0]], axes_list[axes[1]]) = (axes_list[axes[1]],
236 axes_list[axes[0]])
238 if k == 1:
239 return transpose(flip(m, axes[1]), axes_list)
240 else:
241 # k == 3
242 return flip(transpose(m, axes_list), axes[1])
245def _flip_dispatcher(m, axis=None):
246 return (m,)
249@array_function_dispatch(_flip_dispatcher)
250def flip(m, axis=None):
251 """
252 Reverse the order of elements in an array along the given axis.
254 The shape of the array is preserved, but the elements are reordered.
256 .. versionadded:: 1.12.0
258 Parameters
259 ----------
260 m : array_like
261 Input array.
262 axis : None or int or tuple of ints, optional
263 Axis or axes along which to flip over. The default,
264 axis=None, will flip over all of the axes of the input array.
265 If axis is negative it counts from the last to the first axis.
267 If axis is a tuple of ints, flipping is performed on all of the axes
268 specified in the tuple.
270 .. versionchanged:: 1.15.0
271 None and tuples of axes are supported
273 Returns
274 -------
275 out : array_like
276 A view of `m` with the entries of axis reversed. Since a view is
277 returned, this operation is done in constant time.
279 See Also
280 --------
281 flipud : Flip an array vertically (axis=0).
282 fliplr : Flip an array horizontally (axis=1).
284 Notes
285 -----
286 flip(m, 0) is equivalent to flipud(m).
288 flip(m, 1) is equivalent to fliplr(m).
290 flip(m, n) corresponds to ``m[...,::-1,...]`` with ``::-1`` at position n.
292 flip(m) corresponds to ``m[::-1,::-1,...,::-1]`` with ``::-1`` at all
293 positions.
295 flip(m, (0, 1)) corresponds to ``m[::-1,::-1,...]`` with ``::-1`` at
296 position 0 and position 1.
298 Examples
299 --------
300 >>> A = np.arange(8).reshape((2,2,2))
301 >>> A
302 array([[[0, 1],
303 [2, 3]],
304 [[4, 5],
305 [6, 7]]])
306 >>> np.flip(A, 0)
307 array([[[4, 5],
308 [6, 7]],
309 [[0, 1],
310 [2, 3]]])
311 >>> np.flip(A, 1)
312 array([[[2, 3],
313 [0, 1]],
314 [[6, 7],
315 [4, 5]]])
316 >>> np.flip(A)
317 array([[[7, 6],
318 [5, 4]],
319 [[3, 2],
320 [1, 0]]])
321 >>> np.flip(A, (0, 2))
322 array([[[5, 4],
323 [7, 6]],
324 [[1, 0],
325 [3, 2]]])
326 >>> A = np.random.randn(3,4,5)
327 >>> np.all(np.flip(A,2) == A[:,:,::-1,...])
328 True
329 """
330 if not hasattr(m, 'ndim'):
331 m = asarray(m)
332 if axis is None:
333 indexer = (np.s_[::-1],) * m.ndim
334 else:
335 axis = _nx.normalize_axis_tuple(axis, m.ndim)
336 indexer = [np.s_[:]] * m.ndim
337 for ax in axis:
338 indexer[ax] = np.s_[::-1]
339 indexer = tuple(indexer)
340 return m[indexer]
343@set_module('numpy')
344def iterable(y):
345 """
346 Check whether or not an object can be iterated over.
348 Parameters
349 ----------
350 y : object
351 Input object.
353 Returns
354 -------
355 b : bool
356 Return ``True`` if the object has an iterator method or is a
357 sequence and ``False`` otherwise.
360 Examples
361 --------
362 >>> np.iterable([1, 2, 3])
363 True
364 >>> np.iterable(2)
365 False
367 Notes
368 -----
369 In most cases, the results of ``np.iterable(obj)`` are consistent with
370 ``isinstance(obj, collections.abc.Iterable)``. One notable exception is
371 the treatment of 0-dimensional arrays::
373 >>> from collections.abc import Iterable
374 >>> a = np.array(1.0) # 0-dimensional numpy array
375 >>> isinstance(a, Iterable)
376 True
377 >>> np.iterable(a)
378 False
380 """
381 try:
382 iter(y)
383 except TypeError:
384 return False
385 return True
388def _weights_are_valid(weights, a, axis):
389 """Validate weights array.
391 We assume, weights is not None.
392 """
393 wgt = np.asanyarray(weights)
395 # Sanity checks
396 if a.shape != wgt.shape:
397 if axis is None:
398 raise TypeError(
399 "Axis must be specified when shapes of a and weights "
400 "differ.")
401 if wgt.shape != tuple(a.shape[ax] for ax in axis):
402 raise ValueError(
403 "Shape of weights must be consistent with "
404 "shape of a along specified axis.")
406 # setup wgt to broadcast along axis
407 wgt = wgt.transpose(np.argsort(axis))
408 wgt = wgt.reshape(tuple((s if ax in axis else 1)
409 for ax, s in enumerate(a.shape)))
410 return wgt
413def _average_dispatcher(a, axis=None, weights=None, returned=None, *,
414 keepdims=None):
415 return (a, weights)
418@array_function_dispatch(_average_dispatcher)
419def average(a, axis=None, weights=None, returned=False, *,
420 keepdims=np._NoValue):
421 """
422 Compute the weighted average along the specified axis.
424 Parameters
425 ----------
426 a : array_like
427 Array containing data to be averaged. If `a` is not an array, a
428 conversion is attempted.
429 axis : None or int or tuple of ints, optional
430 Axis or axes along which to average `a`. The default,
431 `axis=None`, will average over all of the elements of the input array.
432 If axis is negative it counts from the last to the first axis.
434 .. versionadded:: 1.7.0
436 If axis is a tuple of ints, averaging is performed on all of the axes
437 specified in the tuple instead of a single axis or all the axes as
438 before.
439 weights : array_like, optional
440 An array of weights associated with the values in `a`. Each value in
441 `a` contributes to the average according to its associated weight.
442 The array of weights must be the same shape as `a` if no axis is
443 specified, otherwise the weights must have dimensions and shape
444 consistent with `a` along the specified axis.
445 If `weights=None`, then all data in `a` are assumed to have a
446 weight equal to one.
447 The calculation is::
449 avg = sum(a * weights) / sum(weights)
451 where the sum is over all included elements.
452 The only constraint on the values of `weights` is that `sum(weights)`
453 must not be 0.
454 returned : bool, optional
455 Default is `False`. If `True`, the tuple (`average`, `sum_of_weights`)
456 is returned, otherwise only the average is returned.
457 If `weights=None`, `sum_of_weights` is equivalent to the number of
458 elements over which the average is taken.
459 keepdims : bool, optional
460 If this is set to True, the axes which are reduced are left
461 in the result as dimensions with size one. With this option,
462 the result will broadcast correctly against the original `a`.
463 *Note:* `keepdims` will not work with instances of `numpy.matrix`
464 or other classes whose methods do not support `keepdims`.
466 .. versionadded:: 1.23.0
468 Returns
469 -------
470 retval, [sum_of_weights] : array_type or double
471 Return the average along the specified axis. When `returned` is `True`,
472 return a tuple with the average as the first element and the sum
473 of the weights as the second element. `sum_of_weights` is of the
474 same type as `retval`. The result dtype follows a general pattern.
475 If `weights` is None, the result dtype will be that of `a` , or ``float64``
476 if `a` is integral. Otherwise, if `weights` is not None and `a` is non-
477 integral, the result type will be the type of lowest precision capable of
478 representing values of both `a` and `weights`. If `a` happens to be
479 integral, the previous rules still applies but the result dtype will
480 at least be ``float64``.
482 Raises
483 ------
484 ZeroDivisionError
485 When all weights along axis are zero. See `numpy.ma.average` for a
486 version robust to this type of error.
487 TypeError
488 When `weights` does not have the same shape as `a`, and `axis=None`.
489 ValueError
490 When `weights` does not have dimensions and shape consistent with `a`
491 along specified `axis`.
493 See Also
494 --------
495 mean
497 ma.average : average for masked arrays -- useful if your data contains
498 "missing" values
499 numpy.result_type : Returns the type that results from applying the
500 numpy type promotion rules to the arguments.
502 Examples
503 --------
504 >>> data = np.arange(1, 5)
505 >>> data
506 array([1, 2, 3, 4])
507 >>> np.average(data)
508 2.5
509 >>> np.average(np.arange(1, 11), weights=np.arange(10, 0, -1))
510 4.0
512 >>> data = np.arange(6).reshape((3, 2))
513 >>> data
514 array([[0, 1],
515 [2, 3],
516 [4, 5]])
517 >>> np.average(data, axis=1, weights=[1./4, 3./4])
518 array([0.75, 2.75, 4.75])
519 >>> np.average(data, weights=[1./4, 3./4])
520 Traceback (most recent call last):
521 ...
522 TypeError: Axis must be specified when shapes of a and weights differ.
524 With ``keepdims=True``, the following result has shape (3, 1).
526 >>> np.average(data, axis=1, keepdims=True)
527 array([[0.5],
528 [2.5],
529 [4.5]])
531 >>> data = np.arange(8).reshape((2, 2, 2))
532 >>> data
533 array([[[0, 1],
534 [2, 3]],
535 [[4, 5],
536 [6, 7]]])
537 >>> np.average(data, axis=(0, 1), weights=[[1./4, 3./4], [1., 1./2]])
538 array([3.4, 4.4])
539 >>> np.average(data, axis=0, weights=[[1./4, 3./4], [1., 1./2]])
540 Traceback (most recent call last):
541 ...
542 ValueError: Shape of weights must be consistent
543 with shape of a along specified axis.
544 """
545 a = np.asanyarray(a)
547 if axis is not None:
548 axis = _nx.normalize_axis_tuple(axis, a.ndim, argname="axis")
550 if keepdims is np._NoValue:
551 # Don't pass on the keepdims argument if one wasn't given.
552 keepdims_kw = {}
553 else:
554 keepdims_kw = {'keepdims': keepdims}
556 if weights is None:
557 avg = a.mean(axis, **keepdims_kw)
558 avg_as_array = np.asanyarray(avg)
559 scl = avg_as_array.dtype.type(a.size/avg_as_array.size)
560 else:
561 wgt = _weights_are_valid(weights=weights, a=a, axis=axis)
563 if issubclass(a.dtype.type, (np.integer, np.bool)):
564 result_dtype = np.result_type(a.dtype, wgt.dtype, 'f8')
565 else:
566 result_dtype = np.result_type(a.dtype, wgt.dtype)
568 scl = wgt.sum(axis=axis, dtype=result_dtype, **keepdims_kw)
569 if np.any(scl == 0.0):
570 raise ZeroDivisionError(
571 "Weights sum to zero, can't be normalized")
573 avg = avg_as_array = np.multiply(a, wgt,
574 dtype=result_dtype).sum(axis, **keepdims_kw) / scl
576 if returned:
577 if scl.shape != avg_as_array.shape:
578 scl = np.broadcast_to(scl, avg_as_array.shape).copy()
579 return avg, scl
580 else:
581 return avg
584@set_module('numpy')
585def asarray_chkfinite(a, dtype=None, order=None):
586 """Convert the input to an array, checking for NaNs or Infs.
588 Parameters
589 ----------
590 a : array_like
591 Input data, in any form that can be converted to an array. This
592 includes lists, lists of tuples, tuples, tuples of tuples, tuples
593 of lists and ndarrays. Success requires no NaNs or Infs.
594 dtype : data-type, optional
595 By default, the data-type is inferred from the input data.
596 order : {'C', 'F', 'A', 'K'}, optional
597 Memory layout. 'A' and 'K' depend on the order of input array a.
598 'C' row-major (C-style),
599 'F' column-major (Fortran-style) memory representation.
600 'A' (any) means 'F' if `a` is Fortran contiguous, 'C' otherwise
601 'K' (keep) preserve input order
602 Defaults to 'C'.
604 Returns
605 -------
606 out : ndarray
607 Array interpretation of `a`. No copy is performed if the input
608 is already an ndarray. If `a` is a subclass of ndarray, a base
609 class ndarray is returned.
611 Raises
612 ------
613 ValueError
614 Raises ValueError if `a` contains NaN (Not a Number) or Inf (Infinity).
616 See Also
617 --------
618 asarray : Create and array.
619 asanyarray : Similar function which passes through subclasses.
620 ascontiguousarray : Convert input to a contiguous array.
621 asfortranarray : Convert input to an ndarray with column-major
622 memory order.
623 fromiter : Create an array from an iterator.
624 fromfunction : Construct an array by executing a function on grid
625 positions.
627 Examples
628 --------
629 Convert a list into an array. If all elements are finite
630 ``asarray_chkfinite`` is identical to ``asarray``.
632 >>> a = [1, 2]
633 >>> np.asarray_chkfinite(a, dtype=float)
634 array([1., 2.])
636 Raises ValueError if array_like contains Nans or Infs.
638 >>> a = [1, 2, np.inf]
639 >>> try:
640 ... np.asarray_chkfinite(a)
641 ... except ValueError:
642 ... print('ValueError')
643 ...
644 ValueError
646 """
647 a = asarray(a, dtype=dtype, order=order)
648 if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
649 raise ValueError(
650 "array must not contain infs or NaNs")
651 return a
654def _piecewise_dispatcher(x, condlist, funclist, *args, **kw):
655 yield x
656 # support the undocumented behavior of allowing scalars
657 if np.iterable(condlist):
658 yield from condlist
661@array_function_dispatch(_piecewise_dispatcher)
662def piecewise(x, condlist, funclist, *args, **kw):
663 """
664 Evaluate a piecewise-defined function.
666 Given a set of conditions and corresponding functions, evaluate each
667 function on the input data wherever its condition is true.
669 Parameters
670 ----------
671 x : ndarray or scalar
672 The input domain.
673 condlist : list of bool arrays or bool scalars
674 Each boolean array corresponds to a function in `funclist`. Wherever
675 `condlist[i]` is True, `funclist[i](x)` is used as the output value.
677 Each boolean array in `condlist` selects a piece of `x`,
678 and should therefore be of the same shape as `x`.
680 The length of `condlist` must correspond to that of `funclist`.
681 If one extra function is given, i.e. if
682 ``len(funclist) == len(condlist) + 1``, then that extra function
683 is the default value, used wherever all conditions are false.
684 funclist : list of callables, f(x,*args,**kw), or scalars
685 Each function is evaluated over `x` wherever its corresponding
686 condition is True. It should take a 1d array as input and give an 1d
687 array or a scalar value as output. If, instead of a callable,
688 a scalar is provided then a constant function (``lambda x: scalar``) is
689 assumed.
690 args : tuple, optional
691 Any further arguments given to `piecewise` are passed to the functions
692 upon execution, i.e., if called ``piecewise(..., ..., 1, 'a')``, then
693 each function is called as ``f(x, 1, 'a')``.
694 kw : dict, optional
695 Keyword arguments used in calling `piecewise` are passed to the
696 functions upon execution, i.e., if called
697 ``piecewise(..., ..., alpha=1)``, then each function is called as
698 ``f(x, alpha=1)``.
700 Returns
701 -------
702 out : ndarray
703 The output is the same shape and type as x and is found by
704 calling the functions in `funclist` on the appropriate portions of `x`,
705 as defined by the boolean arrays in `condlist`. Portions not covered
706 by any condition have a default value of 0.
709 See Also
710 --------
711 choose, select, where
713 Notes
714 -----
715 This is similar to choose or select, except that functions are
716 evaluated on elements of `x` that satisfy the corresponding condition from
717 `condlist`.
719 The result is::
721 |--
722 |funclist[0](x[condlist[0]])
723 out = |funclist[1](x[condlist[1]])
724 |...
725 |funclist[n2](x[condlist[n2]])
726 |--
728 Examples
729 --------
730 Define the signum function, which is -1 for ``x < 0`` and +1 for ``x >= 0``.
732 >>> x = np.linspace(-2.5, 2.5, 6)
733 >>> np.piecewise(x, [x < 0, x >= 0], [-1, 1])
734 array([-1., -1., -1., 1., 1., 1.])
736 Define the absolute value, which is ``-x`` for ``x <0`` and ``x`` for
737 ``x >= 0``.
739 >>> np.piecewise(x, [x < 0, x >= 0], [lambda x: -x, lambda x: x])
740 array([2.5, 1.5, 0.5, 0.5, 1.5, 2.5])
742 Apply the same function to a scalar value.
744 >>> y = -2
745 >>> np.piecewise(y, [y < 0, y >= 0], [lambda x: -x, lambda x: x])
746 array(2)
748 """
749 x = asanyarray(x)
750 n2 = len(funclist)
752 # undocumented: single condition is promoted to a list of one condition
753 if isscalar(condlist) or (
754 not isinstance(condlist[0], (list, ndarray)) and x.ndim != 0):
755 condlist = [condlist]
757 condlist = asarray(condlist, dtype=bool)
758 n = len(condlist)
760 if n == n2 - 1: # compute the "otherwise" condition.
761 condelse = ~np.any(condlist, axis=0, keepdims=True)
762 condlist = np.concatenate([condlist, condelse], axis=0)
763 n += 1
764 elif n != n2:
765 raise ValueError(
766 "with {} condition(s), either {} or {} functions are expected"
767 .format(n, n, n+1)
768 )
770 y = zeros_like(x)
771 for cond, func in zip(condlist, funclist):
772 if not isinstance(func, collections.abc.Callable):
773 y[cond] = func
774 else:
775 vals = x[cond]
776 if vals.size > 0:
777 y[cond] = func(vals, *args, **kw)
779 return y
782def _select_dispatcher(condlist, choicelist, default=None):
783 yield from condlist
784 yield from choicelist
787@array_function_dispatch(_select_dispatcher)
788def select(condlist, choicelist, default=0):
789 """
790 Return an array drawn from elements in choicelist, depending on conditions.
792 Parameters
793 ----------
794 condlist : list of bool ndarrays
795 The list of conditions which determine from which array in `choicelist`
796 the output elements are taken. When multiple conditions are satisfied,
797 the first one encountered in `condlist` is used.
798 choicelist : list of ndarrays
799 The list of arrays from which the output elements are taken. It has
800 to be of the same length as `condlist`.
801 default : scalar, optional
802 The element inserted in `output` when all conditions evaluate to False.
804 Returns
805 -------
806 output : ndarray
807 The output at position m is the m-th element of the array in
808 `choicelist` where the m-th element of the corresponding array in
809 `condlist` is True.
811 See Also
812 --------
813 where : Return elements from one of two arrays depending on condition.
814 take, choose, compress, diag, diagonal
816 Examples
817 --------
818 Beginning with an array of integers from 0 to 5 (inclusive),
819 elements less than ``3`` are negated, elements greater than ``3``
820 are squared, and elements not meeting either of these conditions
821 (exactly ``3``) are replaced with a `default` value of ``42``.
823 >>> x = np.arange(6)
824 >>> condlist = [x<3, x>3]
825 >>> choicelist = [x, x**2]
826 >>> np.select(condlist, choicelist, 42)
827 array([ 0, 1, 2, 42, 16, 25])
829 When multiple conditions are satisfied, the first one encountered in
830 `condlist` is used.
832 >>> condlist = [x<=4, x>3]
833 >>> choicelist = [x, x**2]
834 >>> np.select(condlist, choicelist, 55)
835 array([ 0, 1, 2, 3, 4, 25])
837 """
838 # Check the size of condlist and choicelist are the same, or abort.
839 if len(condlist) != len(choicelist):
840 raise ValueError(
841 'list of cases must be same length as list of conditions')
843 # Now that the dtype is known, handle the deprecated select([], []) case
844 if len(condlist) == 0:
845 raise ValueError("select with an empty condition list is not possible")
847 # TODO: This preserves the Python int, float, complex manually to get the
848 # right `result_type` with NEP 50. Most likely we will grow a better
849 # way to spell this (and this can be replaced).
850 choicelist = [
851 choice if type(choice) in (int, float, complex) else np.asarray(choice)
852 for choice in choicelist]
853 choicelist.append(default if type(default) in (int, float, complex)
854 else np.asarray(default))
856 try:
857 dtype = np.result_type(*choicelist)
858 except TypeError as e:
859 msg = f'Choicelist and default value do not have a common dtype: {e}'
860 raise TypeError(msg) from None
862 # Convert conditions to arrays and broadcast conditions and choices
863 # as the shape is needed for the result. Doing it separately optimizes
864 # for example when all choices are scalars.
865 condlist = np.broadcast_arrays(*condlist)
866 choicelist = np.broadcast_arrays(*choicelist)
868 # If cond array is not an ndarray in boolean format or scalar bool, abort.
869 for i, cond in enumerate(condlist):
870 if cond.dtype.type is not np.bool:
871 raise TypeError(
872 'invalid entry {} in condlist: should be boolean ndarray'.format(i))
874 if choicelist[0].ndim == 0:
875 # This may be common, so avoid the call.
876 result_shape = condlist[0].shape
877 else:
878 result_shape = np.broadcast_arrays(condlist[0], choicelist[0])[0].shape
880 result = np.full(result_shape, choicelist[-1], dtype)
882 # Use np.copyto to burn each choicelist array onto result, using the
883 # corresponding condlist as a boolean mask. This is done in reverse
884 # order since the first choice should take precedence.
885 choicelist = choicelist[-2::-1]
886 condlist = condlist[::-1]
887 for choice, cond in zip(choicelist, condlist):
888 np.copyto(result, choice, where=cond)
890 return result
893def _copy_dispatcher(a, order=None, subok=None):
894 return (a,)
897@array_function_dispatch(_copy_dispatcher)
898def copy(a, order='K', subok=False):
899 """
900 Return an array copy of the given object.
902 Parameters
903 ----------
904 a : array_like
905 Input data.
906 order : {'C', 'F', 'A', 'K'}, optional
907 Controls the memory layout of the copy. 'C' means C-order,
908 'F' means F-order, 'A' means 'F' if `a` is Fortran contiguous,
909 'C' otherwise. 'K' means match the layout of `a` as closely
910 as possible. (Note that this function and :meth:`ndarray.copy` are very
911 similar, but have different default values for their order=
912 arguments.)
913 subok : bool, optional
914 If True, then sub-classes will be passed-through, otherwise the
915 returned array will be forced to be a base-class array (defaults to False).
917 .. versionadded:: 1.19.0
919 Returns
920 -------
921 arr : ndarray
922 Array interpretation of `a`.
924 See Also
925 --------
926 ndarray.copy : Preferred method for creating an array copy
928 Notes
929 -----
930 This is equivalent to:
932 >>> np.array(a, copy=True) #doctest: +SKIP
934 The copy made of the data is shallow, i.e., for arrays with object dtype,
935 the new array will point to the same objects.
936 See Examples from `ndarray.copy`.
938 Examples
939 --------
940 Create an array x, with a reference y and a copy z:
942 >>> x = np.array([1, 2, 3])
943 >>> y = x
944 >>> z = np.copy(x)
946 Note that, when we modify x, y changes, but not z:
948 >>> x[0] = 10
949 >>> x[0] == y[0]
950 True
951 >>> x[0] == z[0]
952 False
954 Note that, np.copy clears previously set WRITEABLE=False flag.
956 >>> a = np.array([1, 2, 3])
957 >>> a.flags["WRITEABLE"] = False
958 >>> b = np.copy(a)
959 >>> b.flags["WRITEABLE"]
960 True
961 >>> b[0] = 3
962 >>> b
963 array([3, 2, 3])
964 """
965 return array(a, order=order, subok=subok, copy=True)
967# Basic operations
970def _gradient_dispatcher(f, *varargs, axis=None, edge_order=None):
971 yield f
972 yield from varargs
975@array_function_dispatch(_gradient_dispatcher)
976def gradient(f, *varargs, axis=None, edge_order=1):
977 """
978 Return the gradient of an N-dimensional array.
980 The gradient is computed using second order accurate central differences
981 in the interior points and either first or second order accurate one-sides
982 (forward or backwards) differences at the boundaries.
983 The returned gradient hence has the same shape as the input array.
985 Parameters
986 ----------
987 f : array_like
988 An N-dimensional array containing samples of a scalar function.
989 varargs : list of scalar or array, optional
990 Spacing between f values. Default unitary spacing for all dimensions.
991 Spacing can be specified using:
993 1. single scalar to specify a sample distance for all dimensions.
994 2. N scalars to specify a constant sample distance for each dimension.
995 i.e. `dx`, `dy`, `dz`, ...
996 3. N arrays to specify the coordinates of the values along each
997 dimension of F. The length of the array must match the size of
998 the corresponding dimension
999 4. Any combination of N scalars/arrays with the meaning of 2. and 3.
1001 If `axis` is given, the number of varargs must equal the number of axes.
1002 Default: 1.
1004 edge_order : {1, 2}, optional
1005 Gradient is calculated using N-th order accurate differences
1006 at the boundaries. Default: 1.
1008 .. versionadded:: 1.9.1
1010 axis : None or int or tuple of ints, optional
1011 Gradient is calculated only along the given axis or axes
1012 The default (axis = None) is to calculate the gradient for all the axes
1013 of the input array. axis may be negative, in which case it counts from
1014 the last to the first axis.
1016 .. versionadded:: 1.11.0
1018 Returns
1019 -------
1020 gradient : ndarray or list of ndarray
1021 A list of ndarrays (or a single ndarray if there is only one dimension)
1022 corresponding to the derivatives of f with respect to each dimension.
1023 Each derivative has the same shape as f.
1025 Examples
1026 --------
1027 >>> f = np.array([1, 2, 4, 7, 11, 16], dtype=float)
1028 >>> np.gradient(f)
1029 array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
1030 >>> np.gradient(f, 2)
1031 array([0.5 , 0.75, 1.25, 1.75, 2.25, 2.5 ])
1033 Spacing can be also specified with an array that represents the coordinates
1034 of the values F along the dimensions.
1035 For instance a uniform spacing:
1037 >>> x = np.arange(f.size)
1038 >>> np.gradient(f, x)
1039 array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
1041 Or a non uniform one:
1043 >>> x = np.array([0., 1., 1.5, 3.5, 4., 6.], dtype=float)
1044 >>> np.gradient(f, x)
1045 array([1. , 3. , 3.5, 6.7, 6.9, 2.5])
1047 For two dimensional arrays, the return will be two arrays ordered by
1048 axis. In this example the first array stands for the gradient in
1049 rows and the second one in columns direction:
1051 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float))
1052 [array([[ 2., 2., -1.],
1053 [ 2., 2., -1.]]), array([[1. , 2.5, 4. ],
1054 [1. , 1. , 1. ]])]
1056 In this example the spacing is also specified:
1057 uniform for axis=0 and non uniform for axis=1
1059 >>> dx = 2.
1060 >>> y = [1., 1.5, 3.5]
1061 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), dx, y)
1062 [array([[ 1. , 1. , -0.5],
1063 [ 1. , 1. , -0.5]]), array([[2. , 2. , 2. ],
1064 [2. , 1.7, 0.5]])]
1066 It is possible to specify how boundaries are treated using `edge_order`
1068 >>> x = np.array([0, 1, 2, 3, 4])
1069 >>> f = x**2
1070 >>> np.gradient(f, edge_order=1)
1071 array([1., 2., 4., 6., 7.])
1072 >>> np.gradient(f, edge_order=2)
1073 array([0., 2., 4., 6., 8.])
1075 The `axis` keyword can be used to specify a subset of axes of which the
1076 gradient is calculated
1078 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), axis=0)
1079 array([[ 2., 2., -1.],
1080 [ 2., 2., -1.]])
1082 Notes
1083 -----
1084 Assuming that :math:`f\\in C^{3}` (i.e., :math:`f` has at least 3 continuous
1085 derivatives) and let :math:`h_{*}` be a non-homogeneous stepsize, we
1086 minimize the "consistency error" :math:`\\eta_{i}` between the true gradient
1087 and its estimate from a linear combination of the neighboring grid-points:
1089 .. math::
1091 \\eta_{i} = f_{i}^{\\left(1\\right)} -
1092 \\left[ \\alpha f\\left(x_{i}\\right) +
1093 \\beta f\\left(x_{i} + h_{d}\\right) +
1094 \\gamma f\\left(x_{i}-h_{s}\\right)
1095 \\right]
1097 By substituting :math:`f(x_{i} + h_{d})` and :math:`f(x_{i} - h_{s})`
1098 with their Taylor series expansion, this translates into solving
1099 the following the linear system:
1101 .. math::
1103 \\left\\{
1104 \\begin{array}{r}
1105 \\alpha+\\beta+\\gamma=0 \\\\
1106 \\beta h_{d}-\\gamma h_{s}=1 \\\\
1107 \\beta h_{d}^{2}+\\gamma h_{s}^{2}=0
1108 \\end{array}
1109 \\right.
1111 The resulting approximation of :math:`f_{i}^{(1)}` is the following:
1113 .. math::
1115 \\hat f_{i}^{(1)} =
1116 \\frac{
1117 h_{s}^{2}f\\left(x_{i} + h_{d}\\right)
1118 + \\left(h_{d}^{2} - h_{s}^{2}\\right)f\\left(x_{i}\\right)
1119 - h_{d}^{2}f\\left(x_{i}-h_{s}\\right)}
1120 { h_{s}h_{d}\\left(h_{d} + h_{s}\\right)}
1121 + \\mathcal{O}\\left(\\frac{h_{d}h_{s}^{2}
1122 + h_{s}h_{d}^{2}}{h_{d}
1123 + h_{s}}\\right)
1125 It is worth noting that if :math:`h_{s}=h_{d}`
1126 (i.e., data are evenly spaced)
1127 we find the standard second order approximation:
1129 .. math::
1131 \\hat f_{i}^{(1)}=
1132 \\frac{f\\left(x_{i+1}\\right) - f\\left(x_{i-1}\\right)}{2h}
1133 + \\mathcal{O}\\left(h^{2}\\right)
1135 With a similar procedure the forward/backward approximations used for
1136 boundaries can be derived.
1138 References
1139 ----------
1140 .. [1] Quarteroni A., Sacco R., Saleri F. (2007) Numerical Mathematics
1141 (Texts in Applied Mathematics). New York: Springer.
1142 .. [2] Durran D. R. (1999) Numerical Methods for Wave Equations
1143 in Geophysical Fluid Dynamics. New York: Springer.
1144 .. [3] Fornberg B. (1988) Generation of Finite Difference Formulas on
1145 Arbitrarily Spaced Grids,
1146 Mathematics of Computation 51, no. 184 : 699-706.
1147 `PDF <https://www.ams.org/journals/mcom/1988-51-184/
1148 S0025-5718-1988-0935077-0/S0025-5718-1988-0935077-0.pdf>`_.
1149 """
1150 f = np.asanyarray(f)
1151 N = f.ndim # number of dimensions
1153 if axis is None:
1154 axes = tuple(range(N))
1155 else:
1156 axes = _nx.normalize_axis_tuple(axis, N)
1158 len_axes = len(axes)
1159 n = len(varargs)
1160 if n == 0:
1161 # no spacing argument - use 1 in all axes
1162 dx = [1.0] * len_axes
1163 elif n == 1 and np.ndim(varargs[0]) == 0:
1164 # single scalar for all axes
1165 dx = varargs * len_axes
1166 elif n == len_axes:
1167 # scalar or 1d array for each axis
1168 dx = list(varargs)
1169 for i, distances in enumerate(dx):
1170 distances = np.asanyarray(distances)
1171 if distances.ndim == 0:
1172 continue
1173 elif distances.ndim != 1:
1174 raise ValueError("distances must be either scalars or 1d")
1175 if len(distances) != f.shape[axes[i]]:
1176 raise ValueError("when 1d, distances must match "
1177 "the length of the corresponding dimension")
1178 if np.issubdtype(distances.dtype, np.integer):
1179 # Convert numpy integer types to float64 to avoid modular
1180 # arithmetic in np.diff(distances).
1181 distances = distances.astype(np.float64)
1182 diffx = np.diff(distances)
1183 # if distances are constant reduce to the scalar case
1184 # since it brings a consistent speedup
1185 if (diffx == diffx[0]).all():
1186 diffx = diffx[0]
1187 dx[i] = diffx
1188 else:
1189 raise TypeError("invalid number of arguments")
1191 if edge_order > 2:
1192 raise ValueError("'edge_order' greater than 2 not supported")
1194 # use central differences on interior and one-sided differences on the
1195 # endpoints. This preserves second order-accuracy over the full domain.
1197 outvals = []
1199 # create slice objects --- initially all are [:, :, ..., :]
1200 slice1 = [slice(None)]*N
1201 slice2 = [slice(None)]*N
1202 slice3 = [slice(None)]*N
1203 slice4 = [slice(None)]*N
1205 otype = f.dtype
1206 if otype.type is np.datetime64:
1207 # the timedelta dtype with the same unit information
1208 otype = np.dtype(otype.name.replace('datetime', 'timedelta'))
1209 # view as timedelta to allow addition
1210 f = f.view(otype)
1211 elif otype.type is np.timedelta64:
1212 pass
1213 elif np.issubdtype(otype, np.inexact):
1214 pass
1215 else:
1216 # All other types convert to floating point.
1217 # First check if f is a numpy integer type; if so, convert f to float64
1218 # to avoid modular arithmetic when computing the changes in f.
1219 if np.issubdtype(otype, np.integer):
1220 f = f.astype(np.float64)
1221 otype = np.float64
1223 for axis, ax_dx in zip(axes, dx):
1224 if f.shape[axis] < edge_order + 1:
1225 raise ValueError(
1226 "Shape of array too small to calculate a numerical gradient, "
1227 "at least (edge_order + 1) elements are required.")
1228 # result allocation
1229 out = np.empty_like(f, dtype=otype)
1231 # spacing for the current axis
1232 uniform_spacing = np.ndim(ax_dx) == 0
1234 # Numerical differentiation: 2nd order interior
1235 slice1[axis] = slice(1, -1)
1236 slice2[axis] = slice(None, -2)
1237 slice3[axis] = slice(1, -1)
1238 slice4[axis] = slice(2, None)
1240 if uniform_spacing:
1241 out[tuple(slice1)] = (f[tuple(slice4)] - f[tuple(slice2)]) / (2. * ax_dx)
1242 else:
1243 dx1 = ax_dx[0:-1]
1244 dx2 = ax_dx[1:]
1245 a = -(dx2)/(dx1 * (dx1 + dx2))
1246 b = (dx2 - dx1) / (dx1 * dx2)
1247 c = dx1 / (dx2 * (dx1 + dx2))
1248 # fix the shape for broadcasting
1249 shape = np.ones(N, dtype=int)
1250 shape[axis] = -1
1251 a.shape = b.shape = c.shape = shape
1252 # 1D equivalent -- out[1:-1] = a * f[:-2] + b * f[1:-1] + c * f[2:]
1253 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
1255 # Numerical differentiation: 1st order edges
1256 if edge_order == 1:
1257 slice1[axis] = 0
1258 slice2[axis] = 1
1259 slice3[axis] = 0
1260 dx_0 = ax_dx if uniform_spacing else ax_dx[0]
1261 # 1D equivalent -- out[0] = (f[1] - f[0]) / (x[1] - x[0])
1262 out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0
1264 slice1[axis] = -1
1265 slice2[axis] = -1
1266 slice3[axis] = -2
1267 dx_n = ax_dx if uniform_spacing else ax_dx[-1]
1268 # 1D equivalent -- out[-1] = (f[-1] - f[-2]) / (x[-1] - x[-2])
1269 out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_n
1271 # Numerical differentiation: 2nd order edges
1272 else:
1273 slice1[axis] = 0
1274 slice2[axis] = 0
1275 slice3[axis] = 1
1276 slice4[axis] = 2
1277 if uniform_spacing:
1278 a = -1.5 / ax_dx
1279 b = 2. / ax_dx
1280 c = -0.5 / ax_dx
1281 else:
1282 dx1 = ax_dx[0]
1283 dx2 = ax_dx[1]
1284 a = -(2. * dx1 + dx2)/(dx1 * (dx1 + dx2))
1285 b = (dx1 + dx2) / (dx1 * dx2)
1286 c = - dx1 / (dx2 * (dx1 + dx2))
1287 # 1D equivalent -- out[0] = a * f[0] + b * f[1] + c * f[2]
1288 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
1290 slice1[axis] = -1
1291 slice2[axis] = -3
1292 slice3[axis] = -2
1293 slice4[axis] = -1
1294 if uniform_spacing:
1295 a = 0.5 / ax_dx
1296 b = -2. / ax_dx
1297 c = 1.5 / ax_dx
1298 else:
1299 dx1 = ax_dx[-2]
1300 dx2 = ax_dx[-1]
1301 a = (dx2) / (dx1 * (dx1 + dx2))
1302 b = - (dx2 + dx1) / (dx1 * dx2)
1303 c = (2. * dx2 + dx1) / (dx2 * (dx1 + dx2))
1304 # 1D equivalent -- out[-1] = a * f[-3] + b * f[-2] + c * f[-1]
1305 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
1307 outvals.append(out)
1309 # reset the slice object in this dimension to ":"
1310 slice1[axis] = slice(None)
1311 slice2[axis] = slice(None)
1312 slice3[axis] = slice(None)
1313 slice4[axis] = slice(None)
1315 if len_axes == 1:
1316 return outvals[0]
1317 return tuple(outvals)
1320def _diff_dispatcher(a, n=None, axis=None, prepend=None, append=None):
1321 return (a, prepend, append)
1324@array_function_dispatch(_diff_dispatcher)
1325def diff(a, n=1, axis=-1, prepend=np._NoValue, append=np._NoValue):
1326 """
1327 Calculate the n-th discrete difference along the given axis.
1329 The first difference is given by ``out[i] = a[i+1] - a[i]`` along
1330 the given axis, higher differences are calculated by using `diff`
1331 recursively.
1333 Parameters
1334 ----------
1335 a : array_like
1336 Input array
1337 n : int, optional
1338 The number of times values are differenced. If zero, the input
1339 is returned as-is.
1340 axis : int, optional
1341 The axis along which the difference is taken, default is the
1342 last axis.
1343 prepend, append : array_like, optional
1344 Values to prepend or append to `a` along axis prior to
1345 performing the difference. Scalar values are expanded to
1346 arrays with length 1 in the direction of axis and the shape
1347 of the input array in along all other axes. Otherwise the
1348 dimension and shape must match `a` except along axis.
1350 .. versionadded:: 1.16.0
1352 Returns
1353 -------
1354 diff : ndarray
1355 The n-th differences. The shape of the output is the same as `a`
1356 except along `axis` where the dimension is smaller by `n`. The
1357 type of the output is the same as the type of the difference
1358 between any two elements of `a`. This is the same as the type of
1359 `a` in most cases. A notable exception is `datetime64`, which
1360 results in a `timedelta64` output array.
1362 See Also
1363 --------
1364 gradient, ediff1d, cumsum
1366 Notes
1367 -----
1368 Type is preserved for boolean arrays, so the result will contain
1369 `False` when consecutive elements are the same and `True` when they
1370 differ.
1372 For unsigned integer arrays, the results will also be unsigned. This
1373 should not be surprising, as the result is consistent with
1374 calculating the difference directly:
1376 >>> u8_arr = np.array([1, 0], dtype=np.uint8)
1377 >>> np.diff(u8_arr)
1378 array([255], dtype=uint8)
1379 >>> u8_arr[1,...] - u8_arr[0,...]
1380 255
1382 If this is not desirable, then the array should be cast to a larger
1383 integer type first:
1385 >>> i16_arr = u8_arr.astype(np.int16)
1386 >>> np.diff(i16_arr)
1387 array([-1], dtype=int16)
1389 Examples
1390 --------
1391 >>> x = np.array([1, 2, 4, 7, 0])
1392 >>> np.diff(x)
1393 array([ 1, 2, 3, -7])
1394 >>> np.diff(x, n=2)
1395 array([ 1, 1, -10])
1397 >>> x = np.array([[1, 3, 6, 10], [0, 5, 6, 8]])
1398 >>> np.diff(x)
1399 array([[2, 3, 4],
1400 [5, 1, 2]])
1401 >>> np.diff(x, axis=0)
1402 array([[-1, 2, 0, -2]])
1404 >>> x = np.arange('1066-10-13', '1066-10-16', dtype=np.datetime64)
1405 >>> np.diff(x)
1406 array([1, 1], dtype='timedelta64[D]')
1408 """
1409 if n == 0:
1410 return a
1411 if n < 0:
1412 raise ValueError(
1413 "order must be non-negative but got " + repr(n))
1415 a = asanyarray(a)
1416 nd = a.ndim
1417 if nd == 0:
1418 raise ValueError("diff requires input that is at least one dimensional")
1419 axis = normalize_axis_index(axis, nd)
1421 combined = []
1422 if prepend is not np._NoValue:
1423 prepend = np.asanyarray(prepend)
1424 if prepend.ndim == 0:
1425 shape = list(a.shape)
1426 shape[axis] = 1
1427 prepend = np.broadcast_to(prepend, tuple(shape))
1428 combined.append(prepend)
1430 combined.append(a)
1432 if append is not np._NoValue:
1433 append = np.asanyarray(append)
1434 if append.ndim == 0:
1435 shape = list(a.shape)
1436 shape[axis] = 1
1437 append = np.broadcast_to(append, tuple(shape))
1438 combined.append(append)
1440 if len(combined) > 1:
1441 a = np.concatenate(combined, axis)
1443 slice1 = [slice(None)] * nd
1444 slice2 = [slice(None)] * nd
1445 slice1[axis] = slice(1, None)
1446 slice2[axis] = slice(None, -1)
1447 slice1 = tuple(slice1)
1448 slice2 = tuple(slice2)
1450 op = not_equal if a.dtype == np.bool else subtract
1451 for _ in range(n):
1452 a = op(a[slice1], a[slice2])
1454 return a
1457def _interp_dispatcher(x, xp, fp, left=None, right=None, period=None):
1458 return (x, xp, fp)
1461@array_function_dispatch(_interp_dispatcher)
1462def interp(x, xp, fp, left=None, right=None, period=None):
1463 """
1464 One-dimensional linear interpolation for monotonically increasing sample points.
1466 Returns the one-dimensional piecewise linear interpolant to a function
1467 with given discrete data points (`xp`, `fp`), evaluated at `x`.
1469 Parameters
1470 ----------
1471 x : array_like
1472 The x-coordinates at which to evaluate the interpolated values.
1474 xp : 1-D sequence of floats
1475 The x-coordinates of the data points, must be increasing if argument
1476 `period` is not specified. Otherwise, `xp` is internally sorted after
1477 normalizing the periodic boundaries with ``xp = xp % period``.
1479 fp : 1-D sequence of float or complex
1480 The y-coordinates of the data points, same length as `xp`.
1482 left : optional float or complex corresponding to fp
1483 Value to return for `x < xp[0]`, default is `fp[0]`.
1485 right : optional float or complex corresponding to fp
1486 Value to return for `x > xp[-1]`, default is `fp[-1]`.
1488 period : None or float, optional
1489 A period for the x-coordinates. This parameter allows the proper
1490 interpolation of angular x-coordinates. Parameters `left` and `right`
1491 are ignored if `period` is specified.
1493 .. versionadded:: 1.10.0
1495 Returns
1496 -------
1497 y : float or complex (corresponding to fp) or ndarray
1498 The interpolated values, same shape as `x`.
1500 Raises
1501 ------
1502 ValueError
1503 If `xp` and `fp` have different length
1504 If `xp` or `fp` are not 1-D sequences
1505 If `period == 0`
1507 See Also
1508 --------
1509 scipy.interpolate
1511 Warnings
1512 --------
1513 The x-coordinate sequence is expected to be increasing, but this is not
1514 explicitly enforced. However, if the sequence `xp` is non-increasing,
1515 interpolation results are meaningless.
1517 Note that, since NaN is unsortable, `xp` also cannot contain NaNs.
1519 A simple check for `xp` being strictly increasing is::
1521 np.all(np.diff(xp) > 0)
1523 Examples
1524 --------
1525 >>> xp = [1, 2, 3]
1526 >>> fp = [3, 2, 0]
1527 >>> np.interp(2.5, xp, fp)
1528 1.0
1529 >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
1530 array([3. , 3. , 2.5 , 0.56, 0. ])
1531 >>> UNDEF = -99.0
1532 >>> np.interp(3.14, xp, fp, right=UNDEF)
1533 -99.0
1535 Plot an interpolant to the sine function:
1537 >>> x = np.linspace(0, 2*np.pi, 10)
1538 >>> y = np.sin(x)
1539 >>> xvals = np.linspace(0, 2*np.pi, 50)
1540 >>> yinterp = np.interp(xvals, x, y)
1541 >>> import matplotlib.pyplot as plt
1542 >>> plt.plot(x, y, 'o')
1543 [<matplotlib.lines.Line2D object at 0x...>]
1544 >>> plt.plot(xvals, yinterp, '-x')
1545 [<matplotlib.lines.Line2D object at 0x...>]
1546 >>> plt.show()
1548 Interpolation with periodic x-coordinates:
1550 >>> x = [-180, -170, -185, 185, -10, -5, 0, 365]
1551 >>> xp = [190, -190, 350, -350]
1552 >>> fp = [5, 10, 3, 4]
1553 >>> np.interp(x, xp, fp, period=360)
1554 array([7.5 , 5. , 8.75, 6.25, 3. , 3.25, 3.5 , 3.75])
1556 Complex interpolation:
1558 >>> x = [1.5, 4.0]
1559 >>> xp = [2,3,5]
1560 >>> fp = [1.0j, 0, 2+3j]
1561 >>> np.interp(x, xp, fp)
1562 array([0.+1.j , 1.+1.5j])
1564 """
1566 fp = np.asarray(fp)
1568 if np.iscomplexobj(fp):
1569 interp_func = compiled_interp_complex
1570 input_dtype = np.complex128
1571 else:
1572 interp_func = compiled_interp
1573 input_dtype = np.float64
1575 if period is not None:
1576 if period == 0:
1577 raise ValueError("period must be a non-zero value")
1578 period = abs(period)
1579 left = None
1580 right = None
1582 x = np.asarray(x, dtype=np.float64)
1583 xp = np.asarray(xp, dtype=np.float64)
1584 fp = np.asarray(fp, dtype=input_dtype)
1586 if xp.ndim != 1 or fp.ndim != 1:
1587 raise ValueError("Data points must be 1-D sequences")
1588 if xp.shape[0] != fp.shape[0]:
1589 raise ValueError("fp and xp are not of the same length")
1590 # normalizing periodic boundaries
1591 x = x % period
1592 xp = xp % period
1593 asort_xp = np.argsort(xp)
1594 xp = xp[asort_xp]
1595 fp = fp[asort_xp]
1596 xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period))
1597 fp = np.concatenate((fp[-1:], fp, fp[0:1]))
1599 return interp_func(x, xp, fp, left, right)
1602def _angle_dispatcher(z, deg=None):
1603 return (z,)
1606@array_function_dispatch(_angle_dispatcher)
1607def angle(z, deg=False):
1608 """
1609 Return the angle of the complex argument.
1611 Parameters
1612 ----------
1613 z : array_like
1614 A complex number or sequence of complex numbers.
1615 deg : bool, optional
1616 Return angle in degrees if True, radians if False (default).
1618 Returns
1619 -------
1620 angle : ndarray or scalar
1621 The counterclockwise angle from the positive real axis on the complex
1622 plane in the range ``(-pi, pi]``, with dtype as numpy.float64.
1624 .. versionchanged:: 1.16.0
1625 This function works on subclasses of ndarray like `ma.array`.
1627 See Also
1628 --------
1629 arctan2
1630 absolute
1632 Notes
1633 -----
1634 This function passes the imaginary and real parts of the argument to
1635 `arctan2` to compute the result; consequently, it follows the convention
1636 of `arctan2` when the magnitude of the argument is zero. See example.
1638 Examples
1639 --------
1640 >>> np.angle([1.0, 1.0j, 1+1j]) # in radians
1641 array([ 0. , 1.57079633, 0.78539816]) # may vary
1642 >>> np.angle(1+1j, deg=True) # in degrees
1643 45.0
1644 >>> np.angle([0., -0., complex(0., -0.), complex(-0., -0.)]) # convention
1645 array([ 0. , 3.14159265, -0. , -3.14159265])
1647 """
1648 z = asanyarray(z)
1649 if issubclass(z.dtype.type, _nx.complexfloating):
1650 zimag = z.imag
1651 zreal = z.real
1652 else:
1653 zimag = 0
1654 zreal = z
1656 a = arctan2(zimag, zreal)
1657 if deg:
1658 a *= 180/pi
1659 return a
1662def _unwrap_dispatcher(p, discont=None, axis=None, *, period=None):
1663 return (p,)
1666@array_function_dispatch(_unwrap_dispatcher)
1667def unwrap(p, discont=None, axis=-1, *, period=2*pi):
1668 r"""
1669 Unwrap by taking the complement of large deltas with respect to the period.
1671 This unwraps a signal `p` by changing elements which have an absolute
1672 difference from their predecessor of more than ``max(discont, period/2)``
1673 to their `period`-complementary values.
1675 For the default case where `period` is :math:`2\pi` and `discont` is
1676 :math:`\pi`, this unwraps a radian phase `p` such that adjacent differences
1677 are never greater than :math:`\pi` by adding :math:`2k\pi` for some
1678 integer :math:`k`.
1680 Parameters
1681 ----------
1682 p : array_like
1683 Input array.
1684 discont : float, optional
1685 Maximum discontinuity between values, default is ``period/2``.
1686 Values below ``period/2`` are treated as if they were ``period/2``.
1687 To have an effect different from the default, `discont` should be
1688 larger than ``period/2``.
1689 axis : int, optional
1690 Axis along which unwrap will operate, default is the last axis.
1691 period : float, optional
1692 Size of the range over which the input wraps. By default, it is
1693 ``2 pi``.
1695 .. versionadded:: 1.21.0
1697 Returns
1698 -------
1699 out : ndarray
1700 Output array.
1702 See Also
1703 --------
1704 rad2deg, deg2rad
1706 Notes
1707 -----
1708 If the discontinuity in `p` is smaller than ``period/2``,
1709 but larger than `discont`, no unwrapping is done because taking
1710 the complement would only make the discontinuity larger.
1712 Examples
1713 --------
1714 >>> phase = np.linspace(0, np.pi, num=5)
1715 >>> phase[3:] += np.pi
1716 >>> phase
1717 array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531]) # may vary
1718 >>> np.unwrap(phase)
1719 array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ]) # may vary
1720 >>> np.unwrap([0, 1, 2, -1, 0], period=4)
1721 array([0, 1, 2, 3, 4])
1722 >>> np.unwrap([ 1, 2, 3, 4, 5, 6, 1, 2, 3], period=6)
1723 array([1, 2, 3, 4, 5, 6, 7, 8, 9])
1724 >>> np.unwrap([2, 3, 4, 5, 2, 3, 4, 5], period=4)
1725 array([2, 3, 4, 5, 6, 7, 8, 9])
1726 >>> phase_deg = np.mod(np.linspace(0 ,720, 19), 360) - 180
1727 >>> np.unwrap(phase_deg, period=360)
1728 array([-180., -140., -100., -60., -20., 20., 60., 100., 140.,
1729 180., 220., 260., 300., 340., 380., 420., 460., 500.,
1730 540.])
1731 """
1732 p = asarray(p)
1733 nd = p.ndim
1734 dd = diff(p, axis=axis)
1735 if discont is None:
1736 discont = period/2
1737 slice1 = [slice(None, None)]*nd # full slices
1738 slice1[axis] = slice(1, None)
1739 slice1 = tuple(slice1)
1740 dtype = np.result_type(dd, period)
1741 if _nx.issubdtype(dtype, _nx.integer):
1742 interval_high, rem = divmod(period, 2)
1743 boundary_ambiguous = rem == 0
1744 else:
1745 interval_high = period / 2
1746 boundary_ambiguous = True
1747 interval_low = -interval_high
1748 ddmod = mod(dd - interval_low, period) + interval_low
1749 if boundary_ambiguous:
1750 # for `mask = (abs(dd) == period/2)`, the above line made
1751 # `ddmod[mask] == -period/2`. correct these such that
1752 # `ddmod[mask] == sign(dd[mask])*period/2`.
1753 _nx.copyto(ddmod, interval_high,
1754 where=(ddmod == interval_low) & (dd > 0))
1755 ph_correct = ddmod - dd
1756 _nx.copyto(ph_correct, 0, where=abs(dd) < discont)
1757 up = array(p, copy=True, dtype=dtype)
1758 up[slice1] = p[slice1] + ph_correct.cumsum(axis)
1759 return up
1762def _sort_complex(a):
1763 return (a,)
1766@array_function_dispatch(_sort_complex)
1767def sort_complex(a):
1768 """
1769 Sort a complex array using the real part first, then the imaginary part.
1771 Parameters
1772 ----------
1773 a : array_like
1774 Input array
1776 Returns
1777 -------
1778 out : complex ndarray
1779 Always returns a sorted complex array.
1781 Examples
1782 --------
1783 >>> np.sort_complex([5, 3, 6, 2, 1])
1784 array([1.+0.j, 2.+0.j, 3.+0.j, 5.+0.j, 6.+0.j])
1786 >>> np.sort_complex([1 + 2j, 2 - 1j, 3 - 2j, 3 - 3j, 3 + 5j])
1787 array([1.+2.j, 2.-1.j, 3.-3.j, 3.-2.j, 3.+5.j])
1789 """
1790 b = array(a, copy=True)
1791 b.sort()
1792 if not issubclass(b.dtype.type, _nx.complexfloating):
1793 if b.dtype.char in 'bhBH':
1794 return b.astype('F')
1795 elif b.dtype.char == 'g':
1796 return b.astype('G')
1797 else:
1798 return b.astype('D')
1799 else:
1800 return b
1803def _trim_zeros(filt, trim=None):
1804 return (filt,)
1807@array_function_dispatch(_trim_zeros)
1808def trim_zeros(filt, trim='fb'):
1809 """
1810 Trim the leading and/or trailing zeros from a 1-D array or sequence.
1812 Parameters
1813 ----------
1814 filt : 1-D array or sequence
1815 Input array.
1816 trim : str, optional
1817 A string with 'f' representing trim from front and 'b' to trim from
1818 back. Default is 'fb', trim zeros from both front and back of the
1819 array.
1821 Returns
1822 -------
1823 trimmed : 1-D array or sequence
1824 The result of trimming the input. The input data type is preserved.
1826 Examples
1827 --------
1828 >>> a = np.array((0, 0, 0, 1, 2, 3, 0, 2, 1, 0))
1829 >>> np.trim_zeros(a)
1830 array([1, 2, 3, 0, 2, 1])
1832 >>> np.trim_zeros(a, 'b')
1833 array([0, 0, 0, ..., 0, 2, 1])
1835 The input data type is preserved, list/tuple in means list/tuple out.
1837 >>> np.trim_zeros([0, 1, 2, 0])
1838 [1, 2]
1840 """
1842 first = 0
1843 trim = trim.upper()
1844 if 'F' in trim:
1845 for i in filt:
1846 if i != 0.:
1847 break
1848 else:
1849 first = first + 1
1850 last = len(filt)
1851 if 'B' in trim:
1852 for i in filt[::-1]:
1853 if i != 0.:
1854 break
1855 else:
1856 last = last - 1
1857 return filt[first:last]
1860def _extract_dispatcher(condition, arr):
1861 return (condition, arr)
1864@array_function_dispatch(_extract_dispatcher)
1865def extract(condition, arr):
1866 """
1867 Return the elements of an array that satisfy some condition.
1869 This is equivalent to ``np.compress(ravel(condition), ravel(arr))``. If
1870 `condition` is boolean ``np.extract`` is equivalent to ``arr[condition]``.
1872 Note that `place` does the exact opposite of `extract`.
1874 Parameters
1875 ----------
1876 condition : array_like
1877 An array whose nonzero or True entries indicate the elements of `arr`
1878 to extract.
1879 arr : array_like
1880 Input array of the same size as `condition`.
1882 Returns
1883 -------
1884 extract : ndarray
1885 Rank 1 array of values from `arr` where `condition` is True.
1887 See Also
1888 --------
1889 take, put, copyto, compress, place
1891 Examples
1892 --------
1893 >>> arr = np.arange(12).reshape((3, 4))
1894 >>> arr
1895 array([[ 0, 1, 2, 3],
1896 [ 4, 5, 6, 7],
1897 [ 8, 9, 10, 11]])
1898 >>> condition = np.mod(arr, 3)==0
1899 >>> condition
1900 array([[ True, False, False, True],
1901 [False, False, True, False],
1902 [False, True, False, False]])
1903 >>> np.extract(condition, arr)
1904 array([0, 3, 6, 9])
1907 If `condition` is boolean:
1909 >>> arr[condition]
1910 array([0, 3, 6, 9])
1912 """
1913 return _nx.take(ravel(arr), nonzero(ravel(condition))[0])
1916def _place_dispatcher(arr, mask, vals):
1917 return (arr, mask, vals)
1920@array_function_dispatch(_place_dispatcher)
1921def place(arr, mask, vals):
1922 """
1923 Change elements of an array based on conditional and input values.
1925 Similar to ``np.copyto(arr, vals, where=mask)``, the difference is that
1926 `place` uses the first N elements of `vals`, where N is the number of
1927 True values in `mask`, while `copyto` uses the elements where `mask`
1928 is True.
1930 Note that `extract` does the exact opposite of `place`.
1932 Parameters
1933 ----------
1934 arr : ndarray
1935 Array to put data into.
1936 mask : array_like
1937 Boolean mask array. Must have the same size as `a`.
1938 vals : 1-D sequence
1939 Values to put into `a`. Only the first N elements are used, where
1940 N is the number of True values in `mask`. If `vals` is smaller
1941 than N, it will be repeated, and if elements of `a` are to be masked,
1942 this sequence must be non-empty.
1944 See Also
1945 --------
1946 copyto, put, take, extract
1948 Examples
1949 --------
1950 >>> arr = np.arange(6).reshape(2, 3)
1951 >>> np.place(arr, arr>2, [44, 55])
1952 >>> arr
1953 array([[ 0, 1, 2],
1954 [44, 55, 44]])
1956 """
1957 return _place(arr, mask, vals)
1960def disp(mesg, device=None, linefeed=True):
1961 """
1962 Display a message on a device.
1964 .. deprecated:: 2.0
1965 Use your own printing function instead.
1967 Parameters
1968 ----------
1969 mesg : str
1970 Message to display.
1971 device : object
1972 Device to write message. If None, defaults to ``sys.stdout`` which is
1973 very similar to ``print``. `device` needs to have ``write()`` and
1974 ``flush()`` methods.
1975 linefeed : bool, optional
1976 Option whether to print a line feed or not. Defaults to True.
1978 Raises
1979 ------
1980 AttributeError
1981 If `device` does not have a ``write()`` or ``flush()`` method.
1983 Examples
1984 --------
1985 Besides ``sys.stdout``, a file-like object can also be used as it has
1986 both required methods:
1988 >>> from io import StringIO
1989 >>> buf = StringIO()
1990 >>> np.disp('"Display" in a file', device=buf)
1991 >>> buf.getvalue()
1992 '"Display" in a file\\n'
1994 """
1996 # Deprecated in NumPy 2.0, 2023-07-11
1997 warnings.warn(
1998 "`disp` is deprecated, "
1999 "use your own printing function instead. "
2000 "(deprecated in NumPy 2.0)",
2001 DeprecationWarning,
2002 stacklevel=2
2003 )
2005 if device is None:
2006 device = sys.stdout
2007 if linefeed:
2008 device.write('%s\n' % mesg)
2009 else:
2010 device.write('%s' % mesg)
2011 device.flush()
2012 return
2015# See https://docs.scipy.org/doc/numpy/reference/c-api.generalized-ufuncs.html
2016_DIMENSION_NAME = r'\w+'
2017_CORE_DIMENSION_LIST = '(?:{0:}(?:,{0:})*)?'.format(_DIMENSION_NAME)
2018_ARGUMENT = r'\({}\)'.format(_CORE_DIMENSION_LIST)
2019_ARGUMENT_LIST = '{0:}(?:,{0:})*'.format(_ARGUMENT)
2020_SIGNATURE = '^{0:}->{0:}$'.format(_ARGUMENT_LIST)
2023def _parse_gufunc_signature(signature):
2024 """
2025 Parse string signatures for a generalized universal function.
2027 Arguments
2028 ---------
2029 signature : string
2030 Generalized universal function signature, e.g., ``(m,n),(n,p)->(m,p)``
2031 for ``np.matmul``.
2033 Returns
2034 -------
2035 Tuple of input and output core dimensions parsed from the signature, each
2036 of the form List[Tuple[str, ...]].
2037 """
2038 signature = re.sub(r'\s+', '', signature)
2040 if not re.match(_SIGNATURE, signature):
2041 raise ValueError(
2042 'not a valid gufunc signature: {}'.format(signature))
2043 return tuple([tuple(re.findall(_DIMENSION_NAME, arg))
2044 for arg in re.findall(_ARGUMENT, arg_list)]
2045 for arg_list in signature.split('->'))
2048def _update_dim_sizes(dim_sizes, arg, core_dims):
2049 """
2050 Incrementally check and update core dimension sizes for a single argument.
2052 Arguments
2053 ---------
2054 dim_sizes : Dict[str, int]
2055 Sizes of existing core dimensions. Will be updated in-place.
2056 arg : ndarray
2057 Argument to examine.
2058 core_dims : Tuple[str, ...]
2059 Core dimensions for this argument.
2060 """
2061 if not core_dims:
2062 return
2064 num_core_dims = len(core_dims)
2065 if arg.ndim < num_core_dims:
2066 raise ValueError(
2067 '%d-dimensional argument does not have enough '
2068 'dimensions for all core dimensions %r'
2069 % (arg.ndim, core_dims))
2071 core_shape = arg.shape[-num_core_dims:]
2072 for dim, size in zip(core_dims, core_shape):
2073 if dim in dim_sizes:
2074 if size != dim_sizes[dim]:
2075 raise ValueError(
2076 'inconsistent size for core dimension %r: %r vs %r'
2077 % (dim, size, dim_sizes[dim]))
2078 else:
2079 dim_sizes[dim] = size
2082def _parse_input_dimensions(args, input_core_dims):
2083 """
2084 Parse broadcast and core dimensions for vectorize with a signature.
2086 Arguments
2087 ---------
2088 args : Tuple[ndarray, ...]
2089 Tuple of input arguments to examine.
2090 input_core_dims : List[Tuple[str, ...]]
2091 List of core dimensions corresponding to each input.
2093 Returns
2094 -------
2095 broadcast_shape : Tuple[int, ...]
2096 Common shape to broadcast all non-core dimensions to.
2097 dim_sizes : Dict[str, int]
2098 Common sizes for named core dimensions.
2099 """
2100 broadcast_args = []
2101 dim_sizes = {}
2102 for arg, core_dims in zip(args, input_core_dims):
2103 _update_dim_sizes(dim_sizes, arg, core_dims)
2104 ndim = arg.ndim - len(core_dims)
2105 dummy_array = np.lib.stride_tricks.as_strided(0, arg.shape[:ndim])
2106 broadcast_args.append(dummy_array)
2107 broadcast_shape = np.lib._stride_tricks_impl._broadcast_shape(
2108 *broadcast_args
2109 )
2110 return broadcast_shape, dim_sizes
2113def _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims):
2114 """Helper for calculating broadcast shapes with core dimensions."""
2115 return [broadcast_shape + tuple(dim_sizes[dim] for dim in core_dims)
2116 for core_dims in list_of_core_dims]
2119def _create_arrays(broadcast_shape, dim_sizes, list_of_core_dims, dtypes,
2120 results=None):
2121 """Helper for creating output arrays in vectorize."""
2122 shapes = _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims)
2123 if dtypes is None:
2124 dtypes = [None] * len(shapes)
2125 if results is None:
2126 arrays = tuple(np.empty(shape=shape, dtype=dtype)
2127 for shape, dtype in zip(shapes, dtypes))
2128 else:
2129 arrays = tuple(np.empty_like(result, shape=shape, dtype=dtype)
2130 for result, shape, dtype
2131 in zip(results, shapes, dtypes))
2132 return arrays
2135@set_module('numpy')
2136class vectorize:
2137 """
2138 vectorize(pyfunc=np._NoValue, otypes=None, doc=None, excluded=None,
2139 cache=False, signature=None)
2141 Returns an object that acts like pyfunc, but takes arrays as input.
2143 Define a vectorized function which takes a nested sequence of objects or
2144 numpy arrays as inputs and returns a single numpy array or a tuple of numpy
2145 arrays. The vectorized function evaluates `pyfunc` over successive tuples
2146 of the input arrays like the python map function, except it uses the
2147 broadcasting rules of numpy.
2149 The data type of the output of `vectorized` is determined by calling
2150 the function with the first element of the input. This can be avoided
2151 by specifying the `otypes` argument.
2153 Parameters
2154 ----------
2155 pyfunc : callable, optional
2156 A python function or method.
2157 Can be omitted to produce a decorator with keyword arguments.
2158 otypes : str or list of dtypes, optional
2159 The output data type. It must be specified as either a string of
2160 typecode characters or a list of data type specifiers. There should
2161 be one data type specifier for each output.
2162 doc : str, optional
2163 The docstring for the function. If None, the docstring will be the
2164 ``pyfunc.__doc__``.
2165 excluded : set, optional
2166 Set of strings or integers representing the positional or keyword
2167 arguments for which the function will not be vectorized. These will be
2168 passed directly to `pyfunc` unmodified.
2170 .. versionadded:: 1.7.0
2172 cache : bool, optional
2173 If `True`, then cache the first function call that determines the number
2174 of outputs if `otypes` is not provided.
2176 .. versionadded:: 1.7.0
2178 signature : string, optional
2179 Generalized universal function signature, e.g., ``(m,n),(n)->(m)`` for
2180 vectorized matrix-vector multiplication. If provided, ``pyfunc`` will
2181 be called with (and expected to return) arrays with shapes given by the
2182 size of corresponding core dimensions. By default, ``pyfunc`` is
2183 assumed to take scalars as input and output.
2185 .. versionadded:: 1.12.0
2187 Returns
2188 -------
2189 out : callable
2190 A vectorized function if ``pyfunc`` was provided,
2191 a decorator otherwise.
2193 See Also
2194 --------
2195 frompyfunc : Takes an arbitrary Python function and returns a ufunc
2197 Notes
2198 -----
2199 The `vectorize` function is provided primarily for convenience, not for
2200 performance. The implementation is essentially a for loop.
2202 If `otypes` is not specified, then a call to the function with the
2203 first argument will be used to determine the number of outputs. The
2204 results of this call will be cached if `cache` is `True` to prevent
2205 calling the function twice. However, to implement the cache, the
2206 original function must be wrapped which will slow down subsequent
2207 calls, so only do this if your function is expensive.
2209 The new keyword argument interface and `excluded` argument support
2210 further degrades performance.
2212 References
2213 ----------
2214 .. [1] :doc:`/reference/c-api/generalized-ufuncs`
2216 Examples
2217 --------
2218 >>> def myfunc(a, b):
2219 ... "Return a-b if a>b, otherwise return a+b"
2220 ... if a > b:
2221 ... return a - b
2222 ... else:
2223 ... return a + b
2225 >>> vfunc = np.vectorize(myfunc)
2226 >>> vfunc([1, 2, 3, 4], 2)
2227 array([3, 4, 1, 2])
2229 The docstring is taken from the input function to `vectorize` unless it
2230 is specified:
2232 >>> vfunc.__doc__
2233 'Return a-b if a>b, otherwise return a+b'
2234 >>> vfunc = np.vectorize(myfunc, doc='Vectorized `myfunc`')
2235 >>> vfunc.__doc__
2236 'Vectorized `myfunc`'
2238 The output type is determined by evaluating the first element of the input,
2239 unless it is specified:
2241 >>> out = vfunc([1, 2, 3, 4], 2)
2242 >>> type(out[0])
2243 <class 'numpy.int64'>
2244 >>> vfunc = np.vectorize(myfunc, otypes=[float])
2245 >>> out = vfunc([1, 2, 3, 4], 2)
2246 >>> type(out[0])
2247 <class 'numpy.float64'>
2249 The `excluded` argument can be used to prevent vectorizing over certain
2250 arguments. This can be useful for array-like arguments of a fixed length
2251 such as the coefficients for a polynomial as in `polyval`:
2253 >>> def mypolyval(p, x):
2254 ... _p = list(p)
2255 ... res = _p.pop(0)
2256 ... while _p:
2257 ... res = res*x + _p.pop(0)
2258 ... return res
2259 >>> vpolyval = np.vectorize(mypolyval, excluded=['p'])
2260 >>> vpolyval(p=[1, 2, 3], x=[0, 1])
2261 array([3, 6])
2263 Positional arguments may also be excluded by specifying their position:
2265 >>> vpolyval.excluded.add(0)
2266 >>> vpolyval([1, 2, 3], x=[0, 1])
2267 array([3, 6])
2269 The `signature` argument allows for vectorizing functions that act on
2270 non-scalar arrays of fixed length. For example, you can use it for a
2271 vectorized calculation of Pearson correlation coefficient and its p-value:
2273 >>> import scipy.stats
2274 >>> pearsonr = np.vectorize(scipy.stats.pearsonr,
2275 ... signature='(n),(n)->(),()')
2276 >>> pearsonr([[0, 1, 2, 3]], [[1, 2, 3, 4], [4, 3, 2, 1]])
2277 (array([ 1., -1.]), array([ 0., 0.]))
2279 Or for a vectorized convolution:
2281 >>> convolve = np.vectorize(np.convolve, signature='(n),(m)->(k)')
2282 >>> convolve(np.eye(4), [1, 2, 1])
2283 array([[1., 2., 1., 0., 0., 0.],
2284 [0., 1., 2., 1., 0., 0.],
2285 [0., 0., 1., 2., 1., 0.],
2286 [0., 0., 0., 1., 2., 1.]])
2288 Decorator syntax is supported. The decorator can be called as
2289 a function to provide keyword arguments:
2291 >>> @np.vectorize
2292 ... def identity(x):
2293 ... return x
2294 ...
2295 >>> identity([0, 1, 2])
2296 array([0, 1, 2])
2297 >>> @np.vectorize(otypes=[float])
2298 ... def as_float(x):
2299 ... return x
2300 ...
2301 >>> as_float([0, 1, 2])
2302 array([0., 1., 2.])
2303 """
2304 def __init__(self, pyfunc=np._NoValue, otypes=None, doc=None,
2305 excluded=None, cache=False, signature=None):
2307 if (pyfunc != np._NoValue) and (not callable(pyfunc)):
2308 #Splitting the error message to keep
2309 #the length below 79 characters.
2310 part1 = "When used as a decorator, "
2311 part2 = "only accepts keyword arguments."
2312 raise TypeError(part1 + part2)
2314 self.pyfunc = pyfunc
2315 self.cache = cache
2316 self.signature = signature
2317 if pyfunc != np._NoValue and hasattr(pyfunc, '__name__'):
2318 self.__name__ = pyfunc.__name__
2320 self._ufunc = {} # Caching to improve default performance
2321 self._doc = None
2322 self.__doc__ = doc
2323 if doc is None and hasattr(pyfunc, '__doc__'):
2324 self.__doc__ = pyfunc.__doc__
2325 else:
2326 self._doc = doc
2328 if isinstance(otypes, str):
2329 for char in otypes:
2330 if char not in typecodes['All']:
2331 raise ValueError("Invalid otype specified: %s" % (char,))
2332 elif iterable(otypes):
2333 otypes = [_nx.dtype(x) for x in otypes]
2334 elif otypes is not None:
2335 raise ValueError("Invalid otype specification")
2336 self.otypes = otypes
2338 # Excluded variable support
2339 if excluded is None:
2340 excluded = set()
2341 self.excluded = set(excluded)
2343 if signature is not None:
2344 self._in_and_out_core_dims = _parse_gufunc_signature(signature)
2345 else:
2346 self._in_and_out_core_dims = None
2348 def _init_stage_2(self, pyfunc, *args, **kwargs):
2349 self.__name__ = pyfunc.__name__
2350 self.pyfunc = pyfunc
2351 if self._doc is None:
2352 self.__doc__ = pyfunc.__doc__
2353 else:
2354 self.__doc__ = self._doc
2356 def _call_as_normal(self, *args, **kwargs):
2357 """
2358 Return arrays with the results of `pyfunc` broadcast (vectorized) over
2359 `args` and `kwargs` not in `excluded`.
2360 """
2361 excluded = self.excluded
2362 if not kwargs and not excluded:
2363 func = self.pyfunc
2364 vargs = args
2365 else:
2366 # The wrapper accepts only positional arguments: we use `names` and
2367 # `inds` to mutate `the_args` and `kwargs` to pass to the original
2368 # function.
2369 nargs = len(args)
2371 names = [_n for _n in kwargs if _n not in excluded]
2372 inds = [_i for _i in range(nargs) if _i not in excluded]
2373 the_args = list(args)
2375 def func(*vargs):
2376 for _n, _i in enumerate(inds):
2377 the_args[_i] = vargs[_n]
2378 kwargs.update(zip(names, vargs[len(inds):]))
2379 return self.pyfunc(*the_args, **kwargs)
2381 vargs = [args[_i] for _i in inds]
2382 vargs.extend([kwargs[_n] for _n in names])
2384 return self._vectorize_call(func=func, args=vargs)
2386 def __call__(self, *args, **kwargs):
2387 if self.pyfunc is np._NoValue:
2388 self._init_stage_2(*args, **kwargs)
2389 return self
2391 return self._call_as_normal(*args, **kwargs)
2393 def _get_ufunc_and_otypes(self, func, args):
2394 """Return (ufunc, otypes)."""
2395 # frompyfunc will fail if args is empty
2396 if not args:
2397 raise ValueError('args can not be empty')
2399 if self.otypes is not None:
2400 otypes = self.otypes
2402 # self._ufunc is a dictionary whose keys are the number of
2403 # arguments (i.e. len(args)) and whose values are ufuncs created
2404 # by frompyfunc. len(args) can be different for different calls if
2405 # self.pyfunc has parameters with default values. We only use the
2406 # cache when func is self.pyfunc, which occurs when the call uses
2407 # only positional arguments and no arguments are excluded.
2409 nin = len(args)
2410 nout = len(self.otypes)
2411 if func is not self.pyfunc or nin not in self._ufunc:
2412 ufunc = frompyfunc(func, nin, nout)
2413 else:
2414 ufunc = None # We'll get it from self._ufunc
2415 if func is self.pyfunc:
2416 ufunc = self._ufunc.setdefault(nin, ufunc)
2417 else:
2418 # Get number of outputs and output types by calling the function on
2419 # the first entries of args. We also cache the result to prevent
2420 # the subsequent call when the ufunc is evaluated.
2421 # Assumes that ufunc first evaluates the 0th elements in the input
2422 # arrays (the input values are not checked to ensure this)
2423 args = [asarray(arg) for arg in args]
2424 if builtins.any(arg.size == 0 for arg in args):
2425 raise ValueError('cannot call `vectorize` on size 0 inputs '
2426 'unless `otypes` is set')
2428 inputs = [arg.flat[0] for arg in args]
2429 outputs = func(*inputs)
2431 # Performance note: profiling indicates that -- for simple
2432 # functions at least -- this wrapping can almost double the
2433 # execution time.
2434 # Hence we make it optional.
2435 if self.cache:
2436 _cache = [outputs]
2438 def _func(*vargs):
2439 if _cache:
2440 return _cache.pop()
2441 else:
2442 return func(*vargs)
2443 else:
2444 _func = func
2446 if isinstance(outputs, tuple):
2447 nout = len(outputs)
2448 else:
2449 nout = 1
2450 outputs = (outputs,)
2452 otypes = ''.join([asarray(outputs[_k]).dtype.char
2453 for _k in range(nout)])
2455 # Performance note: profiling indicates that creating the ufunc is
2456 # not a significant cost compared with wrapping so it seems not
2457 # worth trying to cache this.
2458 ufunc = frompyfunc(_func, len(args), nout)
2460 return ufunc, otypes
2462 def _vectorize_call(self, func, args):
2463 """Vectorized call to `func` over positional `args`."""
2464 if self.signature is not None:
2465 res = self._vectorize_call_with_signature(func, args)
2466 elif not args:
2467 res = func()
2468 else:
2469 ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args)
2471 # Convert args to object arrays first
2472 inputs = [asanyarray(a, dtype=object) for a in args]
2474 outputs = ufunc(*inputs)
2476 if ufunc.nout == 1:
2477 res = asanyarray(outputs, dtype=otypes[0])
2478 else:
2479 res = tuple([asanyarray(x, dtype=t)
2480 for x, t in zip(outputs, otypes)])
2481 return res
2483 def _vectorize_call_with_signature(self, func, args):
2484 """Vectorized call over positional arguments with a signature."""
2485 input_core_dims, output_core_dims = self._in_and_out_core_dims
2487 if len(args) != len(input_core_dims):
2488 raise TypeError('wrong number of positional arguments: '
2489 'expected %r, got %r'
2490 % (len(input_core_dims), len(args)))
2491 args = tuple(asanyarray(arg) for arg in args)
2493 broadcast_shape, dim_sizes = _parse_input_dimensions(
2494 args, input_core_dims)
2495 input_shapes = _calculate_shapes(broadcast_shape, dim_sizes,
2496 input_core_dims)
2497 args = [np.broadcast_to(arg, shape, subok=True)
2498 for arg, shape in zip(args, input_shapes)]
2500 outputs = None
2501 otypes = self.otypes
2502 nout = len(output_core_dims)
2504 for index in np.ndindex(*broadcast_shape):
2505 results = func(*(arg[index] for arg in args))
2507 n_results = len(results) if isinstance(results, tuple) else 1
2509 if nout != n_results:
2510 raise ValueError(
2511 'wrong number of outputs from pyfunc: expected %r, got %r'
2512 % (nout, n_results))
2514 if nout == 1:
2515 results = (results,)
2517 if outputs is None:
2518 for result, core_dims in zip(results, output_core_dims):
2519 _update_dim_sizes(dim_sizes, result, core_dims)
2521 outputs = _create_arrays(broadcast_shape, dim_sizes,
2522 output_core_dims, otypes, results)
2524 for output, result in zip(outputs, results):
2525 output[index] = result
2527 if outputs is None:
2528 # did not call the function even once
2529 if otypes is None:
2530 raise ValueError('cannot call `vectorize` on size 0 inputs '
2531 'unless `otypes` is set')
2532 if builtins.any(dim not in dim_sizes
2533 for dims in output_core_dims
2534 for dim in dims):
2535 raise ValueError('cannot call `vectorize` with a signature '
2536 'including new output dimensions on size 0 '
2537 'inputs')
2538 outputs = _create_arrays(broadcast_shape, dim_sizes,
2539 output_core_dims, otypes)
2541 return outputs[0] if nout == 1 else outputs
2544def _cov_dispatcher(m, y=None, rowvar=None, bias=None, ddof=None,
2545 fweights=None, aweights=None, *, dtype=None):
2546 return (m, y, fweights, aweights)
2549@array_function_dispatch(_cov_dispatcher)
2550def cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None,
2551 aweights=None, *, dtype=None):
2552 """
2553 Estimate a covariance matrix, given data and weights.
2555 Covariance indicates the level to which two variables vary together.
2556 If we examine N-dimensional samples, :math:`X = [x_1, x_2, ... x_N]^T`,
2557 then the covariance matrix element :math:`C_{ij}` is the covariance of
2558 :math:`x_i` and :math:`x_j`. The element :math:`C_{ii}` is the variance
2559 of :math:`x_i`.
2561 See the notes for an outline of the algorithm.
2563 Parameters
2564 ----------
2565 m : array_like
2566 A 1-D or 2-D array containing multiple variables and observations.
2567 Each row of `m` represents a variable, and each column a single
2568 observation of all those variables. Also see `rowvar` below.
2569 y : array_like, optional
2570 An additional set of variables and observations. `y` has the same form
2571 as that of `m`.
2572 rowvar : bool, optional
2573 If `rowvar` is True (default), then each row represents a
2574 variable, with observations in the columns. Otherwise, the relationship
2575 is transposed: each column represents a variable, while the rows
2576 contain observations.
2577 bias : bool, optional
2578 Default normalization (False) is by ``(N - 1)``, where ``N`` is the
2579 number of observations given (unbiased estimate). If `bias` is True,
2580 then normalization is by ``N``. These values can be overridden by using
2581 the keyword ``ddof`` in numpy versions >= 1.5.
2582 ddof : int, optional
2583 If not ``None`` the default value implied by `bias` is overridden.
2584 Note that ``ddof=1`` will return the unbiased estimate, even if both
2585 `fweights` and `aweights` are specified, and ``ddof=0`` will return
2586 the simple average. See the notes for the details. The default value
2587 is ``None``.
2589 .. versionadded:: 1.5
2590 fweights : array_like, int, optional
2591 1-D array of integer frequency weights; the number of times each
2592 observation vector should be repeated.
2594 .. versionadded:: 1.10
2595 aweights : array_like, optional
2596 1-D array of observation vector weights. These relative weights are
2597 typically large for observations considered "important" and smaller for
2598 observations considered less "important". If ``ddof=0`` the array of
2599 weights can be used to assign probabilities to observation vectors.
2601 .. versionadded:: 1.10
2602 dtype : data-type, optional
2603 Data-type of the result. By default, the return data-type will have
2604 at least `numpy.float64` precision.
2606 .. versionadded:: 1.20
2608 Returns
2609 -------
2610 out : ndarray
2611 The covariance matrix of the variables.
2613 See Also
2614 --------
2615 corrcoef : Normalized covariance matrix
2617 Notes
2618 -----
2619 Assume that the observations are in the columns of the observation
2620 array `m` and let ``f = fweights`` and ``a = aweights`` for brevity. The
2621 steps to compute the weighted covariance are as follows::
2623 >>> m = np.arange(10, dtype=np.float64)
2624 >>> f = np.arange(10) * 2
2625 >>> a = np.arange(10) ** 2.
2626 >>> ddof = 1
2627 >>> w = f * a
2628 >>> v1 = np.sum(w)
2629 >>> v2 = np.sum(w * a)
2630 >>> m -= np.sum(m * w, axis=None, keepdims=True) / v1
2631 >>> cov = np.dot(m * w, m.T) * v1 / (v1**2 - ddof * v2)
2633 Note that when ``a == 1``, the normalization factor
2634 ``v1 / (v1**2 - ddof * v2)`` goes over to ``1 / (np.sum(f) - ddof)``
2635 as it should.
2637 Examples
2638 --------
2639 Consider two variables, :math:`x_0` and :math:`x_1`, which
2640 correlate perfectly, but in opposite directions:
2642 >>> x = np.array([[0, 2], [1, 1], [2, 0]]).T
2643 >>> x
2644 array([[0, 1, 2],
2645 [2, 1, 0]])
2647 Note how :math:`x_0` increases while :math:`x_1` decreases. The covariance
2648 matrix shows this clearly:
2650 >>> np.cov(x)
2651 array([[ 1., -1.],
2652 [-1., 1.]])
2654 Note that element :math:`C_{0,1}`, which shows the correlation between
2655 :math:`x_0` and :math:`x_1`, is negative.
2657 Further, note how `x` and `y` are combined:
2659 >>> x = [-2.1, -1, 4.3]
2660 >>> y = [3, 1.1, 0.12]
2661 >>> X = np.stack((x, y), axis=0)
2662 >>> np.cov(X)
2663 array([[11.71 , -4.286 ], # may vary
2664 [-4.286 , 2.144133]])
2665 >>> np.cov(x, y)
2666 array([[11.71 , -4.286 ], # may vary
2667 [-4.286 , 2.144133]])
2668 >>> np.cov(x)
2669 array(11.71)
2671 """
2672 # Check inputs
2673 if ddof is not None and ddof != int(ddof):
2674 raise ValueError(
2675 "ddof must be integer")
2677 # Handles complex arrays too
2678 m = np.asarray(m)
2679 if m.ndim > 2:
2680 raise ValueError("m has more than 2 dimensions")
2682 if y is not None:
2683 y = np.asarray(y)
2684 if y.ndim > 2:
2685 raise ValueError("y has more than 2 dimensions")
2687 if dtype is None:
2688 if y is None:
2689 dtype = np.result_type(m, np.float64)
2690 else:
2691 dtype = np.result_type(m, y, np.float64)
2693 X = array(m, ndmin=2, dtype=dtype)
2694 if not rowvar and X.shape[0] != 1:
2695 X = X.T
2696 if X.shape[0] == 0:
2697 return np.array([]).reshape(0, 0)
2698 if y is not None:
2699 y = array(y, copy=None, ndmin=2, dtype=dtype)
2700 if not rowvar and y.shape[0] != 1:
2701 y = y.T
2702 X = np.concatenate((X, y), axis=0)
2704 if ddof is None:
2705 if bias == 0:
2706 ddof = 1
2707 else:
2708 ddof = 0
2710 # Get the product of frequencies and weights
2711 w = None
2712 if fweights is not None:
2713 fweights = np.asarray(fweights, dtype=float)
2714 if not np.all(fweights == np.around(fweights)):
2715 raise TypeError(
2716 "fweights must be integer")
2717 if fweights.ndim > 1:
2718 raise RuntimeError(
2719 "cannot handle multidimensional fweights")
2720 if fweights.shape[0] != X.shape[1]:
2721 raise RuntimeError(
2722 "incompatible numbers of samples and fweights")
2723 if any(fweights < 0):
2724 raise ValueError(
2725 "fweights cannot be negative")
2726 w = fweights
2727 if aweights is not None:
2728 aweights = np.asarray(aweights, dtype=float)
2729 if aweights.ndim > 1:
2730 raise RuntimeError(
2731 "cannot handle multidimensional aweights")
2732 if aweights.shape[0] != X.shape[1]:
2733 raise RuntimeError(
2734 "incompatible numbers of samples and aweights")
2735 if any(aweights < 0):
2736 raise ValueError(
2737 "aweights cannot be negative")
2738 if w is None:
2739 w = aweights
2740 else:
2741 w *= aweights
2743 avg, w_sum = average(X, axis=1, weights=w, returned=True)
2744 w_sum = w_sum[0]
2746 # Determine the normalization
2747 if w is None:
2748 fact = X.shape[1] - ddof
2749 elif ddof == 0:
2750 fact = w_sum
2751 elif aweights is None:
2752 fact = w_sum - ddof
2753 else:
2754 fact = w_sum - ddof*sum(w*aweights)/w_sum
2756 if fact <= 0:
2757 warnings.warn("Degrees of freedom <= 0 for slice",
2758 RuntimeWarning, stacklevel=2)
2759 fact = 0.0
2761 X -= avg[:, None]
2762 if w is None:
2763 X_T = X.T
2764 else:
2765 X_T = (X*w).T
2766 c = dot(X, X_T.conj())
2767 c *= np.true_divide(1, fact)
2768 return c.squeeze()
2771def _corrcoef_dispatcher(x, y=None, rowvar=None, bias=None, ddof=None, *,
2772 dtype=None):
2773 return (x, y)
2776@array_function_dispatch(_corrcoef_dispatcher)
2777def corrcoef(x, y=None, rowvar=True, bias=np._NoValue, ddof=np._NoValue, *,
2778 dtype=None):
2779 """
2780 Return Pearson product-moment correlation coefficients.
2782 Please refer to the documentation for `cov` for more detail. The
2783 relationship between the correlation coefficient matrix, `R`, and the
2784 covariance matrix, `C`, is
2786 .. math:: R_{ij} = \\frac{ C_{ij} } { \\sqrt{ C_{ii} C_{jj} } }
2788 The values of `R` are between -1 and 1, inclusive.
2790 Parameters
2791 ----------
2792 x : array_like
2793 A 1-D or 2-D array containing multiple variables and observations.
2794 Each row of `x` represents a variable, and each column a single
2795 observation of all those variables. Also see `rowvar` below.
2796 y : array_like, optional
2797 An additional set of variables and observations. `y` has the same
2798 shape as `x`.
2799 rowvar : bool, optional
2800 If `rowvar` is True (default), then each row represents a
2801 variable, with observations in the columns. Otherwise, the relationship
2802 is transposed: each column represents a variable, while the rows
2803 contain observations.
2804 bias : _NoValue, optional
2805 Has no effect, do not use.
2807 .. deprecated:: 1.10.0
2808 ddof : _NoValue, optional
2809 Has no effect, do not use.
2811 .. deprecated:: 1.10.0
2812 dtype : data-type, optional
2813 Data-type of the result. By default, the return data-type will have
2814 at least `numpy.float64` precision.
2816 .. versionadded:: 1.20
2818 Returns
2819 -------
2820 R : ndarray
2821 The correlation coefficient matrix of the variables.
2823 See Also
2824 --------
2825 cov : Covariance matrix
2827 Notes
2828 -----
2829 Due to floating point rounding the resulting array may not be Hermitian,
2830 the diagonal elements may not be 1, and the elements may not satisfy the
2831 inequality abs(a) <= 1. The real and imaginary parts are clipped to the
2832 interval [-1, 1] in an attempt to improve on that situation but is not
2833 much help in the complex case.
2835 This function accepts but discards arguments `bias` and `ddof`. This is
2836 for backwards compatibility with previous versions of this function. These
2837 arguments had no effect on the return values of the function and can be
2838 safely ignored in this and previous versions of numpy.
2840 Examples
2841 --------
2842 In this example we generate two random arrays, ``xarr`` and ``yarr``, and
2843 compute the row-wise and column-wise Pearson correlation coefficients,
2844 ``R``. Since ``rowvar`` is true by default, we first find the row-wise
2845 Pearson correlation coefficients between the variables of ``xarr``.
2847 >>> import numpy as np
2848 >>> rng = np.random.default_rng(seed=42)
2849 >>> xarr = rng.random((3, 3))
2850 >>> xarr
2851 array([[0.77395605, 0.43887844, 0.85859792],
2852 [0.69736803, 0.09417735, 0.97562235],
2853 [0.7611397 , 0.78606431, 0.12811363]])
2854 >>> R1 = np.corrcoef(xarr)
2855 >>> R1
2856 array([[ 1. , 0.99256089, -0.68080986],
2857 [ 0.99256089, 1. , -0.76492172],
2858 [-0.68080986, -0.76492172, 1. ]])
2860 If we add another set of variables and observations ``yarr``, we can
2861 compute the row-wise Pearson correlation coefficients between the
2862 variables in ``xarr`` and ``yarr``.
2864 >>> yarr = rng.random((3, 3))
2865 >>> yarr
2866 array([[0.45038594, 0.37079802, 0.92676499],
2867 [0.64386512, 0.82276161, 0.4434142 ],
2868 [0.22723872, 0.55458479, 0.06381726]])
2869 >>> R2 = np.corrcoef(xarr, yarr)
2870 >>> R2
2871 array([[ 1. , 0.99256089, -0.68080986, 0.75008178, -0.934284 ,
2872 -0.99004057],
2873 [ 0.99256089, 1. , -0.76492172, 0.82502011, -0.97074098,
2874 -0.99981569],
2875 [-0.68080986, -0.76492172, 1. , -0.99507202, 0.89721355,
2876 0.77714685],
2877 [ 0.75008178, 0.82502011, -0.99507202, 1. , -0.93657855,
2878 -0.83571711],
2879 [-0.934284 , -0.97074098, 0.89721355, -0.93657855, 1. ,
2880 0.97517215],
2881 [-0.99004057, -0.99981569, 0.77714685, -0.83571711, 0.97517215,
2882 1. ]])
2884 Finally if we use the option ``rowvar=False``, the columns are now
2885 being treated as the variables and we will find the column-wise Pearson
2886 correlation coefficients between variables in ``xarr`` and ``yarr``.
2888 >>> R3 = np.corrcoef(xarr, yarr, rowvar=False)
2889 >>> R3
2890 array([[ 1. , 0.77598074, -0.47458546, -0.75078643, -0.9665554 ,
2891 0.22423734],
2892 [ 0.77598074, 1. , -0.92346708, -0.99923895, -0.58826587,
2893 -0.44069024],
2894 [-0.47458546, -0.92346708, 1. , 0.93773029, 0.23297648,
2895 0.75137473],
2896 [-0.75078643, -0.99923895, 0.93773029, 1. , 0.55627469,
2897 0.47536961],
2898 [-0.9665554 , -0.58826587, 0.23297648, 0.55627469, 1. ,
2899 -0.46666491],
2900 [ 0.22423734, -0.44069024, 0.75137473, 0.47536961, -0.46666491,
2901 1. ]])
2903 """
2904 if bias is not np._NoValue or ddof is not np._NoValue:
2905 # 2015-03-15, 1.10
2906 warnings.warn('bias and ddof have no effect and are deprecated',
2907 DeprecationWarning, stacklevel=2)
2908 c = cov(x, y, rowvar, dtype=dtype)
2909 try:
2910 d = diag(c)
2911 except ValueError:
2912 # scalar covariance
2913 # nan if incorrect value (nan, inf, 0), 1 otherwise
2914 return c / c
2915 stddev = sqrt(d.real)
2916 c /= stddev[:, None]
2917 c /= stddev[None, :]
2919 # Clip real and imaginary parts to [-1, 1]. This does not guarantee
2920 # abs(a[i,j]) <= 1 for complex arrays, but is the best we can do without
2921 # excessive work.
2922 np.clip(c.real, -1, 1, out=c.real)
2923 if np.iscomplexobj(c):
2924 np.clip(c.imag, -1, 1, out=c.imag)
2926 return c
2929@set_module('numpy')
2930def blackman(M):
2931 """
2932 Return the Blackman window.
2934 The Blackman window is a taper formed by using the first three
2935 terms of a summation of cosines. It was designed to have close to the
2936 minimal leakage possible. It is close to optimal, only slightly worse
2937 than a Kaiser window.
2939 Parameters
2940 ----------
2941 M : int
2942 Number of points in the output window. If zero or less, an empty
2943 array is returned.
2945 Returns
2946 -------
2947 out : ndarray
2948 The window, with the maximum value normalized to one (the value one
2949 appears only if the number of samples is odd).
2951 See Also
2952 --------
2953 bartlett, hamming, hanning, kaiser
2955 Notes
2956 -----
2957 The Blackman window is defined as
2959 .. math:: w(n) = 0.42 - 0.5 \\cos(2\\pi n/M) + 0.08 \\cos(4\\pi n/M)
2961 Most references to the Blackman window come from the signal processing
2962 literature, where it is used as one of many windowing functions for
2963 smoothing values. It is also known as an apodization (which means
2964 "removing the foot", i.e. smoothing discontinuities at the beginning
2965 and end of the sampled signal) or tapering function. It is known as a
2966 "near optimal" tapering function, almost as good (by some measures)
2967 as the kaiser window.
2969 References
2970 ----------
2971 Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra,
2972 Dover Publications, New York.
2974 Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing.
2975 Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471.
2977 Examples
2978 --------
2979 >>> import matplotlib.pyplot as plt
2980 >>> np.blackman(12)
2981 array([-1.38777878e-17, 3.26064346e-02, 1.59903635e-01, # may vary
2982 4.14397981e-01, 7.36045180e-01, 9.67046769e-01,
2983 9.67046769e-01, 7.36045180e-01, 4.14397981e-01,
2984 1.59903635e-01, 3.26064346e-02, -1.38777878e-17])
2986 Plot the window and the frequency response.
2988 .. plot::
2989 :include-source:
2991 import matplotlib.pyplot as plt
2992 from numpy.fft import fft, fftshift
2993 window = np.blackman(51)
2994 plt.plot(window)
2995 plt.title("Blackman window")
2996 plt.ylabel("Amplitude")
2997 plt.xlabel("Sample")
2998 plt.show() # doctest: +SKIP
3000 plt.figure()
3001 A = fft(window, 2048) / 25.5
3002 mag = np.abs(fftshift(A))
3003 freq = np.linspace(-0.5, 0.5, len(A))
3004 with np.errstate(divide='ignore', invalid='ignore'):
3005 response = 20 * np.log10(mag)
3006 response = np.clip(response, -100, 100)
3007 plt.plot(freq, response)
3008 plt.title("Frequency response of Blackman window")
3009 plt.ylabel("Magnitude [dB]")
3010 plt.xlabel("Normalized frequency [cycles per sample]")
3011 plt.axis('tight')
3012 plt.show()
3014 """
3015 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3016 # to double is safe for a range.
3017 values = np.array([0.0, M])
3018 M = values[1]
3020 if M < 1:
3021 return array([], dtype=values.dtype)
3022 if M == 1:
3023 return ones(1, dtype=values.dtype)
3024 n = arange(1-M, M, 2)
3025 return 0.42 + 0.5*cos(pi*n/(M-1)) + 0.08*cos(2.0*pi*n/(M-1))
3028@set_module('numpy')
3029def bartlett(M):
3030 """
3031 Return the Bartlett window.
3033 The Bartlett window is very similar to a triangular window, except
3034 that the end points are at zero. It is often used in signal
3035 processing for tapering a signal, without generating too much
3036 ripple in the frequency domain.
3038 Parameters
3039 ----------
3040 M : int
3041 Number of points in the output window. If zero or less, an
3042 empty array is returned.
3044 Returns
3045 -------
3046 out : array
3047 The triangular window, with the maximum value normalized to one
3048 (the value one appears only if the number of samples is odd), with
3049 the first and last samples equal to zero.
3051 See Also
3052 --------
3053 blackman, hamming, hanning, kaiser
3055 Notes
3056 -----
3057 The Bartlett window is defined as
3059 .. math:: w(n) = \\frac{2}{M-1} \\left(
3060 \\frac{M-1}{2} - \\left|n - \\frac{M-1}{2}\\right|
3061 \\right)
3063 Most references to the Bartlett window come from the signal processing
3064 literature, where it is used as one of many windowing functions for
3065 smoothing values. Note that convolution with this window produces linear
3066 interpolation. It is also known as an apodization (which means "removing
3067 the foot", i.e. smoothing discontinuities at the beginning and end of the
3068 sampled signal) or tapering function. The Fourier transform of the
3069 Bartlett window is the product of two sinc functions. Note the excellent
3070 discussion in Kanasewich [2]_.
3072 References
3073 ----------
3074 .. [1] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra",
3075 Biometrika 37, 1-16, 1950.
3076 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
3077 The University of Alberta Press, 1975, pp. 109-110.
3078 .. [3] A.V. Oppenheim and R.W. Schafer, "Discrete-Time Signal
3079 Processing", Prentice-Hall, 1999, pp. 468-471.
3080 .. [4] Wikipedia, "Window function",
3081 https://en.wikipedia.org/wiki/Window_function
3082 .. [5] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3083 "Numerical Recipes", Cambridge University Press, 1986, page 429.
3085 Examples
3086 --------
3087 >>> import matplotlib.pyplot as plt
3088 >>> np.bartlett(12)
3089 array([ 0. , 0.18181818, 0.36363636, 0.54545455, 0.72727273, # may vary
3090 0.90909091, 0.90909091, 0.72727273, 0.54545455, 0.36363636,
3091 0.18181818, 0. ])
3093 Plot the window and its frequency response (requires SciPy and matplotlib).
3095 .. plot::
3096 :include-source:
3098 import matplotlib.pyplot as plt
3099 from numpy.fft import fft, fftshift
3100 window = np.bartlett(51)
3101 plt.plot(window)
3102 plt.title("Bartlett window")
3103 plt.ylabel("Amplitude")
3104 plt.xlabel("Sample")
3105 plt.show()
3106 plt.figure()
3107 A = fft(window, 2048) / 25.5
3108 mag = np.abs(fftshift(A))
3109 freq = np.linspace(-0.5, 0.5, len(A))
3110 with np.errstate(divide='ignore', invalid='ignore'):
3111 response = 20 * np.log10(mag)
3112 response = np.clip(response, -100, 100)
3113 plt.plot(freq, response)
3114 plt.title("Frequency response of Bartlett window")
3115 plt.ylabel("Magnitude [dB]")
3116 plt.xlabel("Normalized frequency [cycles per sample]")
3117 plt.axis('tight')
3118 plt.show()
3120 """
3121 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3122 # to double is safe for a range.
3123 values = np.array([0.0, M])
3124 M = values[1]
3126 if M < 1:
3127 return array([], dtype=values.dtype)
3128 if M == 1:
3129 return ones(1, dtype=values.dtype)
3130 n = arange(1-M, M, 2)
3131 return where(less_equal(n, 0), 1 + n/(M-1), 1 - n/(M-1))
3134@set_module('numpy')
3135def hanning(M):
3136 """
3137 Return the Hanning window.
3139 The Hanning window is a taper formed by using a weighted cosine.
3141 Parameters
3142 ----------
3143 M : int
3144 Number of points in the output window. If zero or less, an
3145 empty array is returned.
3147 Returns
3148 -------
3149 out : ndarray, shape(M,)
3150 The window, with the maximum value normalized to one (the value
3151 one appears only if `M` is odd).
3153 See Also
3154 --------
3155 bartlett, blackman, hamming, kaiser
3157 Notes
3158 -----
3159 The Hanning window is defined as
3161 .. math:: w(n) = 0.5 - 0.5\\cos\\left(\\frac{2\\pi{n}}{M-1}\\right)
3162 \\qquad 0 \\leq n \\leq M-1
3164 The Hanning was named for Julius von Hann, an Austrian meteorologist.
3165 It is also known as the Cosine Bell. Some authors prefer that it be
3166 called a Hann window, to help avoid confusion with the very similar
3167 Hamming window.
3169 Most references to the Hanning window come from the signal processing
3170 literature, where it is used as one of many windowing functions for
3171 smoothing values. It is also known as an apodization (which means
3172 "removing the foot", i.e. smoothing discontinuities at the beginning
3173 and end of the sampled signal) or tapering function.
3175 References
3176 ----------
3177 .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
3178 spectra, Dover Publications, New York.
3179 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
3180 The University of Alberta Press, 1975, pp. 106-108.
3181 .. [3] Wikipedia, "Window function",
3182 https://en.wikipedia.org/wiki/Window_function
3183 .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3184 "Numerical Recipes", Cambridge University Press, 1986, page 425.
3186 Examples
3187 --------
3188 >>> np.hanning(12)
3189 array([0. , 0.07937323, 0.29229249, 0.57115742, 0.82743037,
3190 0.97974649, 0.97974649, 0.82743037, 0.57115742, 0.29229249,
3191 0.07937323, 0. ])
3193 Plot the window and its frequency response.
3195 .. plot::
3196 :include-source:
3198 import matplotlib.pyplot as plt
3199 from numpy.fft import fft, fftshift
3200 window = np.hanning(51)
3201 plt.plot(window)
3202 plt.title("Hann window")
3203 plt.ylabel("Amplitude")
3204 plt.xlabel("Sample")
3205 plt.show()
3207 plt.figure()
3208 A = fft(window, 2048) / 25.5
3209 mag = np.abs(fftshift(A))
3210 freq = np.linspace(-0.5, 0.5, len(A))
3211 with np.errstate(divide='ignore', invalid='ignore'):
3212 response = 20 * np.log10(mag)
3213 response = np.clip(response, -100, 100)
3214 plt.plot(freq, response)
3215 plt.title("Frequency response of the Hann window")
3216 plt.ylabel("Magnitude [dB]")
3217 plt.xlabel("Normalized frequency [cycles per sample]")
3218 plt.axis('tight')
3219 plt.show()
3221 """
3222 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3223 # to double is safe for a range.
3224 values = np.array([0.0, M])
3225 M = values[1]
3227 if M < 1:
3228 return array([], dtype=values.dtype)
3229 if M == 1:
3230 return ones(1, dtype=values.dtype)
3231 n = arange(1-M, M, 2)
3232 return 0.5 + 0.5*cos(pi*n/(M-1))
3235@set_module('numpy')
3236def hamming(M):
3237 """
3238 Return the Hamming window.
3240 The Hamming window is a taper formed by using a weighted cosine.
3242 Parameters
3243 ----------
3244 M : int
3245 Number of points in the output window. If zero or less, an
3246 empty array is returned.
3248 Returns
3249 -------
3250 out : ndarray
3251 The window, with the maximum value normalized to one (the value
3252 one appears only if the number of samples is odd).
3254 See Also
3255 --------
3256 bartlett, blackman, hanning, kaiser
3258 Notes
3259 -----
3260 The Hamming window is defined as
3262 .. math:: w(n) = 0.54 - 0.46\\cos\\left(\\frac{2\\pi{n}}{M-1}\\right)
3263 \\qquad 0 \\leq n \\leq M-1
3265 The Hamming was named for R. W. Hamming, an associate of J. W. Tukey
3266 and is described in Blackman and Tukey. It was recommended for
3267 smoothing the truncated autocovariance function in the time domain.
3268 Most references to the Hamming window come from the signal processing
3269 literature, where it is used as one of many windowing functions for
3270 smoothing values. It is also known as an apodization (which means
3271 "removing the foot", i.e. smoothing discontinuities at the beginning
3272 and end of the sampled signal) or tapering function.
3274 References
3275 ----------
3276 .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
3277 spectra, Dover Publications, New York.
3278 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
3279 University of Alberta Press, 1975, pp. 109-110.
3280 .. [3] Wikipedia, "Window function",
3281 https://en.wikipedia.org/wiki/Window_function
3282 .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
3283 "Numerical Recipes", Cambridge University Press, 1986, page 425.
3285 Examples
3286 --------
3287 >>> np.hamming(12)
3288 array([ 0.08 , 0.15302337, 0.34890909, 0.60546483, 0.84123594, # may vary
3289 0.98136677, 0.98136677, 0.84123594, 0.60546483, 0.34890909,
3290 0.15302337, 0.08 ])
3292 Plot the window and the frequency response.
3294 .. plot::
3295 :include-source:
3297 import matplotlib.pyplot as plt
3298 from numpy.fft import fft, fftshift
3299 window = np.hamming(51)
3300 plt.plot(window)
3301 plt.title("Hamming window")
3302 plt.ylabel("Amplitude")
3303 plt.xlabel("Sample")
3304 plt.show()
3306 plt.figure()
3307 A = fft(window, 2048) / 25.5
3308 mag = np.abs(fftshift(A))
3309 freq = np.linspace(-0.5, 0.5, len(A))
3310 response = 20 * np.log10(mag)
3311 response = np.clip(response, -100, 100)
3312 plt.plot(freq, response)
3313 plt.title("Frequency response of Hamming window")
3314 plt.ylabel("Magnitude [dB]")
3315 plt.xlabel("Normalized frequency [cycles per sample]")
3316 plt.axis('tight')
3317 plt.show()
3319 """
3320 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3321 # to double is safe for a range.
3322 values = np.array([0.0, M])
3323 M = values[1]
3325 if M < 1:
3326 return array([], dtype=values.dtype)
3327 if M == 1:
3328 return ones(1, dtype=values.dtype)
3329 n = arange(1-M, M, 2)
3330 return 0.54 + 0.46*cos(pi*n/(M-1))
3333## Code from cephes for i0
3335_i0A = [
3336 -4.41534164647933937950E-18,
3337 3.33079451882223809783E-17,
3338 -2.43127984654795469359E-16,
3339 1.71539128555513303061E-15,
3340 -1.16853328779934516808E-14,
3341 7.67618549860493561688E-14,
3342 -4.85644678311192946090E-13,
3343 2.95505266312963983461E-12,
3344 -1.72682629144155570723E-11,
3345 9.67580903537323691224E-11,
3346 -5.18979560163526290666E-10,
3347 2.65982372468238665035E-9,
3348 -1.30002500998624804212E-8,
3349 6.04699502254191894932E-8,
3350 -2.67079385394061173391E-7,
3351 1.11738753912010371815E-6,
3352 -4.41673835845875056359E-6,
3353 1.64484480707288970893E-5,
3354 -5.75419501008210370398E-5,
3355 1.88502885095841655729E-4,
3356 -5.76375574538582365885E-4,
3357 1.63947561694133579842E-3,
3358 -4.32430999505057594430E-3,
3359 1.05464603945949983183E-2,
3360 -2.37374148058994688156E-2,
3361 4.93052842396707084878E-2,
3362 -9.49010970480476444210E-2,
3363 1.71620901522208775349E-1,
3364 -3.04682672343198398683E-1,
3365 6.76795274409476084995E-1
3366 ]
3368_i0B = [
3369 -7.23318048787475395456E-18,
3370 -4.83050448594418207126E-18,
3371 4.46562142029675999901E-17,
3372 3.46122286769746109310E-17,
3373 -2.82762398051658348494E-16,
3374 -3.42548561967721913462E-16,
3375 1.77256013305652638360E-15,
3376 3.81168066935262242075E-15,
3377 -9.55484669882830764870E-15,
3378 -4.15056934728722208663E-14,
3379 1.54008621752140982691E-14,
3380 3.85277838274214270114E-13,
3381 7.18012445138366623367E-13,
3382 -1.79417853150680611778E-12,
3383 -1.32158118404477131188E-11,
3384 -3.14991652796324136454E-11,
3385 1.18891471078464383424E-11,
3386 4.94060238822496958910E-10,
3387 3.39623202570838634515E-9,
3388 2.26666899049817806459E-8,
3389 2.04891858946906374183E-7,
3390 2.89137052083475648297E-6,
3391 6.88975834691682398426E-5,
3392 3.36911647825569408990E-3,
3393 8.04490411014108831608E-1
3394 ]
3397def _chbevl(x, vals):
3398 b0 = vals[0]
3399 b1 = 0.0
3401 for i in range(1, len(vals)):
3402 b2 = b1
3403 b1 = b0
3404 b0 = x*b1 - b2 + vals[i]
3406 return 0.5*(b0 - b2)
3409def _i0_1(x):
3410 return exp(x) * _chbevl(x/2.0-2, _i0A)
3413def _i0_2(x):
3414 return exp(x) * _chbevl(32.0/x - 2.0, _i0B) / sqrt(x)
3417def _i0_dispatcher(x):
3418 return (x,)
3421@array_function_dispatch(_i0_dispatcher)
3422def i0(x):
3423 """
3424 Modified Bessel function of the first kind, order 0.
3426 Usually denoted :math:`I_0`.
3428 Parameters
3429 ----------
3430 x : array_like of float
3431 Argument of the Bessel function.
3433 Returns
3434 -------
3435 out : ndarray, shape = x.shape, dtype = float
3436 The modified Bessel function evaluated at each of the elements of `x`.
3438 See Also
3439 --------
3440 scipy.special.i0, scipy.special.iv, scipy.special.ive
3442 Notes
3443 -----
3444 The scipy implementation is recommended over this function: it is a
3445 proper ufunc written in C, and more than an order of magnitude faster.
3447 We use the algorithm published by Clenshaw [1]_ and referenced by
3448 Abramowitz and Stegun [2]_, for which the function domain is
3449 partitioned into the two intervals [0,8] and (8,inf), and Chebyshev
3450 polynomial expansions are employed in each interval. Relative error on
3451 the domain [0,30] using IEEE arithmetic is documented [3]_ as having a
3452 peak of 5.8e-16 with an rms of 1.4e-16 (n = 30000).
3454 References
3455 ----------
3456 .. [1] C. W. Clenshaw, "Chebyshev series for mathematical functions", in
3457 *National Physical Laboratory Mathematical Tables*, vol. 5, London:
3458 Her Majesty's Stationery Office, 1962.
3459 .. [2] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical
3460 Functions*, 10th printing, New York: Dover, 1964, pp. 379.
3461 https://personal.math.ubc.ca/~cbm/aands/page_379.htm
3462 .. [3] https://metacpan.org/pod/distribution/Math-Cephes/lib/Math/Cephes.pod#i0:-Modified-Bessel-function-of-order-zero
3464 Examples
3465 --------
3466 >>> np.i0(0.)
3467 array(1.0)
3468 >>> np.i0([0, 1, 2, 3])
3469 array([1. , 1.26606588, 2.2795853 , 4.88079259])
3471 """
3472 x = np.asanyarray(x)
3473 if x.dtype.kind == 'c':
3474 raise TypeError("i0 not supported for complex values")
3475 if x.dtype.kind != 'f':
3476 x = x.astype(float)
3477 x = np.abs(x)
3478 return piecewise(x, [x <= 8.0], [_i0_1, _i0_2])
3480## End of cephes code for i0
3483@set_module('numpy')
3484def kaiser(M, beta):
3485 """
3486 Return the Kaiser window.
3488 The Kaiser window is a taper formed by using a Bessel function.
3490 Parameters
3491 ----------
3492 M : int
3493 Number of points in the output window. If zero or less, an
3494 empty array is returned.
3495 beta : float
3496 Shape parameter for window.
3498 Returns
3499 -------
3500 out : array
3501 The window, with the maximum value normalized to one (the value
3502 one appears only if the number of samples is odd).
3504 See Also
3505 --------
3506 bartlett, blackman, hamming, hanning
3508 Notes
3509 -----
3510 The Kaiser window is defined as
3512 .. math:: w(n) = I_0\\left( \\beta \\sqrt{1-\\frac{4n^2}{(M-1)^2}}
3513 \\right)/I_0(\\beta)
3515 with
3517 .. math:: \\quad -\\frac{M-1}{2} \\leq n \\leq \\frac{M-1}{2},
3519 where :math:`I_0` is the modified zeroth-order Bessel function.
3521 The Kaiser was named for Jim Kaiser, who discovered a simple
3522 approximation to the DPSS window based on Bessel functions. The Kaiser
3523 window is a very good approximation to the Digital Prolate Spheroidal
3524 Sequence, or Slepian window, which is the transform which maximizes the
3525 energy in the main lobe of the window relative to total energy.
3527 The Kaiser can approximate many other windows by varying the beta
3528 parameter.
3530 ==== =======================
3531 beta Window shape
3532 ==== =======================
3533 0 Rectangular
3534 5 Similar to a Hamming
3535 6 Similar to a Hanning
3536 8.6 Similar to a Blackman
3537 ==== =======================
3539 A beta value of 14 is probably a good starting point. Note that as beta
3540 gets large, the window narrows, and so the number of samples needs to be
3541 large enough to sample the increasingly narrow spike, otherwise NaNs will
3542 get returned.
3544 Most references to the Kaiser window come from the signal processing
3545 literature, where it is used as one of many windowing functions for
3546 smoothing values. It is also known as an apodization (which means
3547 "removing the foot", i.e. smoothing discontinuities at the beginning
3548 and end of the sampled signal) or tapering function.
3550 References
3551 ----------
3552 .. [1] J. F. Kaiser, "Digital Filters" - Ch 7 in "Systems analysis by
3553 digital computer", Editors: F.F. Kuo and J.F. Kaiser, p 218-285.
3554 John Wiley and Sons, New York, (1966).
3555 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
3556 University of Alberta Press, 1975, pp. 177-178.
3557 .. [3] Wikipedia, "Window function",
3558 https://en.wikipedia.org/wiki/Window_function
3560 Examples
3561 --------
3562 >>> import matplotlib.pyplot as plt
3563 >>> np.kaiser(12, 14)
3564 array([7.72686684e-06, 3.46009194e-03, 4.65200189e-02, # may vary
3565 2.29737120e-01, 5.99885316e-01, 9.45674898e-01,
3566 9.45674898e-01, 5.99885316e-01, 2.29737120e-01,
3567 4.65200189e-02, 3.46009194e-03, 7.72686684e-06])
3570 Plot the window and the frequency response.
3572 .. plot::
3573 :include-source:
3575 import matplotlib.pyplot as plt
3576 from numpy.fft import fft, fftshift
3577 window = np.kaiser(51, 14)
3578 plt.plot(window)
3579 plt.title("Kaiser window")
3580 plt.ylabel("Amplitude")
3581 plt.xlabel("Sample")
3582 plt.show()
3584 plt.figure()
3585 A = fft(window, 2048) / 25.5
3586 mag = np.abs(fftshift(A))
3587 freq = np.linspace(-0.5, 0.5, len(A))
3588 response = 20 * np.log10(mag)
3589 response = np.clip(response, -100, 100)
3590 plt.plot(freq, response)
3591 plt.title("Frequency response of Kaiser window")
3592 plt.ylabel("Magnitude [dB]")
3593 plt.xlabel("Normalized frequency [cycles per sample]")
3594 plt.axis('tight')
3595 plt.show()
3597 """
3598 # Ensures at least float64 via 0.0. M should be an integer, but conversion
3599 # to double is safe for a range. (Simplified result_type with 0.0
3600 # strongly typed. result-type is not/less order sensitive, but that mainly
3601 # matters for integers anyway.)
3602 values = np.array([0.0, M, beta])
3603 M = values[1]
3604 beta = values[2]
3606 if M == 1:
3607 return np.ones(1, dtype=values.dtype)
3608 n = arange(0, M)
3609 alpha = (M-1)/2.0
3610 return i0(beta * sqrt(1-((n-alpha)/alpha)**2.0))/i0(beta)
3613def _sinc_dispatcher(x):
3614 return (x,)
3617@array_function_dispatch(_sinc_dispatcher)
3618def sinc(x):
3619 r"""
3620 Return the normalized sinc function.
3622 The sinc function is equal to :math:`\sin(\pi x)/(\pi x)` for any argument
3623 :math:`x\ne 0`. ``sinc(0)`` takes the limit value 1, making ``sinc`` not
3624 only everywhere continuous but also infinitely differentiable.
3626 .. note::
3628 Note the normalization factor of ``pi`` used in the definition.
3629 This is the most commonly used definition in signal processing.
3630 Use ``sinc(x / np.pi)`` to obtain the unnormalized sinc function
3631 :math:`\sin(x)/x` that is more common in mathematics.
3633 Parameters
3634 ----------
3635 x : ndarray
3636 Array (possibly multi-dimensional) of values for which to calculate
3637 ``sinc(x)``.
3639 Returns
3640 -------
3641 out : ndarray
3642 ``sinc(x)``, which has the same shape as the input.
3644 Notes
3645 -----
3646 The name sinc is short for "sine cardinal" or "sinus cardinalis".
3648 The sinc function is used in various signal processing applications,
3649 including in anti-aliasing, in the construction of a Lanczos resampling
3650 filter, and in interpolation.
3652 For bandlimited interpolation of discrete-time signals, the ideal
3653 interpolation kernel is proportional to the sinc function.
3655 References
3656 ----------
3657 .. [1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram Web
3658 Resource. https://mathworld.wolfram.com/SincFunction.html
3659 .. [2] Wikipedia, "Sinc function",
3660 https://en.wikipedia.org/wiki/Sinc_function
3662 Examples
3663 --------
3664 >>> import matplotlib.pyplot as plt
3665 >>> x = np.linspace(-4, 4, 41)
3666 >>> np.sinc(x)
3667 array([-3.89804309e-17, -4.92362781e-02, -8.40918587e-02, # may vary
3668 -8.90384387e-02, -5.84680802e-02, 3.89804309e-17,
3669 6.68206631e-02, 1.16434881e-01, 1.26137788e-01,
3670 8.50444803e-02, -3.89804309e-17, -1.03943254e-01,
3671 -1.89206682e-01, -2.16236208e-01, -1.55914881e-01,
3672 3.89804309e-17, 2.33872321e-01, 5.04551152e-01,
3673 7.56826729e-01, 9.35489284e-01, 1.00000000e+00,
3674 9.35489284e-01, 7.56826729e-01, 5.04551152e-01,
3675 2.33872321e-01, 3.89804309e-17, -1.55914881e-01,
3676 -2.16236208e-01, -1.89206682e-01, -1.03943254e-01,
3677 -3.89804309e-17, 8.50444803e-02, 1.26137788e-01,
3678 1.16434881e-01, 6.68206631e-02, 3.89804309e-17,
3679 -5.84680802e-02, -8.90384387e-02, -8.40918587e-02,
3680 -4.92362781e-02, -3.89804309e-17])
3682 >>> plt.plot(x, np.sinc(x))
3683 [<matplotlib.lines.Line2D object at 0x...>]
3684 >>> plt.title("Sinc Function")
3685 Text(0.5, 1.0, 'Sinc Function')
3686 >>> plt.ylabel("Amplitude")
3687 Text(0, 0.5, 'Amplitude')
3688 >>> plt.xlabel("X")
3689 Text(0.5, 0, 'X')
3690 >>> plt.show()
3692 """
3693 x = np.asanyarray(x)
3694 y = pi * where(x == 0, 1.0e-20, x)
3695 return sin(y)/y
3698def _ureduce(a, func, keepdims=False, **kwargs):
3699 """
3700 Internal Function.
3701 Call `func` with `a` as first argument swapping the axes to use extended
3702 axis on functions that don't support it natively.
3704 Returns result and a.shape with axis dims set to 1.
3706 Parameters
3707 ----------
3708 a : array_like
3709 Input array or object that can be converted to an array.
3710 func : callable
3711 Reduction function capable of receiving a single axis argument.
3712 It is called with `a` as first argument followed by `kwargs`.
3713 kwargs : keyword arguments
3714 additional keyword arguments to pass to `func`.
3716 Returns
3717 -------
3718 result : tuple
3719 Result of func(a, **kwargs) and a.shape with axis dims set to 1
3720 which can be used to reshape the result to the same shape a ufunc with
3721 keepdims=True would produce.
3723 """
3724 a = np.asanyarray(a)
3725 axis = kwargs.get('axis', None)
3726 out = kwargs.get('out', None)
3728 if keepdims is np._NoValue:
3729 keepdims = False
3731 nd = a.ndim
3732 if axis is not None:
3733 axis = _nx.normalize_axis_tuple(axis, nd)
3735 if keepdims:
3736 if out is not None:
3737 index_out = tuple(
3738 0 if i in axis else slice(None) for i in range(nd))
3739 kwargs['out'] = out[(Ellipsis, ) + index_out]
3741 if len(axis) == 1:
3742 kwargs['axis'] = axis[0]
3743 else:
3744 keep = set(range(nd)) - set(axis)
3745 nkeep = len(keep)
3746 # swap axis that should not be reduced to front
3747 for i, s in enumerate(sorted(keep)):
3748 a = a.swapaxes(i, s)
3749 # merge reduced axis
3750 a = a.reshape(a.shape[:nkeep] + (-1,))
3751 kwargs['axis'] = -1
3752 else:
3753 if keepdims:
3754 if out is not None:
3755 index_out = (0, ) * nd
3756 kwargs['out'] = out[(Ellipsis, ) + index_out]
3758 r = func(a, **kwargs)
3760 if out is not None:
3761 return out
3763 if keepdims:
3764 if axis is None:
3765 index_r = (np.newaxis, ) * nd
3766 else:
3767 index_r = tuple(
3768 np.newaxis if i in axis else slice(None)
3769 for i in range(nd))
3770 r = r[(Ellipsis, ) + index_r]
3772 return r
3775def _median_dispatcher(
3776 a, axis=None, out=None, overwrite_input=None, keepdims=None):
3777 return (a, out)
3780@array_function_dispatch(_median_dispatcher)
3781def median(a, axis=None, out=None, overwrite_input=False, keepdims=False):
3782 """
3783 Compute the median along the specified axis.
3785 Returns the median of the array elements.
3787 Parameters
3788 ----------
3789 a : array_like
3790 Input array or object that can be converted to an array.
3791 axis : {int, sequence of int, None}, optional
3792 Axis or axes along which the medians are computed. The default,
3793 axis=None, will compute the median along a flattened version of
3794 the array.
3796 .. versionadded:: 1.9.0
3798 If a sequence of axes, the array is first flattened along the
3799 given axes, then the median is computed along the resulting
3800 flattened axis.
3801 out : ndarray, optional
3802 Alternative output array in which to place the result. It must
3803 have the same shape and buffer length as the expected output,
3804 but the type (of the output) will be cast if necessary.
3805 overwrite_input : bool, optional
3806 If True, then allow use of memory of input array `a` for
3807 calculations. The input array will be modified by the call to
3808 `median`. This will save memory when you do not need to preserve
3809 the contents of the input array. Treat the input as undefined,
3810 but it will probably be fully or partially sorted. Default is
3811 False. If `overwrite_input` is ``True`` and `a` is not already an
3812 `ndarray`, an error will be raised.
3813 keepdims : bool, optional
3814 If this is set to True, the axes which are reduced are left
3815 in the result as dimensions with size one. With this option,
3816 the result will broadcast correctly against the original `arr`.
3818 .. versionadded:: 1.9.0
3820 Returns
3821 -------
3822 median : ndarray
3823 A new array holding the result. If the input contains integers
3824 or floats smaller than ``float64``, then the output data-type is
3825 ``np.float64``. Otherwise, the data-type of the output is the
3826 same as that of the input. If `out` is specified, that array is
3827 returned instead.
3829 See Also
3830 --------
3831 mean, percentile
3833 Notes
3834 -----
3835 Given a vector ``V`` of length ``N``, the median of ``V`` is the
3836 middle value of a sorted copy of ``V``, ``V_sorted`` - i
3837 e., ``V_sorted[(N-1)/2]``, when ``N`` is odd, and the average of the
3838 two middle values of ``V_sorted`` when ``N`` is even.
3840 Examples
3841 --------
3842 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
3843 >>> a
3844 array([[10, 7, 4],
3845 [ 3, 2, 1]])
3846 >>> np.median(a)
3847 np.float64(3.5)
3848 >>> np.median(a, axis=0)
3849 array([6.5, 4.5, 2.5])
3850 >>> np.median(a, axis=1)
3851 array([7., 2.])
3852 >>> np.median(a, axis=(0, 1))
3853 np.float64(3.5)
3854 >>> m = np.median(a, axis=0)
3855 >>> out = np.zeros_like(m)
3856 >>> np.median(a, axis=0, out=m)
3857 array([6.5, 4.5, 2.5])
3858 >>> m
3859 array([6.5, 4.5, 2.5])
3860 >>> b = a.copy()
3861 >>> np.median(b, axis=1, overwrite_input=True)
3862 array([7., 2.])
3863 >>> assert not np.all(a==b)
3864 >>> b = a.copy()
3865 >>> np.median(b, axis=None, overwrite_input=True)
3866 np.float64(3.5)
3867 >>> assert not np.all(a==b)
3869 """
3870 return _ureduce(a, func=_median, keepdims=keepdims, axis=axis, out=out,
3871 overwrite_input=overwrite_input)
3874def _median(a, axis=None, out=None, overwrite_input=False):
3875 # can't be reasonably be implemented in terms of percentile as we have to
3876 # call mean to not break astropy
3877 a = np.asanyarray(a)
3879 # Set the partition indexes
3880 if axis is None:
3881 sz = a.size
3882 else:
3883 sz = a.shape[axis]
3884 if sz % 2 == 0:
3885 szh = sz // 2
3886 kth = [szh - 1, szh]
3887 else:
3888 kth = [(sz - 1) // 2]
3890 # We have to check for NaNs (as of writing 'M' doesn't actually work).
3891 supports_nans = np.issubdtype(a.dtype, np.inexact) or a.dtype.kind in 'Mm'
3892 if supports_nans:
3893 kth.append(-1)
3895 if overwrite_input:
3896 if axis is None:
3897 part = a.ravel()
3898 part.partition(kth)
3899 else:
3900 a.partition(kth, axis=axis)
3901 part = a
3902 else:
3903 part = partition(a, kth, axis=axis)
3905 if part.shape == ():
3906 # make 0-D arrays work
3907 return part.item()
3908 if axis is None:
3909 axis = 0
3911 indexer = [slice(None)] * part.ndim
3912 index = part.shape[axis] // 2
3913 if part.shape[axis] % 2 == 1:
3914 # index with slice to allow mean (below) to work
3915 indexer[axis] = slice(index, index+1)
3916 else:
3917 indexer[axis] = slice(index-1, index+1)
3918 indexer = tuple(indexer)
3920 # Use mean in both odd and even case to coerce data type,
3921 # using out array if needed.
3922 rout = mean(part[indexer], axis=axis, out=out)
3923 if supports_nans and sz > 0:
3924 # If nans are possible, warn and replace by nans like mean would.
3925 rout = np.lib._utils_impl._median_nancheck(part, rout, axis)
3927 return rout
3930def _percentile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
3931 method=None, keepdims=None, *, weights=None,
3932 interpolation=None):
3933 return (a, q, out, weights)
3936@array_function_dispatch(_percentile_dispatcher)
3937def percentile(a,
3938 q,
3939 axis=None,
3940 out=None,
3941 overwrite_input=False,
3942 method="linear",
3943 keepdims=False,
3944 *,
3945 weights=None,
3946 interpolation=None):
3947 """
3948 Compute the q-th percentile of the data along the specified axis.
3950 Returns the q-th percentile(s) of the array elements.
3952 Parameters
3953 ----------
3954 a : array_like of real numbers
3955 Input array or object that can be converted to an array.
3956 q : array_like of float
3957 Percentage or sequence of percentages for the percentiles to compute.
3958 Values must be between 0 and 100 inclusive.
3959 axis : {int, tuple of int, None}, optional
3960 Axis or axes along which the percentiles are computed. The
3961 default is to compute the percentile(s) along a flattened
3962 version of the array.
3964 .. versionchanged:: 1.9.0
3965 A tuple of axes is supported
3966 out : ndarray, optional
3967 Alternative output array in which to place the result. It must
3968 have the same shape and buffer length as the expected output,
3969 but the type (of the output) will be cast if necessary.
3970 overwrite_input : bool, optional
3971 If True, then allow the input array `a` to be modified by intermediate
3972 calculations, to save memory. In this case, the contents of the input
3973 `a` after this function completes is undefined.
3974 method : str, optional
3975 This parameter specifies the method to use for estimating the
3976 percentile. There are many different methods, some unique to NumPy.
3977 See the notes for explanation. The options sorted by their R type
3978 as summarized in the H&F paper [1]_ are:
3980 1. 'inverted_cdf'
3981 2. 'averaged_inverted_cdf'
3982 3. 'closest_observation'
3983 4. 'interpolated_inverted_cdf'
3984 5. 'hazen'
3985 6. 'weibull'
3986 7. 'linear' (default)
3987 8. 'median_unbiased'
3988 9. 'normal_unbiased'
3990 The first three methods are discontinuous. NumPy further defines the
3991 following discontinuous variations of the default 'linear' (7.) option:
3993 * 'lower'
3994 * 'higher',
3995 * 'midpoint'
3996 * 'nearest'
3998 .. versionchanged:: 1.22.0
3999 This argument was previously called "interpolation" and only
4000 offered the "linear" default and last four options.
4002 keepdims : bool, optional
4003 If this is set to True, the axes which are reduced are left in
4004 the result as dimensions with size one. With this option, the
4005 result will broadcast correctly against the original array `a`.
4007 .. versionadded:: 1.9.0
4009 weights : array_like, optional
4010 An array of weights associated with the values in `a`. Each value in
4011 `a` contributes to the percentile according to its associated weight.
4012 The weights array can either be 1-D (in which case its length must be
4013 the size of `a` along the given axis) or of the same shape as `a`.
4014 If `weights=None`, then all data in `a` are assumed to have a
4015 weight equal to one.
4016 Only `method="inverted_cdf"` supports weights.
4017 See the notes for more details.
4019 .. versionadded:: 2.0.0
4021 interpolation : str, optional
4022 Deprecated name for the method keyword argument.
4024 .. deprecated:: 1.22.0
4026 Returns
4027 -------
4028 percentile : scalar or ndarray
4029 If `q` is a single percentile and `axis=None`, then the result
4030 is a scalar. If multiple percentiles are given, first axis of
4031 the result corresponds to the percentiles. The other axes are
4032 the axes that remain after the reduction of `a`. If the input
4033 contains integers or floats smaller than ``float64``, the output
4034 data-type is ``float64``. Otherwise, the output data-type is the
4035 same as that of the input. If `out` is specified, that array is
4036 returned instead.
4038 See Also
4039 --------
4040 mean
4041 median : equivalent to ``percentile(..., 50)``
4042 nanpercentile
4043 quantile : equivalent to percentile, except q in the range [0, 1].
4045 Notes
4046 -----
4047 In general, the percentile at percentage level :math:`q` of a cumulative
4048 distribution function :math:`F(y)=P(Y \\leq y)` with probability measure
4049 :math:`P` is defined as any number :math:`x` that fulfills the
4050 *coverage conditions*
4052 .. math:: P(Y < x) \\leq q/100 \\quad\\text{and}
4053 \\quad P(Y \\leq x) \\geq q/100
4055 with random variable :math:`Y\\sim P`.
4056 Sample percentiles, the result of ``percentile``, provide nonparametric
4057 estimation of the underlying population counterparts, represented by the
4058 unknown :math:`F`, given a data vector ``a`` of length ``n``.
4060 One type of estimators arises when one considers :math:`F` as the empirical
4061 distribution function of the data, i.e.
4062 :math:`F(y) = \\frac{1}{n} \\sum_i 1_{a_i \\leq y}`.
4063 Then, different methods correspond to different choices of :math:`x` that
4064 fulfill the above inequalities. Methods that follow this approach are
4065 ``inverted_cdf`` and ``averaged_inverted_cdf``.
4067 A more general way to define sample percentile estimators is as follows.
4068 The empirical q-percentile of ``a`` is the ``n * q/100``-th value of the
4069 way from the minimum to the maximum in a sorted copy of ``a``. The values
4070 and distances of the two nearest neighbors as well as the `method`
4071 parameter will determine the percentile if the normalized ranking does not
4072 match the location of ``n * q/100`` exactly. This function is the same as
4073 the median if ``q=50``, the same as the minimum if ``q=0`` and the same
4074 as the maximum if ``q=100``.
4076 The optional `method` parameter specifies the method to use when the
4077 desired percentile lies between two indexes ``i`` and ``j = i + 1``.
4078 In that case, we first determine ``i + g``, a virtual index that lies
4079 between ``i`` and ``j``, where ``i`` is the floor and ``g`` is the
4080 fractional part of the index. The final result is, then, an interpolation
4081 of ``a[i]`` and ``a[j]`` based on ``g``. During the computation of ``g``,
4082 ``i`` and ``j`` are modified using correction constants ``alpha`` and
4083 ``beta`` whose choices depend on the ``method`` used. Finally, note that
4084 since Python uses 0-based indexing, the code subtracts another 1 from the
4085 index internally.
4087 The following formula determines the virtual index ``i + g``, the location
4088 of the percentile in the sorted sample:
4090 .. math::
4091 i + g = (q / 100) * ( n - alpha - beta + 1 ) + alpha
4093 The different methods then work as follows
4095 inverted_cdf:
4096 method 1 of H&F [1]_.
4097 This method gives discontinuous results:
4099 * if g > 0 ; then take j
4100 * if g = 0 ; then take i
4102 averaged_inverted_cdf:
4103 method 2 of H&F [1]_.
4104 This method gives discontinuous results:
4106 * if g > 0 ; then take j
4107 * if g = 0 ; then average between bounds
4109 closest_observation:
4110 method 3 of H&F [1]_.
4111 This method gives discontinuous results:
4113 * if g > 0 ; then take j
4114 * if g = 0 and index is odd ; then take j
4115 * if g = 0 and index is even ; then take i
4117 interpolated_inverted_cdf:
4118 method 4 of H&F [1]_.
4119 This method gives continuous results using:
4121 * alpha = 0
4122 * beta = 1
4124 hazen:
4125 method 5 of H&F [1]_.
4126 This method gives continuous results using:
4128 * alpha = 1/2
4129 * beta = 1/2
4131 weibull:
4132 method 6 of H&F [1]_.
4133 This method gives continuous results using:
4135 * alpha = 0
4136 * beta = 0
4138 linear:
4139 method 7 of H&F [1]_.
4140 This method gives continuous results using:
4142 * alpha = 1
4143 * beta = 1
4145 median_unbiased:
4146 method 8 of H&F [1]_.
4147 This method is probably the best method if the sample
4148 distribution function is unknown (see reference).
4149 This method gives continuous results using:
4151 * alpha = 1/3
4152 * beta = 1/3
4154 normal_unbiased:
4155 method 9 of H&F [1]_.
4156 This method is probably the best method if the sample
4157 distribution function is known to be normal.
4158 This method gives continuous results using:
4160 * alpha = 3/8
4161 * beta = 3/8
4163 lower:
4164 NumPy method kept for backwards compatibility.
4165 Takes ``i`` as the interpolation point.
4167 higher:
4168 NumPy method kept for backwards compatibility.
4169 Takes ``j`` as the interpolation point.
4171 nearest:
4172 NumPy method kept for backwards compatibility.
4173 Takes ``i`` or ``j``, whichever is nearest.
4175 midpoint:
4176 NumPy method kept for backwards compatibility.
4177 Uses ``(i + j) / 2``.
4179 For weighted percentiles, the above coverage conditions still hold. The
4180 empirical cumulative distribution is simply replaced by its weighted
4181 version, i.e.
4182 :math:`P(Y \\leq t) = \\frac{1}{\\sum_i w_i} \\sum_i w_i 1_{x_i \\leq t}`.
4183 Only ``method="inverted_cdf"`` supports weights.
4186 Examples
4187 --------
4188 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
4189 >>> a
4190 array([[10, 7, 4],
4191 [ 3, 2, 1]])
4192 >>> np.percentile(a, 50)
4193 3.5
4194 >>> np.percentile(a, 50, axis=0)
4195 array([6.5, 4.5, 2.5])
4196 >>> np.percentile(a, 50, axis=1)
4197 array([7., 2.])
4198 >>> np.percentile(a, 50, axis=1, keepdims=True)
4199 array([[7.],
4200 [2.]])
4202 >>> m = np.percentile(a, 50, axis=0)
4203 >>> out = np.zeros_like(m)
4204 >>> np.percentile(a, 50, axis=0, out=out)
4205 array([6.5, 4.5, 2.5])
4206 >>> m
4207 array([6.5, 4.5, 2.5])
4209 >>> b = a.copy()
4210 >>> np.percentile(b, 50, axis=1, overwrite_input=True)
4211 array([7., 2.])
4212 >>> assert not np.all(a == b)
4214 The different methods can be visualized graphically:
4216 .. plot::
4218 import matplotlib.pyplot as plt
4220 a = np.arange(4)
4221 p = np.linspace(0, 100, 6001)
4222 ax = plt.gca()
4223 lines = [
4224 ('linear', '-', 'C0'),
4225 ('inverted_cdf', ':', 'C1'),
4226 # Almost the same as `inverted_cdf`:
4227 ('averaged_inverted_cdf', '-.', 'C1'),
4228 ('closest_observation', ':', 'C2'),
4229 ('interpolated_inverted_cdf', '--', 'C1'),
4230 ('hazen', '--', 'C3'),
4231 ('weibull', '-.', 'C4'),
4232 ('median_unbiased', '--', 'C5'),
4233 ('normal_unbiased', '-.', 'C6'),
4234 ]
4235 for method, style, color in lines:
4236 ax.plot(
4237 p, np.percentile(a, p, method=method),
4238 label=method, linestyle=style, color=color)
4239 ax.set(
4240 title='Percentiles for different methods and data: ' + str(a),
4241 xlabel='Percentile',
4242 ylabel='Estimated percentile value',
4243 yticks=a)
4244 ax.legend(bbox_to_anchor=(1.03, 1))
4245 plt.tight_layout()
4246 plt.show()
4248 References
4249 ----------
4250 .. [1] R. J. Hyndman and Y. Fan,
4251 "Sample quantiles in statistical packages,"
4252 The American Statistician, 50(4), pp. 361-365, 1996
4254 """
4255 if interpolation is not None:
4256 method = _check_interpolation_as_method(
4257 method, interpolation, "percentile")
4259 a = np.asanyarray(a)
4260 if a.dtype.kind == "c":
4261 raise TypeError("a must be an array of real numbers")
4263 # Use dtype of array if possible (e.g., if q is a python int or float)
4264 # by making the divisor have the dtype of the data array.
4265 q = np.true_divide(q, a.dtype.type(100) if a.dtype.kind == "f" else 100)
4266 q = asanyarray(q) # undo any decay that the ufunc performed (see gh-13105)
4267 if not _quantile_is_valid(q):
4268 raise ValueError("Percentiles must be in the range [0, 100]")
4270 if weights is not None:
4271 if method != "inverted_cdf":
4272 msg = ("Only method 'inverted_cdf' supports weights. "
4273 f"Got: {method}.")
4274 raise ValueError(msg)
4275 if axis is not None:
4276 axis = _nx.normalize_axis_tuple(axis, a.ndim, argname="axis")
4277 weights = _weights_are_valid(weights=weights, a=a, axis=axis)
4278 if np.any(weights < 0):
4279 raise ValueError("Weights must be non-negative.")
4281 return _quantile_unchecked(
4282 a, q, axis, out, overwrite_input, method, keepdims, weights)
4285def _quantile_dispatcher(a, q, axis=None, out=None, overwrite_input=None,
4286 method=None, keepdims=None, *, weights=None,
4287 interpolation=None):
4288 return (a, q, out, weights)
4291@array_function_dispatch(_quantile_dispatcher)
4292def quantile(a,
4293 q,
4294 axis=None,
4295 out=None,
4296 overwrite_input=False,
4297 method="linear",
4298 keepdims=False,
4299 *,
4300 weights=None,
4301 interpolation=None):
4302 """
4303 Compute the q-th quantile of the data along the specified axis.
4305 .. versionadded:: 1.15.0
4307 Parameters
4308 ----------
4309 a : array_like of real numbers
4310 Input array or object that can be converted to an array.
4311 q : array_like of float
4312 Probability or sequence of probabilities for the quantiles to compute.
4313 Values must be between 0 and 1 inclusive.
4314 axis : {int, tuple of int, None}, optional
4315 Axis or axes along which the quantiles are computed. The default is
4316 to compute the quantile(s) along a flattened version of the array.
4317 out : ndarray, optional
4318 Alternative output array in which to place the result. It must have
4319 the same shape and buffer length as the expected output, but the
4320 type (of the output) will be cast if necessary.
4321 overwrite_input : bool, optional
4322 If True, then allow the input array `a` to be modified by
4323 intermediate calculations, to save memory. In this case, the
4324 contents of the input `a` after this function completes is
4325 undefined.
4326 method : str, optional
4327 This parameter specifies the method to use for estimating the
4328 quantile. There are many different methods, some unique to NumPy.
4329 See the notes for explanation. The options sorted by their R type
4330 as summarized in the H&F paper [1]_ are:
4332 1. 'inverted_cdf'
4333 2. 'averaged_inverted_cdf'
4334 3. 'closest_observation'
4335 4. 'interpolated_inverted_cdf'
4336 5. 'hazen'
4337 6. 'weibull'
4338 7. 'linear' (default)
4339 8. 'median_unbiased'
4340 9. 'normal_unbiased'
4342 The first three methods are discontinuous. NumPy further defines the
4343 following discontinuous variations of the default 'linear' (7.) option:
4345 * 'lower'
4346 * 'higher',
4347 * 'midpoint'
4348 * 'nearest'
4350 .. versionchanged:: 1.22.0
4351 This argument was previously called "interpolation" and only
4352 offered the "linear" default and last four options.
4354 keepdims : bool, optional
4355 If this is set to True, the axes which are reduced are left in
4356 the result as dimensions with size one. With this option, the
4357 result will broadcast correctly against the original array `a`.
4359 weights : array_like, optional
4360 An array of weights associated with the values in `a`. Each value in
4361 `a` contributes to the quantile according to its associated weight.
4362 The weights array can either be 1-D (in which case its length must be
4363 the size of `a` along the given axis) or of the same shape as `a`.
4364 If `weights=None`, then all data in `a` are assumed to have a
4365 weight equal to one.
4366 Only `method="inverted_cdf"` supports weights.
4367 See the notes for more details.
4369 .. versionadded:: 2.0.0
4371 interpolation : str, optional
4372 Deprecated name for the method keyword argument.
4374 .. deprecated:: 1.22.0
4376 Returns
4377 -------
4378 quantile : scalar or ndarray
4379 If `q` is a single probability and `axis=None`, then the result
4380 is a scalar. If multiple probability levels are given, first axis
4381 of the result corresponds to the quantiles. The other axes are
4382 the axes that remain after the reduction of `a`. If the input
4383 contains integers or floats smaller than ``float64``, the output
4384 data-type is ``float64``. Otherwise, the output data-type is the
4385 same as that of the input. If `out` is specified, that array is
4386 returned instead.
4388 See Also
4389 --------
4390 mean
4391 percentile : equivalent to quantile, but with q in the range [0, 100].
4392 median : equivalent to ``quantile(..., 0.5)``
4393 nanquantile
4395 Notes
4396 -----
4397 In general, the quantile at probability level :math:`q` of a cumulative
4398 distribution function :math:`F(y)=P(Y \\leq y)` with probability measure
4399 :math:`P` is defined as any number :math:`x` that fulfills the
4400 *coverage conditions*
4402 .. math:: P(Y < x) \\leq q \\quad\\text{and}\\quad P(Y \\leq x) \\geq q
4404 with random variable :math:`Y\\sim P`.
4405 Sample quantiles, the result of ``quantile``, provide nonparametric
4406 estimation of the underlying population counterparts, represented by the
4407 unknown :math:`F`, given a data vector ``a`` of length ``n``.
4409 One type of estimators arises when one considers :math:`F` as the empirical
4410 distribution function of the data, i.e.
4411 :math:`F(y) = \\frac{1}{n} \\sum_i 1_{a_i \\leq y}`.
4412 Then, different methods correspond to different choices of :math:`x` that
4413 fulfill the above inequalities. Methods that follow this approach are
4414 ``inverted_cdf`` and ``averaged_inverted_cdf``.
4416 A more general way to define sample quantile estimators is as follows.
4417 The empirical q-quantile of ``a`` is the ``n * q``-th value of the
4418 way from the minimum to the maximum in a sorted copy of ``a``. The values
4419 and distances of the two nearest neighbors as well as the `method`
4420 parameter will determine the quantile if the normalized ranking does not
4421 match the location of ``n * q`` exactly. This function is the same as
4422 the median if ``q=0.5``, the same as the minimum if ``q=0.0`` and the same
4423 as the maximum if ``q=1.0``.
4425 The optional `method` parameter specifies the method to use when the
4426 desired quantile lies between two indexes ``i`` and ``j = i + 1``.
4427 In that case, we first determine ``i + g``, a virtual index that lies
4428 between ``i`` and ``j``, where ``i`` is the floor and ``g`` is the
4429 fractional part of the index. The final result is, then, an interpolation
4430 of ``a[i]`` and ``a[j]`` based on ``g``. During the computation of ``g``,
4431 ``i`` and ``j`` are modified using correction constants ``alpha`` and
4432 ``beta`` whose choices depend on the ``method`` used. Finally, note that
4433 since Python uses 0-based indexing, the code subtracts another 1 from the
4434 index internally.
4436 The following formula determines the virtual index ``i + g``, the location
4437 of the quantile in the sorted sample:
4439 .. math::
4440 i + g = q * ( n - alpha - beta + 1 ) + alpha
4442 The different methods then work as follows
4444 inverted_cdf:
4445 method 1 of H&F [1]_.
4446 This method gives discontinuous results:
4448 * if g > 0 ; then take j
4449 * if g = 0 ; then take i
4451 averaged_inverted_cdf:
4452 method 2 of H&F [1]_.
4453 This method gives discontinuous results:
4455 * if g > 0 ; then take j
4456 * if g = 0 ; then average between bounds
4458 closest_observation:
4459 method 3 of H&F [1]_.
4460 This method gives discontinuous results:
4462 * if g > 0 ; then take j
4463 * if g = 0 and index is odd ; then take j
4464 * if g = 0 and index is even ; then take i
4466 interpolated_inverted_cdf:
4467 method 4 of H&F [1]_.
4468 This method gives continuous results using:
4470 * alpha = 0
4471 * beta = 1
4473 hazen:
4474 method 5 of H&F [1]_.
4475 This method gives continuous results using:
4477 * alpha = 1/2
4478 * beta = 1/2
4480 weibull:
4481 method 6 of H&F [1]_.
4482 This method gives continuous results using:
4484 * alpha = 0
4485 * beta = 0
4487 linear:
4488 method 7 of H&F [1]_.
4489 This method gives continuous results using:
4491 * alpha = 1
4492 * beta = 1
4494 median_unbiased:
4495 method 8 of H&F [1]_.
4496 This method is probably the best method if the sample
4497 distribution function is unknown (see reference).
4498 This method gives continuous results using:
4500 * alpha = 1/3
4501 * beta = 1/3
4503 normal_unbiased:
4504 method 9 of H&F [1]_.
4505 This method is probably the best method if the sample
4506 distribution function is known to be normal.
4507 This method gives continuous results using:
4509 * alpha = 3/8
4510 * beta = 3/8
4512 lower:
4513 NumPy method kept for backwards compatibility.
4514 Takes ``i`` as the interpolation point.
4516 higher:
4517 NumPy method kept for backwards compatibility.
4518 Takes ``j`` as the interpolation point.
4520 nearest:
4521 NumPy method kept for backwards compatibility.
4522 Takes ``i`` or ``j``, whichever is nearest.
4524 midpoint:
4525 NumPy method kept for backwards compatibility.
4526 Uses ``(i + j) / 2``.
4528 **Weighted quantiles:**
4529 For weighted quantiles, the above coverage conditions still hold. The
4530 empirical cumulative distribution is simply replaced by its weighted
4531 version, i.e.
4532 :math:`P(Y \\leq t) = \\frac{1}{\\sum_i w_i} \\sum_i w_i 1_{x_i \\leq t}`.
4533 Only ``method="inverted_cdf"`` supports weights.
4535 Examples
4536 --------
4537 >>> a = np.array([[10, 7, 4], [3, 2, 1]])
4538 >>> a
4539 array([[10, 7, 4],
4540 [ 3, 2, 1]])
4541 >>> np.quantile(a, 0.5)
4542 3.5
4543 >>> np.quantile(a, 0.5, axis=0)
4544 array([6.5, 4.5, 2.5])
4545 >>> np.quantile(a, 0.5, axis=1)
4546 array([7., 2.])
4547 >>> np.quantile(a, 0.5, axis=1, keepdims=True)
4548 array([[7.],
4549 [2.]])
4550 >>> m = np.quantile(a, 0.5, axis=0)
4551 >>> out = np.zeros_like(m)
4552 >>> np.quantile(a, 0.5, axis=0, out=out)
4553 array([6.5, 4.5, 2.5])
4554 >>> m
4555 array([6.5, 4.5, 2.5])
4556 >>> b = a.copy()
4557 >>> np.quantile(b, 0.5, axis=1, overwrite_input=True)
4558 array([7., 2.])
4559 >>> assert not np.all(a == b)
4561 See also `numpy.percentile` for a visualization of most methods.
4563 References
4564 ----------
4565 .. [1] R. J. Hyndman and Y. Fan,
4566 "Sample quantiles in statistical packages,"
4567 The American Statistician, 50(4), pp. 361-365, 1996
4569 """
4570 if interpolation is not None:
4571 method = _check_interpolation_as_method(
4572 method, interpolation, "quantile")
4574 a = np.asanyarray(a)
4575 if a.dtype.kind == "c":
4576 raise TypeError("a must be an array of real numbers")
4578 # Use dtype of array if possible (e.g., if q is a python int or float).
4579 if isinstance(q, (int, float)) and a.dtype.kind == "f":
4580 q = np.asanyarray(q, dtype=a.dtype)
4581 else:
4582 q = np.asanyarray(q)
4584 if not _quantile_is_valid(q):
4585 raise ValueError("Quantiles must be in the range [0, 1]")
4587 if weights is not None:
4588 if method != "inverted_cdf":
4589 msg = ("Only method 'inverted_cdf' supports weights. "
4590 f"Got: {method}.")
4591 raise ValueError(msg)
4592 if axis is not None:
4593 axis = _nx.normalize_axis_tuple(axis, a.ndim, argname="axis")
4594 weights = _weights_are_valid(weights=weights, a=a, axis=axis)
4595 if np.any(weights < 0):
4596 raise ValueError("Weights must be non-negative.")
4598 return _quantile_unchecked(
4599 a, q, axis, out, overwrite_input, method, keepdims, weights)
4602def _quantile_unchecked(a,
4603 q,
4604 axis=None,
4605 out=None,
4606 overwrite_input=False,
4607 method="linear",
4608 keepdims=False,
4609 weights=None):
4610 """Assumes that q is in [0, 1], and is an ndarray"""
4611 return _ureduce(a,
4612 func=_quantile_ureduce_func,
4613 q=q,
4614 weights=weights,
4615 keepdims=keepdims,
4616 axis=axis,
4617 out=out,
4618 overwrite_input=overwrite_input,
4619 method=method)
4622def _quantile_is_valid(q):
4623 # avoid expensive reductions, relevant for arrays with < O(1000) elements
4624 if q.ndim == 1 and q.size < 10:
4625 for i in range(q.size):
4626 if not (0.0 <= q[i] <= 1.0):
4627 return False
4628 else:
4629 if not (q.min() >= 0 and q.max() <= 1):
4630 return False
4631 return True
4634def _check_interpolation_as_method(method, interpolation, fname):
4635 # Deprecated NumPy 1.22, 2021-11-08
4636 warnings.warn(
4637 f"the `interpolation=` argument to {fname} was renamed to "
4638 "`method=`, which has additional options.\n"
4639 "Users of the modes 'nearest', 'lower', 'higher', or "
4640 "'midpoint' are encouraged to review the method they used. "
4641 "(Deprecated NumPy 1.22)",
4642 DeprecationWarning, stacklevel=4)
4643 if method != "linear":
4644 # sanity check, we assume this basically never happens
4645 raise TypeError(
4646 "You shall not pass both `method` and `interpolation`!\n"
4647 "(`interpolation` is Deprecated in favor of `method`)")
4648 return interpolation
4651def _compute_virtual_index(n, quantiles, alpha: float, beta: float):
4652 """
4653 Compute the floating point indexes of an array for the linear
4654 interpolation of quantiles.
4655 n : array_like
4656 The sample sizes.
4657 quantiles : array_like
4658 The quantiles values.
4659 alpha : float
4660 A constant used to correct the index computed.
4661 beta : float
4662 A constant used to correct the index computed.
4664 alpha and beta values depend on the chosen method
4665 (see quantile documentation)
4667 Reference:
4668 Hyndman&Fan paper "Sample Quantiles in Statistical Packages",
4669 DOI: 10.1080/00031305.1996.10473566
4670 """
4671 return n * quantiles + (
4672 alpha + quantiles * (1 - alpha - beta)
4673 ) - 1
4676def _get_gamma(virtual_indexes, previous_indexes, method):
4677 """
4678 Compute gamma (a.k.a 'm' or 'weight') for the linear interpolation
4679 of quantiles.
4681 virtual_indexes : array_like
4682 The indexes where the percentile is supposed to be found in the sorted
4683 sample.
4684 previous_indexes : array_like
4685 The floor values of virtual_indexes.
4686 interpolation : dict
4687 The interpolation method chosen, which may have a specific rule
4688 modifying gamma.
4690 gamma is usually the fractional part of virtual_indexes but can be modified
4691 by the interpolation method.
4692 """
4693 gamma = np.asanyarray(virtual_indexes - previous_indexes)
4694 gamma = method["fix_gamma"](gamma, virtual_indexes)
4695 # Ensure both that we have an array, and that we keep the dtype
4696 # (which may have been matched to the input array).
4697 return np.asanyarray(gamma, dtype=virtual_indexes.dtype)
4700def _lerp(a, b, t, out=None):
4701 """
4702 Compute the linear interpolation weighted by gamma on each point of
4703 two same shape array.
4705 a : array_like
4706 Left bound.
4707 b : array_like
4708 Right bound.
4709 t : array_like
4710 The interpolation weight.
4711 out : array_like
4712 Output array.
4713 """
4714 diff_b_a = subtract(b, a)
4715 # asanyarray is a stop-gap until gh-13105
4716 lerp_interpolation = asanyarray(add(a, diff_b_a * t, out=out))
4717 subtract(b, diff_b_a * (1 - t), out=lerp_interpolation, where=t >= 0.5,
4718 casting='unsafe', dtype=type(lerp_interpolation.dtype))
4719 if lerp_interpolation.ndim == 0 and out is None:
4720 lerp_interpolation = lerp_interpolation[()] # unpack 0d arrays
4721 return lerp_interpolation
4724def _get_gamma_mask(shape, default_value, conditioned_value, where):
4725 out = np.full(shape, default_value)
4726 np.copyto(out, conditioned_value, where=where, casting="unsafe")
4727 return out
4730def _discret_interpolation_to_boundaries(index, gamma_condition_fun):
4731 previous = np.floor(index)
4732 next = previous + 1
4733 gamma = index - previous
4734 res = _get_gamma_mask(shape=index.shape,
4735 default_value=next,
4736 conditioned_value=previous,
4737 where=gamma_condition_fun(gamma, index)
4738 ).astype(np.intp)
4739 # Some methods can lead to out-of-bound integers, clip them:
4740 res[res < 0] = 0
4741 return res
4744def _closest_observation(n, quantiles):
4745 gamma_fun = lambda gamma, index: (gamma == 0) & (np.floor(index) % 2 == 0)
4746 return _discret_interpolation_to_boundaries((n * quantiles) - 1 - 0.5,
4747 gamma_fun)
4750def _inverted_cdf(n, quantiles):
4751 gamma_fun = lambda gamma, _: (gamma == 0)
4752 return _discret_interpolation_to_boundaries((n * quantiles) - 1,
4753 gamma_fun)
4756def _quantile_ureduce_func(
4757 a: np.array,
4758 q: np.array,
4759 weights: np.array,
4760 axis: int = None,
4761 out=None,
4762 overwrite_input: bool = False,
4763 method="linear",
4764) -> np.array:
4765 if q.ndim > 2:
4766 # The code below works fine for nd, but it might not have useful
4767 # semantics. For now, keep the supported dimensions the same as it was
4768 # before.
4769 raise ValueError("q must be a scalar or 1d")
4770 if overwrite_input:
4771 if axis is None:
4772 axis = 0
4773 arr = a.ravel()
4774 wgt = None if weights is None else weights.ravel()
4775 else:
4776 arr = a
4777 wgt = weights
4778 else:
4779 if axis is None:
4780 axis = 0
4781 arr = a.flatten()
4782 wgt = None if weights is None else weights.flatten()
4783 else:
4784 arr = a.copy()
4785 wgt = weights
4786 result = _quantile(arr,
4787 quantiles=q,
4788 axis=axis,
4789 method=method,
4790 out=out,
4791 weights=wgt)
4792 return result
4795def _get_indexes(arr, virtual_indexes, valid_values_count):
4796 """
4797 Get the valid indexes of arr neighbouring virtual_indexes.
4798 Note
4799 This is a companion function to linear interpolation of
4800 Quantiles
4802 Returns
4803 -------
4804 (previous_indexes, next_indexes): Tuple
4805 A Tuple of virtual_indexes neighbouring indexes
4806 """
4807 previous_indexes = np.asanyarray(np.floor(virtual_indexes))
4808 next_indexes = np.asanyarray(previous_indexes + 1)
4809 indexes_above_bounds = virtual_indexes >= valid_values_count - 1
4810 # When indexes is above max index, take the max value of the array
4811 if indexes_above_bounds.any():
4812 previous_indexes[indexes_above_bounds] = -1
4813 next_indexes[indexes_above_bounds] = -1
4814 # When indexes is below min index, take the min value of the array
4815 indexes_below_bounds = virtual_indexes < 0
4816 if indexes_below_bounds.any():
4817 previous_indexes[indexes_below_bounds] = 0
4818 next_indexes[indexes_below_bounds] = 0
4819 if np.issubdtype(arr.dtype, np.inexact):
4820 # After the sort, slices having NaNs will have for last element a NaN
4821 virtual_indexes_nans = np.isnan(virtual_indexes)
4822 if virtual_indexes_nans.any():
4823 previous_indexes[virtual_indexes_nans] = -1
4824 next_indexes[virtual_indexes_nans] = -1
4825 previous_indexes = previous_indexes.astype(np.intp)
4826 next_indexes = next_indexes.astype(np.intp)
4827 return previous_indexes, next_indexes
4830def _quantile(
4831 arr: np.array,
4832 quantiles: np.array,
4833 axis: int = -1,
4834 method="linear",
4835 out=None,
4836 weights=None,
4837):
4838 """
4839 Private function that doesn't support extended axis or keepdims.
4840 These methods are extended to this function using _ureduce
4841 See nanpercentile for parameter usage
4842 It computes the quantiles of the array for the given axis.
4843 A linear interpolation is performed based on the `interpolation`.
4845 By default, the method is "linear" where alpha == beta == 1 which
4846 performs the 7th method of Hyndman&Fan.
4847 With "median_unbiased" we get alpha == beta == 1/3
4848 thus the 8th method of Hyndman&Fan.
4849 """
4850 # --- Setup
4851 arr = np.asanyarray(arr)
4852 values_count = arr.shape[axis]
4853 # The dimensions of `q` are prepended to the output shape, so we need the
4854 # axis being sampled from `arr` to be last.
4855 if axis != 0: # But moveaxis is slow, so only call it if necessary.
4856 arr = np.moveaxis(arr, axis, destination=0)
4857 supports_nans = (
4858 np.issubdtype(arr.dtype, np.inexact) or arr.dtype.kind in 'Mm'
4859 )
4861 if weights is None:
4862 # --- Computation of indexes
4863 # Index where to find the value in the sorted array.
4864 # Virtual because it is a floating point value, not an valid index.
4865 # The nearest neighbours are used for interpolation
4866 try:
4867 method_props = _QuantileMethods[method]
4868 except KeyError:
4869 raise ValueError(
4870 f"{method!r} is not a valid method. Use one of: "
4871 f"{_QuantileMethods.keys()}") from None
4872 virtual_indexes = method_props["get_virtual_index"](values_count,
4873 quantiles)
4874 virtual_indexes = np.asanyarray(virtual_indexes)
4876 if method_props["fix_gamma"] is None:
4877 supports_integers = True
4878 else:
4879 int_virtual_indices = np.issubdtype(virtual_indexes.dtype,
4880 np.integer)
4881 supports_integers = method == 'linear' and int_virtual_indices
4883 if supports_integers:
4884 # No interpolation needed, take the points along axis
4885 if supports_nans:
4886 # may contain nan, which would sort to the end
4887 arr.partition(
4888 concatenate((virtual_indexes.ravel(), [-1])), axis=0,
4889 )
4890 slices_having_nans = np.isnan(arr[-1, ...])
4891 else:
4892 # cannot contain nan
4893 arr.partition(virtual_indexes.ravel(), axis=0)
4894 slices_having_nans = np.array(False, dtype=bool)
4895 result = take(arr, virtual_indexes, axis=0, out=out)
4896 else:
4897 previous_indexes, next_indexes = _get_indexes(arr,
4898 virtual_indexes,
4899 values_count)
4900 # --- Sorting
4901 arr.partition(
4902 np.unique(np.concatenate(([0, -1],
4903 previous_indexes.ravel(),
4904 next_indexes.ravel(),
4905 ))),
4906 axis=0)
4907 if supports_nans:
4908 slices_having_nans = np.isnan(arr[-1, ...])
4909 else:
4910 slices_having_nans = None
4911 # --- Get values from indexes
4912 previous = arr[previous_indexes]
4913 next = arr[next_indexes]
4914 # --- Linear interpolation
4915 gamma = _get_gamma(virtual_indexes, previous_indexes, method_props)
4916 result_shape = virtual_indexes.shape + (1,) * (arr.ndim - 1)
4917 gamma = gamma.reshape(result_shape)
4918 result = _lerp(previous,
4919 next,
4920 gamma,
4921 out=out)
4922 else:
4923 # Weighted case
4924 # This implements method="inverted_cdf", the only supported weighted
4925 # method, which needs to sort anyway.
4926 weights = np.asanyarray(weights)
4927 if axis != 0:
4928 weights = np.moveaxis(weights, axis, destination=0)
4929 index_array = np.argsort(arr, axis=0, kind="stable")
4931 # arr = arr[index_array, ...] # but this adds trailing dimensions of
4932 # 1.
4933 arr = np.take_along_axis(arr, index_array, axis=0)
4934 if weights.shape == arr.shape:
4935 weights = np.take_along_axis(weights, index_array, axis=0)
4936 else:
4937 # weights is 1d
4938 weights = weights.reshape(-1)[index_array, ...]
4940 if supports_nans:
4941 # may contain nan, which would sort to the end
4942 slices_having_nans = np.isnan(arr[-1, ...])
4943 else:
4944 # cannot contain nan
4945 slices_having_nans = np.array(False, dtype=bool)
4947 # We use the weights to calculate the empirical cumulative
4948 # distribution function cdf
4949 cdf = weights.cumsum(axis=0, dtype=np.float64)
4950 cdf /= cdf[-1, ...] # normalization to 1
4951 # Search index i such that
4952 # sum(weights[j], j=0..i-1) < quantile <= sum(weights[j], j=0..i)
4953 # is then equivalent to
4954 # cdf[i-1] < quantile <= cdf[i]
4955 # Unfortunately, searchsorted only accepts 1-d arrays as first
4956 # argument, so we will need to iterate over dimensions.
4958 # Without the following cast, searchsorted can return surprising
4959 # results, e.g.
4960 # np.searchsorted(np.array([0.2, 0.4, 0.6, 0.8, 1.]),
4961 # np.array(0.4, dtype=np.float32), side="left")
4962 # returns 2 instead of 1 because 0.4 is not binary representable.
4963 if quantiles.dtype.kind == "f":
4964 cdf = cdf.astype(quantiles.dtype)
4966 def find_cdf_1d(arr, cdf):
4967 indices = np.searchsorted(cdf, quantiles, side="left")
4968 # We might have reached the maximum with i = len(arr), e.g. for
4969 # quantiles = 1, and need to cut it to len(arr) - 1.
4970 indices = minimum(indices, values_count - 1)
4971 result = take(arr, indices, axis=0)
4972 return result
4974 r_shape = arr.shape[1:]
4975 if quantiles.ndim > 0:
4976 r_shape = quantiles.shape + r_shape
4977 if out is None:
4978 result = np.empty_like(arr, shape=r_shape)
4979 else:
4980 if out.shape != r_shape:
4981 msg = (f"Wrong shape of argument 'out', shape={r_shape} is "
4982 f"required; got shape={out.shape}.")
4983 raise ValueError(msg)
4984 result = out
4986 # See apply_along_axis, which we do for axis=0. Note that Ni = (,)
4987 # always, so we remove it here.
4988 Nk = arr.shape[1:]
4989 for kk in np.ndindex(Nk):
4990 result[(...,) + kk] = find_cdf_1d(
4991 arr[np.s_[:, ] + kk], cdf[np.s_[:, ] + kk]
4992 )
4994 # Make result the same as in unweighted inverted_cdf.
4995 if result.shape == () and result.dtype == np.dtype("O"):
4996 result = result.item()
4998 if np.any(slices_having_nans):
4999 if result.ndim == 0 and out is None:
5000 # can't write to a scalar, but indexing will be correct
5001 result = arr[-1]
5002 else:
5003 np.copyto(result, arr[-1, ...], where=slices_having_nans)
5004 return result
5007def _trapezoid_dispatcher(y, x=None, dx=None, axis=None):
5008 return (y, x)
5011@array_function_dispatch(_trapezoid_dispatcher)
5012def trapezoid(y, x=None, dx=1.0, axis=-1):
5013 r"""
5014 Integrate along the given axis using the composite trapezoidal rule.
5016 If `x` is provided, the integration happens in sequence along its
5017 elements - they are not sorted.
5019 Integrate `y` (`x`) along each 1d slice on the given axis, compute
5020 :math:`\int y(x) dx`.
5021 When `x` is specified, this integrates along the parametric curve,
5022 computing :math:`\int_t y(t) dt =
5023 \int_t y(t) \left.\frac{dx}{dt}\right|_{x=x(t)} dt`.
5025 .. versionadded:: 2.0.0
5027 Parameters
5028 ----------
5029 y : array_like
5030 Input array to integrate.
5031 x : array_like, optional
5032 The sample points corresponding to the `y` values. If `x` is None,
5033 the sample points are assumed to be evenly spaced `dx` apart. The
5034 default is None.
5035 dx : scalar, optional
5036 The spacing between sample points when `x` is None. The default is 1.
5037 axis : int, optional
5038 The axis along which to integrate.
5040 Returns
5041 -------
5042 trapezoid : float or ndarray
5043 Definite integral of `y` = n-dimensional array as approximated along
5044 a single axis by the trapezoidal rule. If `y` is a 1-dimensional array,
5045 then the result is a float. If `n` is greater than 1, then the result
5046 is an `n`-1 dimensional array.
5048 See Also
5049 --------
5050 sum, cumsum
5052 Notes
5053 -----
5054 Image [2]_ illustrates trapezoidal rule -- y-axis locations of points
5055 will be taken from `y` array, by default x-axis distances between
5056 points will be 1.0, alternatively they can be provided with `x` array
5057 or with `dx` scalar. Return value will be equal to combined area under
5058 the red lines.
5061 References
5062 ----------
5063 .. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule
5065 .. [2] Illustration image:
5066 https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
5068 Examples
5069 --------
5070 Use the trapezoidal rule on evenly spaced points:
5072 >>> np.trapezoid([1, 2, 3])
5073 4.0
5075 The spacing between sample points can be selected by either the
5076 ``x`` or ``dx`` arguments:
5078 >>> np.trapezoid([1, 2, 3], x=[4, 6, 8])
5079 8.0
5080 >>> np.trapezoid([1, 2, 3], dx=2)
5081 8.0
5083 Using a decreasing ``x`` corresponds to integrating in reverse:
5085 >>> np.trapezoid([1, 2, 3], x=[8, 6, 4])
5086 -8.0
5088 More generally ``x`` is used to integrate along a parametric curve. We can
5089 estimate the integral :math:`\int_0^1 x^2 = 1/3` using:
5091 >>> x = np.linspace(0, 1, num=50)
5092 >>> y = x**2
5093 >>> np.trapezoid(y, x)
5094 0.33340274885464394
5096 Or estimate the area of a circle, noting we repeat the sample which closes
5097 the curve:
5099 >>> theta = np.linspace(0, 2 * np.pi, num=1000, endpoint=True)
5100 >>> np.trapezoid(np.cos(theta), x=np.sin(theta))
5101 3.141571941375841
5103 ``np.trapezoid`` can be applied along a specified axis to do multiple
5104 computations in one call:
5106 >>> a = np.arange(6).reshape(2, 3)
5107 >>> a
5108 array([[0, 1, 2],
5109 [3, 4, 5]])
5110 >>> np.trapezoid(a, axis=0)
5111 array([1.5, 2.5, 3.5])
5112 >>> np.trapezoid(a, axis=1)
5113 array([2., 8.])
5114 """
5116 y = asanyarray(y)
5117 if x is None:
5118 d = dx
5119 else:
5120 x = asanyarray(x)
5121 if x.ndim == 1:
5122 d = diff(x)
5123 # reshape to correct shape
5124 shape = [1]*y.ndim
5125 shape[axis] = d.shape[0]
5126 d = d.reshape(shape)
5127 else:
5128 d = diff(x, axis=axis)
5129 nd = y.ndim
5130 slice1 = [slice(None)]*nd
5131 slice2 = [slice(None)]*nd
5132 slice1[axis] = slice(1, None)
5133 slice2[axis] = slice(None, -1)
5134 try:
5135 ret = (d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0).sum(axis)
5136 except ValueError:
5137 # Operations didn't work, cast to ndarray
5138 d = np.asarray(d)
5139 y = np.asarray(y)
5140 ret = add.reduce(d * (y[tuple(slice1)]+y[tuple(slice2)])/2.0, axis)
5141 return ret
5144@set_module('numpy')
5145def trapz(y, x=None, dx=1.0, axis=-1):
5146 """
5147 `trapz` is deprecated in NumPy 2.0.
5149 Please use `trapezoid` instead, or one of the numerical integration
5150 functions in `scipy.integrate`.
5151 """
5152 # Deprecated in NumPy 2.0, 2023-08-18
5153 warnings.warn(
5154 "`trapz` is deprecated. Use `trapezoid` instead, or one of the "
5155 "numerical integration functions in `scipy.integrate`.",
5156 DeprecationWarning,
5157 stacklevel=2
5158 )
5159 return trapezoid(y, x=x, dx=dx, axis=axis)
5162def _meshgrid_dispatcher(*xi, copy=None, sparse=None, indexing=None):
5163 return xi
5166# Based on scitools meshgrid
5167@array_function_dispatch(_meshgrid_dispatcher)
5168def meshgrid(*xi, copy=True, sparse=False, indexing='xy'):
5169 """
5170 Return a tuple of coordinate matrices from coordinate vectors.
5172 Make N-D coordinate arrays for vectorized evaluations of
5173 N-D scalar/vector fields over N-D grids, given
5174 one-dimensional coordinate arrays x1, x2,..., xn.
5176 .. versionchanged:: 1.9
5177 1-D and 0-D cases are allowed.
5179 Parameters
5180 ----------
5181 x1, x2,..., xn : array_like
5182 1-D arrays representing the coordinates of a grid.
5183 indexing : {'xy', 'ij'}, optional
5184 Cartesian ('xy', default) or matrix ('ij') indexing of output.
5185 See Notes for more details.
5187 .. versionadded:: 1.7.0
5188 sparse : bool, optional
5189 If True the shape of the returned coordinate array for dimension *i*
5190 is reduced from ``(N1, ..., Ni, ... Nn)`` to
5191 ``(1, ..., 1, Ni, 1, ..., 1)``. These sparse coordinate grids are
5192 intended to be use with :ref:`basics.broadcasting`. When all
5193 coordinates are used in an expression, broadcasting still leads to a
5194 fully-dimensonal result array.
5196 Default is False.
5198 .. versionadded:: 1.7.0
5199 copy : bool, optional
5200 If False, a view into the original arrays are returned in order to
5201 conserve memory. Default is True. Please note that
5202 ``sparse=False, copy=False`` will likely return non-contiguous
5203 arrays. Furthermore, more than one element of a broadcast array
5204 may refer to a single memory location. If you need to write to the
5205 arrays, make copies first.
5207 .. versionadded:: 1.7.0
5209 Returns
5210 -------
5211 X1, X2,..., XN : tuple of ndarrays
5212 For vectors `x1`, `x2`,..., `xn` with lengths ``Ni=len(xi)``,
5213 returns ``(N1, N2, N3,..., Nn)`` shaped arrays if indexing='ij'
5214 or ``(N2, N1, N3,..., Nn)`` shaped arrays if indexing='xy'
5215 with the elements of `xi` repeated to fill the matrix along
5216 the first dimension for `x1`, the second for `x2` and so on.
5218 Notes
5219 -----
5220 This function supports both indexing conventions through the indexing
5221 keyword argument. Giving the string 'ij' returns a meshgrid with
5222 matrix indexing, while 'xy' returns a meshgrid with Cartesian indexing.
5223 In the 2-D case with inputs of length M and N, the outputs are of shape
5224 (N, M) for 'xy' indexing and (M, N) for 'ij' indexing. In the 3-D case
5225 with inputs of length M, N and P, outputs are of shape (N, M, P) for
5226 'xy' indexing and (M, N, P) for 'ij' indexing. The difference is
5227 illustrated by the following code snippet::
5229 xv, yv = np.meshgrid(x, y, indexing='ij')
5230 for i in range(nx):
5231 for j in range(ny):
5232 # treat xv[i,j], yv[i,j]
5234 xv, yv = np.meshgrid(x, y, indexing='xy')
5235 for i in range(nx):
5236 for j in range(ny):
5237 # treat xv[j,i], yv[j,i]
5239 In the 1-D and 0-D case, the indexing and sparse keywords have no effect.
5241 See Also
5242 --------
5243 mgrid : Construct a multi-dimensional "meshgrid" using indexing notation.
5244 ogrid : Construct an open multi-dimensional "meshgrid" using indexing
5245 notation.
5246 :ref:`how-to-index`
5248 Examples
5249 --------
5250 >>> nx, ny = (3, 2)
5251 >>> x = np.linspace(0, 1, nx)
5252 >>> y = np.linspace(0, 1, ny)
5253 >>> xv, yv = np.meshgrid(x, y)
5254 >>> xv
5255 array([[0. , 0.5, 1. ],
5256 [0. , 0.5, 1. ]])
5257 >>> yv
5258 array([[0., 0., 0.],
5259 [1., 1., 1.]])
5261 The result of `meshgrid` is a coordinate grid:
5263 >>> import matplotlib.pyplot as plt
5264 >>> plt.plot(xv, yv, marker='o', color='k', linestyle='none')
5265 >>> plt.show()
5267 You can create sparse output arrays to save memory and computation time.
5269 >>> xv, yv = np.meshgrid(x, y, sparse=True)
5270 >>> xv
5271 array([[0. , 0.5, 1. ]])
5272 >>> yv
5273 array([[0.],
5274 [1.]])
5276 `meshgrid` is very useful to evaluate functions on a grid. If the
5277 function depends on all coordinates, both dense and sparse outputs can be
5278 used.
5280 >>> x = np.linspace(-5, 5, 101)
5281 >>> y = np.linspace(-5, 5, 101)
5282 >>> # full coordinate arrays
5283 >>> xx, yy = np.meshgrid(x, y)
5284 >>> zz = np.sqrt(xx**2 + yy**2)
5285 >>> xx.shape, yy.shape, zz.shape
5286 ((101, 101), (101, 101), (101, 101))
5287 >>> # sparse coordinate arrays
5288 >>> xs, ys = np.meshgrid(x, y, sparse=True)
5289 >>> zs = np.sqrt(xs**2 + ys**2)
5290 >>> xs.shape, ys.shape, zs.shape
5291 ((1, 101), (101, 1), (101, 101))
5292 >>> np.array_equal(zz, zs)
5293 True
5295 >>> h = plt.contourf(x, y, zs)
5296 >>> plt.axis('scaled')
5297 >>> plt.colorbar()
5298 >>> plt.show()
5299 """
5300 ndim = len(xi)
5302 if indexing not in ['xy', 'ij']:
5303 raise ValueError(
5304 "Valid values for `indexing` are 'xy' and 'ij'.")
5306 s0 = (1,) * ndim
5307 output = [np.asanyarray(x).reshape(s0[:i] + (-1,) + s0[i + 1:])
5308 for i, x in enumerate(xi)]
5310 if indexing == 'xy' and ndim > 1:
5311 # switch first and second axis
5312 output[0].shape = (1, -1) + s0[2:]
5313 output[1].shape = (-1, 1) + s0[2:]
5315 if not sparse:
5316 # Return the full N-D matrix (not only the 1-D vector)
5317 output = np.broadcast_arrays(*output, subok=True)
5319 if copy:
5320 output = tuple(x.copy() for x in output)
5322 return output
5325def _delete_dispatcher(arr, obj, axis=None):
5326 return (arr, obj)
5329@array_function_dispatch(_delete_dispatcher)
5330def delete(arr, obj, axis=None):
5331 """
5332 Return a new array with sub-arrays along an axis deleted. For a one
5333 dimensional array, this returns those entries not returned by
5334 `arr[obj]`.
5336 Parameters
5337 ----------
5338 arr : array_like
5339 Input array.
5340 obj : slice, int or array of ints
5341 Indicate indices of sub-arrays to remove along the specified axis.
5343 .. versionchanged:: 1.19.0
5344 Boolean indices are now treated as a mask of elements to remove,
5345 rather than being cast to the integers 0 and 1.
5347 axis : int, optional
5348 The axis along which to delete the subarray defined by `obj`.
5349 If `axis` is None, `obj` is applied to the flattened array.
5351 Returns
5352 -------
5353 out : ndarray
5354 A copy of `arr` with the elements specified by `obj` removed. Note
5355 that `delete` does not occur in-place. If `axis` is None, `out` is
5356 a flattened array.
5358 See Also
5359 --------
5360 insert : Insert elements into an array.
5361 append : Append elements at the end of an array.
5363 Notes
5364 -----
5365 Often it is preferable to use a boolean mask. For example:
5367 >>> arr = np.arange(12) + 1
5368 >>> mask = np.ones(len(arr), dtype=bool)
5369 >>> mask[[0,2,4]] = False
5370 >>> result = arr[mask,...]
5372 Is equivalent to ``np.delete(arr, [0,2,4], axis=0)``, but allows further
5373 use of `mask`.
5375 Examples
5376 --------
5377 >>> arr = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]])
5378 >>> arr
5379 array([[ 1, 2, 3, 4],
5380 [ 5, 6, 7, 8],
5381 [ 9, 10, 11, 12]])
5382 >>> np.delete(arr, 1, 0)
5383 array([[ 1, 2, 3, 4],
5384 [ 9, 10, 11, 12]])
5386 >>> np.delete(arr, np.s_[::2], 1)
5387 array([[ 2, 4],
5388 [ 6, 8],
5389 [10, 12]])
5390 >>> np.delete(arr, [1,3,5], None)
5391 array([ 1, 3, 5, 7, 8, 9, 10, 11, 12])
5393 """
5394 conv = _array_converter(arr)
5395 arr, = conv.as_arrays(subok=False)
5397 ndim = arr.ndim
5398 arrorder = 'F' if arr.flags.fnc else 'C'
5399 if axis is None:
5400 if ndim != 1:
5401 arr = arr.ravel()
5402 # needed for np.matrix, which is still not 1d after being ravelled
5403 ndim = arr.ndim
5404 axis = ndim - 1
5405 else:
5406 axis = normalize_axis_index(axis, ndim)
5408 slobj = [slice(None)]*ndim
5409 N = arr.shape[axis]
5410 newshape = list(arr.shape)
5412 if isinstance(obj, slice):
5413 start, stop, step = obj.indices(N)
5414 xr = range(start, stop, step)
5415 numtodel = len(xr)
5417 if numtodel <= 0:
5418 return conv.wrap(arr.copy(order=arrorder), to_scalar=False)
5420 # Invert if step is negative:
5421 if step < 0:
5422 step = -step
5423 start = xr[-1]
5424 stop = xr[0] + 1
5426 newshape[axis] -= numtodel
5427 new = empty(newshape, arr.dtype, arrorder)
5428 # copy initial chunk
5429 if start == 0:
5430 pass
5431 else:
5432 slobj[axis] = slice(None, start)
5433 new[tuple(slobj)] = arr[tuple(slobj)]
5434 # copy end chunk
5435 if stop == N:
5436 pass
5437 else:
5438 slobj[axis] = slice(stop-numtodel, None)
5439 slobj2 = [slice(None)]*ndim
5440 slobj2[axis] = slice(stop, None)
5441 new[tuple(slobj)] = arr[tuple(slobj2)]
5442 # copy middle pieces
5443 if step == 1:
5444 pass
5445 else: # use array indexing.
5446 keep = ones(stop-start, dtype=bool)
5447 keep[:stop-start:step] = False
5448 slobj[axis] = slice(start, stop-numtodel)
5449 slobj2 = [slice(None)]*ndim
5450 slobj2[axis] = slice(start, stop)
5451 arr = arr[tuple(slobj2)]
5452 slobj2[axis] = keep
5453 new[tuple(slobj)] = arr[tuple(slobj2)]
5455 return conv.wrap(new, to_scalar=False)
5457 if isinstance(obj, (int, integer)) and not isinstance(obj, bool):
5458 single_value = True
5459 else:
5460 single_value = False
5461 _obj = obj
5462 obj = np.asarray(obj)
5463 # `size == 0` to allow empty lists similar to indexing, but (as there)
5464 # is really too generic:
5465 if obj.size == 0 and not isinstance(_obj, np.ndarray):
5466 obj = obj.astype(intp)
5467 elif obj.size == 1 and obj.dtype.kind in "ui":
5468 # For a size 1 integer array we can use the single-value path
5469 # (most dtypes, except boolean, should just fail later).
5470 obj = obj.item()
5471 single_value = True
5473 if single_value:
5474 # optimization for a single value
5475 if (obj < -N or obj >= N):
5476 raise IndexError(
5477 "index %i is out of bounds for axis %i with "
5478 "size %i" % (obj, axis, N))
5479 if (obj < 0):
5480 obj += N
5481 newshape[axis] -= 1
5482 new = empty(newshape, arr.dtype, arrorder)
5483 slobj[axis] = slice(None, obj)
5484 new[tuple(slobj)] = arr[tuple(slobj)]
5485 slobj[axis] = slice(obj, None)
5486 slobj2 = [slice(None)]*ndim
5487 slobj2[axis] = slice(obj+1, None)
5488 new[tuple(slobj)] = arr[tuple(slobj2)]
5489 else:
5490 if obj.dtype == bool:
5491 if obj.shape != (N,):
5492 raise ValueError('boolean array argument obj to delete '
5493 'must be one dimensional and match the axis '
5494 'length of {}'.format(N))
5496 # optimization, the other branch is slower
5497 keep = ~obj
5498 else:
5499 keep = ones(N, dtype=bool)
5500 keep[obj,] = False
5502 slobj[axis] = keep
5503 new = arr[tuple(slobj)]
5505 return conv.wrap(new, to_scalar=False)
5508def _insert_dispatcher(arr, obj, values, axis=None):
5509 return (arr, obj, values)
5512@array_function_dispatch(_insert_dispatcher)
5513def insert(arr, obj, values, axis=None):
5514 """
5515 Insert values along the given axis before the given indices.
5517 Parameters
5518 ----------
5519 arr : array_like
5520 Input array.
5521 obj : int, slice or sequence of ints
5522 Object that defines the index or indices before which `values` is
5523 inserted.
5525 .. versionadded:: 1.8.0
5527 Support for multiple insertions when `obj` is a single scalar or a
5528 sequence with one element (similar to calling insert multiple
5529 times).
5530 values : array_like
5531 Values to insert into `arr`. If the type of `values` is different
5532 from that of `arr`, `values` is converted to the type of `arr`.
5533 `values` should be shaped so that ``arr[...,obj,...] = values``
5534 is legal.
5535 axis : int, optional
5536 Axis along which to insert `values`. If `axis` is None then `arr`
5537 is flattened first.
5539 Returns
5540 -------
5541 out : ndarray
5542 A copy of `arr` with `values` inserted. Note that `insert`
5543 does not occur in-place: a new array is returned. If
5544 `axis` is None, `out` is a flattened array.
5546 See Also
5547 --------
5548 append : Append elements at the end of an array.
5549 concatenate : Join a sequence of arrays along an existing axis.
5550 delete : Delete elements from an array.
5552 Notes
5553 -----
5554 Note that for higher dimensional inserts ``obj=0`` behaves very different
5555 from ``obj=[0]`` just like ``arr[:,0,:] = values`` is different from
5556 ``arr[:,[0],:] = values``.
5558 Examples
5559 --------
5560 >>> a = np.array([[1, 1], [2, 2], [3, 3]])
5561 >>> a
5562 array([[1, 1],
5563 [2, 2],
5564 [3, 3]])
5565 >>> np.insert(a, 1, 5)
5566 array([1, 5, 1, ..., 2, 3, 3])
5567 >>> np.insert(a, 1, 5, axis=1)
5568 array([[1, 5, 1],
5569 [2, 5, 2],
5570 [3, 5, 3]])
5572 Difference between sequence and scalars:
5574 >>> np.insert(a, [1], [[1],[2],[3]], axis=1)
5575 array([[1, 1, 1],
5576 [2, 2, 2],
5577 [3, 3, 3]])
5578 >>> np.array_equal(np.insert(a, 1, [1, 2, 3], axis=1),
5579 ... np.insert(a, [1], [[1],[2],[3]], axis=1))
5580 True
5582 >>> b = a.flatten()
5583 >>> b
5584 array([1, 1, 2, 2, 3, 3])
5585 >>> np.insert(b, [2, 2], [5, 6])
5586 array([1, 1, 5, ..., 2, 3, 3])
5588 >>> np.insert(b, slice(2, 4), [5, 6])
5589 array([1, 1, 5, ..., 2, 3, 3])
5591 >>> np.insert(b, [2, 2], [7.13, False]) # type casting
5592 array([1, 1, 7, ..., 2, 3, 3])
5594 >>> x = np.arange(8).reshape(2, 4)
5595 >>> idx = (1, 3)
5596 >>> np.insert(x, idx, 999, axis=1)
5597 array([[ 0, 999, 1, 2, 999, 3],
5598 [ 4, 999, 5, 6, 999, 7]])
5600 """
5601 conv = _array_converter(arr)
5602 arr, = conv.as_arrays(subok=False)
5604 ndim = arr.ndim
5605 arrorder = 'F' if arr.flags.fnc else 'C'
5606 if axis is None:
5607 if ndim != 1:
5608 arr = arr.ravel()
5609 # needed for np.matrix, which is still not 1d after being ravelled
5610 ndim = arr.ndim
5611 axis = ndim - 1
5612 else:
5613 axis = normalize_axis_index(axis, ndim)
5614 slobj = [slice(None)]*ndim
5615 N = arr.shape[axis]
5616 newshape = list(arr.shape)
5618 if isinstance(obj, slice):
5619 # turn it into a range object
5620 indices = arange(*obj.indices(N), dtype=intp)
5621 else:
5622 # need to copy obj, because indices will be changed in-place
5623 indices = np.array(obj)
5624 if indices.dtype == bool:
5625 # See also delete
5626 # 2012-10-11, NumPy 1.8
5627 warnings.warn(
5628 "in the future insert will treat boolean arrays and "
5629 "array-likes as a boolean index instead of casting it to "
5630 "integer", FutureWarning, stacklevel=2)
5631 indices = indices.astype(intp)
5632 # Code after warning period:
5633 #if obj.ndim != 1:
5634 # raise ValueError('boolean array argument obj to insert '
5635 # 'must be one dimensional')
5636 #indices = np.flatnonzero(obj)
5637 elif indices.ndim > 1:
5638 raise ValueError(
5639 "index array argument obj to insert must be one dimensional "
5640 "or scalar")
5641 if indices.size == 1:
5642 index = indices.item()
5643 if index < -N or index > N:
5644 raise IndexError(f"index {obj} is out of bounds for axis {axis} "
5645 f"with size {N}")
5646 if (index < 0):
5647 index += N
5649 # There are some object array corner cases here, but we cannot avoid
5650 # that:
5651 values = array(values, copy=None, ndmin=arr.ndim, dtype=arr.dtype)
5652 if indices.ndim == 0:
5653 # broadcasting is very different here, since a[:,0,:] = ... behaves
5654 # very different from a[:,[0],:] = ...! This changes values so that
5655 # it works likes the second case. (here a[:,0:1,:])
5656 values = np.moveaxis(values, 0, axis)
5657 numnew = values.shape[axis]
5658 newshape[axis] += numnew
5659 new = empty(newshape, arr.dtype, arrorder)
5660 slobj[axis] = slice(None, index)
5661 new[tuple(slobj)] = arr[tuple(slobj)]
5662 slobj[axis] = slice(index, index+numnew)
5663 new[tuple(slobj)] = values
5664 slobj[axis] = slice(index+numnew, None)
5665 slobj2 = [slice(None)] * ndim
5666 slobj2[axis] = slice(index, None)
5667 new[tuple(slobj)] = arr[tuple(slobj2)]
5669 return conv.wrap(new, to_scalar=False)
5671 elif indices.size == 0 and not isinstance(obj, np.ndarray):
5672 # Can safely cast the empty list to intp
5673 indices = indices.astype(intp)
5675 indices[indices < 0] += N
5677 numnew = len(indices)
5678 order = indices.argsort(kind='mergesort') # stable sort
5679 indices[order] += np.arange(numnew)
5681 newshape[axis] += numnew
5682 old_mask = ones(newshape[axis], dtype=bool)
5683 old_mask[indices] = False
5685 new = empty(newshape, arr.dtype, arrorder)
5686 slobj2 = [slice(None)]*ndim
5687 slobj[axis] = indices
5688 slobj2[axis] = old_mask
5689 new[tuple(slobj)] = values
5690 new[tuple(slobj2)] = arr
5692 return conv.wrap(new, to_scalar=False)
5695def _append_dispatcher(arr, values, axis=None):
5696 return (arr, values)
5699@array_function_dispatch(_append_dispatcher)
5700def append(arr, values, axis=None):
5701 """
5702 Append values to the end of an array.
5704 Parameters
5705 ----------
5706 arr : array_like
5707 Values are appended to a copy of this array.
5708 values : array_like
5709 These values are appended to a copy of `arr`. It must be of the
5710 correct shape (the same shape as `arr`, excluding `axis`). If
5711 `axis` is not specified, `values` can be any shape and will be
5712 flattened before use.
5713 axis : int, optional
5714 The axis along which `values` are appended. If `axis` is not
5715 given, both `arr` and `values` are flattened before use.
5717 Returns
5718 -------
5719 append : ndarray
5720 A copy of `arr` with `values` appended to `axis`. Note that
5721 `append` does not occur in-place: a new array is allocated and
5722 filled. If `axis` is None, `out` is a flattened array.
5724 See Also
5725 --------
5726 insert : Insert elements into an array.
5727 delete : Delete elements from an array.
5729 Examples
5730 --------
5731 >>> np.append([1, 2, 3], [[4, 5, 6], [7, 8, 9]])
5732 array([1, 2, 3, ..., 7, 8, 9])
5734 When `axis` is specified, `values` must have the correct shape.
5736 >>> np.append([[1, 2, 3], [4, 5, 6]], [[7, 8, 9]], axis=0)
5737 array([[1, 2, 3],
5738 [4, 5, 6],
5739 [7, 8, 9]])
5740 >>> np.append([[1, 2, 3], [4, 5, 6]], [7, 8, 9], axis=0)
5741 Traceback (most recent call last):
5742 ...
5743 ValueError: all the input arrays must have same number of dimensions, but
5744 the array at index 0 has 2 dimension(s) and the array at index 1 has 1
5745 dimension(s)
5747 """
5748 arr = asanyarray(arr)
5749 if axis is None:
5750 if arr.ndim != 1:
5751 arr = arr.ravel()
5752 values = ravel(values)
5753 axis = arr.ndim-1
5754 return concatenate((arr, values), axis=axis)
5757def _digitize_dispatcher(x, bins, right=None):
5758 return (x, bins)
5761@array_function_dispatch(_digitize_dispatcher)
5762def digitize(x, bins, right=False):
5763 """
5764 Return the indices of the bins to which each value in input array belongs.
5766 ========= ============= ============================
5767 `right` order of bins returned index `i` satisfies
5768 ========= ============= ============================
5769 ``False`` increasing ``bins[i-1] <= x < bins[i]``
5770 ``True`` increasing ``bins[i-1] < x <= bins[i]``
5771 ``False`` decreasing ``bins[i-1] > x >= bins[i]``
5772 ``True`` decreasing ``bins[i-1] >= x > bins[i]``
5773 ========= ============= ============================
5775 If values in `x` are beyond the bounds of `bins`, 0 or ``len(bins)`` is
5776 returned as appropriate.
5778 Parameters
5779 ----------
5780 x : array_like
5781 Input array to be binned. Prior to NumPy 1.10.0, this array had to
5782 be 1-dimensional, but can now have any shape.
5783 bins : array_like
5784 Array of bins. It has to be 1-dimensional and monotonic.
5785 right : bool, optional
5786 Indicating whether the intervals include the right or the left bin
5787 edge. Default behavior is (right==False) indicating that the interval
5788 does not include the right edge. The left bin end is open in this
5789 case, i.e., bins[i-1] <= x < bins[i] is the default behavior for
5790 monotonically increasing bins.
5792 Returns
5793 -------
5794 indices : ndarray of ints
5795 Output array of indices, of same shape as `x`.
5797 Raises
5798 ------
5799 ValueError
5800 If `bins` is not monotonic.
5801 TypeError
5802 If the type of the input is complex.
5804 See Also
5805 --------
5806 bincount, histogram, unique, searchsorted
5808 Notes
5809 -----
5810 If values in `x` are such that they fall outside the bin range,
5811 attempting to index `bins` with the indices that `digitize` returns
5812 will result in an IndexError.
5814 .. versionadded:: 1.10.0
5816 `numpy.digitize` is implemented in terms of `numpy.searchsorted`.
5817 This means that a binary search is used to bin the values, which scales
5818 much better for larger number of bins than the previous linear search.
5819 It also removes the requirement for the input array to be 1-dimensional.
5821 For monotonically *increasing* `bins`, the following are equivalent::
5823 np.digitize(x, bins, right=True)
5824 np.searchsorted(bins, x, side='left')
5826 Note that as the order of the arguments are reversed, the side must be too.
5827 The `searchsorted` call is marginally faster, as it does not do any
5828 monotonicity checks. Perhaps more importantly, it supports all dtypes.
5830 Examples
5831 --------
5832 >>> x = np.array([0.2, 6.4, 3.0, 1.6])
5833 >>> bins = np.array([0.0, 1.0, 2.5, 4.0, 10.0])
5834 >>> inds = np.digitize(x, bins)
5835 >>> inds
5836 array([1, 4, 3, 2])
5837 >>> for n in range(x.size):
5838 ... print(bins[inds[n]-1], "<=", x[n], "<", bins[inds[n]])
5839 ...
5840 0.0 <= 0.2 < 1.0
5841 4.0 <= 6.4 < 10.0
5842 2.5 <= 3.0 < 4.0
5843 1.0 <= 1.6 < 2.5
5845 >>> x = np.array([1.2, 10.0, 12.4, 15.5, 20.])
5846 >>> bins = np.array([0, 5, 10, 15, 20])
5847 >>> np.digitize(x,bins,right=True)
5848 array([1, 2, 3, 4, 4])
5849 >>> np.digitize(x,bins,right=False)
5850 array([1, 3, 3, 4, 5])
5851 """
5852 x = _nx.asarray(x)
5853 bins = _nx.asarray(bins)
5855 # here for compatibility, searchsorted below is happy to take this
5856 if np.issubdtype(x.dtype, _nx.complexfloating):
5857 raise TypeError("x may not be complex")
5859 mono = _monotonicity(bins)
5860 if mono == 0:
5861 raise ValueError("bins must be monotonically increasing or decreasing")
5863 # this is backwards because the arguments below are swapped
5864 side = 'left' if right else 'right'
5865 if mono == -1:
5866 # reverse the bins, and invert the results
5867 return len(bins) - _nx.searchsorted(bins[::-1], x, side=side)
5868 else:
5869 return _nx.searchsorted(bins, x, side=side)