codekingpro/portable-devtools
115k
1import builtins
2import collections.abc
3import functools
4import re
5import warnings
6
7import numpy as np
8import numpy._core.numeric as _nx
9from numpy._core import overrides, transpose
10from numpy._core._multiarray_umath import _array_converter
11from numpy._core.fromnumeric import any, mean, nonzero, partition, ravel, sum
12from numpy._core.multiarray import (
13 _monotonicity,
14 _place,
15 bincount,
16 interp as compiled_interp,
17 interp_complex as compiled_interp_complex,
18 normalize_axis_index,
19)
20from numpy._core.numeric import (
21 absolute,
22 arange,
23 array,
24 asanyarray,
25 asarray,
26 concatenate,
27 dot,
28 empty,
29 integer,
30 intp,
31 isscalar,
32 ndarray,
33 ones,
34 take,
35 where,
36 zeros_like,
37)
38from numpy._core.numerictypes import typecodes
39from numpy._core.umath import (
40 add,
41 arctan2,
42 cos,
43 exp,
44 floor,
45 frompyfunc,
46 less_equal,
47 minimum,
48 mod,
49 not_equal,
50 pi,
51 sin,
52 sqrt,
53 subtract,
54)
55from numpy._utils import set_module
56
57# needed in this module for compatibility
58from numpy.lib._histograms_impl import histogram, histogramdd # noqa: F401
59from numpy.lib._twodim_base_impl import diag
60
61array_function_dispatch = functools.partial(
62 overrides.array_function_dispatch, module='numpy')
63
64
65__all__ = [
66 'select', 'piecewise', 'trim_zeros', 'copy', 'iterable', 'percentile',
67 'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'flip',
68 'rot90', 'extract', 'place', 'vectorize', 'asarray_chkfinite', 'average',
69 'bincount', 'digitize', 'cov', 'corrcoef',
70 'median', 'sinc', 'hamming', 'hanning', 'bartlett',
71 'blackman', 'kaiser', 'trapezoid', 'i0',
72 'meshgrid', 'delete', 'insert', 'append', 'interp',
73 'quantile'
74 ]
75
76# _QuantileMethods is a dictionary listing all the supported methods to
77# compute quantile/percentile.
78#
79# Below virtual_index refers to the index of the element where the percentile
80# would be found in the sorted sample.
81# When the sample contains exactly the percentile wanted, the virtual_index is
82# an integer to the index of this element.
83# When the percentile wanted is in between two elements, the virtual_index
84# is made of a integer part (a.k.a 'i' or 'left') and a fractional part
85# (a.k.a 'g' or 'gamma')
86#
87# Each method in _QuantileMethods has two properties
88# get_virtual_index : Callable
89# The function used to compute the virtual_index.
90# fix_gamma : Callable
91# A function used for discrete methods to force the index to a specific value.
92_QuantileMethods = {
93 # --- HYNDMAN and FAN METHODS
94 # Discrete methods
95 'inverted_cdf': {
96 'get_virtual_index': lambda n, quantiles: _inverted_cdf(n, quantiles), # noqa: PLW0108
97 'fix_gamma': None, # should never be called
98 },
99 'averaged_inverted_cdf': {
100 'get_virtual_index': lambda n, quantiles: (n * quantiles) - 1,
101 'fix_gamma': lambda gamma, _: _get_gamma_mask(
102 shape=gamma.shape,
103 default_value=1.,
104 conditioned_value=0.5,
105 where=gamma == 0),
106 },
107 'closest_observation': {
108 'get_virtual_index': lambda n, quantiles: _closest_observation(n, quantiles), # noqa: PLW0108
109 'fix_gamma': None, # should never be called
110 },
111 # Continuous methods
112 'interpolated_inverted_cdf': {
113 'get_virtual_index': lambda n, quantiles:
114 _compute_virtual_index(n, quantiles, 0, 1),
115 'fix_gamma': lambda gamma, _: gamma,
116 },
117 'hazen': {
118 'get_virtual_index': lambda n, quantiles:
119 _compute_virtual_index(n, quantiles, 0.5, 0.5),
120 'fix_gamma': lambda gamma, _: gamma,
121 },
122 'weibull': {
123 'get_virtual_index': lambda n, quantiles:
124 _compute_virtual_index(n, quantiles, 0, 0),
125 'fix_gamma': lambda gamma, _: gamma,
126 },
127 # Default method.
128 # To avoid some rounding issues, `(n-1) * quantiles` is preferred to
129 # `_compute_virtual_index(n, quantiles, 1, 1)`.
130 # They are mathematically equivalent.
131 'linear': {
132 'get_virtual_index': lambda n, quantiles: (n - 1) * quantiles,
133 'fix_gamma': lambda gamma, _: gamma,
134 },
135 'median_unbiased': {
136 'get_virtual_index': lambda n, quantiles:
137 _compute_virtual_index(n, quantiles, 1 / 3.0, 1 / 3.0),
138 'fix_gamma': lambda gamma, _: gamma,
139 },
140 'normal_unbiased': {
141 'get_virtual_index': lambda n, quantiles:
142 _compute_virtual_index(n, quantiles, 3 / 8.0, 3 / 8.0),
143 'fix_gamma': lambda gamma, _: gamma,
144 },
145 # --- OTHER METHODS
146 'lower': {
147 'get_virtual_index': lambda n, quantiles: np.floor(
148 (n - 1) * quantiles).astype(np.intp),
149 'fix_gamma': None, # should never be called, index dtype is int
150 },
151 'higher': {
152 'get_virtual_index': lambda n, quantiles: np.ceil(
153 (n - 1) * quantiles).astype(np.intp),
154 'fix_gamma': None, # should never be called, index dtype is int
155 },
156 'midpoint': {
157 'get_virtual_index': lambda n, quantiles: 0.5 * (
158 np.floor((n - 1) * quantiles)
159 + np.ceil((n - 1) * quantiles)),
160 'fix_gamma': lambda gamma, index: _get_gamma_mask(
161 shape=gamma.shape,
162 default_value=0.5,
163 conditioned_value=0.,
164 where=index % 1 == 0),
165 },
166 'nearest': {
167 'get_virtual_index': lambda n, quantiles: np.around(
168 (n - 1) * quantiles).astype(np.intp),
169 'fix_gamma': None,
170 # should never be called, index dtype is int
171 }}
172
173
174def _rot90_dispatcher(m, k=None, axes=None):
175 return (m,)
176
177
178@array_function_dispatch(_rot90_dispatcher)
179def rot90(m, k=1, axes=(0, 1)):
180 """
181 Rotate an array by 90 degrees in the plane specified by axes.
182
183 Rotation direction is from the first towards the second axis.
184 This means for a 2D array with the default `k` and `axes`, the
185 rotation will be counterclockwise.
186
187 Parameters
188 ----------
189 m : array_like
190 Array of two or more dimensions.
191 k : integer
192 Number of times the array is rotated by 90 degrees.
193 axes : (2,) array_like
194 The array is rotated in the plane defined by the axes.
195 Axes must be different.
196
197 Returns
198 -------
199 y : ndarray
200 A rotated view of `m`.
201
202 See Also
203 --------
204 flip : Reverse the order of elements in an array along the given axis.
205 fliplr : Flip an array horizontally.
206 flipud : Flip an array vertically.
207
208 Notes
209 -----
210 ``rot90(m, k=1, axes=(1,0))`` is the reverse of
211 ``rot90(m, k=1, axes=(0,1))``
212
213 ``rot90(m, k=1, axes=(1,0))`` is equivalent to
214 ``rot90(m, k=-1, axes=(0,1))``
215
216 Examples
217 --------
218 >>> import numpy as np
219 >>> m = np.array([[1,2],[3,4]], int)
220 >>> m
221 array([[1, 2],
222 [3, 4]])
223 >>> np.rot90(m)
224 array([[2, 4],
225 [1, 3]])
226 >>> np.rot90(m, 2)
227 array([[4, 3],
228 [2, 1]])
229 >>> m = np.arange(8).reshape((2,2,2))
230 >>> np.rot90(m, 1, (1,2))
231 array([[[1, 3],
232 [0, 2]],
233 [[5, 7],
234 [4, 6]]])
235
236 """
237 axes = tuple(axes)
238 if len(axes) != 2:
239 raise ValueError("len(axes) must be 2.")
240
241 m = asanyarray(m)
242
243 if axes[0] == axes[1] or absolute(axes[0] - axes[1]) == m.ndim:
244 raise ValueError("Axes must be different.")
245
246 if (axes[0] >= m.ndim or axes[0] < -m.ndim
247 or axes[1] >= m.ndim or axes[1] < -m.ndim):
248 raise ValueError(f"Axes={axes} out of range for array of ndim={m.ndim}.")
249
250 k %= 4
251
252 if k == 0:
253 return m[:]
254 if k == 2:
255 return flip(flip(m, axes[0]), axes[1])
256
257 axes_list = arange(0, m.ndim)
258 (axes_list[axes[0]], axes_list[axes[1]]) = (axes_list[axes[1]],
259 axes_list[axes[0]])
260
261 if k == 1:
262 return transpose(flip(m, axes[1]), axes_list)
263 else:
264 # k == 3
265 return flip(transpose(m, axes_list), axes[1])
266
267
268def _flip_dispatcher(m, axis=None):
269 return (m,)
270
271
272@array_function_dispatch(_flip_dispatcher)
273def flip(m, axis=None):
274 """
275 Reverse the order of elements in an array along the given axis.
276
277 The shape of the array is preserved, but the elements are reordered.
278
279 Parameters
280 ----------
281 m : array_like
282 Input array.
283 axis : None or int or tuple of ints, optional
284 Axis or axes along which to flip over. The default,
285 axis=None, will flip over all of the axes of the input array.
286 If axis is negative it counts from the last to the first axis.
287
288 If axis is a tuple of ints, flipping is performed on all of the axes
289 specified in the tuple.
290
291 Returns
292 -------
293 out : array_like
294 A view of `m` with the entries of axis reversed. Since a view is
295 returned, this operation is done in constant time.
296
297 See Also
298 --------
299 flipud : Flip an array vertically (axis=0).
300 fliplr : Flip an array horizontally (axis=1).
301
302 Notes
303 -----
304 flip(m, 0) is equivalent to flipud(m).
305
306 flip(m, 1) is equivalent to fliplr(m).
307
308 flip(m, n) corresponds to ``m[...,::-1,...]`` with ``::-1`` at position n.
309
310 flip(m) corresponds to ``m[::-1,::-1,...,::-1]`` with ``::-1`` at all
311 positions.
312
313 flip(m, (0, 1)) corresponds to ``m[::-1,::-1,...]`` with ``::-1`` at
314 position 0 and position 1.
315
316 Examples
317 --------
318 >>> import numpy as np
319 >>> A = np.arange(8).reshape((2,2,2))
320 >>> A
321 array([[[0, 1],
322 [2, 3]],
323 [[4, 5],
324 [6, 7]]])
325 >>> np.flip(A, 0)
326 array([[[4, 5],
327 [6, 7]],
328 [[0, 1],
329 [2, 3]]])
330 >>> np.flip(A, 1)
331 array([[[2, 3],
332 [0, 1]],
333 [[6, 7],
334 [4, 5]]])
335 >>> np.flip(A)
336 array([[[7, 6],
337 [5, 4]],
338 [[3, 2],
339 [1, 0]]])
340 >>> np.flip(A, (0, 2))
341 array([[[5, 4],
342 [7, 6]],
343 [[1, 0],
344 [3, 2]]])
345 >>> rng = np.random.default_rng()
346 >>> A = rng.normal(size=(3,4,5))
347 >>> np.all(np.flip(A,2) == A[:,:,::-1,...])
348 True
349 """
350 if not hasattr(m, 'ndim'):
351 m = asarray(m)
352 if axis is None:
353 indexer = (np.s_[::-1],) * m.ndim
354 else:
355 axis = _nx.normalize_axis_tuple(axis, m.ndim)
356 indexer = [np.s_[:]] * m.ndim
357 for ax in axis:
358 indexer[ax] = np.s_[::-1]
359 indexer = tuple(indexer)
360 return m[indexer]
361
362
363@set_module('numpy')
364def iterable(y):
365 """
366 Check whether or not an object can be iterated over.
367
368 Parameters
369 ----------
370 y : object
371 Input object.
372
373 Returns
374 -------
375 b : bool
376 Return ``True`` if the object has an iterator method or is a
377 sequence and ``False`` otherwise.
378
379
380 Examples
381 --------
382 >>> import numpy as np
383 >>> np.iterable([1, 2, 3])
384 True
385 >>> np.iterable(2)
386 False
387
388 Notes
389 -----
390 In most cases, the results of ``np.iterable(obj)`` are consistent with
391 ``isinstance(obj, collections.abc.Iterable)``. One notable exception is
392 the treatment of 0-dimensional arrays::
393
394 >>> from collections.abc import Iterable
395 >>> a = np.array(1.0) # 0-dimensional numpy array
396 >>> isinstance(a, Iterable)
397 True
398 >>> np.iterable(a)
399 False
400
401 """
402 try:
403 iter(y)
404 except TypeError:
405 return False
406 return True
407
408
409def _weights_are_valid(weights, a, axis):
410 """Validate weights array.
411
412 We assume, weights is not None.
413 """
414 wgt = np.asanyarray(weights)
415
416 # Sanity checks
417 if a.shape != wgt.shape:
418 if axis is None:
419 raise TypeError(
420 "Axis must be specified when shapes of a and weights "
421 "differ.")
422 if wgt.shape != tuple(a.shape[ax] for ax in axis):
423 raise ValueError(
424 "Shape of weights must be consistent with "
425 "shape of a along specified axis.")
426
427 # setup wgt to broadcast along axis
428 wgt = wgt.transpose(np.argsort(axis))
429 wgt = wgt.reshape(tuple((s if ax in axis else 1)
430 for ax, s in enumerate(a.shape)))
431 return wgt
432
433
434def _average_dispatcher(a, axis=None, weights=None, returned=None, *,
435 keepdims=None):
436 return (a, weights)
437
438
439@array_function_dispatch(_average_dispatcher)
440def average(a, axis=None, weights=None, returned=False, *,
441 keepdims=np._NoValue):
442 """
443 Compute the weighted average along the specified axis.
444
445 Parameters
446 ----------
447 a : array_like
448 Array containing data to be averaged. If `a` is not an array, a
449 conversion is attempted.
450 axis : None or int or tuple of ints, optional
451 Axis or axes along which to average `a`. The default,
452 `axis=None`, will average over all of the elements of the input array.
453 If axis is negative it counts from the last to the first axis.
454 If axis is a tuple of ints, averaging is performed on all of the axes
455 specified in the tuple instead of a single axis or all the axes as
456 before.
457 weights : array_like, optional
458 An array of weights associated with the values in `a`. Each value in
459 `a` contributes to the average according to its associated weight.
460 The array of weights must be the same shape as `a` if no axis is
461 specified, otherwise the weights must have dimensions and shape
462 consistent with `a` along the specified axis.
463 If `weights=None`, then all data in `a` are assumed to have a
464 weight equal to one.
465 The calculation is::
466
467 avg = sum(a * weights) / sum(weights)
468
469 where the sum is over all included elements.
470 The only constraint on the values of `weights` is that `sum(weights)`
471 must not be 0.
472 returned : bool, optional
473 Default is `False`. If `True`, the tuple (`average`, `sum_of_weights`)
474 is returned, otherwise only the average is returned.
475 If `weights=None`, `sum_of_weights` is equivalent to the number of
476 elements over which the average is taken.
477 keepdims : bool, optional
478 If this is set to True, the axes which are reduced are left
479 in the result as dimensions with size one. With this option,
480 the result will broadcast correctly against the original `a`.
481 *Note:* `keepdims` will not work with instances of `numpy.matrix`
482 or other classes whose methods do not support `keepdims`.
483
484 .. versionadded:: 1.23.0
485
486 Returns
487 -------
488 retval, [sum_of_weights] : array_type or double
489 Return the average along the specified axis. When `returned` is `True`,
490 return a tuple with the average as the first element and the sum
491 of the weights as the second element. `sum_of_weights` is of the
492 same type as `retval`. The result dtype follows a general pattern.
493 If `weights` is None, the result dtype will be that of `a` , or ``float64``
494 if `a` is integral. Otherwise, if `weights` is not None and `a` is non-
495 integral, the result type will be the type of lowest precision capable of
496 representing values of both `a` and `weights`. If `a` happens to be
497 integral, the previous rules still applies but the result dtype will
498 at least be ``float64``.
499
500 Raises
501 ------
502 ZeroDivisionError
503 When all weights along axis are zero. See `numpy.ma.average` for a
504 version robust to this type of error.
505 TypeError
506 When `weights` does not have the same shape as `a`, and `axis=None`.
507 ValueError
508 When `weights` does not have dimensions and shape consistent with `a`
509 along specified `axis`.
510
511 See Also
512 --------
513 mean
514
515 ma.average : average for masked arrays -- useful if your data contains
516 "missing" values
517 numpy.result_type : Returns the type that results from applying the
518 numpy type promotion rules to the arguments.
519
520 Examples
521 --------
522 >>> import numpy as np
523 >>> data = np.arange(1, 5)
524 >>> data
525 array([1, 2, 3, 4])
526 >>> np.average(data)
527 2.5
528 >>> np.average(np.arange(1, 11), weights=np.arange(10, 0, -1))
529 4.0
530
531 >>> data = np.arange(6).reshape((3, 2))
532 >>> data
533 array([[0, 1],
534 [2, 3],
535 [4, 5]])
536 >>> np.average(data, axis=1, weights=[1./4, 3./4])
537 array([0.75, 2.75, 4.75])
538 >>> np.average(data, weights=[1./4, 3./4])
539 Traceback (most recent call last):
540 ...
541 TypeError: Axis must be specified when shapes of a and weights differ.
542
543 With ``keepdims=True``, the following result has shape (3, 1).
544
545 >>> np.average(data, axis=1, keepdims=True)
546 array([[0.5],
547 [2.5],
548 [4.5]])
549
550 >>> data = np.arange(8).reshape((2, 2, 2))
551 >>> data
552 array([[[0, 1],
553 [2, 3]],
554 [[4, 5],
555 [6, 7]]])
556 >>> np.average(data, axis=(0, 1), weights=[[1./4, 3./4], [1., 1./2]])
557 array([3.4, 4.4])
558 >>> np.average(data, axis=0, weights=[[1./4, 3./4], [1., 1./2]])
559 Traceback (most recent call last):
560 ...
561 ValueError: Shape of weights must be consistent
562 with shape of a along specified axis.
563 """
564 a = np.asanyarray(a)
565
566 if axis is not None:
567 axis = _nx.normalize_axis_tuple(axis, a.ndim, argname="axis")
568
569 if keepdims is np._NoValue:
570 # Don't pass on the keepdims argument if one wasn't given.
571 keepdims_kw = {}
572 else:
573 keepdims_kw = {'keepdims': keepdims}
574
575 if weights is None:
576 avg = a.mean(axis, **keepdims_kw)
577 avg_as_array = np.asanyarray(avg)
578 scl = avg_as_array.dtype.type(a.size / avg_as_array.size)
579 else:
580 wgt = _weights_are_valid(weights=weights, a=a, axis=axis)
581
582 if issubclass(a.dtype.type, (np.integer, np.bool)):
583 result_dtype = np.result_type(a.dtype, wgt.dtype, 'f8')
584 else:
585 result_dtype = np.result_type(a.dtype, wgt.dtype)
586
587 scl = wgt.sum(axis=axis, dtype=result_dtype, **keepdims_kw)
588 if np.any(scl == 0.0):
589 raise ZeroDivisionError(
590 "Weights sum to zero, can't be normalized")
591
592 avg = avg_as_array = np.multiply(a, wgt,
593 dtype=result_dtype).sum(axis, **keepdims_kw) / scl
594
595 if returned:
596 if scl.shape != avg_as_array.shape:
597 scl = np.broadcast_to(scl, avg_as_array.shape, subok=True).copy()
598 return avg, scl
599 else:
600 return avg
601
602
603@set_module('numpy')
604def asarray_chkfinite(a, dtype=None, order=None):
605 """Convert the input to an array, checking for NaNs or Infs.
606
607 Parameters
608 ----------
609 a : array_like
610 Input data, in any form that can be converted to an array. This
611 includes lists, lists of tuples, tuples, tuples of tuples, tuples
612 of lists and ndarrays. Success requires no NaNs or Infs.
613 dtype : data-type, optional
614 By default, the data-type is inferred from the input data.
615 order : {'C', 'F', 'A', 'K'}, optional
616 The memory layout of the output.
617 'C' gives a row-major layout (C-style),
618 'F' gives a column-major layout (Fortran-style).
619 'C' and 'F' will copy if needed to ensure the output format.
620 'A' (any) is equivalent to 'F' if input a is non-contiguous or
621 Fortran-contiguous, otherwise, it is equivalent to 'C'.
622 Unlike 'C' or 'F', 'A' does not ensure that the result is contiguous.
623 'K' (keep) preserves the input order for the output.
624 'C' is the default.
625
626 Returns
627 -------
628 out : ndarray
629 Array interpretation of `a`. No copy is performed if the input
630 is already an ndarray. If `a` is a subclass of ndarray, a base
631 class ndarray is returned.
632
633 Raises
634 ------
635 ValueError
636 Raises ValueError if `a` contains NaN (Not a Number) or Inf (Infinity).
637
638 See Also
639 --------
640 asarray : Create and array.
641 asanyarray : Similar function which passes through subclasses.
642 ascontiguousarray : Convert input to a contiguous array.
643 asfortranarray : Convert input to an ndarray with column-major
644 memory order.
645 fromiter : Create an array from an iterator.
646 fromfunction : Construct an array by executing a function on grid
647 positions.
648
649 Examples
650 --------
651 >>> import numpy as np
652
653 Convert a list into an array. If all elements are finite, then
654 ``asarray_chkfinite`` is identical to ``asarray``.
655
656 >>> a = [1, 2]
657 >>> np.asarray_chkfinite(a, dtype=float)
658 array([1., 2.])
659
660 Raises ValueError if array_like contains Nans or Infs.
661
662 >>> a = [1, 2, np.inf]
663 >>> try:
664 ... np.asarray_chkfinite(a)
665 ... except ValueError:
666 ... print('ValueError')
667 ...
668 ValueError
669
670 """
671 a = asarray(a, dtype=dtype, order=order)
672 if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
673 raise ValueError(
674 "array must not contain infs or NaNs")
675 return a
676
677
678def _piecewise_dispatcher(x, condlist, funclist, *args, **kw):
679 yield x
680 # support the undocumented behavior of allowing scalars
681 if np.iterable(condlist):
682 yield from condlist
683
684
685@array_function_dispatch(_piecewise_dispatcher)
686def piecewise(x, condlist, funclist, *args, **kw):
687 """
688 Evaluate a piecewise-defined function.
689
690 Given a set of conditions and corresponding functions, evaluate each
691 function on the input data wherever its condition is true.
692
693 Parameters
694 ----------
695 x : ndarray or scalar
696 The input domain.
697 condlist : list of bool arrays or bool scalars
698 Each boolean array corresponds to a function in `funclist`. Wherever
699 `condlist[i]` is True, `funclist[i](x)` is used as the output value.
700
701 Each boolean array in `condlist` selects a piece of `x`,
702 and should therefore be of the same shape as `x`.
703
704 The length of `condlist` must correspond to that of `funclist`.
705 If one extra function is given, i.e. if
706 ``len(funclist) == len(condlist) + 1``, then that extra function
707 is the default value, used wherever all conditions are false.
708 funclist : list of callables, f(x,*args,**kw), or scalars
709 Each function is evaluated over `x` wherever its corresponding
710 condition is True. It should take a 1d array as input and give a 1d
711 array or a scalar value as output. If, instead of a callable,
712 a scalar is provided then a constant function (``lambda x: scalar``) is
713 assumed.
714 args : tuple, optional
715 Any further arguments given to `piecewise` are passed to the functions
716 upon execution, i.e., if called ``piecewise(..., ..., 1, 'a')``, then
717 each function is called as ``f(x, 1, 'a')``.
718 kw : dict, optional
719 Keyword arguments used in calling `piecewise` are passed to the
720 functions upon execution, i.e., if called
721 ``piecewise(..., ..., alpha=1)``, then each function is called as
722 ``f(x, alpha=1)``.
723
724 Returns
725 -------
726 out : ndarray
727 The output is the same shape and type as x and is found by
728 calling the functions in `funclist` on the appropriate portions of `x`,
729 as defined by the boolean arrays in `condlist`. Portions not covered
730 by any condition have a default value of 0.
731
732
733 See Also
734 --------
735 choose, select, where
736
737 Notes
738 -----
739 This is similar to choose or select, except that functions are
740 evaluated on elements of `x` that satisfy the corresponding condition from
741 `condlist`.
742
743 The result is::
744
745 |--
746 |funclist[0](x[condlist[0]])
747 out = |funclist[1](x[condlist[1]])
748 |...
749 |funclist[n2](x[condlist[n2]])
750 |--
751
752 Examples
753 --------
754 >>> import numpy as np
755
756 Define the signum function, which is -1 for ``x < 0`` and +1 for ``x >= 0``.
757
758 >>> x = np.linspace(-2.5, 2.5, 6)
759 >>> np.piecewise(x, [x < 0, x >= 0], [-1, 1])
760 array([-1., -1., -1., 1., 1., 1.])
761
762 Define the absolute value, which is ``-x`` for ``x <0`` and ``x`` for
763 ``x >= 0``.
764
765 >>> np.piecewise(x, [x < 0, x >= 0], [lambda x: -x, lambda x: x])
766 array([2.5, 1.5, 0.5, 0.5, 1.5, 2.5])
767
768 Apply the same function to a scalar value.
769
770 >>> y = -2
771 >>> np.piecewise(y, [y < 0, y >= 0], [lambda x: -x, lambda x: x])
772 array(2)
773
774 """
775 x = asanyarray(x)
776 n2 = len(funclist)
777
778 # undocumented: single condition is promoted to a list of one condition
779 if isscalar(condlist) or (
780 not isinstance(condlist[0], (list, ndarray)) and x.ndim != 0):
781 condlist = [condlist]
782
783 condlist = asarray(condlist, dtype=bool)
784 n = len(condlist)
785
786 if n == n2 - 1: # compute the "otherwise" condition.
787 condelse = ~np.any(condlist, axis=0, keepdims=True)
788 condlist = np.concatenate([condlist, condelse], axis=0)
789 n += 1
790 elif n != n2:
791 raise ValueError(
792 f"with {n} condition(s), either {n} or {n + 1} functions are expected"
793 )
794
795 y = zeros_like(x)
796 for cond, func in zip(condlist, funclist):
797 if not isinstance(func, collections.abc.Callable):
798 y[cond] = func
799 else:
800 vals = x[cond]
801 if vals.size > 0:
802 y[cond] = func(vals, *args, **kw)
803
804 return y
805
806
807def _select_dispatcher(condlist, choicelist, default=None):
808 yield from condlist
809 yield from choicelist
810
811
812@array_function_dispatch(_select_dispatcher)
813def select(condlist, choicelist, default=0):
814 """
815 Return an array drawn from elements in choicelist, depending on conditions.
816
817 Parameters
818 ----------
819 condlist : list of bool ndarrays
820 The list of conditions which determine from which array in `choicelist`
821 the output elements are taken. When multiple conditions are satisfied,
822 the first one encountered in `condlist` is used.
823 choicelist : list of ndarrays
824 The list of arrays from which the output elements are taken. It has
825 to be of the same length as `condlist`.
826 default : array_like, optional
827 The element inserted in `output` when all conditions evaluate to False.
828
829 Returns
830 -------
831 output : ndarray
832 The output at position m is the m-th element of the array in
833 `choicelist` where the m-th element of the corresponding array in
834 `condlist` is True.
835
836 See Also
837 --------
838 where : Return elements from one of two arrays depending on condition.
839 take, choose, compress, diag, diagonal
840
841 Examples
842 --------
843 >>> import numpy as np
844
845 Beginning with an array of integers from 0 to 5 (inclusive),
846 elements less than ``3`` are negated, elements greater than ``3``
847 are squared, and elements not meeting either of these conditions
848 (exactly ``3``) are replaced with a `default` value of ``42``.
849
850 >>> x = np.arange(6)
851 >>> condlist = [x<3, x>3]
852 >>> choicelist = [-x, x**2]
853 >>> np.select(condlist, choicelist, 42)
854 array([ 0, -1, -2, 42, 16, 25])
855
856 When multiple conditions are satisfied, the first one encountered in
857 `condlist` is used.
858
859 >>> condlist = [x<=4, x>3]
860 >>> choicelist = [x, x**2]
861 >>> np.select(condlist, choicelist, 55)
862 array([ 0, 1, 2, 3, 4, 25])
863
864 """
865 # Check the size of condlist and choicelist are the same, or abort.
866 if len(condlist) != len(choicelist):
867 raise ValueError(
868 'list of cases must be same length as list of conditions')
869
870 # Now that the dtype is known, handle the deprecated select([], []) case
871 if len(condlist) == 0:
872 raise ValueError("select with an empty condition list is not possible")
873
874 # TODO: This preserves the Python int, float, complex manually to get the
875 # right `result_type` with NEP 50. Most likely we will grow a better
876 # way to spell this (and this can be replaced).
877 choicelist = [
878 choice if type(choice) in (int, float, complex) else np.asarray(choice)
879 for choice in choicelist]
880 choicelist.append(default if type(default) in (int, float, complex)
881 else np.asarray(default))
882
883 try:
884 dtype = np.result_type(*choicelist)
885 except TypeError as e:
886 msg = f'Choicelist and default value do not have a common dtype: {e}'
887 raise TypeError(msg) from None
888
889 # Convert conditions to arrays and broadcast conditions and choices
890 # as the shape is needed for the result. Doing it separately optimizes
891 # for example when all choices are scalars.
892 condlist = np.broadcast_arrays(*condlist)
893 choicelist = np.broadcast_arrays(*choicelist)
894
895 # If cond array is not an ndarray in boolean format or scalar bool, abort.
896 for i, cond in enumerate(condlist):
897 if cond.dtype.type is not np.bool:
898 raise TypeError(
899 f'invalid entry {i} in condlist: should be boolean ndarray')
900
901 if choicelist[0].ndim == 0:
902 # This may be common, so avoid the call.
903 result_shape = condlist[0].shape
904 else:
905 result_shape = np.broadcast_arrays(condlist[0], choicelist[0])[0].shape
906
907 result = np.full(result_shape, choicelist[-1], dtype)
908
909 # Use np.copyto to burn each choicelist array onto result, using the
910 # corresponding condlist as a boolean mask. This is done in reverse
911 # order since the first choice should take precedence.
912 choicelist = choicelist[-2::-1]
913 condlist = condlist[::-1]
914 for choice, cond in zip(choicelist, condlist):
915 np.copyto(result, choice, where=cond)
916
917 return result
918
919
920def _copy_dispatcher(a, order=None, subok=None):
921 return (a,)
922
923
924@array_function_dispatch(_copy_dispatcher)
925def copy(a, order='K', subok=False):
926 """
927 Return an array copy of the given object.
928
929 Parameters
930 ----------
931 a : array_like
932 Input data.
933 order : {'C', 'F', 'A', 'K'}, optional
934 Controls the memory layout of the copy. 'C' means C-order,
935 'F' means F-order, 'A' means 'F' if `a` is Fortran contiguous,
936 'C' otherwise. 'K' means match the layout of `a` as closely
937 as possible. (Note that this function and :meth:`ndarray.copy` are very
938 similar, but have different default values for their order=
939 arguments.)
940 subok : bool, optional
941 If True, then sub-classes will be passed-through, otherwise the
942 returned array will be forced to be a base-class array (defaults to False).
943
944 Returns
945 -------
946 arr : ndarray
947 Array interpretation of `a`.
948
949 See Also
950 --------
951 ndarray.copy : Preferred method for creating an array copy
952
953 Notes
954 -----
955 This is equivalent to:
956
957 >>> np.array(a, copy=True) #doctest: +SKIP
958
959 The copy made of the data is shallow, i.e., for arrays with object dtype,
960 the new array will point to the same objects.
961 See Examples from `ndarray.copy`.
962
963 Examples
964 --------
965 >>> import numpy as np
966
967 Create an array x, with a reference y and a copy z:
968
969 >>> x = np.array([1, 2, 3])
970 >>> y = x
971 >>> z = np.copy(x)
972
973 Note that, when we modify x, y changes, but not z:
974
975 >>> x[0] = 10
976 >>> x[0] == y[0]
977 True
978 >>> x[0] == z[0]
979 False
980
981 Note that, np.copy clears previously set WRITEABLE=False flag.
982
983 >>> a = np.array([1, 2, 3])
984 >>> a.flags["WRITEABLE"] = False
985 >>> b = np.copy(a)
986 >>> b.flags["WRITEABLE"]
987 True
988 >>> b[0] = 3
989 >>> b
990 array([3, 2, 3])
991 """
992 return array(a, order=order, subok=subok, copy=True)
993
994# Basic operations
995
996
997def _gradient_dispatcher(f, *varargs, axis=None, edge_order=None):
998 yield f
999 yield from varargs
1000
1001
1002@array_function_dispatch(_gradient_dispatcher)
1003def gradient(f, *varargs, axis=None, edge_order=1):
1004 """
1005 Return the gradient of an N-dimensional array.
1006
1007 The gradient is computed using second order accurate central differences
1008 in the interior points and either first or second order accurate one-sides
1009 (forward or backwards) differences at the boundaries.
1010 The returned gradient hence has the same shape as the input array.
1011
1012 Parameters
1013 ----------
1014 f : array_like
1015 An N-dimensional array containing samples of a scalar function.
1016 varargs : list of scalar or array, optional
1017 Spacing between f values. Default unitary spacing for all dimensions.
1018 Spacing can be specified using:
1019
1020 1. single scalar to specify a sample distance for all dimensions.
1021 2. N scalars to specify a constant sample distance for each dimension.
1022 i.e. `dx`, `dy`, `dz`, ...
1023 3. N arrays to specify the coordinates of the values along each
1024 dimension of F. The length of the array must match the size of
1025 the corresponding dimension
1026 4. Any combination of N scalars/arrays with the meaning of 2. and 3.
1027
1028 If `axis` is given, the number of varargs must equal the number of axes
1029 specified in the axis parameter.
1030 Default: 1. (see Examples below).
1031
1032 edge_order : {1, 2}, optional
1033 Gradient is calculated using N-th order accurate differences
1034 at the boundaries. Default: 1.
1035 axis : None or int or tuple of ints, optional
1036 Gradient is calculated only along the given axis or axes
1037 The default (axis = None) is to calculate the gradient for all the axes
1038 of the input array. axis may be negative, in which case it counts from
1039 the last to the first axis.
1040
1041 Returns
1042 -------
1043 gradient : ndarray or tuple of ndarray
1044 A tuple of ndarrays (or a single ndarray if there is only one
1045 dimension) corresponding to the derivatives of f with respect
1046 to each dimension. Each derivative has the same shape as f.
1047
1048 Examples
1049 --------
1050 >>> import numpy as np
1051 >>> f = np.array([1, 2, 4, 7, 11, 16])
1052 >>> np.gradient(f)
1053 array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
1054 >>> np.gradient(f, 2)
1055 array([0.5 , 0.75, 1.25, 1.75, 2.25, 2.5 ])
1056
1057 Spacing can be also specified with an array that represents the coordinates
1058 of the values F along the dimensions.
1059 For instance a uniform spacing:
1060
1061 >>> x = np.arange(f.size)
1062 >>> np.gradient(f, x)
1063 array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
1064
1065 Or a non uniform one:
1066
1067 >>> x = np.array([0., 1., 1.5, 3.5, 4., 6.])
1068 >>> np.gradient(f, x)
1069 array([1. , 3. , 3.5, 6.7, 6.9, 2.5])
1070
1071 For two dimensional arrays, the return will be two arrays ordered by
1072 axis. In this example the first array stands for the gradient in
1073 rows and the second one in columns direction:
1074
1075 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]]))
1076 (array([[ 2., 2., -1.],
1077 [ 2., 2., -1.]]),
1078 array([[1. , 2.5, 4. ],
1079 [1. , 1. , 1. ]]))
1080
1081 In this example the spacing is also specified:
1082 uniform for axis=0 and non uniform for axis=1
1083
1084 >>> dx = 2.
1085 >>> y = [1., 1.5, 3.5]
1086 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]]), dx, y)
1087 (array([[ 1. , 1. , -0.5],
1088 [ 1. , 1. , -0.5]]),
1089 array([[2. , 2. , 2. ],
1090 [2. , 1.7, 0.5]]))
1091
1092 It is possible to specify how boundaries are treated using `edge_order`
1093
1094 >>> x = np.array([0, 1, 2, 3, 4])
1095 >>> f = x**2
1096 >>> np.gradient(f, edge_order=1)
1097 array([1., 2., 4., 6., 7.])
1098 >>> np.gradient(f, edge_order=2)
1099 array([0., 2., 4., 6., 8.])
1100
1101 The `axis` keyword can be used to specify a subset of axes of which the
1102 gradient is calculated
1103
1104 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]]), axis=0)
1105 array([[ 2., 2., -1.],
1106 [ 2., 2., -1.]])
1107
1108 The `varargs` argument defines the spacing between sample points in the
1109 input array. It can take two forms:
1110
1111 1. An array, specifying coordinates, which may be unevenly spaced:
1112
1113 >>> x = np.array([0., 2., 3., 6., 8.])
1114 >>> y = x ** 2
1115 >>> np.gradient(y, x, edge_order=2)
1116 array([ 0., 4., 6., 12., 16.])
1117
1118 2. A scalar, representing the fixed sample distance:
1119
1120 >>> dx = 2
1121 >>> x = np.array([0., 2., 4., 6., 8.])
1122 >>> y = x ** 2
1123 >>> np.gradient(y, dx, edge_order=2)
1124 array([ 0., 4., 8., 12., 16.])
1125
1126 It's possible to provide different data for spacing along each dimension.
1127 The number of arguments must match the number of dimensions in the input
1128 data.
1129
1130 >>> dx = 2
1131 >>> dy = 3
1132 >>> x = np.arange(0, 6, dx)
1133 >>> y = np.arange(0, 9, dy)
1134 >>> xs, ys = np.meshgrid(x, y)
1135 >>> zs = xs + 2 * ys
1136 >>> np.gradient(zs, dy, dx) # Passing two scalars
1137 (array([[2., 2., 2.],
1138 [2., 2., 2.],
1139 [2., 2., 2.]]),
1140 array([[1., 1., 1.],
1141 [1., 1., 1.],
1142 [1., 1., 1.]]))
1143
1144 Mixing scalars and arrays is also allowed:
1145
1146 >>> np.gradient(zs, y, dx) # Passing one array and one scalar
1147 (array([[2., 2., 2.],
1148 [2., 2., 2.],
1149 [2., 2., 2.]]),
1150 array([[1., 1., 1.],
1151 [1., 1., 1.],
1152 [1., 1., 1.]]))
1153
1154 Notes
1155 -----
1156 Assuming that :math:`f\\in C^{3}` (i.e., :math:`f` has at least 3 continuous
1157 derivatives) and let :math:`h_{*}` be a non-homogeneous stepsize, we
1158 minimize the "consistency error" :math:`\\eta_{i}` between the true gradient
1159 and its estimate from a linear combination of the neighboring grid-points:
1160
1161 .. math::
1162
1163 \\eta_{i} = f_{i}^{\\left(1\\right)} -
1164 \\left[ \\alpha f\\left(x_{i}\\right) +
1165 \\beta f\\left(x_{i} + h_{d}\\right) +
1166 \\gamma f\\left(x_{i}-h_{s}\\right)
1167 \\right]
1168
1169 By substituting :math:`f(x_{i} + h_{d})` and :math:`f(x_{i} - h_{s})`
1170 with their Taylor series expansion, this translates into solving
1171 the following the linear system:
1172
1173 .. math::
1174
1175 \\left\\{
1176 \\begin{array}{r}
1177 \\alpha+\\beta+\\gamma=0 \\\\
1178 \\beta h_{d}-\\gamma h_{s}=1 \\\\
1179 \\beta h_{d}^{2}+\\gamma h_{s}^{2}=0
1180 \\end{array}
1181 \\right.
1182
1183 The resulting approximation of :math:`f_{i}^{(1)}` is the following:
1184
1185 .. math::
1186
1187 \\hat f_{i}^{(1)} =
1188 \\frac{
1189 h_{s}^{2}f\\left(x_{i} + h_{d}\\right)
1190 + \\left(h_{d}^{2} - h_{s}^{2}\\right)f\\left(x_{i}\\right)
1191 - h_{d}^{2}f\\left(x_{i}-h_{s}\\right)}
1192 { h_{s}h_{d}\\left(h_{d} + h_{s}\\right)}
1193 + \\mathcal{O}\\left(\\frac{h_{d}h_{s}^{2}
1194 + h_{s}h_{d}^{2}}{h_{d}
1195 + h_{s}}\\right)
1196
1197 It is worth noting that if :math:`h_{s}=h_{d}`
1198 (i.e., data are evenly spaced)
1199 we find the standard second order approximation:
1200