1import collections.abc 2import functools 3import re 4import sys 5import warnings 6 7import numpy as np 8import numpy.core.numeric as _nx 9from numpy.core import transpose 10from numpy.core.numeric import ( 11 ones, zeros, arange, concatenate, array, asarray, asanyarray, empty, 12 ndarray, around, floor, ceil, take, dot, where, intp, 13 integer, isscalar, absolute 14 ) 15from numpy.core.umath import ( 16 pi, add, arctan2, frompyfunc, cos, less_equal, sqrt, sin, 17 mod, exp, not_equal, subtract 18 ) 19from numpy.core.fromnumeric import ( 20 ravel, nonzero, partition, mean, any, sum 21 ) 22from numpy.core.numerictypes import typecodes 23from numpy.core.overrides import set_module 24from numpy.core import overrides 25from numpy.core.function_base import add_newdoc 26from numpy.lib.twodim_base import diag 27from numpy.core.multiarray import ( 28 _insert, add_docstring, bincount, normalize_axis_index, _monotonicity, 29 interp as compiled_interp, interp_complex as compiled_interp_complex 30 ) 31from numpy.core.umath import _add_newdoc_ufunc as add_newdoc_ufunc 32 33import builtins 34 35# needed in this module for compatibility 36from numpy.lib.histograms import histogram, histogramdd 37 38 39array_function_dispatch = functools.partial( 40 overrides.array_function_dispatch, module='numpy') 41 42 43__all__ = [ 44 'select', 'piecewise', 'trim_zeros', 'copy', 'iterable', 'percentile', 45 'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'disp', 'flip', 46 'rot90', 'extract', 'place', 'vectorize', 'asarray_chkfinite', 'average', 47 'bincount', 'digitize', 'cov', 'corrcoef', 48 'msort', 'median', 'sinc', 'hamming', 'hanning', 'bartlett', 49 'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc', 'add_docstring', 50 'meshgrid', 'delete', 'insert', 'append', 'interp', 'add_newdoc_ufunc', 51 'quantile' 52 ] 53 54 55def _rot90_dispatcher(m, k=None, axes=None): 56 return (m,) 57 58 59@array_function_dispatch(_rot90_dispatcher) 60def rot90(m, k=1, axes=(0, 1)): 61 """ 62 Rotate an array by 90 degrees in the plane specified by axes. 63 64 Rotation direction is from the first towards the second axis. 65 66 Parameters 67 ---------- 68 m : array_like 69 Array of two or more dimensions. 70 k : integer 71 Number of times the array is rotated by 90 degrees. 72 axes: (2,) array_like 73 The array is rotated in the plane defined by the axes. 74 Axes must be different. 75 76 .. versionadded:: 1.12.0 77 78 Returns 79 ------- 80 y : ndarray 81 A rotated view of `m`. 82 83 See Also 84 -------- 85 flip : Reverse the order of elements in an array along the given axis. 86 fliplr : Flip an array horizontally. 87 flipud : Flip an array vertically. 88 89 Notes 90 ----- 91 rot90(m, k=1, axes=(1,0)) is the reverse of rot90(m, k=1, axes=(0,1)) 92 rot90(m, k=1, axes=(1,0)) is equivalent to rot90(m, k=-1, axes=(0,1)) 93 94 Examples 95 -------- 96 >>> m = np.array([[1,2],[3,4]], int) 97 >>> m 98 array([[1, 2], 99 [3, 4]]) 100 >>> np.rot90(m) 101 array([[2, 4], 102 [1, 3]]) 103 >>> np.rot90(m, 2) 104 array([[4, 3], 105 [2, 1]]) 106 >>> m = np.arange(8).reshape((2,2,2)) 107 >>> np.rot90(m, 1, (1,2)) 108 array([[[1, 3], 109 [0, 2]], 110 [[5, 7], 111 [4, 6]]]) 112 113 """ 114 axes = tuple(axes) 115 if len(axes) != 2: 116 raise ValueError("len(axes) must be 2.") 117 118 m = asanyarray(m) 119 120 if axes[0] == axes[1] or absolute(axes[0] - axes[1]) == m.ndim: 121 raise ValueError("Axes must be different.") 122 123 if (axes[0] >= m.ndim or axes[0] < -m.ndim 124 or axes[1] >= m.ndim or axes[1] < -m.ndim): 125 raise ValueError("Axes={} out of range for array of ndim={}." 126 .format(axes, m.ndim)) 127 128 k %= 4 129 130 if k == 0: 131 return m[:] 132 if k == 2: 133 return flip(flip(m, axes[0]), axes[1]) 134 135 axes_list = arange(0, m.ndim) 136 (axes_list[axes[0]], axes_list[axes[1]]) = (axes_list[axes[1]], 137 axes_list[axes[0]]) 138 139 if k == 1: 140 return transpose(flip(m, axes[1]), axes_list) 141 else: 142 # k == 3 143 return flip(transpose(m, axes_list), axes[1]) 144 145 146def _flip_dispatcher(m, axis=None): 147 return (m,) 148 149 150@array_function_dispatch(_flip_dispatcher) 151def flip(m, axis=None): 152 """ 153 Reverse the order of elements in an array along the given axis. 154 155 The shape of the array is preserved, but the elements are reordered. 156 157 .. versionadded:: 1.12.0 158 159 Parameters 160 ---------- 161 m : array_like 162 Input array. 163 axis : None or int or tuple of ints, optional 164 Axis or axes along which to flip over. The default, 165 axis=None, will flip over all of the axes of the input array. 166 If axis is negative it counts from the last to the first axis. 167 168 If axis is a tuple of ints, flipping is performed on all of the axes 169 specified in the tuple. 170 171 .. versionchanged:: 1.15.0 172 None and tuples of axes are supported 173 174 Returns 175 ------- 176 out : array_like 177 A view of `m` with the entries of axis reversed. Since a view is 178 returned, this operation is done in constant time. 179 180 See Also 181 -------- 182 flipud : Flip an array vertically (axis=0). 183 fliplr : Flip an array horizontally (axis=1). 184 185 Notes 186 ----- 187 flip(m, 0) is equivalent to flipud(m). 188 189 flip(m, 1) is equivalent to fliplr(m). 190 191 flip(m, n) corresponds to ``m[...,::-1,...]`` with ``::-1`` at position n. 192 193 flip(m) corresponds to ``m[::-1,::-1,...,::-1]`` with ``::-1`` at all 194 positions. 195 196 flip(m, (0, 1)) corresponds to ``m[::-1,::-1,...]`` with ``::-1`` at 197 position 0 and position 1. 198 199 Examples 200 -------- 201 >>> A = np.arange(8).reshape((2,2,2)) 202 >>> A 203 array([[[0, 1], 204 [2, 3]], 205 [[4, 5], 206 [6, 7]]]) 207 >>> np.flip(A, 0) 208 array([[[4, 5], 209 [6, 7]], 210 [[0, 1], 211 [2, 3]]]) 212 >>> np.flip(A, 1) 213 array([[[2, 3], 214 [0, 1]], 215 [[6, 7], 216 [4, 5]]]) 217 >>> np.flip(A) 218 array([[[7, 6], 219 [5, 4]], 220 [[3, 2], 221 [1, 0]]]) 222 >>> np.flip(A, (0, 2)) 223 array([[[5, 4], 224 [7, 6]], 225 [[1, 0], 226 [3, 2]]]) 227 >>> A = np.random.randn(3,4,5) 228 >>> np.all(np.flip(A,2) == A[:,:,::-1,...]) 229 True 230 """ 231 if not hasattr(m, 'ndim'): 232 m = asarray(m) 233 if axis is None: 234 indexer = (np.s_[::-1],) * m.ndim 235 else: 236 axis = _nx.normalize_axis_tuple(axis, m.ndim) 237 indexer = [np.s_[:]] * m.ndim 238 for ax in axis: 239 indexer[ax] = np.s_[::-1] 240 indexer = tuple(indexer) 241 return m[indexer] 242 243 244@set_module('numpy') 245def iterable(y): 246 """ 247 Check whether or not an object can be iterated over. 248 249 Parameters 250 ---------- 251 y : object 252 Input object. 253 254 Returns 255 ------- 256 b : bool 257 Return ``True`` if the object has an iterator method or is a 258 sequence and ``False`` otherwise. 259 260 261 Examples 262 -------- 263 >>> np.iterable([1, 2, 3]) 264 True 265 >>> np.iterable(2) 266 False 267 268 """ 269 try: 270 iter(y) 271 except TypeError: 272 return False 273 return True 274 275 276def _average_dispatcher(a, axis=None, weights=None, returned=None): 277 return (a, weights) 278 279 280@array_function_dispatch(_average_dispatcher) 281def average(a, axis=None, weights=None, returned=False): 282 """ 283 Compute the weighted average along the specified axis. 284 285 Parameters 286 ---------- 287 a : array_like 288 Array containing data to be averaged. If `a` is not an array, a 289 conversion is attempted. 290 axis : None or int or tuple of ints, optional 291 Axis or axes along which to average `a`. The default, 292 axis=None, will average over all of the elements of the input array. 293 If axis is negative it counts from the last to the first axis. 294 295 .. versionadded:: 1.7.0 296 297 If axis is a tuple of ints, averaging is performed on all of the axes 298 specified in the tuple instead of a single axis or all the axes as 299 before. 300 weights : array_like, optional 301 An array of weights associated with the values in `a`. Each value in 302 `a` contributes to the average according to its associated weight. 303 The weights array can either be 1-D (in which case its length must be 304 the size of `a` along the given axis) or of the same shape as `a`. 305 If `weights=None`, then all data in `a` are assumed to have a 306 weight equal to one. The 1-D calculation is:: 307 308 avg = sum(a * weights) / sum(weights) 309 310 The only constraint on `weights` is that `sum(weights)` must not be 0. 311 returned : bool, optional 312 Default is `False`. If `True`, the tuple (`average`, `sum_of_weights`) 313 is returned, otherwise only the average is returned. 314 If `weights=None`, `sum_of_weights` is equivalent to the number of 315 elements over which the average is taken. 316 317 Returns 318 ------- 319 retval, [sum_of_weights] : array_type or double 320 Return the average along the specified axis. When `returned` is `True`, 321 return a tuple with the average as the first element and the sum 322 of the weights as the second element. `sum_of_weights` is of the 323 same type as `retval`. The result dtype follows a genereal pattern. 324 If `weights` is None, the result dtype will be that of `a` , or ``float64`` 325 if `a` is integral. Otherwise, if `weights` is not None and `a` is non- 326 integral, the result type will be the type of lowest precision capable of 327 representing values of both `a` and `weights`. If `a` happens to be 328 integral, the previous rules still applies but the result dtype will 329 at least be ``float64``. 330 331 Raises 332 ------ 333 ZeroDivisionError 334 When all weights along axis are zero. See `numpy.ma.average` for a 335 version robust to this type of error. 336 TypeError 337 When the length of 1D `weights` is not the same as the shape of `a` 338 along axis. 339 340 See Also 341 -------- 342 mean 343 344 ma.average : average for masked arrays -- useful if your data contains 345 "missing" values 346 numpy.result_type : Returns the type that results from applying the 347 numpy type promotion rules to the arguments. 348 349 Examples 350 -------- 351 >>> data = np.arange(1, 5) 352 >>> data 353 array([1, 2, 3, 4]) 354 >>> np.average(data) 355 2.5 356 >>> np.average(np.arange(1, 11), weights=np.arange(10, 0, -1)) 357 4.0 358 359 >>> data = np.arange(6).reshape((3,2)) 360 >>> data 361 array([[0, 1], 362 [2, 3], 363 [4, 5]]) 364 >>> np.average(data, axis=1, weights=[1./4, 3./4]) 365 array([0.75, 2.75, 4.75]) 366 >>> np.average(data, weights=[1./4, 3./4]) 367 Traceback (most recent call last): 368 ... 369 TypeError: Axis must be specified when shapes of a and weights differ. 370 371 >>> a = np.ones(5, dtype=np.float128) 372 >>> w = np.ones(5, dtype=np.complex64) 373 >>> avg = np.average(a, weights=w) 374 >>> print(avg.dtype) 375 complex256 376 """ 377 a = np.asanyarray(a) 378 379 if weights is None: 380 avg = a.mean(axis) 381 scl = avg.dtype.type(a.size/avg.size) 382 else: 383 wgt = np.asanyarray(weights) 384 385 if issubclass(a.dtype.type, (np.integer, np.bool_)): 386 result_dtype = np.result_type(a.dtype, wgt.dtype, 'f8') 387 else: 388 result_dtype = np.result_type(a.dtype, wgt.dtype) 389 390 # Sanity checks 391 if a.shape != wgt.shape: 392 if axis is None: 393 raise TypeError( 394 "Axis must be specified when shapes of a and weights " 395 "differ.") 396 if wgt.ndim != 1: 397 raise TypeError( 398 "1D weights expected when shapes of a and weights differ.") 399 if wgt.shape[0] != a.shape[axis]: 400 raise ValueError( 401 "Length of weights not compatible with specified axis.") 402 403 # setup wgt to broadcast along axis 404 wgt = np.broadcast_to(wgt, (a.ndim-1)*(1,) + wgt.shape) 405 wgt = wgt.swapaxes(-1, axis) 406 407 scl = wgt.sum(axis=axis, dtype=result_dtype) 408 if np.any(scl == 0.0): 409 raise ZeroDivisionError( 410 "Weights sum to zero, can't be normalized") 411 412 avg = np.multiply(a, wgt, dtype=result_dtype).sum(axis)/scl 413 414 if returned: 415 if scl.shape != avg.shape: 416 scl = np.broadcast_to(scl, avg.shape).copy() 417 return avg, scl 418 else: 419 return avg 420 421 422@set_module('numpy') 423def asarray_chkfinite(a, dtype=None, order=None): 424 """Convert the input to an array, checking for NaNs or Infs. 425 426 Parameters 427 ---------- 428 a : array_like 429 Input data, in any form that can be converted to an array. This 430 includes lists, lists of tuples, tuples, tuples of tuples, tuples 431 of lists and ndarrays. Success requires no NaNs or Infs. 432 dtype : data-type, optional 433 By default, the data-type is inferred from the input data. 434 order : {'C', 'F', 'A', 'K'}, optional 435 Memory layout. 'A' and 'K' depend on the order of input array a. 436 'C' row-major (C-style), 437 'F' column-major (Fortran-style) memory representation. 438 'A' (any) means 'F' if `a` is Fortran contiguous, 'C' otherwise 439 'K' (keep) preserve input order 440 Defaults to 'C'. 441 442 Returns 443 ------- 444 out : ndarray 445 Array interpretation of `a`. No copy is performed if the input 446 is already an ndarray. If `a` is a subclass of ndarray, a base 447 class ndarray is returned. 448 449 Raises 450 ------ 451 ValueError 452 Raises ValueError if `a` contains NaN (Not a Number) or Inf (Infinity). 453 454 See Also 455 -------- 456 asarray : Create and array. 457 asanyarray : Similar function which passes through subclasses. 458 ascontiguousarray : Convert input to a contiguous array. 459 asfarray : Convert input to a floating point ndarray. 460 asfortranarray : Convert input to an ndarray with column-major 461 memory order. 462 fromiter : Create an array from an iterator. 463 fromfunction : Construct an array by executing a function on grid 464 positions. 465 466 Examples 467 -------- 468 Convert a list into an array. If all elements are finite 469 ``asarray_chkfinite`` is identical to ``asarray``. 470 471 >>> a = [1, 2] 472 >>> np.asarray_chkfinite(a, dtype=float) 473 array([1., 2.]) 474 475 Raises ValueError if array_like contains Nans or Infs. 476 477 >>> a = [1, 2, np.inf] 478 >>> try: 479 ... np.asarray_chkfinite(a) 480 ... except ValueError: 481 ... print('ValueError') 482 ... 483 ValueError 484 485 """ 486 a = asarray(a, dtype=dtype, order=order) 487 if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all(): 488 raise ValueError( 489 "array must not contain infs or NaNs") 490 return a 491 492 493def _piecewise_dispatcher(x, condlist, funclist, *args, **kw): 494 yield x 495 # support the undocumented behavior of allowing scalars 496 if np.iterable(condlist): 497 yield from condlist 498 499 500@array_function_dispatch(_piecewise_dispatcher) 501def piecewise(x, condlist, funclist, *args, **kw): 502 """ 503 Evaluate a piecewise-defined function. 504 505 Given a set of conditions and corresponding functions, evaluate each 506 function on the input data wherever its condition is true. 507 508 Parameters 509 ---------- 510 x : ndarray or scalar 511 The input domain. 512 condlist : list of bool arrays or bool scalars 513 Each boolean array corresponds to a function in `funclist`. Wherever 514 `condlist[i]` is True, `funclist[i](x)` is used as the output value. 515 516 Each boolean array in `condlist` selects a piece of `x`, 517 and should therefore be of the same shape as `x`. 518 519 The length of `condlist` must correspond to that of `funclist`. 520 If one extra function is given, i.e. if 521 ``len(funclist) == len(condlist) + 1``, then that extra function 522 is the default value, used wherever all conditions are false. 523 funclist : list of callables, f(x,*args,**kw), or scalars 524 Each function is evaluated over `x` wherever its corresponding 525 condition is True. It should take a 1d array as input and give an 1d 526 array or a scalar value as output. If, instead of a callable, 527 a scalar is provided then a constant function (``lambda x: scalar``) is 528 assumed. 529 args : tuple, optional 530 Any further arguments given to `piecewise` are passed to the functions 531 upon execution, i.e., if called ``piecewise(..., ..., 1, 'a')``, then 532 each function is called as ``f(x, 1, 'a')``. 533 kw : dict, optional 534 Keyword arguments used in calling `piecewise` are passed to the 535 functions upon execution, i.e., if called 536 ``piecewise(..., ..., alpha=1)``, then each function is called as 537 ``f(x, alpha=1)``. 538 539 Returns 540 ------- 541 out : ndarray 542 The output is the same shape and type as x and is found by 543 calling the functions in `funclist` on the appropriate portions of `x`, 544 as defined by the boolean arrays in `condlist`. Portions not covered 545 by any condition have a default value of 0. 546 547 548 See Also 549 -------- 550 choose, select, where 551 552 Notes 553 ----- 554 This is similar to choose or select, except that functions are 555 evaluated on elements of `x` that satisfy the corresponding condition from 556 `condlist`. 557 558 The result is:: 559 560 |-- 561 |funclist[0](x[condlist[0]]) 562 out = |funclist[1](x[condlist[1]]) 563 |... 564 |funclist[n2](x[condlist[n2]]) 565 |-- 566 567 Examples 568 -------- 569 Define the sigma function, which is -1 for ``x < 0`` and +1 for ``x >= 0``. 570 571 >>> x = np.linspace(-2.5, 2.5, 6) 572 >>> np.piecewise(x, [x < 0, x >= 0], [-1, 1]) 573 array([-1., -1., -1., 1., 1., 1.]) 574 575 Define the absolute value, which is ``-x`` for ``x <0`` and ``x`` for 576 ``x >= 0``. 577 578 >>> np.piecewise(x, [x < 0, x >= 0], [lambda x: -x, lambda x: x]) 579 array([2.5, 1.5, 0.5, 0.5, 1.5, 2.5]) 580 581 Apply the same function to a scalar value. 582 583 >>> y = -2 584 >>> np.piecewise(y, [y < 0, y >= 0], [lambda x: -x, lambda x: x]) 585 array(2) 586 587 """ 588 x = asanyarray(x) 589 n2 = len(funclist) 590 591 # undocumented: single condition is promoted to a list of one condition 592 if isscalar(condlist) or ( 593 not isinstance(condlist[0], (list, ndarray)) and x.ndim != 0): 594 condlist = [condlist] 595 596 condlist = array(condlist, dtype=bool) 597 n = len(condlist) 598 599 if n == n2 - 1: # compute the "otherwise" condition. 600 condelse = ~np.any(condlist, axis=0, keepdims=True) 601 condlist = np.concatenate([condlist, condelse], axis=0) 602 n += 1 603 elif n != n2: 604 raise ValueError( 605 "with {} condition(s), either {} or {} functions are expected" 606 .format(n, n, n+1) 607 ) 608 609 y = zeros(x.shape, x.dtype) 610 for cond, func in zip(condlist, funclist): 611 if not isinstance(func, collections.abc.Callable): 612 y[cond] = func 613 else: 614 vals = x[cond] 615 if vals.size > 0: 616 y[cond] = func(vals, *args, **kw) 617 618 return y 619 620 621def _select_dispatcher(condlist, choicelist, default=None): 622 yield from condlist 623 yield from choicelist 624 625 626@array_function_dispatch(_select_dispatcher) 627def select(condlist, choicelist, default=0): 628 """ 629 Return an array drawn from elements in choicelist, depending on conditions. 630 631 Parameters 632 ---------- 633 condlist : list of bool ndarrays 634 The list of conditions which determine from which array in `choicelist` 635 the output elements are taken. When multiple conditions are satisfied, 636 the first one encountered in `condlist` is used. 637 choicelist : list of ndarrays 638 The list of arrays from which the output elements are taken. It has 639 to be of the same length as `condlist`. 640 default : scalar, optional 641 The element inserted in `output` when all conditions evaluate to False. 642 643 Returns 644 ------- 645 output : ndarray 646 The output at position m is the m-th element of the array in 647 `choicelist` where the m-th element of the corresponding array in 648 `condlist` is True. 649 650 See Also 651 -------- 652 where : Return elements from one of two arrays depending on condition. 653 take, choose, compress, diag, diagonal 654 655 Examples 656 -------- 657 >>> x = np.arange(10) 658 >>> condlist = [x<3, x>5] 659 >>> choicelist = [x, x**2] 660 >>> np.select(condlist, choicelist) 661 array([ 0, 1, 2, ..., 49, 64, 81]) 662 663 """ 664 # Check the size of condlist and choicelist are the same, or abort. 665 if len(condlist) != len(choicelist): 666 raise ValueError( 667 'list of cases must be same length as list of conditions') 668 669 # Now that the dtype is known, handle the deprecated select([], []) case 670 if len(condlist) == 0: 671 raise ValueError("select with an empty condition list is not possible") 672 673 choicelist = [np.asarray(choice) for choice in choicelist] 674 choicelist.append(np.asarray(default)) 675 676 # need to get the result type before broadcasting for correct scalar 677 # behaviour 678 dtype = np.result_type(*choicelist) 679 680 # Convert conditions to arrays and broadcast conditions and choices 681 # as the shape is needed for the result. Doing it separately optimizes 682 # for example when all choices are scalars. 683 condlist = np.broadcast_arrays(*condlist) 684 choicelist = np.broadcast_arrays(*choicelist) 685 686 # If cond array is not an ndarray in boolean format or scalar bool, abort. 687 for i, cond in enumerate(condlist): 688 if cond.dtype.type is not np.bool_: 689 raise TypeError( 690 'invalid entry {} in condlist: should be boolean ndarray'.format(i)) 691 692 if choicelist[0].ndim == 0: 693 # This may be common, so avoid the call. 694 result_shape = condlist[0].shape 695 else: 696 result_shape = np.broadcast_arrays(condlist[0], choicelist[0])[0].shape 697 698 result = np.full(result_shape, choicelist[-1], dtype) 699 700 # Use np.copyto to burn each choicelist array onto result, using the 701 # corresponding condlist as a boolean mask. This is done in reverse 702 # order since the first choice should take precedence. 703 choicelist = choicelist[-2::-1] 704 condlist = condlist[::-1] 705 for choice, cond in zip(choicelist, condlist): 706 np.copyto(result, choice, where=cond) 707 708 return result 709 710 711def _copy_dispatcher(a, order=None, subok=None): 712 return (a,) 713 714 715@array_function_dispatch(_copy_dispatcher) 716def copy(a, order='K', subok=False): 717 """ 718 Return an array copy of the given object. 719 720 Parameters 721 ---------- 722 a : array_like 723 Input data. 724 order : {'C', 'F', 'A', 'K'}, optional 725 Controls the memory layout of the copy. 'C' means C-order, 726 'F' means F-order, 'A' means 'F' if `a` is Fortran contiguous, 727 'C' otherwise. 'K' means match the layout of `a` as closely 728 as possible. (Note that this function and :meth:`ndarray.copy` are very 729 similar, but have different default values for their order= 730 arguments.) 731 subok : bool, optional 732 If True, then sub-classes will be passed-through, otherwise the 733 returned array will be forced to be a base-class array (defaults to False). 734 735 .. versionadded:: 1.19.0 736 737 Returns 738 ------- 739 arr : ndarray 740 Array interpretation of `a`. 741 742 See Also 743 -------- 744 ndarray.copy : Preferred method for creating an array copy 745 746 Notes 747 ----- 748 This is equivalent to: 749 750 >>> np.array(a, copy=True) #doctest: +SKIP 751 752 Examples 753 -------- 754 Create an array x, with a reference y and a copy z: 755 756 >>> x = np.array([1, 2, 3]) 757 >>> y = x 758 >>> z = np.copy(x) 759 760 Note that, when we modify x, y changes, but not z: 761 762 >>> x[0] = 10 763 >>> x[0] == y[0] 764 True 765 >>> x[0] == z[0] 766 False 767 768 Note that np.copy is a shallow copy and will not copy object 769 elements within arrays. This is mainly important for arrays 770 containing Python objects. The new array will contain the 771 same object which may lead to surprises if that object can 772 be modified (is mutable): 773 774 >>> a = np.array([1, 'm', [2, 3, 4]], dtype=object) 775 >>> b = np.copy(a) 776 >>> b[2][0] = 10 777 >>> a 778 array([1, 'm', list([10, 3, 4])], dtype=object) 779 780 To ensure all elements within an ``object`` array are copied, 781 use `copy.deepcopy`: 782 783 >>> import copy 784 >>> a = np.array([1, 'm', [2, 3, 4]], dtype=object) 785 >>> c = copy.deepcopy(a) 786 >>> c[2][0] = 10 787 >>> c 788 array([1, 'm', list([10, 3, 4])], dtype=object) 789 >>> a 790 array([1, 'm', list([2, 3, 4])], dtype=object) 791 792 """ 793 return array(a, order=order, subok=subok, copy=True) 794 795# Basic operations 796 797 798def _gradient_dispatcher(f, *varargs, axis=None, edge_order=None): 799 yield f 800 yield from varargs 801 802 803@array_function_dispatch(_gradient_dispatcher) 804def gradient(f, *varargs, axis=None, edge_order=1): 805 """ 806 Return the gradient of an N-dimensional array. 807 808 The gradient is computed using second order accurate central differences 809 in the interior points and either first or second order accurate one-sides 810 (forward or backwards) differences at the boundaries. 811 The returned gradient hence has the same shape as the input array. 812 813 Parameters 814 ---------- 815 f : array_like 816 An N-dimensional array containing samples of a scalar function. 817 varargs : list of scalar or array, optional 818 Spacing between f values. Default unitary spacing for all dimensions. 819 Spacing can be specified using: 820 821 1. single scalar to specify a sample distance for all dimensions. 822 2. N scalars to specify a constant sample distance for each dimension. 823 i.e. `dx`, `dy`, `dz`, ... 824 3. N arrays to specify the coordinates of the values along each 825 dimension of F. The length of the array must match the size of 826 the corresponding dimension 827 4. Any combination of N scalars/arrays with the meaning of 2. and 3. 828 829 If `axis` is given, the number of varargs must equal the number of axes. 830 Default: 1. 831 832 edge_order : {1, 2}, optional 833 Gradient is calculated using N-th order accurate differences 834 at the boundaries. Default: 1. 835 836 .. versionadded:: 1.9.1 837 838 axis : None or int or tuple of ints, optional 839 Gradient is calculated only along the given axis or axes 840 The default (axis = None) is to calculate the gradient for all the axes 841 of the input array. axis may be negative, in which case it counts from 842 the last to the first axis. 843 844 .. versionadded:: 1.11.0 845 846 Returns 847 ------- 848 gradient : ndarray or list of ndarray 849 A set of ndarrays (or a single ndarray if there is only one dimension) 850 corresponding to the derivatives of f with respect to each dimension. 851 Each derivative has the same shape as f. 852 853 Examples 854 -------- 855 >>> f = np.array([1, 2, 4, 7, 11, 16], dtype=float) 856 >>> np.gradient(f) 857 array([1. , 1.5, 2.5, 3.5, 4.5, 5. ]) 858 >>> np.gradient(f, 2) 859 array([0.5 , 0.75, 1.25, 1.75, 2.25, 2.5 ]) 860 861 Spacing can be also specified with an array that represents the coordinates 862 of the values F along the dimensions. 863 For instance a uniform spacing: 864 865 >>> x = np.arange(f.size) 866 >>> np.gradient(f, x) 867 array([1. , 1.5, 2.5, 3.5, 4.5, 5. ]) 868 869 Or a non uniform one: 870 871 >>> x = np.array([0., 1., 1.5, 3.5, 4., 6.], dtype=float) 872 >>> np.gradient(f, x) 873 array([1. , 3. , 3.5, 6.7, 6.9, 2.5]) 874 875 For two dimensional arrays, the return will be two arrays ordered by 876 axis. In this example the first array stands for the gradient in 877 rows and the second one in columns direction: 878 879 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float)) 880 [array([[ 2., 2., -1.], 881 [ 2., 2., -1.]]), array([[1. , 2.5, 4. ], 882 [1. , 1. , 1. ]])] 883 884 In this example the spacing is also specified: 885 uniform for axis=0 and non uniform for axis=1 886 887 >>> dx = 2. 888 >>> y = [1., 1.5, 3.5] 889 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), dx, y) 890 [array([[ 1. , 1. , -0.5], 891 [ 1. , 1. , -0.5]]), array([[2. , 2. , 2. ], 892 [2. , 1.7, 0.5]])] 893 894 It is possible to specify how boundaries are treated using `edge_order` 895 896 >>> x = np.array([0, 1, 2, 3, 4]) 897 >>> f = x**2 898 >>> np.gradient(f, edge_order=1) 899 array([1., 2., 4., 6., 7.]) 900 >>> np.gradient(f, edge_order=2) 901 array([0., 2., 4., 6., 8.]) 902 903 The `axis` keyword can be used to specify a subset of axes of which the 904 gradient is calculated 905 906 >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), axis=0) 907 array([[ 2., 2., -1.], 908 [ 2., 2., -1.]]) 909 910 Notes 911 ----- 912 Assuming that :math:`f\\in C^{3}` (i.e., :math:`f` has at least 3 continuous 913 derivatives) and let :math:`h_{*}` be a non-homogeneous stepsize, we 914 minimize the "consistency error" :math:`\\eta_{i}` between the true gradient 915 and its estimate from a linear combination of the neighboring grid-points: 916 917 .. math:: 918 919 \\eta_{i} = f_{i}^{\\left(1\\right)} - 920 \\left[ \\alpha f\\left(x_{i}\\right) + 921 \\beta f\\left(x_{i} + h_{d}\\right) + 922 \\gamma f\\left(x_{i}-h_{s}\\right) 923 \\right] 924 925 By substituting :math:`f(x_{i} + h_{d})` and :math:`f(x_{i} - h_{s})` 926 with their Taylor series expansion, this translates into solving 927 the following the linear system: 928 929 .. math:: 930 931 \\left\\{ 932 \\begin{array}{r} 933 \\alpha+\\beta+\\gamma=0 \\\\ 934 \\beta h_{d}-\\gamma h_{s}=1 \\\\ 935 \\beta h_{d}^{2}+\\gamma h_{s}^{2}=0 936 \\end{array} 937 \\right. 938 939 The resulting approximation of :math:`f_{i}^{(1)}` is the following: 940 941 .. math:: 942 943 \\hat f_{i}^{(1)} = 944 \\frac{ 945 h_{s}^{2}f\\left(x_{i} + h_{d}\\right) 946 + \\left(h_{d}^{2} - h_{s}^{2}\\right)f\\left(x_{i}\\right) 947 - h_{d}^{2}f\\left(x_{i}-h_{s}\\right)} 948 { h_{s}h_{d}\\left(h_{d} + h_{s}\\right)} 949 + \\mathcal{O}\\left(\\frac{h_{d}h_{s}^{2} 950 + h_{s}h_{d}^{2}}{h_{d} 951 + h_{s}}\\right) 952 953 It is worth noting that if :math:`h_{s}=h_{d}` 954 (i.e., data are evenly spaced) 955 we find the standard second order approximation: 956 957 .. math:: 958 959 \\hat f_{i}^{(1)}= 960 \\frac{f\\left(x_{i+1}\\right) - f\\left(x_{i-1}\\right)}{2h} 961 + \\mathcal{O}\\left(h^{2}\\right) 962 963 With a similar procedure the forward/backward approximations used for 964 boundaries can be derived. 965 966 References 967 ---------- 968 .. [1] Quarteroni A., Sacco R., Saleri F. (2007) Numerical Mathematics 969 (Texts in Applied Mathematics). New York: Springer. 970 .. [2] Durran D. R. (1999) Numerical Methods for Wave Equations 971 in Geophysical Fluid Dynamics. New York: Springer. 972 .. [3] Fornberg B. (1988) Generation of Finite Difference Formulas on 973 Arbitrarily Spaced Grids, 974 Mathematics of Computation 51, no. 184 : 699-706. 975 `PDF <http://www.ams.org/journals/mcom/1988-51-184/ 976 S0025-5718-1988-0935077-0/S0025-5718-1988-0935077-0.pdf>`_. 977 """ 978 f = np.asanyarray(f) 979 N = f.ndim # number of dimensions 980 981 if axis is None: 982 axes = tuple(range(N)) 983 else: 984 axes = _nx.normalize_axis_tuple(axis, N) 985 986 len_axes = len(axes) 987 n = len(varargs) 988 if n == 0: 989 # no spacing argument - use 1 in all axes 990 dx = [1.0] * len_axes 991 elif n == 1 and np.ndim(varargs[0]) == 0: 992 # single scalar for all axes 993 dx = varargs * len_axes 994 elif n == len_axes: 995 # scalar or 1d array for each axis 996 dx = list(varargs) 997 for i, distances in enumerate(dx): 998 distances = np.asanyarray(distances) 999 if distances.ndim == 0: 1000 continue 1001 elif distances.ndim != 1: 1002 raise ValueError("distances must be either scalars or 1d") 1003 if len(distances) != f.shape[axes[i]]: 1004 raise ValueError("when 1d, distances must match " 1005 "the length of the corresponding dimension") 1006 if np.issubdtype(distances.dtype, np.integer): 1007 # Convert numpy integer types to float64 to avoid modular 1008 # arithmetic in np.diff(distances). 1009 distances = distances.astype(np.float64) 1010 diffx = np.diff(distances) 1011 # if distances are constant reduce to the scalar case 1012 # since it brings a consistent speedup 1013 if (diffx == diffx[0]).all(): 1014 diffx = diffx[0] 1015 dx[i] = diffx 1016 else: 1017 raise TypeError("invalid number of arguments") 1018 1019 if edge_order > 2: 1020 raise ValueError("'edge_order' greater than 2 not supported") 1021 1022 # use central differences on interior and one-sided differences on the 1023 # endpoints. This preserves second order-accuracy over the full domain. 1024 1025 outvals = [] 1026 1027 # create slice objects --- initially all are [:, :, ..., :] 1028 slice1 = [slice(None)]*N 1029 slice2 = [slice(None)]*N 1030 slice3 = [slice(None)]*N 1031 slice4 = [slice(None)]*N 1032 1033 otype = f.dtype 1034 if otype.type is np.datetime64: 1035 # the timedelta dtype with the same unit information 1036 otype = np.dtype(otype.name.replace('datetime', 'timedelta')) 1037 # view as timedelta to allow addition 1038 f = f.view(otype) 1039 elif otype.type is np.timedelta64: 1040 pass 1041 elif np.issubdtype(otype, np.inexact): 1042 pass 1043 else: 1044 # All other types convert to floating point. 1045 # First check if f is a numpy integer type; if so, convert f to float64 1046 # to avoid modular arithmetic when computing the changes in f. 1047 if np.issubdtype(otype, np.integer): 1048 f = f.astype(np.float64) 1049 otype = np.float64 1050 1051 for axis, ax_dx in zip(axes, dx): 1052 if f.shape[axis] < edge_order + 1: 1053 raise ValueError( 1054 "Shape of array too small to calculate a numerical gradient, " 1055 "at least (edge_order + 1) elements are required.") 1056 # result allocation 1057 out = np.empty_like(f, dtype=otype) 1058 1059 # spacing for the current axis 1060 uniform_spacing = np.ndim(ax_dx) == 0 1061 1062 # Numerical differentiation: 2nd order interior 1063 slice1[axis] = slice(1, -1) 1064 slice2[axis] = slice(None, -2) 1065 slice3[axis] = slice(1, -1) 1066 slice4[axis] = slice(2, None) 1067 1068 if uniform_spacing: 1069 out[tuple(slice1)] = (f[tuple(slice4)] - f[tuple(slice2)]) / (2. * ax_dx) 1070 else: 1071 dx1 = ax_dx[0:-1] 1072 dx2 = ax_dx[1:] 1073 a = -(dx2)/(dx1 * (dx1 + dx2)) 1074 b = (dx2 - dx1) / (dx1 * dx2) 1075 c = dx1 / (dx2 * (dx1 + dx2)) 1076 # fix the shape for broadcasting 1077 shape = np.ones(N, dtype=int) 1078 shape[axis] = -1 1079 a.shape = b.shape = c.shape = shape 1080 # 1D equivalent -- out[1:-1] = a * f[:-2] + b * f[1:-1] + c * f[2:] 1081 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)] 1082 1083 # Numerical differentiation: 1st order edges 1084 if edge_order == 1: 1085 slice1[axis] = 0 1086 slice2[axis] = 1 1087 slice3[axis] = 0 1088 dx_0 = ax_dx if uniform_spacing else ax_dx[0] 1089 # 1D equivalent -- out[0] = (f[1] - f[0]) / (x[1] - x[0]) 1090 out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0 1091 1092 slice1[axis] = -1 1093 slice2[axis] = -1 1094 slice3[axis] = -2 1095 dx_n = ax_dx if uniform_spacing else ax_dx[-1] 1096 # 1D equivalent -- out[-1] = (f[-1] - f[-2]) / (x[-1] - x[-2]) 1097 out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_n 1098 1099 # Numerical differentiation: 2nd order edges 1100 else: 1101 slice1[axis] = 0 1102 slice2[axis] = 0 1103 slice3[axis] = 1 1104 slice4[axis] = 2 1105 if uniform_spacing: 1106 a = -1.5 / ax_dx 1107 b = 2. / ax_dx 1108 c = -0.5 / ax_dx 1109 else: 1110 dx1 = ax_dx[0] 1111 dx2 = ax_dx[1] 1112 a = -(2. * dx1 + dx2)/(dx1 * (dx1 + dx2)) 1113 b = (dx1 + dx2) / (dx1 * dx2) 1114 c = - dx1 / (dx2 * (dx1 + dx2)) 1115 # 1D equivalent -- out[0] = a * f[0] + b * f[1] + c * f[2] 1116 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)] 1117 1118 slice1[axis] = -1 1119 slice2[axis] = -3 1120 slice3[axis] = -2 1121 slice4[axis] = -1 1122 if uniform_spacing: 1123 a = 0.5 / ax_dx 1124 b = -2. / ax_dx 1125 c = 1.5 / ax_dx 1126 else: 1127 dx1 = ax_dx[-2] 1128 dx2 = ax_dx[-1] 1129 a = (dx2) / (dx1 * (dx1 + dx2)) 1130 b = - (dx2 + dx1) / (dx1 * dx2) 1131 c = (2. * dx2 + dx1) / (dx2 * (dx1 + dx2)) 1132 # 1D equivalent -- out[-1] = a * f[-3] + b * f[-2] + c * f[-1] 1133 out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)] 1134 1135 outvals.append(out) 1136 1137 # reset the slice object in this dimension to ":" 1138 slice1[axis] = slice(None) 1139 slice2[axis] = slice(None) 1140 slice3[axis] = slice(None) 1141 slice4[axis] = slice(None) 1142 1143 if len_axes == 1: 1144 return outvals[0] 1145 else: 1146 return outvals 1147 1148 1149def _diff_dispatcher(a, n=None, axis=None, prepend=None, append=None): 1150 return (a, prepend, append) 1151 1152 1153@array_function_dispatch(_diff_dispatcher) 1154def diff(a, n=1, axis=-1, prepend=np._NoValue, append=np._NoValue): 1155 """ 1156 Calculate the n-th discrete difference along the given axis. 1157 1158 The first difference is given by ``out[i] = a[i+1] - a[i]`` along 1159 the given axis, higher differences are calculated by using `diff` 1160 recursively. 1161 1162 Parameters 1163 ---------- 1164 a : array_like 1165 Input array 1166 n : int, optional 1167 The number of times values are differenced. If zero, the input 1168 is returned as-is. 1169 axis : int, optional 1170 The axis along which the difference is taken, default is the 1171 last axis. 1172 prepend, append : array_like, optional 1173 Values to prepend or append to `a` along axis prior to 1174 performing the difference. Scalar values are expanded to 1175 arrays with length 1 in the direction of axis and the shape 1176 of the input array in along all other axes. Otherwise the 1177 dimension and shape must match `a` except along axis. 1178 1179 .. versionadded:: 1.16.0 1180 1181 Returns 1182 ------- 1183 diff : ndarray 1184 The n-th differences. The shape of the output is the same as `a` 1185 except along `axis` where the dimension is smaller by `n`. The 1186 type of the output is the same as the type of the difference 1187 between any two elements of `a`. This is the same as the type of 1188 `a` in most cases. A notable exception is `datetime64`, which 1189 results in a `timedelta64` output array. 1190 1191 See Also 1192 -------- 1193 gradient, ediff1d, cumsum 1194 1195 Notes 1196 ----- 1197 Type is preserved for boolean arrays, so the result will contain 1198 `False` when consecutive elements are the same and `True` when they 1199 differ. 1200 1201 For unsigned integer arrays, the results will also be unsigned. This 1202 should not be surprising, as the result is consistent with 1203 calculating the difference directly: 1204 1205 >>> u8_arr = np.array([1, 0], dtype=np.uint8) 1206 >>> np.diff(u8_arr) 1207 array([255], dtype=uint8) 1208 >>> u8_arr[1,...] - u8_arr[0,...] 1209 255 1210 1211 If this is not desirable, then the array should be cast to a larger 1212 integer type first: 1213 1214 >>> i16_arr = u8_arr.astype(np.int16) 1215 >>> np.diff(i16_arr) 1216 array([-1], dtype=int16) 1217 1218 Examples 1219 -------- 1220 >>> x = np.array([1, 2, 4, 7, 0]) 1221 >>> np.diff(x) 1222 array([ 1, 2, 3, -7]) 1223 >>> np.diff(x, n=2) 1224 array([ 1, 1, -10]) 1225 1226 >>> x = np.array([[1, 3, 6, 10], [0, 5, 6, 8]]) 1227 >>> np.diff(x) 1228 array([[2, 3, 4], 1229 [5, 1, 2]]) 1230 >>> np.diff(x, axis=0) 1231 array([[-1, 2, 0, -2]]) 1232 1233 >>> x = np.arange('1066-10-13', '1066-10-16', dtype=np.datetime64) 1234 >>> np.diff(x) 1235 array([1, 1], dtype='timedelta64[D]') 1236 1237 """ 1238 if n == 0: 1239 return a 1240 if n < 0: 1241 raise ValueError( 1242 "order must be non-negative but got " + repr(n)) 1243 1244 a = asanyarray(a) 1245 nd = a.ndim 1246 if nd == 0: 1247 raise ValueError("diff requires input that is at least one dimensional") 1248 axis = normalize_axis_index(axis, nd) 1249 1250 combined = [] 1251 if prepend is not np._NoValue: 1252 prepend = np.asanyarray(prepend) 1253 if prepend.ndim == 0: 1254 shape = list(a.shape) 1255 shape[axis] = 1 1256 prepend = np.broadcast_to(prepend, tuple(shape)) 1257 combined.append(prepend) 1258 1259 combined.append(a) 1260 1261 if append is not np._NoValue: 1262 append = np.asanyarray(append) 1263 if append.ndim == 0: 1264 shape = list(a.shape) 1265 shape[axis] = 1 1266 append = np.broadcast_to(append, tuple(shape)) 1267 combined.append(append) 1268 1269 if len(combined) > 1: 1270 a = np.concatenate(combined, axis) 1271 1272 slice1 = [slice(None)] * nd 1273 slice2 = [slice(None)] * nd 1274 slice1[axis] = slice(1, None) 1275 slice2[axis] = slice(None, -1) 1276 slice1 = tuple(slice1) 1277 slice2 = tuple(slice2) 1278 1279 op = not_equal if a.dtype == np.bool_ else subtract 1280 for _ in range(n): 1281 a = op(a[slice1], a[slice2]) 1282 1283 return a 1284 1285 1286def _interp_dispatcher(x, xp, fp, left=None, right=None, period=None): 1287 return (x, xp, fp) 1288 1289 1290@array_function_dispatch(_interp_dispatcher) 1291def interp(x, xp, fp, left=None, right=None, period=None): 1292 """ 1293 One-dimensional linear interpolation. 1294 1295 Returns the one-dimensional piecewise linear interpolant to a function 1296 with given discrete data points (`xp`, `fp`), evaluated at `x`. 1297 1298 Parameters 1299 ---------- 1300 x : array_like 1301 The x-coordinates at which to evaluate the interpolated values. 1302 1303 xp : 1-D sequence of floats 1304 The x-coordinates of the data points, must be increasing if argument 1305 `period` is not specified. Otherwise, `xp` is internally sorted after 1306 normalizing the periodic boundaries with ``xp = xp % period``. 1307 1308 fp : 1-D sequence of float or complex 1309 The y-coordinates of the data points, same length as `xp`. 1310 1311 left : optional float or complex corresponding to fp 1312 Value to return for `x < xp[0]`, default is `fp[0]`. 1313 1314 right : optional float or complex corresponding to fp 1315 Value to return for `x > xp[-1]`, default is `fp[-1]`. 1316 1317 period : None or float, optional 1318 A period for the x-coordinates. This parameter allows the proper 1319 interpolation of angular x-coordinates. Parameters `left` and `right` 1320 are ignored if `period` is specified. 1321 1322 .. versionadded:: 1.10.0 1323 1324 Returns 1325 ------- 1326 y : float or complex (corresponding to fp) or ndarray 1327 The interpolated values, same shape as `x`. 1328 1329 Raises 1330 ------ 1331 ValueError 1332 If `xp` and `fp` have different length 1333 If `xp` or `fp` are not 1-D sequences 1334 If `period == 0` 1335 1336 See Also 1337 -------- 1338 scipy.interpolate 1339 1340 Notes 1341 ----- 1342 The x-coordinate sequence is expected to be increasing, but this is not 1343 explicitly enforced. However, if the sequence `xp` is non-increasing, 1344 interpolation results are meaningless. 1345 1346 Note that, since NaN is unsortable, `xp` also cannot contain NaNs. 1347 1348 A simple check for `xp` being strictly increasing is:: 1349 1350 np.all(np.diff(xp) > 0) 1351 1352 Examples 1353 -------- 1354 >>> xp = [1, 2, 3] 1355 >>> fp = [3, 2, 0] 1356 >>> np.interp(2.5, xp, fp) 1357 1.0 1358 >>> np.interp([0, 1, 1.5, 2.72, 3.14], xp, fp) 1359 array([3. , 3. , 2.5 , 0.56, 0. ]) 1360 >>> UNDEF = -99.0 1361 >>> np.interp(3.14, xp, fp, right=UNDEF) 1362 -99.0 1363 1364 Plot an interpolant to the sine function: 1365 1366 >>> x = np.linspace(0, 2*np.pi, 10) 1367 >>> y = np.sin(x) 1368 >>> xvals = np.linspace(0, 2*np.pi, 50) 1369 >>> yinterp = np.interp(xvals, x, y) 1370 >>> import matplotlib.pyplot as plt 1371 >>> plt.plot(x, y, 'o') 1372 [<matplotlib.lines.Line2D object at 0x...>] 1373 >>> plt.plot(xvals, yinterp, '-x') 1374 [<matplotlib.lines.Line2D object at 0x...>] 1375 >>> plt.show() 1376 1377 Interpolation with periodic x-coordinates: 1378 1379 >>> x = [-180, -170, -185, 185, -10, -5, 0, 365] 1380 >>> xp = [190, -190, 350, -350] 1381 >>> fp = [5, 10, 3, 4] 1382 >>> np.interp(x, xp, fp, period=360) 1383 array([7.5 , 5. , 8.75, 6.25, 3. , 3.25, 3.5 , 3.75]) 1384 1385 Complex interpolation: 1386 1387 >>> x = [1.5, 4.0] 1388 >>> xp = [2,3,5] 1389 >>> fp = [1.0j, 0, 2+3j] 1390 >>> np.interp(x, xp, fp) 1391 array([0.+1.j , 1.+1.5j]) 1392 1393 """ 1394 1395 fp = np.asarray(fp) 1396 1397 if np.iscomplexobj(fp): 1398 interp_func = compiled_interp_complex 1399 input_dtype = np.complex128 1400 else: 1401 interp_func = compiled_interp 1402 input_dtype = np.float64 1403 1404 if period is not None: 1405 if period == 0: 1406 raise ValueError("period must be a non-zero value") 1407 period = abs(period) 1408 left = None 1409 right = None 1410 1411 x = np.asarray(x, dtype=np.float64) 1412 xp = np.asarray(xp, dtype=np.float64) 1413 fp = np.asarray(fp, dtype=input_dtype) 1414 1415 if xp.ndim != 1 or fp.ndim != 1: 1416 raise ValueError("Data points must be 1-D sequences") 1417 if xp.shape[0] != fp.shape[0]: 1418 raise ValueError("fp and xp are not of the same length") 1419 # normalizing periodic boundaries 1420 x = x % period 1421 xp = xp % period 1422 asort_xp = np.argsort(xp) 1423 xp = xp[asort_xp] 1424 fp = fp[asort_xp] 1425 xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period)) 1426 fp = np.concatenate((fp[-1:], fp, fp[0:1])) 1427 1428 return interp_func(x, xp, fp, left, right) 1429 1430 1431def _angle_dispatcher(z, deg=None): 1432 return (z,) 1433 1434 1435@array_function_dispatch(_angle_dispatcher) 1436def angle(z, deg=False): 1437 """ 1438 Return the angle of the complex argument. 1439 1440 Parameters 1441 ---------- 1442 z : array_like 1443 A complex number or sequence of complex numbers. 1444 deg : bool, optional 1445 Return angle in degrees if True, radians if False (default). 1446 1447 Returns 1448 ------- 1449 angle : ndarray or scalar 1450 The counterclockwise angle from the positive real axis on the complex 1451 plane in the range ``(-pi, pi]``, with dtype as numpy.float64. 1452 1453 .. versionchanged:: 1.16.0 1454 This function works on subclasses of ndarray like `ma.array`. 1455 1456 See Also 1457 -------- 1458 arctan2 1459 absolute 1460 1461 Notes 1462 ----- 1463 Although the angle of the complex number 0 is undefined, ``numpy.angle(0)`` 1464 returns the value 0. 1465 1466 Examples 1467 -------- 1468 >>> np.angle([1.0, 1.0j, 1+1j]) # in radians 1469 array([ 0. , 1.57079633, 0.78539816]) # may vary 1470 >>> np.angle(1+1j, deg=True) # in degrees 1471 45.0 1472 1473 """ 1474 z = asanyarray(z) 1475 if issubclass(z.dtype.type, _nx.complexfloating): 1476 zimag = z.imag 1477 zreal = z.real 1478 else: 1479 zimag = 0 1480 zreal = z 1481 1482 a = arctan2(zimag, zreal) 1483 if deg: 1484 a *= 180/pi 1485 return a 1486 1487 1488def _unwrap_dispatcher(p, discont=None, axis=None): 1489 return (p,) 1490 1491 1492@array_function_dispatch(_unwrap_dispatcher) 1493def unwrap(p, discont=pi, axis=-1): 1494 """ 1495 Unwrap by changing deltas between values to 2*pi complement. 1496 1497 Unwrap radian phase `p` by changing absolute jumps greater than 1498 `discont` to their 2*pi complement along the given axis. 1499 1500 Parameters 1501 ---------- 1502 p : array_like 1503 Input array. 1504 discont : float, optional 1505 Maximum discontinuity between values, default is ``pi``. 1506 axis : int, optional 1507 Axis along which unwrap will operate, default is the last axis. 1508 1509 Returns 1510 ------- 1511 out : ndarray 1512 Output array. 1513 1514 See Also 1515 -------- 1516 rad2deg, deg2rad 1517 1518 Notes 1519 ----- 1520 If the discontinuity in `p` is smaller than ``pi``, but larger than 1521 `discont`, no unwrapping is done because taking the 2*pi complement 1522 would only make the discontinuity larger. 1523 1524 Examples 1525 -------- 1526 >>> phase = np.linspace(0, np.pi, num=5) 1527 >>> phase[3:] += np.pi 1528 >>> phase 1529 array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531]) # may vary 1530 >>> np.unwrap(phase) 1531 array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ]) # may vary 1532 1533 """ 1534 p = asarray(p) 1535 nd = p.ndim 1536 dd = diff(p, axis=axis) 1537 slice1 = [slice(None, None)]*nd # full slices 1538 slice1[axis] = slice(1, None) 1539 slice1 = tuple(slice1) 1540 ddmod = mod(dd + pi, 2*pi) - pi 1541 _nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0)) 1542 ph_correct = ddmod - dd 1543 _nx.copyto(ph_correct, 0, where=abs(dd) < discont) 1544 up = array(p, copy=True, dtype='d') 1545 up[slice1] = p[slice1] + ph_correct.cumsum(axis) 1546 return up 1547 1548 1549def _sort_complex(a): 1550 return (a,) 1551 1552 1553@array_function_dispatch(_sort_complex) 1554def sort_complex(a): 1555 """ 1556 Sort a complex array using the real part first, then the imaginary part. 1557 1558 Parameters 1559 ---------- 1560 a : array_like 1561 Input array 1562 1563 Returns 1564 ------- 1565 out : complex ndarray 1566 Always returns a sorted complex array. 1567 1568 Examples 1569 -------- 1570 >>> np.sort_complex([5, 3, 6, 2, 1]) 1571 array([1.+0.j, 2.+0.j, 3.+0.j, 5.+0.j, 6.+0.j]) 1572 1573 >>> np.sort_complex([1 + 2j, 2 - 1j, 3 - 2j, 3 - 3j, 3 + 5j]) 1574 array([1.+2.j, 2.-1.j, 3.-3.j, 3.-2.j, 3.+5.j]) 1575 1576 """ 1577 b = array(a, copy=True) 1578 b.sort() 1579 if not issubclass(b.dtype.type, _nx.complexfloating): 1580 if b.dtype.char in 'bhBH': 1581 return b.astype('F') 1582 elif b.dtype.char == 'g': 1583 return b.astype('G') 1584 else: 1585 return b.astype('D') 1586 else: 1587 return b 1588 1589 1590def _trim_zeros(filt, trim=None): 1591 return (filt,) 1592 1593 1594@array_function_dispatch(_trim_zeros) 1595def trim_zeros(filt, trim='fb'): 1596 """ 1597 Trim the leading and/or trailing zeros from a 1-D array or sequence. 1598 1599 Parameters 1600 ---------- 1601 filt : 1-D array or sequence 1602 Input array. 1603 trim : str, optional 1604 A string with 'f' representing trim from front and 'b' to trim from 1605 back. Default is 'fb', trim zeros from both front and back of the 1606 array. 1607 1608 Returns 1609 ------- 1610 trimmed : 1-D array or sequence 1611 The result of trimming the input. The input data type is preserved. 1612 1613 Examples 1614 -------- 1615 >>> a = np.array((0, 0, 0, 1, 2, 3, 0, 2, 1, 0)) 1616 >>> np.trim_zeros(a) 1617 array([1, 2, 3, 0, 2, 1]) 1618 1619 >>> np.trim_zeros(a, 'b') 1620 array([0, 0, 0, ..., 0, 2, 1]) 1621 1622 The input data type is preserved, list/tuple in means list/tuple out. 1623 1624 >>> np.trim_zeros([0, 1, 2, 0]) 1625 [1, 2] 1626 1627 """ 1628 1629 first = 0 1630 trim = trim.upper() 1631 if 'F' in trim: 1632 for i in filt: 1633 if i != 0.: 1634 break 1635 else: 1636 first = first + 1 1637 last = len(filt) 1638 if 'B' in trim: 1639 for i in filt[::-1]: 1640 if i != 0.: 1641 break 1642 else: 1643 last = last - 1 1644 return filt[first:last] 1645 1646 1647def _extract_dispatcher(condition, arr): 1648 return (condition, arr) 1649 1650 1651@array_function_dispatch(_extract_dispatcher) 1652def extract(condition, arr): 1653 """ 1654 Return the elements of an array that satisfy some condition. 1655 1656 This is equivalent to ``np.compress(ravel(condition), ravel(arr))``. If 1657 `condition` is boolean ``np.extract`` is equivalent to ``arr[condition]``. 1658 1659 Note that `place` does the exact opposite of `extract`. 1660 1661 Parameters 1662 ---------- 1663 condition : array_like 1664 An array whose nonzero or True entries indicate the elements of `arr` 1665 to extract. 1666 arr : array_like 1667 Input array of the same size as `condition`. 1668 1669 Returns 1670 ------- 1671 extract : ndarray 1672 Rank 1 array of values from `arr` where `condition` is True. 1673 1674 See Also 1675 -------- 1676 take, put, copyto, compress, place 1677 1678 Examples 1679 -------- 1680 >>> arr = np.arange(12).reshape((3, 4)) 1681 >>> arr 1682 array([[ 0, 1, 2, 3], 1683 [ 4, 5, 6, 7], 1684 [ 8, 9, 10, 11]]) 1685 >>> condition = np.mod(arr, 3)==0 1686 >>> condition 1687 array([[ True, False, False, True], 1688 [False, False, True, False], 1689 [False, True, False, False]]) 1690 >>> np.extract(condition, arr) 1691 array([0, 3, 6, 9]) 1692 1693 1694 If `condition` is boolean: 1695 1696 >>> arr[condition] 1697 array([0, 3, 6, 9]) 1698 1699 """ 1700 return _nx.take(ravel(arr), nonzero(ravel(condition))[0]) 1701 1702 1703def _place_dispatcher(arr, mask, vals): 1704 return (arr, mask, vals) 1705 1706 1707@array_function_dispatch(_place_dispatcher) 1708def place(arr, mask, vals): 1709 """ 1710 Change elements of an array based on conditional and input values. 1711 1712 Similar to ``np.copyto(arr, vals, where=mask)``, the difference is that 1713 `place` uses the first N elements of `vals`, where N is the number of 1714 True values in `mask`, while `copyto` uses the elements where `mask` 1715 is True. 1716 1717 Note that `extract` does the exact opposite of `place`. 1718 1719 Parameters 1720 ---------- 1721 arr : ndarray 1722 Array to put data into. 1723 mask : array_like 1724 Boolean mask array. Must have the same size as `a`. 1725 vals : 1-D sequence 1726 Values to put into `a`. Only the first N elements are used, where 1727 N is the number of True values in `mask`. If `vals` is smaller 1728 than N, it will be repeated, and if elements of `a` are to be masked, 1729 this sequence must be non-empty. 1730 1731 See Also 1732 -------- 1733 copyto, put, take, extract 1734 1735 Examples 1736 -------- 1737 >>> arr = np.arange(6).reshape(2, 3) 1738 >>> np.place(arr, arr>2, [44, 55]) 1739 >>> arr 1740 array([[ 0, 1, 2], 1741 [44, 55, 44]]) 1742 1743 """ 1744 if not isinstance(arr, np.ndarray): 1745 raise TypeError("argument 1 must be numpy.ndarray, " 1746 "not {name}".format(name=type(arr).__name__)) 1747 1748 return _insert(arr, mask, vals) 1749 1750 1751def disp(mesg, device=None, linefeed=True): 1752 """ 1753 Display a message on a device. 1754 1755 Parameters 1756 ---------- 1757 mesg : str 1758 Message to display. 1759 device : object 1760 Device to write message. If None, defaults to ``sys.stdout`` which is 1761 very similar to ``print``. `device` needs to have ``write()`` and 1762 ``flush()`` methods. 1763 linefeed : bool, optional 1764 Option whether to print a line feed or not. Defaults to True. 1765 1766 Raises 1767 ------ 1768 AttributeError 1769 If `device` does not have a ``write()`` or ``flush()`` method. 1770 1771 Examples 1772 -------- 1773 Besides ``sys.stdout``, a file-like object can also be used as it has 1774 both required methods: 1775 1776 >>> from io import StringIO 1777 >>> buf = StringIO() 1778 >>> np.disp(u'"Display" in a file', device=buf) 1779 >>> buf.getvalue() 1780 '"Display" in a file\\n' 1781 1782 """ 1783 if device is None: 1784 device = sys.stdout 1785 if linefeed: 1786 device.write('%s\n' % mesg) 1787 else: 1788 device.write('%s' % mesg) 1789 device.flush() 1790 return 1791 1792 1793# See https://docs.scipy.org/doc/numpy/reference/c-api.generalized-ufuncs.html 1794_DIMENSION_NAME = r'\w+' 1795_CORE_DIMENSION_LIST = '(?:{0:}(?:,{0:})*)?'.format(_DIMENSION_NAME) 1796_ARGUMENT = r'\({}\)'.format(_CORE_DIMENSION_LIST) 1797_ARGUMENT_LIST = '{0:}(?:,{0:})*'.format(_ARGUMENT) 1798_SIGNATURE = '^{0:}->{0:}$'.format(_ARGUMENT_LIST) 1799 1800 1801def _parse_gufunc_signature(signature): 1802 """ 1803 Parse string signatures for a generalized universal function. 1804 1805 Arguments 1806 --------- 1807 signature : string 1808 Generalized universal function signature, e.g., ``(m,n),(n,p)->(m,p)`` 1809 for ``np.matmul``. 1810 1811 Returns 1812 ------- 1813 Tuple of input and output core dimensions parsed from the signature, each 1814 of the form List[Tuple[str, ...]]. 1815 """ 1816 if not re.match(_SIGNATURE, signature): 1817 raise ValueError( 1818 'not a valid gufunc signature: {}'.format(signature)) 1819 return tuple([tuple(re.findall(_DIMENSION_NAME, arg)) 1820 for arg in re.findall(_ARGUMENT, arg_list)] 1821 for arg_list in signature.split('->')) 1822 1823 1824def _update_dim_sizes(dim_sizes, arg, core_dims): 1825 """ 1826 Incrementally check and update core dimension sizes for a single argument. 1827 1828 Arguments 1829 --------- 1830 dim_sizes : Dict[str, int] 1831 Sizes of existing core dimensions. Will be updated in-place. 1832 arg : ndarray 1833 Argument to examine. 1834 core_dims : Tuple[str, ...] 1835 Core dimensions for this argument. 1836 """ 1837 if not core_dims: 1838 return 1839 1840 num_core_dims = len(core_dims) 1841 if arg.ndim < num_core_dims: 1842 raise ValueError( 1843 '%d-dimensional argument does not have enough ' 1844 'dimensions for all core dimensions %r' 1845 % (arg.ndim, core_dims)) 1846 1847 core_shape = arg.shape[-num_core_dims:] 1848 for dim, size in zip(core_dims, core_shape): 1849 if dim in dim_sizes: 1850 if size != dim_sizes[dim]: 1851 raise ValueError( 1852 'inconsistent size for core dimension %r: %r vs %r' 1853 % (dim, size, dim_sizes[dim])) 1854 else: 1855 dim_sizes[dim] = size 1856 1857 1858def _parse_input_dimensions(args, input_core_dims): 1859 """ 1860 Parse broadcast and core dimensions for vectorize with a signature. 1861 1862 Arguments 1863 --------- 1864 args : Tuple[ndarray, ...] 1865 Tuple of input arguments to examine. 1866 input_core_dims : List[Tuple[str, ...]] 1867 List of core dimensions corresponding to each input. 1868 1869 Returns 1870 ------- 1871 broadcast_shape : Tuple[int, ...] 1872 Common shape to broadcast all non-core dimensions to. 1873 dim_sizes : Dict[str, int] 1874 Common sizes for named core dimensions. 1875 """ 1876 broadcast_args = [] 1877 dim_sizes = {} 1878 for arg, core_dims in zip(args, input_core_dims): 1879 _update_dim_sizes(dim_sizes, arg, core_dims) 1880 ndim = arg.ndim - len(core_dims) 1881 dummy_array = np.lib.stride_tricks.as_strided(0, arg.shape[:ndim]) 1882 broadcast_args.append(dummy_array) 1883 broadcast_shape = np.lib.stride_tricks._broadcast_shape(*broadcast_args) 1884 return broadcast_shape, dim_sizes 1885 1886 1887def _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims): 1888 """Helper for calculating broadcast shapes with core dimensions.""" 1889 return [broadcast_shape + tuple(dim_sizes[dim] for dim in core_dims) 1890 for core_dims in list_of_core_dims] 1891 1892 1893def _create_arrays(broadcast_shape, dim_sizes, list_of_core_dims, dtypes): 1894 """Helper for creating output arrays in vectorize.""" 1895 shapes = _calculate_shapes(broadcast_shape, dim_sizes, list_of_core_dims) 1896 arrays = tuple(np.empty(shape, dtype=dtype) 1897 for shape, dtype in zip(shapes, dtypes)) 1898 return arrays 1899 1900 1901@set_module('numpy') 1902class vectorize: 1903 """ 1904 vectorize(pyfunc, otypes=None, doc=None, excluded=None, cache=False, 1905 signature=None) 1906 1907 Generalized function class. 1908 1909 Define a vectorized function which takes a nested sequence of objects or 1910 numpy arrays as inputs and returns a single numpy array or a tuple of numpy 1911 arrays. The vectorized function evaluates `pyfunc` over successive tuples 1912 of the input arrays like the python map function, except it uses the 1913 broadcasting rules of numpy. 1914 1915 The data type of the output of `vectorized` is determined by calling 1916 the function with the first element of the input. This can be avoided 1917 by specifying the `otypes` argument. 1918 1919 Parameters 1920 ---------- 1921 pyfunc : callable 1922 A python function or method. 1923 otypes : str or list of dtypes, optional 1924 The output data type. It must be specified as either a string of 1925 typecode characters or a list of data type specifiers. There should 1926 be one data type specifier for each output. 1927 doc : str, optional 1928 The docstring for the function. If None, the docstring will be the 1929 ``pyfunc.__doc__``. 1930 excluded : set, optional 1931 Set of strings or integers representing the positional or keyword 1932 arguments for which the function will not be vectorized. These will be 1933 passed directly to `pyfunc` unmodified. 1934 1935 .. versionadded:: 1.7.0 1936 1937 cache : bool, optional 1938 If `True`, then cache the first function call that determines the number 1939 of outputs if `otypes` is not provided. 1940 1941 .. versionadded:: 1.7.0 1942 1943 signature : string, optional 1944 Generalized universal function signature, e.g., ``(m,n),(n)->(m)`` for 1945 vectorized matrix-vector multiplication. If provided, ``pyfunc`` will 1946 be called with (and expected to return) arrays with shapes given by the 1947 size of corresponding core dimensions. By default, ``pyfunc`` is 1948 assumed to take scalars as input and output. 1949 1950 .. versionadded:: 1.12.0 1951 1952 Returns 1953 ------- 1954 vectorized : callable 1955 Vectorized function. 1956 1957 See Also 1958 -------- 1959 frompyfunc : Takes an arbitrary Python function and returns a ufunc 1960 1961 Notes 1962 ----- 1963 The `vectorize` function is provided primarily for convenience, not for 1964 performance. The implementation is essentially a for loop. 1965 1966 If `otypes` is not specified, then a call to the function with the 1967 first argument will be used to determine the number of outputs. The 1968 results of this call will be cached if `cache` is `True` to prevent 1969 calling the function twice. However, to implement the cache, the 1970 original function must be wrapped which will slow down subsequent 1971 calls, so only do this if your function is expensive. 1972 1973 The new keyword argument interface and `excluded` argument support 1974 further degrades performance. 1975 1976 References 1977 ---------- 1978 .. [1] :doc:`/reference/c-api/generalized-ufuncs` 1979 1980 Examples 1981 -------- 1982 >>> def myfunc(a, b): 1983 ... "Return a-b if a>b, otherwise return a+b" 1984 ... if a > b: 1985 ... return a - b 1986 ... else: 1987 ... return a + b 1988 1989 >>> vfunc = np.vectorize(myfunc) 1990 >>> vfunc([1, 2, 3, 4], 2) 1991 array([3, 4, 1, 2]) 1992 1993 The docstring is taken from the input function to `vectorize` unless it 1994 is specified: 1995 1996 >>> vfunc.__doc__ 1997 'Return a-b if a>b, otherwise return a+b' 1998 >>> vfunc = np.vectorize(myfunc, doc='Vectorized `myfunc`') 1999 >>> vfunc.__doc__ 2000 'Vectorized `myfunc`' 2001 2002 The output type is determined by evaluating the first element of the input, 2003 unless it is specified: 2004 2005 >>> out = vfunc([1, 2, 3, 4], 2) 2006 >>> type(out[0]) 2007 <class 'numpy.int64'> 2008 >>> vfunc = np.vectorize(myfunc, otypes=[float]) 2009 >>> out = vfunc([1, 2, 3, 4], 2) 2010 >>> type(out[0]) 2011 <class 'numpy.float64'> 2012 2013 The `excluded` argument can be used to prevent vectorizing over certain 2014 arguments. This can be useful for array-like arguments of a fixed length 2015 such as the coefficients for a polynomial as in `polyval`: 2016 2017 >>> def mypolyval(p, x): 2018 ... _p = list(p) 2019 ... res = _p.pop(0) 2020 ... while _p: 2021 ... res = res*x + _p.pop(0) 2022 ... return res 2023 >>> vpolyval = np.vectorize(mypolyval, excluded=['p']) 2024 >>> vpolyval(p=[1, 2, 3], x=[0, 1]) 2025 array([3, 6]) 2026 2027 Positional arguments may also be excluded by specifying their position: 2028 2029 >>> vpolyval.excluded.add(0) 2030 >>> vpolyval([1, 2, 3], x=[0, 1]) 2031 array([3, 6]) 2032 2033 The `signature` argument allows for vectorizing functions that act on 2034 non-scalar arrays of fixed length. For example, you can use it for a 2035 vectorized calculation of Pearson correlation coefficient and its p-value: 2036 2037 >>> import scipy.stats 2038 >>> pearsonr = np.vectorize(scipy.stats.pearsonr, 2039 ... signature='(n),(n)->(),()') 2040 >>> pearsonr([[0, 1, 2, 3]], [[1, 2, 3, 4], [4, 3, 2, 1]]) 2041 (array([ 1., -1.]), array([ 0., 0.])) 2042 2043 Or for a vectorized convolution: 2044 2045 >>> convolve = np.vectorize(np.convolve, signature='(n),(m)->(k)') 2046 >>> convolve(np.eye(4), [1, 2, 1]) 2047 array([[1., 2., 1., 0., 0., 0.], 2048 [0., 1., 2., 1., 0., 0.], 2049 [0., 0., 1., 2., 1., 0.], 2050 [0., 0., 0., 1., 2., 1.]]) 2051 2052 """ 2053 def __init__(self, pyfunc, otypes=None, doc=None, excluded=None, 2054 cache=False, signature=None): 2055 self.pyfunc = pyfunc 2056 self.cache = cache 2057 self.signature = signature 2058 self._ufunc = {} # Caching to improve default performance 2059 2060 if doc is None: 2061 self.__doc__ = pyfunc.__doc__ 2062 else: 2063 self.__doc__ = doc 2064 2065 if isinstance(otypes, str): 2066 for char in otypes: 2067 if char not in typecodes['All']: 2068 raise ValueError("Invalid otype specified: %s" % (char,)) 2069 elif iterable(otypes): 2070 otypes = ''.join([_nx.dtype(x).char for x in otypes]) 2071 elif otypes is not None: 2072 raise ValueError("Invalid otype specification") 2073 self.otypes = otypes 2074 2075 # Excluded variable support 2076 if excluded is None: 2077 excluded = set() 2078 self.excluded = set(excluded) 2079 2080 if signature is not None: 2081 self._in_and_out_core_dims = _parse_gufunc_signature(signature) 2082 else: 2083 self._in_and_out_core_dims = None 2084 2085 def __call__(self, *args, **kwargs): 2086 """ 2087 Return arrays with the results of `pyfunc` broadcast (vectorized) over 2088 `args` and `kwargs` not in `excluded`. 2089 """ 2090 excluded = self.excluded 2091 if not kwargs and not excluded: 2092 func = self.pyfunc 2093 vargs = args 2094 else: 2095 # The wrapper accepts only positional arguments: we use `names` and 2096 # `inds` to mutate `the_args` and `kwargs` to pass to the original 2097 # function. 2098 nargs = len(args) 2099 2100 names = [_n for _n in kwargs if _n not in excluded] 2101 inds = [_i for _i in range(nargs) if _i not in excluded] 2102 the_args = list(args) 2103 2104 def func(*vargs): 2105 for _n, _i in enumerate(inds): 2106 the_args[_i] = vargs[_n] 2107 kwargs.update(zip(names, vargs[len(inds):])) 2108 return self.pyfunc(*the_args, **kwargs) 2109 2110 vargs = [args[_i] for _i in inds] 2111 vargs.extend([kwargs[_n] for _n in names]) 2112 2113 return self._vectorize_call(func=func, args=vargs) 2114 2115 def _get_ufunc_and_otypes(self, func, args): 2116 """Return (ufunc, otypes).""" 2117 # frompyfunc will fail if args is empty 2118 if not args: 2119 raise ValueError('args can not be empty') 2120 2121 if self.otypes is not None: 2122 otypes = self.otypes 2123 2124 # self._ufunc is a dictionary whose keys are the number of 2125 # arguments (i.e. len(args)) and whose values are ufuncs created 2126 # by frompyfunc. len(args) can be different for different calls if 2127 # self.pyfunc has parameters with default values. We only use the 2128 # cache when func is self.pyfunc, which occurs when the call uses 2129 # only positional arguments and no arguments are excluded. 2130 2131 nin = len(args) 2132 nout = len(self.otypes) 2133 if func is not self.pyfunc or nin not in self._ufunc: 2134 ufunc = frompyfunc(func, nin, nout) 2135 else: 2136 ufunc = None # We'll get it from self._ufunc 2137 if func is self.pyfunc: 2138 ufunc = self._ufunc.setdefault(nin, ufunc) 2139 else: 2140 # Get number of outputs and output types by calling the function on 2141 # the first entries of args. We also cache the result to prevent 2142 # the subsequent call when the ufunc is evaluated. 2143 # Assumes that ufunc first evaluates the 0th elements in the input 2144 # arrays (the input values are not checked to ensure this) 2145 args = [asarray(arg) for arg in args] 2146 if builtins.any(arg.size == 0 for arg in args): 2147 raise ValueError('cannot call `vectorize` on size 0 inputs ' 2148 'unless `otypes` is set') 2149 2150 inputs = [arg.flat[0] for arg in args] 2151 outputs = func(*inputs) 2152 2153 # Performance note: profiling indicates that -- for simple 2154 # functions at least -- this wrapping can almost double the 2155 # execution time. 2156 # Hence we make it optional. 2157 if self.cache: 2158 _cache = [outputs] 2159 2160 def _func(*vargs): 2161 if _cache: 2162 return _cache.pop() 2163 else: 2164 return func(*vargs) 2165 else: 2166 _func = func 2167 2168 if isinstance(outputs, tuple): 2169 nout = len(outputs) 2170 else: 2171 nout = 1 2172 outputs = (outputs,) 2173 2174 otypes = ''.join([asarray(outputs[_k]).dtype.char 2175 for _k in range(nout)]) 2176 2177 # Performance note: profiling indicates that creating the ufunc is 2178 # not a significant cost compared with wrapping so it seems not 2179 # worth trying to cache this. 2180 ufunc = frompyfunc(_func, len(args), nout) 2181 2182 return ufunc, otypes 2183 2184 def _vectorize_call(self, func, args): 2185 """Vectorized call to `func` over positional `args`.""" 2186 if self.signature is not None: 2187 res = self._vectorize_call_with_signature(func, args) 2188 elif not args: 2189 res = func() 2190 else: 2191 ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args) 2192 2193 # Convert args to object arrays first 2194 inputs = [array(a, copy=False, subok=True, dtype=object) 2195 for a in args] 2196 2197 outputs = ufunc(*inputs) 2198 2199 if ufunc.nout == 1: 2200 res = array(outputs, copy=False, subok=True, dtype=otypes[0]) 2201 else: 2202 res = tuple([array(x, copy=False, subok=True, dtype=t) 2203 for x, t in zip(outputs, otypes)]) 2204 return res 2205 2206 def _vectorize_call_with_signature(self, func, args): 2207 """Vectorized call over positional arguments with a signature.""" 2208 input_core_dims, output_core_dims = self._in_and_out_core_dims 2209 2210 if len(args) != len(input_core_dims): 2211 raise TypeError('wrong number of positional arguments: ' 2212 'expected %r, got %r' 2213 % (len(input_core_dims), len(args))) 2214 args = tuple(asanyarray(arg) for arg in args) 2215 2216 broadcast_shape, dim_sizes = _parse_input_dimensions( 2217 args, input_core_dims) 2218 input_shapes = _calculate_shapes(broadcast_shape, dim_sizes, 2219 input_core_dims) 2220 args = [np.broadcast_to(arg, shape, subok=True) 2221 for arg, shape in zip(args, input_shapes)] 2222 2223 outputs = None 2224 otypes = self.otypes 2225 nout = len(output_core_dims) 2226 2227 for index in np.ndindex(*broadcast_shape): 2228 results = func(*(arg[index] for arg in args)) 2229 2230 n_results = len(results) if isinstance(results, tuple) else 1 2231 2232 if nout != n_results: 2233 raise ValueError( 2234 'wrong number of outputs from pyfunc: expected %r, got %r' 2235 % (nout, n_results)) 2236 2237 if nout == 1: 2238 results = (results,) 2239 2240 if outputs is None: 2241 for result, core_dims in zip(results, output_core_dims): 2242 _update_dim_sizes(dim_sizes, result, core_dims) 2243 2244 if otypes is None: 2245 otypes = [asarray(result).dtype for result in results] 2246 2247 outputs = _create_arrays(broadcast_shape, dim_sizes, 2248 output_core_dims, otypes) 2249 2250 for output, result in zip(outputs, results): 2251 output[index] = result 2252 2253 if outputs is None: 2254 # did not call the function even once 2255 if otypes is None: 2256 raise ValueError('cannot call `vectorize` on size 0 inputs ' 2257 'unless `otypes` is set') 2258 if builtins.any(dim not in dim_sizes 2259 for dims in output_core_dims 2260 for dim in dims): 2261 raise ValueError('cannot call `vectorize` with a signature ' 2262 'including new output dimensions on size 0 ' 2263 'inputs') 2264 outputs = _create_arrays(broadcast_shape, dim_sizes, 2265 output_core_dims, otypes) 2266 2267 return outputs[0] if nout == 1 else outputs 2268 2269 2270def _cov_dispatcher(m, y=None, rowvar=None, bias=None, ddof=None, 2271 fweights=None, aweights=None, *, dtype=None): 2272 return (m, y, fweights, aweights) 2273 2274 2275@array_function_dispatch(_cov_dispatcher) 2276def cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None, 2277 aweights=None, *, dtype=None): 2278 """ 2279 Estimate a covariance matrix, given data and weights. 2280 2281 Covariance indicates the level to which two variables vary together. 2282 If we examine N-dimensional samples, :math:`X = [x_1, x_2, ... x_N]^T`, 2283 then the covariance matrix element :math:`C_{ij}` is the covariance of 2284 :math:`x_i` and :math:`x_j`. The element :math:`C_{ii}` is the variance 2285 of :math:`x_i`. 2286 2287 See the notes for an outline of the algorithm. 2288 2289 Parameters 2290 ---------- 2291 m : array_like 2292 A 1-D or 2-D array containing multiple variables and observations. 2293 Each row of `m` represents a variable, and each column a single 2294 observation of all those variables. Also see `rowvar` below. 2295 y : array_like, optional 2296 An additional set of variables and observations. `y` has the same form 2297 as that of `m`. 2298 rowvar : bool, optional 2299 If `rowvar` is True (default), then each row represents a 2300 variable, with observations in the columns. Otherwise, the relationship 2301 is transposed: each column represents a variable, while the rows 2302 contain observations. 2303 bias : bool, optional 2304 Default normalization (False) is by ``(N - 1)``, where ``N`` is the 2305 number of observations given (unbiased estimate). If `bias` is True, 2306 then normalization is by ``N``. These values can be overridden by using 2307 the keyword ``ddof`` in numpy versions >= 1.5. 2308 ddof : int, optional 2309 If not ``None`` the default value implied by `bias` is overridden. 2310 Note that ``ddof=1`` will return the unbiased estimate, even if both 2311 `fweights` and `aweights` are specified, and ``ddof=0`` will return 2312 the simple average. See the notes for the details. The default value 2313 is ``None``. 2314 2315 .. versionadded:: 1.5 2316 fweights : array_like, int, optional 2317 1-D array of integer frequency weights; the number of times each 2318 observation vector should be repeated. 2319 2320 .. versionadded:: 1.10 2321 aweights : array_like, optional 2322 1-D array of observation vector weights. These relative weights are 2323 typically large for observations considered "important" and smaller for 2324 observations considered less "important". If ``ddof=0`` the array of 2325 weights can be used to assign probabilities to observation vectors. 2326 2327 .. versionadded:: 1.10 2328 dtype : data-type, optional 2329 Data-type of the result. By default, the return data-type will have 2330 at least `numpy.float64` precision. 2331 2332 .. versionadded:: 1.20 2333 2334 Returns 2335 ------- 2336 out : ndarray 2337 The covariance matrix of the variables. 2338 2339 See Also 2340 -------- 2341 corrcoef : Normalized covariance matrix 2342 2343 Notes 2344 ----- 2345 Assume that the observations are in the columns of the observation 2346 array `m` and let ``f = fweights`` and ``a = aweights`` for brevity. The 2347 steps to compute the weighted covariance are as follows:: 2348 2349 >>> m = np.arange(10, dtype=np.float64) 2350 >>> f = np.arange(10) * 2 2351 >>> a = np.arange(10) ** 2. 2352 >>> ddof = 1 2353 >>> w = f * a 2354 >>> v1 = np.sum(w) 2355 >>> v2 = np.sum(w * a) 2356 >>> m -= np.sum(m * w, axis=None, keepdims=True) / v1 2357 >>> cov = np.dot(m * w, m.T) * v1 / (v1**2 - ddof * v2) 2358 2359 Note that when ``a == 1``, the normalization factor 2360 ``v1 / (v1**2 - ddof * v2)`` goes over to ``1 / (np.sum(f) - ddof)`` 2361 as it should. 2362 2363 Examples 2364 -------- 2365 Consider two variables, :math:`x_0` and :math:`x_1`, which 2366 correlate perfectly, but in opposite directions: 2367 2368 >>> x = np.array([[0, 2], [1, 1], [2, 0]]).T 2369 >>> x 2370 array([[0, 1, 2], 2371 [2, 1, 0]]) 2372 2373 Note how :math:`x_0` increases while :math:`x_1` decreases. The covariance 2374 matrix shows this clearly: 2375 2376 >>> np.cov(x) 2377 array([[ 1., -1.], 2378 [-1., 1.]]) 2379 2380 Note that element :math:`C_{0,1}`, which shows the correlation between 2381 :math:`x_0` and :math:`x_1`, is negative. 2382 2383 Further, note how `x` and `y` are combined: 2384 2385 >>> x = [-2.1, -1, 4.3] 2386 >>> y = [3, 1.1, 0.12] 2387 >>> X = np.stack((x, y), axis=0) 2388 >>> np.cov(X) 2389 array([[11.71 , -4.286 ], # may vary 2390 [-4.286 , 2.144133]]) 2391 >>> np.cov(x, y) 2392 array([[11.71 , -4.286 ], # may vary 2393 [-4.286 , 2.144133]]) 2394 >>> np.cov(x) 2395 array(11.71) 2396 2397 """ 2398 # Check inputs 2399 if ddof is not None and ddof != int(ddof): 2400 raise ValueError( 2401 "ddof must be integer") 2402 2403 # Handles complex arrays too 2404 m = np.asarray(m) 2405 if m.ndim > 2: 2406 raise ValueError("m has more than 2 dimensions") 2407 2408 if y is not None: 2409 y = np.asarray(y) 2410 if y.ndim > 2: 2411 raise ValueError("y has more than 2 dimensions") 2412 2413 if dtype is None: 2414 if y is None: 2415 dtype = np.result_type(m, np.float64) 2416 else: 2417 dtype = np.result_type(m, y, np.float64) 2418 2419 X = array(m, ndmin=2, dtype=dtype) 2420 if not rowvar and X.shape[0] != 1: 2421 X = X.T 2422 if X.shape[0] == 0: 2423 return np.array([]).reshape(0, 0) 2424 if y is not None: 2425 y = array(y, copy=False, ndmin=2, dtype=dtype) 2426 if not rowvar and y.shape[0] != 1: 2427 y = y.T 2428 X = np.concatenate((X, y), axis=0) 2429 2430 if ddof is None: 2431 if bias == 0: 2432 ddof = 1 2433 else: 2434 ddof = 0 2435 2436 # Get the product of frequencies and weights 2437 w = None 2438 if fweights is not None: 2439 fweights = np.asarray(fweights, dtype=float) 2440 if not np.all(fweights == np.around(fweights)): 2441 raise TypeError( 2442 "fweights must be integer") 2443 if fweights.ndim > 1: 2444 raise RuntimeError( 2445 "cannot handle multidimensional fweights") 2446 if fweights.shape[0] != X.shape[1]: 2447 raise RuntimeError( 2448 "incompatible numbers of samples and fweights") 2449 if any(fweights < 0): 2450 raise ValueError( 2451 "fweights cannot be negative") 2452 w = fweights 2453 if aweights is not None: 2454 aweights = np.asarray(aweights, dtype=float) 2455 if aweights.ndim > 1: 2456 raise RuntimeError( 2457 "cannot handle multidimensional aweights") 2458 if aweights.shape[0] != X.shape[1]: 2459 raise RuntimeError( 2460 "incompatible numbers of samples and aweights") 2461 if any(aweights < 0): 2462 raise ValueError( 2463 "aweights cannot be negative") 2464 if w is None: 2465 w = aweights 2466 else: 2467 w *= aweights 2468 2469 avg, w_sum = average(X, axis=1, weights=w, returned=True) 2470 w_sum = w_sum[0] 2471 2472 # Determine the normalization 2473 if w is None: 2474 fact = X.shape[1] - ddof 2475 elif ddof == 0: 2476 fact = w_sum 2477 elif aweights is None: 2478 fact = w_sum - ddof 2479 else: 2480 fact = w_sum - ddof*sum(w*aweights)/w_sum 2481 2482 if fact <= 0: 2483 warnings.warn("Degrees of freedom <= 0 for slice", 2484 RuntimeWarning, stacklevel=3) 2485 fact = 0.0 2486 2487 X -= avg[:, None] 2488 if w is None: 2489 X_T = X.T 2490 else: 2491 X_T = (X*w).T 2492 c = dot(X, X_T.conj()) 2493 c *= np.true_divide(1, fact) 2494 return c.squeeze() 2495 2496 2497def _corrcoef_dispatcher(x, y=None, rowvar=None, bias=None, ddof=None, *, 2498 dtype=None): 2499 return (x, y) 2500 2501 2502@array_function_dispatch(_corrcoef_dispatcher) 2503def corrcoef(x, y=None, rowvar=True, bias=np._NoValue, ddof=np._NoValue, *, 2504 dtype=None): 2505 """ 2506 Return Pearson product-moment correlation coefficients. 2507 2508 Please refer to the documentation for `cov` for more detail. The 2509 relationship between the correlation coefficient matrix, `R`, and the 2510 covariance matrix, `C`, is 2511 2512 .. math:: R_{ij} = \\frac{ C_{ij} } { \\sqrt{ C_{ii} * C_{jj} } } 2513 2514 The values of `R` are between -1 and 1, inclusive. 2515 2516 Parameters 2517 ---------- 2518 x : array_like 2519 A 1-D or 2-D array containing multiple variables and observations. 2520 Each row of `x` represents a variable, and each column a single 2521 observation of all those variables. Also see `rowvar` below. 2522 y : array_like, optional 2523 An additional set of variables and observations. `y` has the same 2524 shape as `x`. 2525 rowvar : bool, optional 2526 If `rowvar` is True (default), then each row represents a 2527 variable, with observations in the columns. Otherwise, the relationship 2528 is transposed: each column represents a variable, while the rows 2529 contain observations. 2530 bias : _NoValue, optional 2531 Has no effect, do not use. 2532 2533 .. deprecated:: 1.10.0 2534 ddof : _NoValue, optional 2535 Has no effect, do not use. 2536 2537 .. deprecated:: 1.10.0 2538 dtype : data-type, optional 2539 Data-type of the result. By default, the return data-type will have 2540 at least `numpy.float64` precision. 2541 2542 .. versionadded:: 1.20 2543 2544 Returns 2545 ------- 2546 R : ndarray 2547 The correlation coefficient matrix of the variables. 2548 2549 See Also 2550 -------- 2551 cov : Covariance matrix 2552 2553 Notes 2554 ----- 2555 Due to floating point rounding the resulting array may not be Hermitian, 2556 the diagonal elements may not be 1, and the elements may not satisfy the 2557 inequality abs(a) <= 1. The real and imaginary parts are clipped to the 2558 interval [-1, 1] in an attempt to improve on that situation but is not 2559 much help in the complex case. 2560 2561 This function accepts but discards arguments `bias` and `ddof`. This is 2562 for backwards compatibility with previous versions of this function. These 2563 arguments had no effect on the return values of the function and can be 2564 safely ignored in this and previous versions of numpy. 2565 2566 Examples 2567 -------- 2568 In this example we generate two random arrays, ``xarr`` and ``yarr``, and 2569 compute the row-wise and column-wise Pearson correlation coefficients, 2570 ``R``. Since ``rowvar`` is true by default, we first find the row-wise 2571 Pearson correlation coefficients between the variables of ``xarr``. 2572 2573 >>> import numpy as np 2574 >>> rng = np.random.default_rng(seed=42) 2575 >>> xarr = rng.random((3, 3)) 2576 >>> xarr 2577 array([[0.77395605, 0.43887844, 0.85859792], 2578 [0.69736803, 0.09417735, 0.97562235], 2579 [0.7611397 , 0.78606431, 0.12811363]]) 2580 >>> R1 = np.corrcoef(xarr) 2581 >>> R1 2582 array([[ 1. , 0.99256089, -0.68080986], 2583 [ 0.99256089, 1. , -0.76492172], 2584 [-0.68080986, -0.76492172, 1. ]]) 2585 2586 If we add another set of variables and observations ``yarr``, we can 2587 compute the row-wise Pearson correlation coefficients between the 2588 variables in ``xarr`` and ``yarr``. 2589 2590 >>> yarr = rng.random((3, 3)) 2591 >>> yarr 2592 array([[0.45038594, 0.37079802, 0.92676499], 2593 [0.64386512, 0.82276161, 0.4434142 ], 2594 [0.22723872, 0.55458479, 0.06381726]]) 2595 >>> R2 = np.corrcoef(xarr, yarr) 2596 >>> R2 2597 array([[ 1. , 0.99256089, -0.68080986, 0.75008178, -0.934284 , 2598 -0.99004057], 2599 [ 0.99256089, 1. , -0.76492172, 0.82502011, -0.97074098, 2600 -0.99981569], 2601 [-0.68080986, -0.76492172, 1. , -0.99507202, 0.89721355, 2602 0.77714685], 2603 [ 0.75008178, 0.82502011, -0.99507202, 1. , -0.93657855, 2604 -0.83571711], 2605 [-0.934284 , -0.97074098, 0.89721355, -0.93657855, 1. , 2606 0.97517215], 2607 [-0.99004057, -0.99981569, 0.77714685, -0.83571711, 0.97517215, 2608 1. ]]) 2609 2610 Finally if we use the option ``rowvar=False``, the columns are now 2611 being treated as the variables and we will find the column-wise Pearson 2612 correlation coefficients between variables in ``xarr`` and ``yarr``. 2613 2614 >>> R3 = np.corrcoef(xarr, yarr, rowvar=False) 2615 >>> R3 2616 array([[ 1. , 0.77598074, -0.47458546, -0.75078643, -0.9665554 , 2617 0.22423734], 2618 [ 0.77598074, 1. , -0.92346708, -0.99923895, -0.58826587, 2619 -0.44069024], 2620 [-0.47458546, -0.92346708, 1. , 0.93773029, 0.23297648, 2621 0.75137473], 2622 [-0.75078643, -0.99923895, 0.93773029, 1. , 0.55627469, 2623 0.47536961], 2624 [-0.9665554 , -0.58826587, 0.23297648, 0.55627469, 1. , 2625 -0.46666491], 2626 [ 0.22423734, -0.44069024, 0.75137473, 0.47536961, -0.46666491, 2627 1. ]]) 2628 2629 """ 2630 if bias is not np._NoValue or ddof is not np._NoValue: 2631 # 2015-03-15, 1.10 2632 warnings.warn('bias and ddof have no effect and are deprecated', 2633 DeprecationWarning, stacklevel=3) 2634 c = cov(x, y, rowvar, dtype=dtype) 2635 try: 2636 d = diag(c) 2637 except ValueError: 2638 # scalar covariance 2639 # nan if incorrect value (nan, inf, 0), 1 otherwise 2640 return c / c 2641 stddev = sqrt(d.real) 2642 c /= stddev[:, None] 2643 c /= stddev[None, :] 2644 2645 # Clip real and imaginary parts to [-1, 1]. This does not guarantee 2646 # abs(a[i,j]) <= 1 for complex arrays, but is the best we can do without 2647 # excessive work. 2648 np.clip(c.real, -1, 1, out=c.real) 2649 if np.iscomplexobj(c): 2650 np.clip(c.imag, -1, 1, out=c.imag) 2651 2652 return c 2653 2654 2655@set_module('numpy') 2656def blackman(M): 2657 """ 2658 Return the Blackman window. 2659 2660 The Blackman window is a taper formed by using the first three 2661 terms of a summation of cosines. It was designed to have close to the 2662 minimal leakage possible. It is close to optimal, only slightly worse 2663 than a Kaiser window. 2664 2665 Parameters 2666 ---------- 2667 M : int 2668 Number of points in the output window. If zero or less, an empty 2669 array is returned. 2670 2671 Returns 2672 ------- 2673 out : ndarray 2674 The window, with the maximum value normalized to one (the value one 2675 appears only if the number of samples is odd). 2676 2677 See Also 2678 -------- 2679 bartlett, hamming, hanning, kaiser 2680 2681 Notes 2682 ----- 2683 The Blackman window is defined as 2684 2685 .. math:: w(n) = 0.42 - 0.5 \\cos(2\\pi n/M) + 0.08 \\cos(4\\pi n/M) 2686 2687 Most references to the Blackman window come from the signal processing 2688 literature, where it is used as one of many windowing functions for 2689 smoothing values. It is also known as an apodization (which means 2690 "removing the foot", i.e. smoothing discontinuities at the beginning 2691 and end of the sampled signal) or tapering function. It is known as a 2692 "near optimal" tapering function, almost as good (by some measures) 2693 as the kaiser window. 2694 2695 References 2696 ---------- 2697 Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra, 2698 Dover Publications, New York. 2699 2700 Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing. 2701 Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471. 2702 2703 Examples 2704 -------- 2705 >>> import matplotlib.pyplot as plt 2706 >>> np.blackman(12) 2707 array([-1.38777878e-17, 3.26064346e-02, 1.59903635e-01, # may vary 2708 4.14397981e-01, 7.36045180e-01, 9.67046769e-01, 2709 9.67046769e-01, 7.36045180e-01, 4.14397981e-01, 2710 1.59903635e-01, 3.26064346e-02, -1.38777878e-17]) 2711 2712 Plot the window and the frequency response: 2713 2714 >>> from numpy.fft import fft, fftshift 2715 >>> window = np.blackman(51) 2716 >>> plt.plot(window) 2717 [<matplotlib.lines.Line2D object at 0x...>] 2718 >>> plt.title("Blackman window") 2719 Text(0.5, 1.0, 'Blackman window') 2720 >>> plt.ylabel("Amplitude") 2721 Text(0, 0.5, 'Amplitude') 2722 >>> plt.xlabel("Sample") 2723 Text(0.5, 0, 'Sample') 2724 >>> plt.show() 2725 2726 >>> plt.figure() 2727 <Figure size 640x480 with 0 Axes> 2728 >>> A = fft(window, 2048) / 25.5 2729 >>> mag = np.abs(fftshift(A)) 2730 >>> freq = np.linspace(-0.5, 0.5, len(A)) 2731 >>> with np.errstate(divide='ignore', invalid='ignore'): 2732 ... response = 20 * np.log10(mag) 2733 ... 2734 >>> response = np.clip(response, -100, 100) 2735 >>> plt.plot(freq, response) 2736 [<matplotlib.lines.Line2D object at 0x...>] 2737 >>> plt.title("Frequency response of Blackman window") 2738 Text(0.5, 1.0, 'Frequency response of Blackman window') 2739 >>> plt.ylabel("Magnitude [dB]") 2740 Text(0, 0.5, 'Magnitude [dB]') 2741 >>> plt.xlabel("Normalized frequency [cycles per sample]") 2742 Text(0.5, 0, 'Normalized frequency [cycles per sample]') 2743 >>> _ = plt.axis('tight') 2744 >>> plt.show() 2745 2746 """ 2747 if M < 1: 2748 return array([]) 2749 if M == 1: 2750 return ones(1, float) 2751 n = arange(1-M, M, 2) 2752 return 0.42 + 0.5*cos(pi*n/(M-1)) + 0.08*cos(2.0*pi*n/(M-1)) 2753 2754 2755@set_module('numpy') 2756def bartlett(M): 2757 """ 2758 Return the Bartlett window. 2759 2760 The Bartlett window is very similar to a triangular window, except 2761 that the end points are at zero. It is often used in signal 2762 processing for tapering a signal, without generating too much 2763 ripple in the frequency domain. 2764 2765 Parameters 2766 ---------- 2767 M : int 2768 Number of points in the output window. If zero or less, an 2769 empty array is returned. 2770 2771 Returns 2772 ------- 2773 out : array 2774 The triangular window, with the maximum value normalized to one 2775 (the value one appears only if the number of samples is odd), with 2776 the first and last samples equal to zero. 2777 2778 See Also 2779 -------- 2780 blackman, hamming, hanning, kaiser 2781 2782 Notes 2783 ----- 2784 The Bartlett window is defined as 2785 2786 .. math:: w(n) = \\frac{2}{M-1} \\left( 2787 \\frac{M-1}{2} - \\left|n - \\frac{M-1}{2}\\right| 2788 \\right) 2789 2790 Most references to the Bartlett window come from the signal 2791 processing literature, where it is used as one of many windowing 2792 functions for smoothing values. Note that convolution with this 2793 window produces linear interpolation. It is also known as an 2794 apodization (which means"removing the foot", i.e. smoothing 2795 discontinuities at the beginning and end of the sampled signal) or 2796 tapering function. The fourier transform of the Bartlett is the product 2797 of two sinc functions. 2798 Note the excellent discussion in Kanasewich. 2799 2800 References 2801 ---------- 2802 .. [1] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra", 2803 Biometrika 37, 1-16, 1950. 2804 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", 2805 The University of Alberta Press, 1975, pp. 109-110. 2806 .. [3] A.V. Oppenheim and R.W. Schafer, "Discrete-Time Signal 2807 Processing", Prentice-Hall, 1999, pp. 468-471. 2808 .. [4] Wikipedia, "Window function", 2809 https://en.wikipedia.org/wiki/Window_function 2810 .. [5] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, 2811 "Numerical Recipes", Cambridge University Press, 1986, page 429. 2812 2813 Examples 2814 -------- 2815 >>> import matplotlib.pyplot as plt 2816 >>> np.bartlett(12) 2817 array([ 0. , 0.18181818, 0.36363636, 0.54545455, 0.72727273, # may vary 2818 0.90909091, 0.90909091, 0.72727273, 0.54545455, 0.36363636, 2819 0.18181818, 0. ]) 2820 2821 Plot the window and its frequency response (requires SciPy and matplotlib): 2822 2823 >>> from numpy.fft import fft, fftshift 2824 >>> window = np.bartlett(51) 2825 >>> plt.plot(window) 2826 [<matplotlib.lines.Line2D object at 0x...>] 2827 >>> plt.title("Bartlett window") 2828 Text(0.5, 1.0, 'Bartlett window') 2829 >>> plt.ylabel("Amplitude") 2830 Text(0, 0.5, 'Amplitude') 2831 >>> plt.xlabel("Sample") 2832 Text(0.5, 0, 'Sample') 2833 >>> plt.show() 2834 2835 >>> plt.figure() 2836 <Figure size 640x480 with 0 Axes> 2837 >>> A = fft(window, 2048) / 25.5 2838 >>> mag = np.abs(fftshift(A)) 2839 >>> freq = np.linspace(-0.5, 0.5, len(A)) 2840 >>> with np.errstate(divide='ignore', invalid='ignore'): 2841 ... response = 20 * np.log10(mag) 2842 ... 2843 >>> response = np.clip(response, -100, 100) 2844 >>> plt.plot(freq, response) 2845 [<matplotlib.lines.Line2D object at 0x...>] 2846 >>> plt.title("Frequency response of Bartlett window") 2847 Text(0.5, 1.0, 'Frequency response of Bartlett window') 2848 >>> plt.ylabel("Magnitude [dB]") 2849 Text(0, 0.5, 'Magnitude [dB]') 2850 >>> plt.xlabel("Normalized frequency [cycles per sample]") 2851 Text(0.5, 0, 'Normalized frequency [cycles per sample]') 2852 >>> _ = plt.axis('tight') 2853 >>> plt.show() 2854 2855 """ 2856 if M < 1: 2857 return array([]) 2858 if M == 1: 2859 return ones(1, float) 2860 n = arange(1-M, M, 2) 2861 return where(less_equal(n, 0), 1 + n/(M-1), 1 - n/(M-1)) 2862 2863 2864@set_module('numpy') 2865def hanning(M): 2866 """ 2867 Return the Hanning window. 2868 2869 The Hanning window is a taper formed by using a weighted cosine. 2870 2871 Parameters 2872 ---------- 2873 M : int 2874 Number of points in the output window. If zero or less, an 2875 empty array is returned. 2876 2877 Returns 2878 ------- 2879 out : ndarray, shape(M,) 2880 The window, with the maximum value normalized to one (the value 2881 one appears only if `M` is odd). 2882 2883 See Also 2884 -------- 2885 bartlett, blackman, hamming, kaiser 2886 2887 Notes 2888 ----- 2889 The Hanning window is defined as 2890 2891 .. math:: w(n) = 0.5 - 0.5cos\\left(\\frac{2\\pi{n}}{M-1}\\right) 2892 \\qquad 0 \\leq n \\leq M-1 2893 2894 The Hanning was named for Julius von Hann, an Austrian meteorologist. 2895 It is also known as the Cosine Bell. Some authors prefer that it be 2896 called a Hann window, to help avoid confusion with the very similar 2897 Hamming window. 2898 2899 Most references to the Hanning window come from the signal processing 2900 literature, where it is used as one of many windowing functions for 2901 smoothing values. It is also known as an apodization (which means 2902 "removing the foot", i.e. smoothing discontinuities at the beginning 2903 and end of the sampled signal) or tapering function. 2904 2905 References 2906 ---------- 2907 .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power 2908 spectra, Dover Publications, New York. 2909 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", 2910 The University of Alberta Press, 1975, pp. 106-108. 2911 .. [3] Wikipedia, "Window function", 2912 https://en.wikipedia.org/wiki/Window_function 2913 .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, 2914 "Numerical Recipes", Cambridge University Press, 1986, page 425. 2915 2916 Examples 2917 -------- 2918 >>> np.hanning(12) 2919 array([0. , 0.07937323, 0.29229249, 0.57115742, 0.82743037, 2920 0.97974649, 0.97974649, 0.82743037, 0.57115742, 0.29229249, 2921 0.07937323, 0. ]) 2922 2923 Plot the window and its frequency response: 2924 2925 >>> import matplotlib.pyplot as plt 2926 >>> from numpy.fft import fft, fftshift 2927 >>> window = np.hanning(51) 2928 >>> plt.plot(window) 2929 [<matplotlib.lines.Line2D object at 0x...>] 2930 >>> plt.title("Hann window") 2931 Text(0.5, 1.0, 'Hann window') 2932 >>> plt.ylabel("Amplitude") 2933 Text(0, 0.5, 'Amplitude') 2934 >>> plt.xlabel("Sample") 2935 Text(0.5, 0, 'Sample') 2936 >>> plt.show() 2937 2938 >>> plt.figure() 2939 <Figure size 640x480 with 0 Axes> 2940 >>> A = fft(window, 2048) / 25.5 2941 >>> mag = np.abs(fftshift(A)) 2942 >>> freq = np.linspace(-0.5, 0.5, len(A)) 2943 >>> with np.errstate(divide='ignore', invalid='ignore'): 2944 ... response = 20 * np.log10(mag) 2945 ... 2946 >>> response = np.clip(response, -100, 100) 2947 >>> plt.plot(freq, response) 2948 [<matplotlib.lines.Line2D object at 0x...>] 2949 >>> plt.title("Frequency response of the Hann window") 2950 Text(0.5, 1.0, 'Frequency response of the Hann window') 2951 >>> plt.ylabel("Magnitude [dB]") 2952 Text(0, 0.5, 'Magnitude [dB]') 2953 >>> plt.xlabel("Normalized frequency [cycles per sample]") 2954 Text(0.5, 0, 'Normalized frequency [cycles per sample]') 2955 >>> plt.axis('tight') 2956 ... 2957 >>> plt.show() 2958 2959 """ 2960 if M < 1: 2961 return array([]) 2962 if M == 1: 2963 return ones(1, float) 2964 n = arange(1-M, M, 2) 2965 return 0.5 + 0.5*cos(pi*n/(M-1)) 2966 2967 2968@set_module('numpy') 2969def hamming(M): 2970 """ 2971 Return the Hamming window. 2972 2973 The Hamming window is a taper formed by using a weighted cosine. 2974 2975 Parameters 2976 ---------- 2977 M : int 2978 Number of points in the output window. If zero or less, an 2979 empty array is returned. 2980 2981 Returns 2982 ------- 2983 out : ndarray 2984 The window, with the maximum value normalized to one (the value 2985 one appears only if the number of samples is odd). 2986 2987 See Also 2988 -------- 2989 bartlett, blackman, hanning, kaiser 2990 2991 Notes 2992 ----- 2993 The Hamming window is defined as 2994 2995 .. math:: w(n) = 0.54 - 0.46cos\\left(\\frac{2\\pi{n}}{M-1}\\right) 2996 \\qquad 0 \\leq n \\leq M-1 2997 2998 The Hamming was named for R. W. Hamming, an associate of J. W. Tukey 2999 and is described in Blackman and Tukey. It was recommended for 3000 smoothing the truncated autocovariance function in the time domain. 3001 Most references to the Hamming window come from the signal processing 3002 literature, where it is used as one of many windowing functions for 3003 smoothing values. It is also known as an apodization (which means 3004 "removing the foot", i.e. smoothing discontinuities at the beginning 3005 and end of the sampled signal) or tapering function. 3006 3007 References 3008 ---------- 3009 .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power 3010 spectra, Dover Publications, New York. 3011 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The 3012 University of Alberta Press, 1975, pp. 109-110. 3013 .. [3] Wikipedia, "Window function", 3014 https://en.wikipedia.org/wiki/Window_function 3015 .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, 3016 "Numerical Recipes", Cambridge University Press, 1986, page 425. 3017 3018 Examples 3019 -------- 3020 >>> np.hamming(12) 3021 array([ 0.08 , 0.15302337, 0.34890909, 0.60546483, 0.84123594, # may vary 3022 0.98136677, 0.98136677, 0.84123594, 0.60546483, 0.34890909, 3023 0.15302337, 0.08 ]) 3024 3025 Plot the window and the frequency response: 3026 3027 >>> import matplotlib.pyplot as plt 3028 >>> from numpy.fft import fft, fftshift 3029 >>> window = np.hamming(51) 3030 >>> plt.plot(window) 3031 [<matplotlib.lines.Line2D object at 0x...>] 3032 >>> plt.title("Hamming window") 3033 Text(0.5, 1.0, 'Hamming window') 3034 >>> plt.ylabel("Amplitude") 3035 Text(0, 0.5, 'Amplitude') 3036 >>> plt.xlabel("Sample") 3037 Text(0.5, 0, 'Sample') 3038 >>> plt.show() 3039 3040 >>> plt.figure() 3041 <Figure size 640x480 with 0 Axes> 3042 >>> A = fft(window, 2048) / 25.5 3043 >>> mag = np.abs(fftshift(A)) 3044 >>> freq = np.linspace(-0.5, 0.5, len(A)) 3045 >>> response = 20 * np.log10(mag) 3046 >>> response = np.clip(response, -100, 100) 3047 >>> plt.plot(freq, response) 3048 [<matplotlib.lines.Line2D object at 0x...>] 3049 >>> plt.title("Frequency response of Hamming window") 3050 Text(0.5, 1.0, 'Frequency response of Hamming window') 3051 >>> plt.ylabel("Magnitude [dB]") 3052 Text(0, 0.5, 'Magnitude [dB]') 3053 >>> plt.xlabel("Normalized frequency [cycles per sample]") 3054 Text(0.5, 0, 'Normalized frequency [cycles per sample]') 3055 >>> plt.axis('tight') 3056 ... 3057 >>> plt.show() 3058 3059 """ 3060 if M < 1: 3061 return array([]) 3062 if M == 1: 3063 return ones(1, float) 3064 n = arange(1-M, M, 2) 3065 return 0.54 + 0.46*cos(pi*n/(M-1)) 3066 3067 3068## Code from cephes for i0 3069 3070_i0A = [ 3071 -4.41534164647933937950E-18, 3072 3.33079451882223809783E-17, 3073 -2.43127984654795469359E-16, 3074 1.71539128555513303061E-15, 3075 -1.16853328779934516808E-14, 3076 7.67618549860493561688E-14, 3077 -4.85644678311192946090E-13, 3078 2.95505266312963983461E-12, 3079 -1.72682629144155570723E-11, 3080 9.67580903537323691224E-11, 3081 -5.18979560163526290666E-10, 3082 2.65982372468238665035E-9, 3083 -1.30002500998624804212E-8, 3084 6.04699502254191894932E-8, 3085 -2.67079385394061173391E-7, 3086 1.11738753912010371815E-6, 3087 -4.41673835845875056359E-6, 3088 1.64484480707288970893E-5, 3089 -5.75419501008210370398E-5, 3090 1.88502885095841655729E-4, 3091 -5.76375574538582365885E-4, 3092 1.63947561694133579842E-3, 3093 -4.32430999505057594430E-3, 3094 1.05464603945949983183E-2, 3095 -2.37374148058994688156E-2, 3096 4.93052842396707084878E-2, 3097 -9.49010970480476444210E-2, 3098 1.71620901522208775349E-1, 3099 -3.04682672343198398683E-1, 3100 6.76795274409476084995E-1 3101 ] 3102 3103_i0B = [ 3104 -7.23318048787475395456E-18, 3105 -4.83050448594418207126E-18, 3106 4.46562142029675999901E-17, 3107 3.46122286769746109310E-17, 3108 -2.82762398051658348494E-16, 3109 -3.42548561967721913462E-16, 3110 1.77256013305652638360E-15, 3111 3.81168066935262242075E-15, 3112 -9.55484669882830764870E-15, 3113 -4.15056934728722208663E-14, 3114 1.54008621752140982691E-14, 3115 3.85277838274214270114E-13, 3116 7.18012445138366623367E-13, 3117 -1.79417853150680611778E-12, 3118 -1.32158118404477131188E-11, 3119 -3.14991652796324136454E-11, 3120 1.18891471078464383424E-11, 3121 4.94060238822496958910E-10, 3122 3.39623202570838634515E-9, 3123 2.26666899049817806459E-8, 3124 2.04891858946906374183E-7, 3125 2.89137052083475648297E-6, 3126 6.88975834691682398426E-5, 3127 3.36911647825569408990E-3, 3128 8.04490411014108831608E-1 3129 ] 3130 3131 3132def _chbevl(x, vals): 3133 b0 = vals[0] 3134 b1 = 0.0 3135 3136 for i in range(1, len(vals)): 3137 b2 = b1 3138 b1 = b0 3139 b0 = x*b1 - b2 + vals[i] 3140 3141 return 0.5*(b0 - b2) 3142 3143 3144def _i0_1(x): 3145 return exp(x) * _chbevl(x/2.0-2, _i0A) 3146 3147 3148def _i0_2(x): 3149 return exp(x) * _chbevl(32.0/x - 2.0, _i0B) / sqrt(x) 3150 3151 3152def _i0_dispatcher(x): 3153 return (x,) 3154 3155 3156@array_function_dispatch(_i0_dispatcher) 3157def i0(x): 3158 """ 3159 Modified Bessel function of the first kind, order 0. 3160 3161 Usually denoted :math:`I_0`. 3162 3163 Parameters 3164 ---------- 3165 x : array_like of float 3166 Argument of the Bessel function. 3167 3168 Returns 3169 ------- 3170 out : ndarray, shape = x.shape, dtype = float 3171 The modified Bessel function evaluated at each of the elements of `x`. 3172 3173 See Also 3174 -------- 3175 scipy.special.i0, scipy.special.iv, scipy.special.ive 3176 3177 Notes 3178 ----- 3179 The scipy implementation is recommended over this function: it is a 3180 proper ufunc written in C, and more than an order of magnitude faster. 3181 3182 We use the algorithm published by Clenshaw [1]_ and referenced by 3183 Abramowitz and Stegun [2]_, for which the function domain is 3184 partitioned into the two intervals [0,8] and (8,inf), and Chebyshev 3185 polynomial expansions are employed in each interval. Relative error on 3186 the domain [0,30] using IEEE arithmetic is documented [3]_ as having a 3187 peak of 5.8e-16 with an rms of 1.4e-16 (n = 30000). 3188 3189 References 3190 ---------- 3191 .. [1] C. W. Clenshaw, "Chebyshev series for mathematical functions", in 3192 *National Physical Laboratory Mathematical Tables*, vol. 5, London: 3193 Her Majesty's Stationery Office, 1962. 3194 .. [2] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical 3195 Functions*, 10th printing, New York: Dover, 1964, pp. 379. 3196 http://www.math.sfu.ca/~cbm/aands/page_379.htm 3197 .. [3] https://metacpan.org/pod/distribution/Math-Cephes/lib/Math/Cephes.pod#i0:-Modified-Bessel-function-of-order-zero 3198 3199 Examples 3200 -------- 3201 >>> np.i0(0.) 3202 array(1.0) 3203 >>> np.i0([0, 1, 2, 3]) 3204 array([1. , 1.26606588, 2.2795853 , 4.88079259]) 3205 3206 """ 3207 x = np.asanyarray(x) 3208 if x.dtype.kind == 'c': 3209 raise TypeError("i0 not supported for complex values") 3210 if x.dtype.kind != 'f': 3211 x = x.astype(float) 3212 x = np.abs(x) 3213 return piecewise(x, [x <= 8.0], [_i0_1, _i0_2]) 3214 3215## End of cephes code for i0 3216 3217 3218@set_module('numpy') 3219def kaiser(M, beta): 3220 """ 3221 Return the Kaiser window. 3222 3223 The Kaiser window is a taper formed by using a Bessel function. 3224 3225 Parameters 3226 ---------- 3227 M : int 3228 Number of points in the output window. If zero or less, an 3229 empty array is returned. 3230 beta : float 3231 Shape parameter for window. 3232 3233 Returns 3234 ------- 3235 out : array 3236 The window, with the maximum value normalized to one (the value 3237 one appears only if the number of samples is odd). 3238 3239 See Also 3240 -------- 3241 bartlett, blackman, hamming, hanning 3242 3243 Notes 3244 ----- 3245 The Kaiser window is defined as 3246 3247 .. math:: w(n) = I_0\\left( \\beta \\sqrt{1-\\frac{4n^2}{(M-1)^2}} 3248 \\right)/I_0(\\beta) 3249 3250 with 3251 3252 .. math:: \\quad -\\frac{M-1}{2} \\leq n \\leq \\frac{M-1}{2}, 3253 3254 where :math:`I_0` is the modified zeroth-order Bessel function. 3255 3256 The Kaiser was named for Jim Kaiser, who discovered a simple 3257 approximation to the DPSS window based on Bessel functions. The Kaiser 3258 window is a very good approximation to the Digital Prolate Spheroidal 3259 Sequence, or Slepian window, which is the transform which maximizes the 3260 energy in the main lobe of the window relative to total energy. 3261 3262 The Kaiser can approximate many other windows by varying the beta 3263 parameter. 3264 3265 ==== ======================= 3266 beta Window shape 3267 ==== ======================= 3268 0 Rectangular 3269 5 Similar to a Hamming 3270 6 Similar to a Hanning 3271 8.6 Similar to a Blackman 3272 ==== ======================= 3273 3274 A beta value of 14 is probably a good starting point. Note that as beta 3275 gets large, the window narrows, and so the number of samples needs to be 3276 large enough to sample the increasingly narrow spike, otherwise NaNs will 3277 get returned. 3278 3279 Most references to the Kaiser window come from the signal processing 3280 literature, where it is used as one of many windowing functions for 3281 smoothing values. It is also known as an apodization (which means 3282 "removing the foot", i.e. smoothing discontinuities at the beginning 3283 and end of the sampled signal) or tapering function. 3284 3285 References 3286 ---------- 3287 .. [1] J. F. Kaiser, "Digital Filters" - Ch 7 in "Systems analysis by 3288 digital computer", Editors: F.F. Kuo and J.F. Kaiser, p 218-285. 3289 John Wiley and Sons, New York, (1966). 3290 .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The 3291 University of Alberta Press, 1975, pp. 177-178. 3292 .. [3] Wikipedia, "Window function", 3293 https://en.wikipedia.org/wiki/Window_function 3294 3295 Examples 3296 -------- 3297 >>> import matplotlib.pyplot as plt 3298 >>> np.kaiser(12, 14) 3299 array([7.72686684e-06, 3.46009194e-03, 4.65200189e-02, # may vary 3300 2.29737120e-01, 5.99885316e-01, 9.45674898e-01, 3301 9.45674898e-01, 5.99885316e-01, 2.29737120e-01, 3302 4.65200189e-02, 3.46009194e-03, 7.72686684e-06]) 3303 3304 3305 Plot the window and the frequency response: 3306 3307 >>> from numpy.fft import fft, fftshift 3308 >>> window = np.kaiser(51, 14) 3309 >>> plt.plot(window) 3310 [<matplotlib.lines.Line2D object at 0x...>] 3311 >>> plt.title("Kaiser window") 3312 Text(0.5, 1.0, 'Kaiser window') 3313 >>> plt.ylabel("Amplitude") 3314 Text(0, 0.5, 'Amplitude') 3315 >>> plt.xlabel("Sample") 3316 Text(0.5, 0, 'Sample') 3317 >>> plt.show() 3318 3319 >>> plt.figure() 3320 <Figure size 640x480 with 0 Axes> 3321 >>> A = fft(window, 2048) / 25.5 3322 >>> mag = np.abs(fftshift(A)) 3323 >>> freq = np.linspace(-0.5, 0.5, len(A)) 3324 >>> response = 20 * np.log10(mag) 3325 >>> response = np.clip(response, -100, 100) 3326 >>> plt.plot(freq, response) 3327 [<matplotlib.lines.Line2D object at 0x...>] 3328 >>> plt.title("Frequency response of Kaiser window") 3329 Text(0.5, 1.0, 'Frequency response of Kaiser window') 3330 >>> plt.ylabel("Magnitude [dB]") 3331 Text(0, 0.5, 'Magnitude [dB]') 3332 >>> plt.xlabel("Normalized frequency [cycles per sample]") 3333 Text(0.5, 0, 'Normalized frequency [cycles per sample]') 3334 >>> plt.axis('tight') 3335 (-0.5, 0.5, -100.0, ...) # may vary 3336 >>> plt.show() 3337 3338 """ 3339 if M == 1: 3340 return np.array([1.]) 3341 n = arange(0, M) 3342 alpha = (M-1)/2.0 3343 return i0(beta * sqrt(1-((n-alpha)/alpha)**2.0))/i0(float(beta)) 3344 3345 3346def _sinc_dispatcher(x): 3347 return (x,) 3348 3349 3350@array_function_dispatch(_sinc_dispatcher) 3351def sinc(x): 3352 r""" 3353 Return the normalized sinc function. 3354 3355 The sinc function is :math:`\sin(\pi x)/(\pi x)`. 3356 3357 .. note:: 3358 3359 Note the normalization factor of ``pi`` used in the definition. 3360 This is the most commonly used definition in signal processing. 3361 Use ``sinc(x / np.pi)`` to obtain the unnormalized sinc function 3362 :math:`\sin(x)/(x)` that is more common in mathematics. 3363 3364 Parameters 3365 ---------- 3366 x : ndarray 3367 Array (possibly multi-dimensional) of values for which to to 3368 calculate ``sinc(x)``. 3369 3370 Returns 3371 ------- 3372 out : ndarray 3373 ``sinc(x)``, which has the same shape as the input. 3374 3375 Notes 3376 ----- 3377 ``sinc(0)`` is the limit value 1. 3378 3379 The name sinc is short for "sine cardinal" or "sinus cardinalis". 3380 3381 The sinc function is used in various signal processing applications, 3382 including in anti-aliasing, in the construction of a Lanczos resampling 3383 filter, and in interpolation. 3384 3385 For bandlimited interpolation of discrete-time signals, the ideal 3386 interpolation kernel is proportional to the sinc function. 3387 3388 References 3389 ---------- 3390 .. [1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram Web 3391 Resource. http://mathworld.wolfram.com/SincFunction.html 3392 .. [2] Wikipedia, "Sinc function", 3393 https://en.wikipedia.org/wiki/Sinc_function 3394 3395 Examples 3396 -------- 3397 >>> import matplotlib.pyplot as plt 3398 >>> x = np.linspace(-4, 4, 41) 3399 >>> np.sinc(x) 3400 array([-3.89804309e-17, -4.92362781e-02, -8.40918587e-02, # may vary 3401 -8.90384387e-02, -5.84680802e-02, 3.89804309e-17, 3402 6.68206631e-02, 1.16434881e-01, 1.26137788e-01, 3403 8.50444803e-02, -3.89804309e-17, -1.03943254e-01, 3404 -1.89206682e-01, -2.16236208e-01, -1.55914881e-01, 3405 3.89804309e-17, 2.33872321e-01, 5.04551152e-01, 3406 7.56826729e-01, 9.35489284e-01, 1.00000000e+00, 3407 9.35489284e-01, 7.56826729e-01, 5.04551152e-01, 3408 2.33872321e-01, 3.89804309e-17, -1.55914881e-01, 3409 -2.16236208e-01, -1.89206682e-01, -1.03943254e-01, 3410 -3.89804309e-17, 8.50444803e-02, 1.26137788e-01, 3411 1.16434881e-01, 6.68206631e-02, 3.89804309e-17, 3412 -5.84680802e-02, -8.90384387e-02, -8.40918587e-02, 3413 -4.92362781e-02, -3.89804309e-17]) 3414 3415 >>> plt.plot(x, np.sinc(x)) 3416 [<matplotlib.lines.Line2D object at 0x...>] 3417 >>> plt.title("Sinc Function") 3418 Text(0.5, 1.0, 'Sinc Function') 3419 >>> plt.ylabel("Amplitude") 3420 Text(0, 0.5, 'Amplitude') 3421 >>> plt.xlabel("X") 3422 Text(0.5, 0, 'X') 3423 >>> plt.show() 3424 3425 """ 3426 x = np.asanyarray(x) 3427 y = pi * where(x == 0, 1.0e-20, x) 3428 return sin(y)/y 3429 3430 3431def _msort_dispatcher(a): 3432 return (a,) 3433 3434 3435@array_function_dispatch(_msort_dispatcher) 3436def msort(a): 3437 """ 3438 Return a copy of an array sorted along the first axis. 3439 3440 Parameters 3441 ---------- 3442 a : array_like 3443 Array to be sorted. 3444 3445 Returns 3446 ------- 3447 sorted_array : ndarray 3448 Array of the same type and shape as `a`. 3449 3450 See Also 3451 -------- 3452 sort 3453 3454 Notes 3455 ----- 3456 ``np.msort(a)`` is equivalent to ``np.sort(a, axis=0)``. 3457 3458 """ 3459 b = array(a, subok=True, copy=True) 3460 b.sort(0) 3461 return b 3462 3463 3464def _ureduce(a, func, **kwargs): 3465 """ 3466 Internal Function. 3467 Call `func` with `a` as first argument swapping the axes to use extended 3468 axis on functions that don't support it natively. 3469 3470 Returns result and a.shape with axis dims set to 1. 3471 3472 Parameters 3473 ---------- 3474 a : array_like 3475 Input array or object that can be converted to an array. 3476 func : callable 3477 Reduction function capable of receiving a single axis argument. 3478 It is called with `a` as first argument followed by `kwargs`. 3479 kwargs : keyword arguments 3480 additional keyword arguments to pass to `func`. 3481 3482 Returns 3483 ------- 3484 result : tuple 3485 Result of func(a, **kwargs) and a.shape with axis dims set to 1 3486 which can be used to reshape the result to the same shape a ufunc with 3487 keepdims=True would produce. 3488 3489 """ 3490 a = np.asanyarray(a) 3491 axis = kwargs.get('axis', None) 3492 if axis is not None: 3493 keepdim = list(a.shape) 3494 nd = a.ndim 3495 axis = _nx.normalize_axis_tuple(axis, nd) 3496 3497 for ax in axis: 3498 keepdim[ax] = 1 3499 3500 if len(axis) == 1: 3501 kwargs['axis'] = axis[0] 3502 else: 3503 keep = set(range(nd)) - set(axis) 3504 nkeep = len(keep) 3505 # swap axis that should not be reduced to front 3506 for i, s in enumerate(sorted(keep)): 3507 a = a.swapaxes(i, s) 3508 # merge reduced axis 3509 a = a.reshape(a.shape[:nkeep] + (-1,)) 3510 kwargs['axis'] = -1 3511 keepdim = tuple(keepdim) 3512 else: 3513 keepdim = (1,) * a.ndim 3514 3515 r = func(a, **kwargs) 3516 return r, keepdim 3517 3518 3519def _median_dispatcher( 3520 a, axis=None, out=None, overwrite_input=None, keepdims=None): 3521 return (a, out) 3522 3523 3524@array_function_dispatch(_median_dispatcher) 3525def median(a, axis=None, out=None, overwrite_input=False, keepdims=False): 3526 """ 3527 Compute the median along the specified axis. 3528 3529 Returns the median of the array elements. 3530 3531 Parameters 3532 ---------- 3533 a : array_like 3534 Input array or object that can be converted to an array. 3535 axis : {int, sequence of int, None}, optional 3536 Axis or axes along which the medians are computed. The default 3537 is to compute the median along a flattened version of the array. 3538 A sequence of axes is supported since version 1.9.0. 3539 out : ndarray, optional 3540 Alternative output array in which to place the result. It must 3541 have the same shape and buffer length as the expected output, 3542 but the type (of the output) will be cast if necessary. 3543 overwrite_input : bool, optional 3544 If True, then allow use of memory of input array `a` for 3545 calculations. The input array will be modified by the call to 3546 `median`. This will save memory when you do not need to preserve 3547 the contents of the input array. Treat the input as undefined, 3548 but it will probably be fully or partially sorted. Default is 3549 False. If `overwrite_input` is ``True`` and `a` is not already an 3550 `ndarray`, an error will be raised. 3551 keepdims : bool, optional 3552 If this is set to True, the axes which are reduced are left 3553 in the result as dimensions with size one. With this option, 3554 the result will broadcast correctly against the original `arr`. 3555 3556 .. versionadded:: 1.9.0 3557 3558 Returns 3559 ------- 3560 median : ndarray 3561 A new array holding the result. If the input contains integers 3562 or floats smaller than ``float64``, then the output data-type is 3563 ``np.float64``. Otherwise, the data-type of the output is the 3564 same as that of the input. If `out` is specified, that array is 3565 returned instead. 3566 3567 See Also 3568 -------- 3569 mean, percentile 3570 3571 Notes 3572 ----- 3573 Given a vector ``V`` of length ``N``, the median of ``V`` is the 3574 middle value of a sorted copy of ``V``, ``V_sorted`` - i 3575 e., ``V_sorted[(N-1)/2]``, when ``N`` is odd, and the average of the 3576 two middle values of ``V_sorted`` when ``N`` is even. 3577 3578 Examples 3579 -------- 3580 >>> a = np.array([[10, 7, 4], [3, 2, 1]]) 3581 >>> a 3582 array([[10, 7, 4], 3583 [ 3, 2, 1]]) 3584 >>> np.median(a) 3585 3.5 3586 >>> np.median(a, axis=0) 3587 array([6.5, 4.5, 2.5]) 3588 >>> np.median(a, axis=1) 3589 array([7., 2.]) 3590 >>> m = np.median(a, axis=0) 3591 >>> out = np.zeros_like(m) 3592 >>> np.median(a, axis=0, out=m) 3593 array([6.5, 4.5, 2.5]) 3594 >>> m 3595 array([6.5, 4.5, 2.5]) 3596 >>> b = a.copy() 3597 >>> np.median(b, axis=1, overwrite_input=True) 3598 array([7., 2.]) 3599 >>> assert not np.all(a==b) 3600 >>> b = a.copy() 3601 >>> np.median(b, axis=None, overwrite_input=True) 3602 3.5 3603 >>> assert not np.all(a==b) 3604 3605 """ 3606 r, k = _ureduce(a, func=_median, axis=axis, out=out, 3607 overwrite_input=overwrite_input) 3608 if keepdims: 3609 return r.reshape(k) 3610 else: 3611 return r 3612 3613 3614def _median(a, axis=None, out=None, overwrite_input=False): 3615 # can't be reasonably be implemented in terms of percentile as we have to 3616 # call mean to not break astropy 3617 a = np.asanyarray(a) 3618 3619 # Set the partition indexes 3620 if axis is None: 3621 sz = a.size 3622 else: 3623 sz = a.shape[axis] 3624 if sz % 2 == 0: 3625 szh = sz // 2 3626 kth = [szh - 1, szh] 3627 else: 3628 kth = [(sz - 1) // 2] 3629 # Check if the array contains any nan's 3630 if np.issubdtype(a.dtype, np.inexact): 3631 kth.append(-1) 3632 3633 if overwrite_input: 3634 if axis is None: 3635 part = a.ravel() 3636 part.partition(kth) 3637 else: 3638 a.partition(kth, axis=axis) 3639 part = a 3640 else: 3641 part = partition(a, kth, axis=axis) 3642 3643 if part.shape == (): 3644 # make 0-D arrays work 3645 return part.item() 3646 if axis is None: 3647 axis = 0 3648 3649 indexer = [slice(None)] * part.ndim 3650 index = part.shape[axis] // 2 3651 if part.shape[axis] % 2 == 1: 3652 # index with slice to allow mean (below) to work 3653 indexer[axis] = slice(index, index+1) 3654 else: 3655 indexer[axis] = slice(index-1, index+1) 3656 indexer = tuple(indexer) 3657 3658 # Check if the array contains any nan's 3659 if np.issubdtype(a.dtype, np.inexact) and sz > 0: 3660 # warn and return nans like mean would 3661 rout = mean(part[indexer], axis=axis, out=out) 3662 return np.lib.utils._median_nancheck(part, rout, axis, out) 3663 else: 3664 # if there are no nans 3665 # Use mean in odd and even case to coerce data type 3666 # and check, use out array. 3667 return mean(part[indexer], axis=axis, out=out) 3668 3669 3670def _percentile_dispatcher(a, q, axis=None, out=None, overwrite_input=None, 3671 interpolation=None, keepdims=None): 3672 return (a, q, out) 3673 3674 3675@array_function_dispatch(_percentile_dispatcher) 3676def percentile(a, q, axis=None, out=None, 3677 overwrite_input=False, interpolation='linear', keepdims=False): 3678 """ 3679 Compute the q-th percentile of the data along the specified axis. 3680 3681 Returns the q-th percentile(s) of the array elements. 3682 3683 Parameters 3684 ---------- 3685 a : array_like 3686 Input array or object that can be converted to an array. 3687 q : array_like of float 3688 Percentile or sequence of percentiles to compute, which must be between 3689 0 and 100 inclusive. 3690 axis : {int, tuple of int, None}, optional 3691 Axis or axes along which the percentiles are computed. The 3692 default is to compute the percentile(s) along a flattened 3693 version of the array. 3694 3695 .. versionchanged:: 1.9.0 3696 A tuple of axes is supported 3697 out : ndarray, optional 3698 Alternative output array in which to place the result. It must 3699 have the same shape and buffer length as the expected output, 3700 but the type (of the output) will be cast if necessary. 3701 overwrite_input : bool, optional 3702 If True, then allow the input array `a` to be modified by intermediate 3703 calculations, to save memory. In this case, the contents of the input 3704 `a` after this function completes is undefined. 3705 3706 interpolation : {'linear', 'lower', 'higher', 'midpoint', 'nearest'} 3707 This optional parameter specifies the interpolation method to 3708 use when the desired percentile lies between two data points 3709 ``i < j``: 3710 3711 * 'linear': ``i + (j - i) * fraction``, where ``fraction`` 3712 is the fractional part of the index surrounded by ``i`` 3713 and ``j``. 3714 * 'lower': ``i``. 3715 * 'higher': ``j``. 3716 * 'nearest': ``i`` or ``j``, whichever is nearest. 3717 * 'midpoint': ``(i + j) / 2``. 3718 3719 .. versionadded:: 1.9.0 3720 keepdims : bool, optional 3721 If this is set to True, the axes which are reduced are left in 3722 the result as dimensions with size one. With this option, the 3723 result will broadcast correctly against the original array `a`. 3724 3725 .. versionadded:: 1.9.0 3726 3727 Returns 3728 ------- 3729 percentile : scalar or ndarray 3730 If `q` is a single percentile and `axis=None`, then the result 3731 is a scalar. If multiple percentiles are given, first axis of 3732 the result corresponds to the percentiles. The other axes are 3733 the axes that remain after the reduction of `a`. If the input 3734 contains integers or floats smaller than ``float64``, the output 3735 data-type is ``float64``. Otherwise, the output data-type is the 3736 same as that of the input. If `out` is specified, that array is 3737 returned instead. 3738 3739 See Also 3740 -------- 3741 mean 3742 median : equivalent to ``percentile(..., 50)`` 3743 nanpercentile 3744 quantile : equivalent to percentile, except with q in the range [0, 1]. 3745 3746 Notes 3747 ----- 3748 Given a vector ``V`` of length ``N``, the q-th percentile of 3749 ``V`` is the value ``q/100`` of the way from the minimum to the 3750 maximum in a sorted copy of ``V``. The values and distances of 3751 the two nearest neighbors as well as the `interpolation` parameter 3752 will determine the percentile if the normalized ranking does not 3753 match the location of ``q`` exactly. This function is the same as 3754 the median if ``q=50``, the same as the minimum if ``q=0`` and the 3755 same as the maximum if ``q=100``. 3756 3757 Examples 3758 -------- 3759 >>> a = np.array([[10, 7, 4], [3, 2, 1]]) 3760 >>> a 3761 array([[10, 7, 4], 3762 [ 3, 2, 1]]) 3763 >>> np.percentile(a, 50) 3764 3.5 3765 >>> np.percentile(a, 50, axis=0) 3766 array([6.5, 4.5, 2.5]) 3767 >>> np.percentile(a, 50, axis=1) 3768 array([7., 2.]) 3769 >>> np.percentile(a, 50, axis=1, keepdims=True) 3770 array([[7.], 3771 [2.]]) 3772 3773 >>> m = np.percentile(a, 50, axis=0) 3774 >>> out = np.zeros_like(m) 3775 >>> np.percentile(a, 50, axis=0, out=out) 3776 array([6.5, 4.5, 2.5]) 3777 >>> m 3778 array([6.5, 4.5, 2.5]) 3779 3780 >>> b = a.copy() 3781 >>> np.percentile(b, 50, axis=1, overwrite_input=True) 3782 array([7., 2.]) 3783 >>> assert not np.all(a == b) 3784 3785 The different types of interpolation can be visualized graphically: 3786 3787 .. plot:: 3788 3789 import matplotlib.pyplot as plt 3790 3791 a = np.arange(4) 3792 p = np.linspace(0, 100, 6001) 3793 ax = plt.gca() 3794 lines = [ 3795 ('linear', None), 3796 ('higher', '--'), 3797 ('lower', '--'), 3798 ('nearest', '-.'), 3799 ('midpoint', '-.'), 3800 ] 3801 for interpolation, style in lines: 3802 ax.plot( 3803 p, np.percentile(a, p, interpolation=interpolation), 3804 label=interpolation, linestyle=style) 3805 ax.set( 3806 title='Interpolation methods for list: ' + str(a), 3807 xlabel='Percentile', 3808 ylabel='List item returned', 3809 yticks=a) 3810 ax.legend() 3811 plt.show() 3812 3813 """ 3814 q = np.true_divide(q, 100) 3815 q = asanyarray(q) # undo any decay that the ufunc performed (see gh-13105) 3816 if not _quantile_is_valid(q): 3817 raise ValueError("Percentiles must be in the range [0, 100]") 3818 return _quantile_unchecked( 3819 a, q, axis, out, overwrite_input, interpolation, keepdims) 3820 3821 3822def _quantile_dispatcher(a, q, axis=None, out=None, overwrite_input=None, 3823 interpolation=None, keepdims=None): 3824 return (a, q, out) 3825 3826 3827@array_function_dispatch(_quantile_dispatcher) 3828def quantile(a, q, axis=None, out=None, 3829 overwrite_input=False, interpolation='linear', keepdims=False): 3830 """ 3831 Compute the q-th quantile of the data along the specified axis. 3832 3833 .. versionadded:: 1.15.0 3834 3835 Parameters 3836 ---------- 3837 a : array_like 3838 Input array or object that can be converted to an array. 3839 q : array_like of float 3840 Quantile or sequence of quantiles to compute, which must be between 3841 0 and 1 inclusive. 3842 axis : {int, tuple of int, None}, optional 3843 Axis or axes along which the quantiles are computed. The 3844 default is to compute the quantile(s) along a flattened 3845 version of the array. 3846 out : ndarray, optional 3847 Alternative output array in which to place the result. It must 3848 have the same shape and buffer length as the expected output, 3849 but the type (of the output) will be cast if necessary. 3850 overwrite_input : bool, optional 3851 If True, then allow the input array `a` to be modified by intermediate 3852 calculations, to save memory. In this case, the contents of the input 3853 `a` after this function completes is undefined. 3854 interpolation : {'linear', 'lower', 'higher', 'midpoint', 'nearest'} 3855 This optional parameter specifies the interpolation method to 3856 use when the desired quantile lies between two data points 3857 ``i < j``: 3858 3859 * linear: ``i + (j - i) * fraction``, where ``fraction`` 3860 is the fractional part of the index surrounded by ``i`` 3861 and ``j``. 3862 * lower: ``i``. 3863 * higher: ``j``. 3864 * nearest: ``i`` or ``j``, whichever is nearest. 3865 * midpoint: ``(i + j) / 2``. 3866 keepdims : bool, optional 3867 If this is set to True, the axes which are reduced are left in 3868 the result as dimensions with size one. With this option, the 3869 result will broadcast correctly against the original array `a`. 3870 3871 Returns 3872 ------- 3873 quantile : scalar or ndarray 3874 If `q` is a single quantile and `axis=None`, then the result 3875 is a scalar. If multiple quantiles are given, first axis of 3876 the result corresponds to the quantiles. The other axes are 3877 the axes that remain after the reduction of `a`. If the input 3878 contains integers or floats smaller than ``float64``, the output 3879 data-type is ``float64``. Otherwise, the output data-type is the 3880 same as that of the input. If `out` is specified, that array is 3881 returned instead. 3882 3883 See Also 3884 -------- 3885 mean 3886 percentile : equivalent to quantile, but with q in the range [0, 100]. 3887 median : equivalent to ``quantile(..., 0.5)`` 3888 nanquantile 3889 3890 Notes 3891 ----- 3892 Given a vector ``V`` of length ``N``, the q-th quantile of 3893 ``V`` is the value ``q`` of the way from the minimum to the 3894 maximum in a sorted copy of ``V``. The values and distances of 3895 the two nearest neighbors as well as the `interpolation` parameter 3896 will determine the quantile if the normalized ranking does not 3897 match the location of ``q`` exactly. This function is the same as 3898 the median if ``q=0.5``, the same as the minimum if ``q=0.0`` and the 3899 same as the maximum if ``q=1.0``. 3900 3901 Examples 3902 -------- 3903 >>> a = np.array([[10, 7, 4], [3, 2, 1]]) 3904 >>> a 3905 array([[10, 7, 4], 3906 [ 3, 2, 1]]) 3907 >>> np.quantile(a, 0.5) 3908 3.5 3909 >>> np.quantile(a, 0.5, axis=0) 3910 array([6.5, 4.5, 2.5]) 3911 >>> np.quantile(a, 0.5, axis=1) 3912 array([7., 2.]) 3913 >>> np.quantile(a, 0.5, axis=1, keepdims=True) 3914 array([[7.], 3915 [2.]]) 3916 >>> m = np.quantile(a, 0.5, axis=0) 3917 >>> out = np.zeros_like(m) 3918 >>> np.quantile(a, 0.5, axis=0, out=out) 3919 array([6.5, 4.5, 2.5]) 3920 >>> m 3921 array([6.5, 4.5, 2.5]) 3922 >>> b = a.copy() 3923 >>> np.quantile(b, 0.5, axis=1, overwrite_input=True) 3924 array([7., 2.]) 3925 >>> assert not np.all(a == b) 3926 """ 3927 q = np.asanyarray(q) 3928 if not _quantile_is_valid(q): 3929 raise ValueError("Quantiles must be in the range [0, 1]") 3930 return _quantile_unchecked( 3931 a, q, axis, out, overwrite_input, interpolation, keepdims) 3932 3933 3934def _quantile_unchecked(a, q, axis=None, out=None, overwrite_input=False, 3935 interpolation='linear', keepdims=False): 3936 """Assumes that q is in [0, 1], and is an ndarray""" 3937 r, k = _ureduce(a, func=_quantile_ureduce_func, q=q, axis=axis, out=out, 3938 overwrite_input=overwrite_input, 3939 interpolation=interpolation) 3940 if keepdims: 3941 return r.reshape(q.shape + k) 3942 else: 3943 return r 3944 3945 3946def _quantile_is_valid(q): 3947 # avoid expensive reductions, relevant for arrays with < O(1000) elements 3948 if q.ndim == 1 and q.size < 10: 3949 for i in range(q.size): 3950 if q[i] < 0.0 or q[i] > 1.0: 3951 return False 3952 else: 3953 # faster than any() 3954 if np.count_nonzero(q < 0.0) or np.count_nonzero(q > 1.0): 3955 return False 3956 return True 3957 3958 3959def _lerp(a, b, t, out=None): 3960 """ Linearly interpolate from a to b by a factor of t """ 3961 diff_b_a = subtract(b, a) 3962 # asanyarray is a stop-gap until gh-13105 3963 lerp_interpolation = asanyarray(add(a, diff_b_a*t, out=out)) 3964 subtract(b, diff_b_a * (1 - t), out=lerp_interpolation, where=t>=0.5) 3965 if lerp_interpolation.ndim == 0 and out is None: 3966 lerp_interpolation = lerp_interpolation[()] # unpack 0d arrays 3967 return lerp_interpolation 3968 3969 3970def _quantile_ureduce_func(a, q, axis=None, out=None, overwrite_input=False, 3971 interpolation='linear', keepdims=False): 3972 a = asarray(a) 3973 3974 # ufuncs cause 0d array results to decay to scalars (see gh-13105), which 3975 # makes them problematic for __setitem__ and attribute access. As a 3976 # workaround, we call this on the result of every ufunc on a possibly-0d 3977 # array. 3978 not_scalar = np.asanyarray 3979 3980 # prepare a for partitioning 3981 if overwrite_input: 3982 if axis is None: 3983 ap = a.ravel() 3984 else: 3985 ap = a 3986 else: 3987 if axis is None: 3988 ap = a.flatten() 3989 else: 3990 ap = a.copy() 3991 3992 if axis is None: 3993 axis = 0 3994 3995 if q.ndim > 2: 3996 # The code below works fine for nd, but it might not have useful 3997 # semantics. For now, keep the supported dimensions the same as it was 3998 # before. 3999 raise ValueError("q must be a scalar or 1d") 4000 4001 Nx = ap.shape[axis] 4002 indices = not_scalar(q * (Nx - 1)) 4003 # round fractional indices according to interpolation method 4004 if interpolation == 'lower': 4005 indices = floor(indices).astype(intp) 4006 elif interpolation == 'higher': 4007 indices = ceil(indices).astype(intp) 4008 elif interpolation == 'midpoint': 4009 indices = 0.5 * (floor(indices) + ceil(indices)) 4010 elif interpolation == 'nearest': 4011 indices = around(indices).astype(intp) 4012 elif interpolation == 'linear': 4013 pass # keep index as fraction and interpolate 4014 else: 4015 raise ValueError( 4016 "interpolation can only be 'linear', 'lower' 'higher', " 4017 "'midpoint', or 'nearest'") 4018 4019 # The dimensions of `q` are prepended to the output shape, so we need the 4020 # axis being sampled from `ap` to be first. 4021 ap = np.moveaxis(ap, axis, 0) 4022 del axis 4023 4024 if np.issubdtype(indices.dtype, np.integer): 4025 # take the points along axis 4026 4027 if np.issubdtype(a.dtype, np.inexact): 4028 # may contain nan, which would sort to the end 4029 ap.partition(concatenate((indices.ravel(), [-1])), axis=0) 4030 n = np.isnan(ap[-1]) 4031 else: 4032 # cannot contain nan 4033 ap.partition(indices.ravel(), axis=0) 4034 n = np.array(False, dtype=bool) 4035 4036 r = take(ap, indices, axis=0, out=out) 4037 4038 else: 4039 # weight the points above and below the indices 4040 4041 indices_below = not_scalar(floor(indices)).astype(intp) 4042 indices_above = not_scalar(indices_below + 1) 4043 indices_above[indices_above > Nx - 1] = Nx - 1 4044 4045 if np.issubdtype(a.dtype, np.inexact): 4046 # may contain nan, which would sort to the end 4047 ap.partition(concatenate(( 4048 indices_below.ravel(), indices_above.ravel(), [-1] 4049 )), axis=0) 4050 n = np.isnan(ap[-1]) 4051 else: 4052 # cannot contain nan 4053 ap.partition(concatenate(( 4054 indices_below.ravel(), indices_above.ravel() 4055 )), axis=0) 4056 n = np.array(False, dtype=bool) 4057 4058 weights_shape = indices.shape + (1,) * (ap.ndim - 1) 4059 weights_above = not_scalar(indices - indices_below).reshape(weights_shape) 4060 4061 x_below = take(ap, indices_below, axis=0) 4062 x_above = take(ap, indices_above, axis=0) 4063 4064 r = _lerp(x_below, x_above, weights_above, out=out) 4065 4066 # if any slice contained a nan, then all results on that slice are also nan 4067 if np.any(n): 4068 if r.ndim == 0 and out is None: 4069 # can't write to a scalar 4070 r = a.dtype.type(np.nan) 4071 else: 4072 r[..., n] = a.dtype.type(np.nan) 4073 4074 return r 4075 4076 4077def _trapz_dispatcher(y, x=None, dx=None, axis=None): 4078 return (y, x) 4079 4080 4081@array_function_dispatch(_trapz_dispatcher) 4082def trapz(y, x=None, dx=1.0, axis=-1): 4083 """ 4084 Integrate along the given axis using the composite trapezoidal rule. 4085 4086 Integrate `y` (`x`) along given axis. 4087 4088 Parameters 4089 ---------- 4090 y : array_like 4091 Input array to integrate. 4092 x : array_like, optional 4093 The sample points corresponding to the `y` values. If `x` is None, 4094 the sample points are assumed to be evenly spaced `dx` apart. The 4095 default is None. 4096 dx : scalar, optional 4097 The spacing between sample points when `x` is None. The default is 1. 4098 axis : int, optional 4099 The axis along which to integrate. 4100 4101 Returns 4102 ------- 4103 trapz : float 4104 Definite integral as approximated by trapezoidal rule. 4105 4106 See Also 4107 -------- 4108 sum, cumsum 4109 4110 Notes 4111 ----- 4112 Image [2]_ illustrates trapezoidal rule -- y-axis locations of points 4113 will be taken from `y` array, by default x-axis distances between 4114 points will be 1.0, alternatively they can be provided with `x` array 4115 or with `dx` scalar. Return value will be equal to combined area under 4116 the red lines. 4117 4118 4119 References 4120 ---------- 4121 .. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule 4122 4123 .. [2] Illustration image: 4124 https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png 4125 4126 Examples 4127 -------- 4128 >>> np.trapz([1,2,3]) 4129 4.0 4130 >>> np.trapz([1,2,3], x=[4,6,8]) 4131 8.0 4132 >>> np.trapz([1,2,3], dx=2) 4133 8.0 4134 >>> a = np.arange(6).reshape(2, 3) 4135 >>> a 4136 array([[0, 1, 2], 4137 [3, 4, 5]]) 4138 >>> np.trapz(a, axis=0) 4139 array([1.5, 2.5, 3.5]) 4140 >>> np.trapz(a, axis=1) 4141 array([2., 8.]) 4142 4143 """ 4144 y = asanyarray(y) 4145 if x is None: 4146 d = dx 4147 else: 4148 x = asanyarray(x) 4149 if x.ndim == 1: 4150 d = diff(x) 4151 # reshape to correct shape 4152 shape = [1]*y.ndim 4153 shape[axis] = d.shape[0] 4154 d = d.reshape(shape) 4155 else: 4156 d = diff(x, axis=axis) 4157 nd = y.ndim 4158 slice1 = [slice(None)]*nd 4159 slice2 = [slice(None)]*nd 4160 slice1[axis] = slice(1, None) 4161 slice2[axis] = slice(None, -1) 4162 try: 4163 ret = (d * (y[tuple(slice1)] + y[tuple(slice2)]) / 2.0).sum(axis) 4164 except ValueError: 4165 # Operations didn't work, cast to ndarray 4166 d = np.asarray(d) 4167 y = np.asarray(y) 4168 ret = add.reduce(d * (y[tuple(slice1)]+y[tuple(slice2)])/2.0, axis) 4169 return ret 4170 4171 4172def _meshgrid_dispatcher(*xi, copy=None, sparse=None, indexing=None): 4173 return xi 4174 4175 4176# Based on scitools meshgrid 4177@array_function_dispatch(_meshgrid_dispatcher) 4178def meshgrid(*xi, copy=True, sparse=False, indexing='xy'): 4179 """ 4180 Return coordinate matrices from coordinate vectors. 4181 4182 Make N-D coordinate arrays for vectorized evaluations of 4183 N-D scalar/vector fields over N-D grids, given 4184 one-dimensional coordinate arrays x1, x2,..., xn. 4185 4186 .. versionchanged:: 1.9 4187 1-D and 0-D cases are allowed. 4188 4189 Parameters 4190 ---------- 4191 x1, x2,..., xn : array_like 4192 1-D arrays representing the coordinates of a grid. 4193 indexing : {'xy', 'ij'}, optional 4194 Cartesian ('xy', default) or matrix ('ij') indexing of output. 4195 See Notes for more details. 4196 4197 .. versionadded:: 1.7.0 4198 sparse : bool, optional 4199 If True a sparse grid is returned in order to conserve memory. 4200 Default is False. 4201 4202 .. versionadded:: 1.7.0 4203 copy : bool, optional 4204 If False, a view into the original arrays are returned in order to 4205 conserve memory. Default is True. Please note that 4206 ``sparse=False, copy=False`` will likely return non-contiguous 4207 arrays. Furthermore, more than one element of a broadcast array 4208 may refer to a single memory location. If you need to write to the 4209 arrays, make copies first. 4210 4211 .. versionadded:: 1.7.0 4212 4213 Returns 4214 ------- 4215 X1, X2,..., XN : ndarray 4216 For vectors `x1`, `x2`,..., 'xn' with lengths ``Ni=len(xi)`` , 4217 return ``(N1, N2, N3,...Nn)`` shaped arrays if indexing='ij' 4218 or ``(N2, N1, N3,...Nn)`` shaped arrays if indexing='xy' 4219 with the elements of `xi` repeated to fill the matrix along 4220 the first dimension for `x1`, the second for `x2` and so on. 4221 4222 Notes 4223 ----- 4224 This function supports both indexing conventions through the indexing 4225 keyword argument. Giving the string 'ij' returns a meshgrid with 4226 matrix indexing, while 'xy' returns a meshgrid with Cartesian indexing. 4227 In the 2-D case with inputs of length M and N, the outputs are of shape 4228 (N, M) for 'xy' indexing and (M, N) for 'ij' indexing. In the 3-D case 4229 with inputs of length M, N and P, outputs are of shape (N, M, P) for 4230 'xy' indexing and (M, N, P) for 'ij' indexing. The difference is 4231 illustrated by the following code snippet:: 4232 4233 xv, yv = np.meshgrid(x, y, sparse=False, indexing='ij') 4234 for i in range(nx): 4235 for j in range(ny): 4236 # treat xv[i,j], yv[i,j] 4237 4238 xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy') 4239 for i in range(nx): 4240 for j in range(ny): 4241 # treat xv[j,i], yv[j,i] 4242 4243 In the 1-D and 0-D case, the indexing and sparse keywords have no effect. 4244 4245 See Also 4246 -------- 4247 mgrid : Construct a multi-dimensional "meshgrid" using indexing notation. 4248 ogrid : Construct an open multi-dimensional "meshgrid" using indexing 4249 notation. 4250 4251 Examples 4252 -------- 4253 >>> nx, ny = (3, 2) 4254 >>> x = np.linspace(0, 1, nx) 4255 >>> y = np.linspace(0, 1, ny) 4256 >>> xv, yv = np.meshgrid(x, y) 4257 >>> xv 4258 array([[0. , 0.5, 1. ], 4259 [0. , 0.5, 1. ]]) 4260 >>> yv 4261 array([[0., 0., 0.], 4262 [1., 1., 1.]]) 4263 >>> xv, yv = np.meshgrid(x, y, sparse=True) # make sparse output arrays 4264 >>> xv 4265 array([[0. , 0.5, 1. ]]) 4266 >>> yv 4267 array([[0.], 4268 [1.]]) 4269 4270 `meshgrid` is very useful to evaluate functions on a grid. 4271 4272 >>> import matplotlib.pyplot as plt 4273 >>> x = np.arange(-5, 5, 0.1) 4274 >>> y = np.arange(-5, 5, 0.1) 4275 >>> xx, yy = np.meshgrid(x, y, sparse=True) 4276 >>> z = np.sin(xx**2 + yy**2) / (xx**2 + yy**2) 4277 >>> h = plt.contourf(x,y,z) 4278 >>> plt.show() 4279 4280 """ 4281 ndim = len(xi) 4282 4283 if indexing not in ['xy', 'ij']: 4284 raise ValueError( 4285 "Valid values for `indexing` are 'xy' and 'ij'.") 4286 4287 s0 = (1,) * ndim 4288 output = [np.asanyarray(x).reshape(s0[:i] + (-1,) + s0[i + 1:]) 4289 for i, x in enumerate(xi)] 4290 4291 if indexing == 'xy' and ndim > 1: 4292 # switch first and second axis 4293 output[0].shape = (1, -1) + s0[2:] 4294 output[1].shape = (-1, 1) + s0[2:] 4295 4296 if not sparse: 4297 # Return the full N-D matrix (not only the 1-D vector) 4298 output = np.broadcast_arrays(*output, subok=True) 4299 4300 if copy: 4301 output = [x.copy() for x in output] 4302 4303 return output 4304 4305 4306def _delete_dispatcher(arr, obj, axis=None): 4307 return (arr, obj) 4308 4309 4310@array_function_dispatch(_delete_dispatcher) 4311def delete(arr, obj, axis=None): 4312 """ 4313 Return a new array with sub-arrays along an axis deleted. For a one 4314 dimensional array, this returns those entries not returned by 4315 `arr[obj]`. 4316 4317 Parameters 4318 ---------- 4319 arr : array_like 4320 Input array. 4321 obj : slice, int or array of ints 4322 Indicate indices of sub-arrays to remove along the specified axis. 4323 4324 .. versionchanged:: 1.19.0 4325 Boolean indices are now treated as a mask of elements to remove, 4326 rather than being cast to the integers 0 and 1. 4327 4328 axis : int, optional 4329 The axis along which to delete the subarray defined by `obj`. 4330 If `axis` is None, `obj` is applied to the flattened array. 4331 4332 Returns 4333 ------- 4334 out : ndarray 4335 A copy of `arr` with the elements specified by `obj` removed. Note 4336 that `delete` does not occur in-place. If `axis` is None, `out` is 4337 a flattened array. 4338 4339 See Also 4340 -------- 4341 insert : Insert elements into an array. 4342 append : Append elements at the end of an array. 4343 4344 Notes 4345 ----- 4346 Often it is preferable to use a boolean mask. For example: 4347 4348 >>> arr = np.arange(12) + 1 4349 >>> mask = np.ones(len(arr), dtype=bool) 4350 >>> mask[[0,2,4]] = False 4351 >>> result = arr[mask,...] 4352 4353 Is equivalent to `np.delete(arr, [0,2,4], axis=0)`, but allows further 4354 use of `mask`. 4355 4356 Examples 4357 -------- 4358 >>> arr = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]]) 4359 >>> arr 4360 array([[ 1, 2, 3, 4], 4361 [ 5, 6, 7, 8], 4362 [ 9, 10, 11, 12]]) 4363 >>> np.delete(arr, 1, 0) 4364 array([[ 1, 2, 3, 4], 4365 [ 9, 10, 11, 12]]) 4366 4367 >>> np.delete(arr, np.s_[::2], 1) 4368 array([[ 2, 4], 4369 [ 6, 8], 4370 [10, 12]]) 4371 >>> np.delete(arr, [1,3,5], None) 4372 array([ 1, 3, 5, 7, 8, 9, 10, 11, 12]) 4373 4374 """ 4375 wrap = None 4376 if type(arr) is not ndarray: 4377 try: 4378 wrap = arr.__array_wrap__ 4379 except AttributeError: 4380 pass 4381 4382 arr = asarray(arr) 4383 ndim = arr.ndim 4384 arrorder = 'F' if arr.flags.fnc else 'C' 4385 if axis is None: 4386 if ndim != 1: 4387 arr = arr.ravel() 4388 # needed for np.matrix, which is still not 1d after being ravelled 4389 ndim = arr.ndim 4390 axis = ndim - 1 4391 else: 4392 axis = normalize_axis_index(axis, ndim) 4393 4394 slobj = [slice(None)]*ndim 4395 N = arr.shape[axis] 4396 newshape = list(arr.shape) 4397 4398 if isinstance(obj, slice): 4399 start, stop, step = obj.indices(N) 4400 xr = range(start, stop, step) 4401 numtodel = len(xr) 4402 4403 if numtodel <= 0: 4404 if wrap: 4405 return wrap(arr.copy(order=arrorder)) 4406 else: 4407 return arr.copy(order=arrorder) 4408 4409 # Invert if step is negative: 4410 if step < 0: 4411 step = -step 4412 start = xr[-1] 4413 stop = xr[0] + 1 4414 4415 newshape[axis] -= numtodel 4416 new = empty(newshape, arr.dtype, arrorder) 4417 # copy initial chunk 4418 if start == 0: 4419 pass 4420 else: 4421 slobj[axis] = slice(None, start) 4422 new[tuple(slobj)] = arr[tuple(slobj)] 4423 # copy end chunk 4424 if stop == N: 4425 pass 4426 else: 4427 slobj[axis] = slice(stop-numtodel, None) 4428 slobj2 = [slice(None)]*ndim 4429 slobj2[axis] = slice(stop, None) 4430 new[tuple(slobj)] = arr[tuple(slobj2)] 4431 # copy middle pieces 4432 if step == 1: 4433 pass 4434 else: # use array indexing. 4435 keep = ones(stop-start, dtype=bool) 4436 keep[:stop-start:step] = False 4437 slobj[axis] = slice(start, stop-numtodel) 4438 slobj2 = [slice(None)]*ndim 4439 slobj2[axis] = slice(start, stop) 4440 arr = arr[tuple(slobj2)] 4441 slobj2[axis] = keep 4442 new[tuple(slobj)] = arr[tuple(slobj2)] 4443 if wrap: 4444 return wrap(new) 4445 else: 4446 return new 4447 4448 if isinstance(obj, (int, integer)) and not isinstance(obj, bool): 4449 # optimization for a single value 4450 if (obj < -N or obj >= N): 4451 raise IndexError( 4452 "index %i is out of bounds for axis %i with " 4453 "size %i" % (obj, axis, N)) 4454 if (obj < 0): 4455 obj += N 4456 newshape[axis] -= 1 4457 new = empty(newshape, arr.dtype, arrorder) 4458 slobj[axis] = slice(None, obj) 4459 new[tuple(slobj)] = arr[tuple(slobj)] 4460 slobj[axis] = slice(obj, None) 4461 slobj2 = [slice(None)]*ndim 4462 slobj2[axis] = slice(obj+1, None) 4463 new[tuple(slobj)] = arr[tuple(slobj2)] 4464 else: 4465 _obj = obj 4466 obj = np.asarray(obj) 4467 if obj.size == 0 and not isinstance(_obj, np.ndarray): 4468 obj = obj.astype(intp) 4469 4470 if obj.dtype == bool: 4471 if obj.shape != (N,): 4472 raise ValueError('boolean array argument obj to delete ' 4473 'must be one dimensional and match the axis ' 4474 'length of {}'.format(N)) 4475 4476 # optimization, the other branch is slower 4477 keep = ~obj 4478 else: 4479 keep = ones(N, dtype=bool) 4480 keep[obj,] = False 4481 4482 slobj[axis] = keep 4483 new = arr[tuple(slobj)] 4484 4485 if wrap: 4486 return wrap(new) 4487 else: 4488 return new 4489 4490 4491def _insert_dispatcher(arr, obj, values, axis=None): 4492 return (arr, obj, values) 4493 4494 4495@array_function_dispatch(_insert_dispatcher) 4496def insert(arr, obj, values, axis=None): 4497 """ 4498 Insert values along the given axis before the given indices. 4499 4500 Parameters 4501 ---------- 4502 arr : array_like 4503 Input array. 4504 obj : int, slice or sequence of ints 4505 Object that defines the index or indices before which `values` is 4506 inserted. 4507 4508 .. versionadded:: 1.8.0 4509 4510 Support for multiple insertions when `obj` is a single scalar or a 4511 sequence with one element (similar to calling insert multiple 4512 times). 4513 values : array_like 4514 Values to insert into `arr`. If the type of `values` is different 4515 from that of `arr`, `values` is converted to the type of `arr`. 4516 `values` should be shaped so that ``arr[...,obj,...] = values`` 4517 is legal. 4518 axis : int, optional 4519 Axis along which to insert `values`. If `axis` is None then `arr` 4520 is flattened first. 4521 4522 Returns 4523 ------- 4524 out : ndarray 4525 A copy of `arr` with `values` inserted. Note that `insert` 4526 does not occur in-place: a new array is returned. If 4527 `axis` is None, `out` is a flattened array. 4528 4529 See Also 4530 -------- 4531 append : Append elements at the end of an array. 4532 concatenate : Join a sequence of arrays along an existing axis. 4533 delete : Delete elements from an array. 4534 4535 Notes 4536 ----- 4537 Note that for higher dimensional inserts `obj=0` behaves very different 4538 from `obj=[0]` just like `arr[:,0,:] = values` is different from 4539 `arr[:,[0],:] = values`. 4540 4541 Examples 4542 -------- 4543 >>> a = np.array([[1, 1], [2, 2], [3, 3]]) 4544 >>> a 4545 array([[1, 1], 4546 [2, 2], 4547 [3, 3]]) 4548 >>> np.insert(a, 1, 5) 4549 array([1, 5, 1, ..., 2, 3, 3]) 4550 >>> np.insert(a, 1, 5, axis=1) 4551 array([[1, 5, 1], 4552 [2, 5, 2], 4553 [3, 5, 3]]) 4554 4555 Difference between sequence and scalars: 4556 4557 >>> np.insert(a, [1], [[1],[2],[3]], axis=1) 4558 array([[1, 1, 1], 4559 [2, 2, 2], 4560 [3, 3, 3]]) 4561 >>> np.array_equal(np.insert(a, 1, [1, 2, 3], axis=1), 4562 ... np.insert(a, [1], [[1],[2],[3]], axis=1)) 4563 True 4564 4565 >>> b = a.flatten() 4566 >>> b 4567 array([1, 1, 2, 2, 3, 3]) 4568 >>> np.insert(b, [2, 2], [5, 6]) 4569 array([1, 1, 5, ..., 2, 3, 3]) 4570 4571 >>> np.insert(b, slice(2, 4), [5, 6]) 4572 array([1, 1, 5, ..., 2, 3, 3]) 4573 4574 >>> np.insert(b, [2, 2], [7.13, False]) # type casting 4575 array([1, 1, 7, ..., 2, 3, 3]) 4576 4577 >>> x = np.arange(8).reshape(2, 4) 4578 >>> idx = (1, 3) 4579 >>> np.insert(x, idx, 999, axis=1) 4580 array([[ 0, 999, 1, 2, 999, 3], 4581 [ 4, 999, 5, 6, 999, 7]]) 4582 4583 """ 4584 wrap = None 4585 if type(arr) is not ndarray: 4586 try: 4587 wrap = arr.__array_wrap__ 4588 except AttributeError: 4589 pass 4590 4591 arr = asarray(arr) 4592 ndim = arr.ndim 4593 arrorder = 'F' if arr.flags.fnc else 'C' 4594 if axis is None: 4595 if ndim != 1: 4596 arr = arr.ravel() 4597 # needed for np.matrix, which is still not 1d after being ravelled 4598 ndim = arr.ndim 4599 axis = ndim - 1 4600 else: 4601 axis = normalize_axis_index(axis, ndim) 4602 slobj = [slice(None)]*ndim 4603 N = arr.shape[axis] 4604 newshape = list(arr.shape) 4605 4606 if isinstance(obj, slice): 4607 # turn it into a range object 4608 indices = arange(*obj.indices(N), dtype=intp) 4609 else: 4610 # need to copy obj, because indices will be changed in-place 4611 indices = np.array(obj) 4612 if indices.dtype == bool: 4613 # See also delete 4614 # 2012-10-11, NumPy 1.8 4615 warnings.warn( 4616 "in the future insert will treat boolean arrays and " 4617 "array-likes as a boolean index instead of casting it to " 4618 "integer", FutureWarning, stacklevel=3) 4619 indices = indices.astype(intp) 4620 # Code after warning period: 4621 #if obj.ndim != 1: 4622 # raise ValueError('boolean array argument obj to insert ' 4623 # 'must be one dimensional') 4624 #indices = np.flatnonzero(obj) 4625 elif indices.ndim > 1: 4626 raise ValueError( 4627 "index array argument obj to insert must be one dimensional " 4628 "or scalar") 4629 if indices.size == 1: 4630 index = indices.item() 4631 if index < -N or index > N: 4632 raise IndexError( 4633 "index %i is out of bounds for axis %i with " 4634 "size %i" % (obj, axis, N)) 4635 if (index < 0): 4636 index += N 4637 4638 # There are some object array corner cases here, but we cannot avoid 4639 # that: 4640 values = array(values, copy=False, ndmin=arr.ndim, dtype=arr.dtype) 4641 if indices.ndim == 0: 4642 # broadcasting is very different here, since a[:,0,:] = ... behaves 4643 # very different from a[:,[0],:] = ...! This changes values so that 4644 # it works likes the second case. (here a[:,0:1,:]) 4645 values = np.moveaxis(values, 0, axis) 4646 numnew = values.shape[axis] 4647 newshape[axis] += numnew 4648 new = empty(newshape, arr.dtype, arrorder) 4649 slobj[axis] = slice(None, index) 4650 new[tuple(slobj)] = arr[tuple(slobj)] 4651 slobj[axis] = slice(index, index+numnew) 4652 new[tuple(slobj)] = values 4653 slobj[axis] = slice(index+numnew, None) 4654 slobj2 = [slice(None)] * ndim 4655 slobj2[axis] = slice(index, None) 4656 new[tuple(slobj)] = arr[tuple(slobj2)] 4657 if wrap: 4658 return wrap(new) 4659 return new 4660 elif indices.size == 0 and not isinstance(obj, np.ndarray): 4661 # Can safely cast the empty list to intp 4662 indices = indices.astype(intp) 4663 4664 indices[indices < 0] += N 4665 4666 numnew = len(indices) 4667 order = indices.argsort(kind='mergesort') # stable sort 4668 indices[order] += np.arange(numnew) 4669 4670 newshape[axis] += numnew 4671 old_mask = ones(newshape[axis], dtype=bool) 4672 old_mask[indices] = False 4673 4674 new = empty(newshape, arr.dtype, arrorder) 4675 slobj2 = [slice(None)]*ndim 4676 slobj[axis] = indices 4677 slobj2[axis] = old_mask 4678 new[tuple(slobj)] = values 4679 new[tuple(slobj2)] = arr 4680 4681 if wrap: 4682 return wrap(new) 4683 return new 4684 4685 4686def _append_dispatcher(arr, values, axis=None): 4687 return (arr, values) 4688 4689 4690@array_function_dispatch(_append_dispatcher) 4691def append(arr, values, axis=None): 4692 """ 4693 Append values to the end of an array. 4694 4695 Parameters 4696 ---------- 4697 arr : array_like 4698 Values are appended to a copy of this array. 4699 values : array_like 4700 These values are appended to a copy of `arr`. It must be of the 4701 correct shape (the same shape as `arr`, excluding `axis`). If 4702 `axis` is not specified, `values` can be any shape and will be 4703 flattened before use. 4704 axis : int, optional 4705 The axis along which `values` are appended. If `axis` is not 4706 given, both `arr` and `values` are flattened before use. 4707 4708 Returns 4709 ------- 4710 append : ndarray 4711 A copy of `arr` with `values` appended to `axis`. Note that 4712 `append` does not occur in-place: a new array is allocated and 4713 filled. If `axis` is None, `out` is a flattened array. 4714 4715 See Also 4716 -------- 4717 insert : Insert elements into an array. 4718 delete : Delete elements from an array. 4719 4720 Examples 4721 -------- 4722 >>> np.append([1, 2, 3], [[4, 5, 6], [7, 8, 9]]) 4723 array([1, 2, 3, ..., 7, 8, 9]) 4724 4725 When `axis` is specified, `values` must have the correct shape. 4726 4727 >>> np.append([[1, 2, 3], [4, 5, 6]], [[7, 8, 9]], axis=0) 4728 array([[1, 2, 3], 4729 [4, 5, 6], 4730 [7, 8, 9]]) 4731 >>> np.append([[1, 2, 3], [4, 5, 6]], [7, 8, 9], axis=0) 4732 Traceback (most recent call last): 4733 ... 4734 ValueError: all the input arrays must have same number of dimensions, but 4735 the array at index 0 has 2 dimension(s) and the array at index 1 has 1 4736 dimension(s) 4737 4738 """ 4739 arr = asanyarray(arr) 4740 if axis is None: 4741 if arr.ndim != 1: 4742 arr = arr.ravel() 4743 values = ravel(values) 4744 axis = arr.ndim-1 4745 return concatenate((arr, values), axis=axis) 4746 4747 4748def _digitize_dispatcher(x, bins, right=None): 4749 return (x, bins) 4750 4751 4752@array_function_dispatch(_digitize_dispatcher) 4753def digitize(x, bins, right=False): 4754 """ 4755 Return the indices of the bins to which each value in input array belongs. 4756 4757 ========= ============= ============================ 4758 `right` order of bins returned index `i` satisfies 4759 ========= ============= ============================ 4760 ``False`` increasing ``bins[i-1] <= x < bins[i]`` 4761 ``True`` increasing ``bins[i-1] < x <= bins[i]`` 4762 ``False`` decreasing ``bins[i-1] > x >= bins[i]`` 4763 ``True`` decreasing ``bins[i-1] >= x > bins[i]`` 4764 ========= ============= ============================ 4765 4766 If values in `x` are beyond the bounds of `bins`, 0 or ``len(bins)`` is 4767 returned as appropriate. 4768 4769 Parameters 4770 ---------- 4771 x : array_like 4772 Input array to be binned. Prior to NumPy 1.10.0, this array had to 4773 be 1-dimensional, but can now have any shape. 4774 bins : array_like 4775 Array of bins. It has to be 1-dimensional and monotonic. 4776 right : bool, optional 4777 Indicating whether the intervals include the right or the left bin 4778 edge. Default behavior is (right==False) indicating that the interval 4779 does not include the right edge. The left bin end is open in this 4780 case, i.e., bins[i-1] <= x < bins[i] is the default behavior for 4781 monotonically increasing bins. 4782 4783 Returns 4784 ------- 4785 indices : ndarray of ints 4786 Output array of indices, of same shape as `x`. 4787 4788 Raises 4789 ------ 4790 ValueError 4791 If `bins` is not monotonic. 4792 TypeError 4793 If the type of the input is complex. 4794 4795 See Also 4796 -------- 4797 bincount, histogram, unique, searchsorted 4798 4799 Notes 4800 ----- 4801 If values in `x` are such that they fall outside the bin range, 4802 attempting to index `bins` with the indices that `digitize` returns 4803 will result in an IndexError. 4804 4805 .. versionadded:: 1.10.0 4806 4807 `np.digitize` is implemented in terms of `np.searchsorted`. This means 4808 that a binary search is used to bin the values, which scales much better 4809 for larger number of bins than the previous linear search. It also removes 4810 the requirement for the input array to be 1-dimensional. 4811 4812 For monotonically _increasing_ `bins`, the following are equivalent:: 4813 4814 np.digitize(x, bins, right=True) 4815 np.searchsorted(bins, x, side='left') 4816 4817 Note that as the order of the arguments are reversed, the side must be too. 4818 The `searchsorted` call is marginally faster, as it does not do any 4819 monotonicity checks. Perhaps more importantly, it supports all dtypes. 4820 4821 Examples 4822 -------- 4823 >>> x = np.array([0.2, 6.4, 3.0, 1.6]) 4824 >>> bins = np.array([0.0, 1.0, 2.5, 4.0, 10.0]) 4825 >>> inds = np.digitize(x, bins) 4826 >>> inds 4827 array([1, 4, 3, 2]) 4828 >>> for n in range(x.size): 4829 ... print(bins[inds[n]-1], "<=", x[n], "<", bins[inds[n]]) 4830 ... 4831 0.0 <= 0.2 < 1.0 4832 4.0 <= 6.4 < 10.0 4833 2.5 <= 3.0 < 4.0 4834 1.0 <= 1.6 < 2.5 4835 4836 >>> x = np.array([1.2, 10.0, 12.4, 15.5, 20.]) 4837 >>> bins = np.array([0, 5, 10, 15, 20]) 4838 >>> np.digitize(x,bins,right=True) 4839 array([1, 2, 3, 4, 4]) 4840 >>> np.digitize(x,bins,right=False) 4841 array([1, 3, 3, 4, 5]) 4842 """ 4843 x = _nx.asarray(x) 4844 bins = _nx.asarray(bins) 4845 4846 # here for compatibility, searchsorted below is happy to take this 4847 if np.issubdtype(x.dtype, _nx.complexfloating): 4848 raise TypeError("x may not be complex") 4849 4850 mono = _monotonicity(bins) 4851 if mono == 0: 4852 raise ValueError("bins must be monotonically increasing or decreasing") 4853 4854 # this is backwards because the arguments below are swapped 4855 side = 'left' if right else 'right' 4856 if mono == -1: 4857 # reverse the bins, and invert the results 4858 return len(bins) - _nx.searchsorted(bins[::-1], x, side=side) 4859 else: 4860 return _nx.searchsorted(bins, x, side=side) 4861