1""" 2Utilities that manipulate strides to achieve desirable effects. 3 4An explanation of strides can be found in the "ndarray.rst" file in the 5NumPy reference guide. 6 7""" 8import numpy as np 9from numpy.core.numeric import normalize_axis_tuple 10from numpy.core.overrides import array_function_dispatch, set_module 11 12__all__ = ['broadcast_to', 'broadcast_arrays', 'broadcast_shapes'] 13 14 15class DummyArray: 16 """Dummy object that just exists to hang __array_interface__ dictionaries 17 and possibly keep alive a reference to a base array. 18 """ 19 20 def __init__(self, interface, base=None): 21 self.__array_interface__ = interface 22 self.base = base 23 24 25def _maybe_view_as_subclass(original_array, new_array): 26 if type(original_array) is not type(new_array): 27 # if input was an ndarray subclass and subclasses were OK, 28 # then view the result as that subclass. 29 new_array = new_array.view(type=type(original_array)) 30 # Since we have done something akin to a view from original_array, we 31 # should let the subclass finalize (if it has it implemented, i.e., is 32 # not None). 33 if new_array.__array_finalize__: 34 new_array.__array_finalize__(original_array) 35 return new_array 36 37 38def as_strided(x, shape=None, strides=None, subok=False, writeable=True): 39 """ 40 Create a view into the array with the given shape and strides. 41 42 .. warning:: This function has to be used with extreme care, see notes. 43 44 Parameters 45 ---------- 46 x : ndarray 47 Array to create a new. 48 shape : sequence of int, optional 49 The shape of the new array. Defaults to ``x.shape``. 50 strides : sequence of int, optional 51 The strides of the new array. Defaults to ``x.strides``. 52 subok : bool, optional 53 .. versionadded:: 1.10 54 55 If True, subclasses are preserved. 56 writeable : bool, optional 57 .. versionadded:: 1.12 58 59 If set to False, the returned array will always be readonly. 60 Otherwise it will be writable if the original array was. It 61 is advisable to set this to False if possible (see Notes). 62 63 Returns 64 ------- 65 view : ndarray 66 67 See also 68 -------- 69 broadcast_to : broadcast an array to a given shape. 70 reshape : reshape an array. 71 lib.stride_tricks.sliding_window_view : 72 userfriendly and safe function for the creation of sliding window views. 73 74 Notes 75 ----- 76 ``as_strided`` creates a view into the array given the exact strides 77 and shape. This means it manipulates the internal data structure of 78 ndarray and, if done incorrectly, the array elements can point to 79 invalid memory and can corrupt results or crash your program. 80 It is advisable to always use the original ``x.strides`` when 81 calculating new strides to avoid reliance on a contiguous memory 82 layout. 83 84 Furthermore, arrays created with this function often contain self 85 overlapping memory, so that two elements are identical. 86 Vectorized write operations on such arrays will typically be 87 unpredictable. They may even give different results for small, large, 88 or transposed arrays. 89 Since writing to these arrays has to be tested and done with great 90 care, you may want to use ``writeable=False`` to avoid accidental write 91 operations. 92 93 For these reasons it is advisable to avoid ``as_strided`` when 94 possible. 95 """ 96 # first convert input to array, possibly keeping subclass 97 x = np.array(x, copy=False, subok=subok) 98 interface = dict(x.__array_interface__) 99 if shape is not None: 100 interface['shape'] = tuple(shape) 101 if strides is not None: 102 interface['strides'] = tuple(strides) 103 104 array = np.asarray(DummyArray(interface, base=x)) 105 # The route via `__interface__` does not preserve structured 106 # dtypes. Since dtype should remain unchanged, we set it explicitly. 107 array.dtype = x.dtype 108 109 view = _maybe_view_as_subclass(x, array) 110 111 if view.flags.writeable and not writeable: 112 view.flags.writeable = False 113 114 return view 115 116 117def _sliding_window_view_dispatcher(x, window_shape, axis=None, *, 118 subok=None, writeable=None): 119 return (x,) 120 121 122@array_function_dispatch(_sliding_window_view_dispatcher) 123def sliding_window_view(x, window_shape, axis=None, *, 124 subok=False, writeable=False): 125 """ 126 Create a sliding window view into the array with the given window shape. 127 128 Also known as rolling or moving window, the window slides across all 129 dimensions of the array and extracts subsets of the array at all window 130 positions. 131 132 .. versionadded:: 1.20.0 133 134 Parameters 135 ---------- 136 x : array_like 137 Array to create the sliding window view from. 138 window_shape : int or tuple of int 139 Size of window over each axis that takes part in the sliding window. 140 If `axis` is not present, must have same length as the number of input 141 array dimensions. Single integers `i` are treated as if they were the 142 tuple `(i,)`. 143 axis : int or tuple of int, optional 144 Axis or axes along which the sliding window is applied. 145 By default, the sliding window is applied to all axes and 146 `window_shape[i]` will refer to axis `i` of `x`. 147 If `axis` is given as a `tuple of int`, `window_shape[i]` will refer to 148 the axis `axis[i]` of `x`. 149 Single integers `i` are treated as if they were the tuple `(i,)`. 150 subok : bool, optional 151 If True, sub-classes will be passed-through, otherwise the returned 152 array will be forced to be a base-class array (default). 153 writeable : bool, optional 154 When true, allow writing to the returned view. The default is false, 155 as this should be used with caution: the returned view contains the 156 same memory location multiple times, so writing to one location will 157 cause others to change. 158 159 Returns 160 ------- 161 view : ndarray 162 Sliding window view of the array. The sliding window dimensions are 163 inserted at the end, and the original dimensions are trimmed as 164 required by the size of the sliding window. 165 That is, ``view.shape = x_shape_trimmed + window_shape``, where 166 ``x_shape_trimmed`` is ``x.shape`` with every entry reduced by one less 167 than the corresponding window size. 168 169 See Also 170 -------- 171 lib.stride_tricks.as_strided: A lower-level and less safe routine for 172 creating arbitrary views from custom shape and strides. 173 broadcast_to: broadcast an array to a given shape. 174 175 Notes 176 ----- 177 For many applications using a sliding window view can be convenient, but 178 potentially very slow. Often specialized solutions exist, for example: 179 180 - `scipy.signal.fftconvolve` 181 182 - filtering functions in `scipy.ndimage` 183 184 - moving window functions provided by 185 `bottleneck <https://github.com/pydata/bottleneck>`_. 186 187 As a rough estimate, a sliding window approach with an input size of `N` 188 and a window size of `W` will scale as `O(N*W)` where frequently a special 189 algorithm can achieve `O(N)`. That means that the sliding window variant 190 for a window size of 100 can be a 100 times slower than a more specialized 191 version. 192 193 Nevertheless, for small window sizes, when no custom algorithm exists, or 194 as a prototyping and developing tool, this function can be a good solution. 195 196 Examples 197 -------- 198 >>> x = np.arange(6) 199 >>> x.shape 200 (6,) 201 >>> v = sliding_window_view(x, 3) 202 >>> v.shape 203 (4, 3) 204 >>> v 205 array([[0, 1, 2], 206 [1, 2, 3], 207 [2, 3, 4], 208 [3, 4, 5]]) 209 210 This also works in more dimensions, e.g. 211 212 >>> i, j = np.ogrid[:3, :4] 213 >>> x = 10*i + j 214 >>> x.shape 215 (3, 4) 216 >>> x 217 array([[ 0, 1, 2, 3], 218 [10, 11, 12, 13], 219 [20, 21, 22, 23]]) 220 >>> shape = (2,2) 221 >>> v = sliding_window_view(x, shape) 222 >>> v.shape 223 (2, 3, 2, 2) 224 >>> v 225 array([[[[ 0, 1], 226 [10, 11]], 227 [[ 1, 2], 228 [11, 12]], 229 [[ 2, 3], 230 [12, 13]]], 231 [[[10, 11], 232 [20, 21]], 233 [[11, 12], 234 [21, 22]], 235 [[12, 13], 236 [22, 23]]]]) 237 238 The axis can be specified explicitly: 239 240 >>> v = sliding_window_view(x, 3, 0) 241 >>> v.shape 242 (1, 4, 3) 243 >>> v 244 array([[[ 0, 10, 20], 245 [ 1, 11, 21], 246 [ 2, 12, 22], 247 [ 3, 13, 23]]]) 248 249 The same axis can be used several times. In that case, every use reduces 250 the corresponding original dimension: 251 252 >>> v = sliding_window_view(x, (2, 3), (1, 1)) 253 >>> v.shape 254 (3, 1, 2, 3) 255 >>> v 256 array([[[[ 0, 1, 2], 257 [ 1, 2, 3]]], 258 [[[10, 11, 12], 259 [11, 12, 13]]], 260 [[[20, 21, 22], 261 [21, 22, 23]]]]) 262 263 Combining with stepped slicing (`::step`), this can be used to take sliding 264 views which skip elements: 265 266 >>> x = np.arange(7) 267 >>> sliding_window_view(x, 5)[:, ::2] 268 array([[0, 2, 4], 269 [1, 3, 5], 270 [2, 4, 6]]) 271 272 or views which move by multiple elements 273 274 >>> x = np.arange(7) 275 >>> sliding_window_view(x, 3)[::2, :] 276 array([[0, 1, 2], 277 [2, 3, 4], 278 [4, 5, 6]]) 279 280 A common application of `sliding_window_view` is the calculation of running 281 statistics. The simplest example is the 282 `moving average <https://en.wikipedia.org/wiki/Moving_average>`_: 283 284 >>> x = np.arange(6) 285 >>> x.shape 286 (6,) 287 >>> v = sliding_window_view(x, 3) 288 >>> v.shape 289 (4, 3) 290 >>> v 291 array([[0, 1, 2], 292 [1, 2, 3], 293 [2, 3, 4], 294 [3, 4, 5]]) 295 >>> moving_average = v.mean(axis=-1) 296 >>> moving_average 297 array([1., 2., 3., 4.]) 298 299 Note that a sliding window approach is often **not** optimal (see Notes). 300 """ 301 window_shape = (tuple(window_shape) 302 if np.iterable(window_shape) 303 else (window_shape,)) 304 # first convert input to array, possibly keeping subclass 305 x = np.array(x, copy=False, subok=subok) 306 307 window_shape_array = np.array(window_shape) 308 if np.any(window_shape_array < 0): 309 raise ValueError('`window_shape` cannot contain negative values') 310 311 if axis is None: 312 axis = tuple(range(x.ndim)) 313 if len(window_shape) != len(axis): 314 raise ValueError(f'Since axis is `None`, must provide ' 315 f'window_shape for all dimensions of `x`; ' 316 f'got {len(window_shape)} window_shape elements ' 317 f'and `x.ndim` is {x.ndim}.') 318 else: 319 axis = normalize_axis_tuple(axis, x.ndim, allow_duplicate=True) 320 if len(window_shape) != len(axis): 321 raise ValueError(f'Must provide matching length window_shape and ' 322 f'axis; got {len(window_shape)} window_shape ' 323 f'elements and {len(axis)} axes elements.') 324 325 out_strides = x.strides + tuple(x.strides[ax] for ax in axis) 326 327 # note: same axis can be windowed repeatedly 328 x_shape_trimmed = list(x.shape) 329 for ax, dim in zip(axis, window_shape): 330 if x_shape_trimmed[ax] < dim: 331 raise ValueError( 332 'window shape cannot be larger than input array shape') 333 x_shape_trimmed[ax] -= dim - 1 334 out_shape = tuple(x_shape_trimmed) + window_shape 335 return as_strided(x, strides=out_strides, shape=out_shape, 336 subok=subok, writeable=writeable) 337 338 339def _broadcast_to(array, shape, subok, readonly): 340 shape = tuple(shape) if np.iterable(shape) else (shape,) 341 array = np.array(array, copy=False, subok=subok) 342 if not shape and array.shape: 343 raise ValueError('cannot broadcast a non-scalar to a scalar array') 344 if any(size < 0 for size in shape): 345 raise ValueError('all elements of broadcast shape must be non-' 346 'negative') 347 extras = [] 348 it = np.nditer( 349 (array,), flags=['multi_index', 'refs_ok', 'zerosize_ok'] + extras, 350 op_flags=['readonly'], itershape=shape, order='C') 351 with it: 352 # never really has writebackifcopy semantics 353 broadcast = it.itviews[0] 354 result = _maybe_view_as_subclass(array, broadcast) 355 # In a future version this will go away 356 if not readonly and array.flags._writeable_no_warn: 357 result.flags.writeable = True 358 result.flags._warn_on_write = True 359 return result 360 361 362def _broadcast_to_dispatcher(array, shape, subok=None): 363 return (array,) 364 365 366@array_function_dispatch(_broadcast_to_dispatcher, module='numpy') 367def broadcast_to(array, shape, subok=False): 368 """Broadcast an array to a new shape. 369 370 Parameters 371 ---------- 372 array : array_like 373 The array to broadcast. 374 shape : tuple 375 The shape of the desired array. 376 subok : bool, optional 377 If True, then sub-classes will be passed-through, otherwise 378 the returned array will be forced to be a base-class array (default). 379 380 Returns 381 ------- 382 broadcast : array 383 A readonly view on the original array with the given shape. It is 384 typically not contiguous. Furthermore, more than one element of a 385 broadcasted array may refer to a single memory location. 386 387 Raises 388 ------ 389 ValueError 390 If the array is not compatible with the new shape according to NumPy's 391 broadcasting rules. 392 393 See Also 394 -------- 395 broadcast 396 broadcast_arrays 397 broadcast_shapes 398 399 Notes 400 ----- 401 .. versionadded:: 1.10.0 402 403 Examples 404 -------- 405 >>> x = np.array([1, 2, 3]) 406 >>> np.broadcast_to(x, (3, 3)) 407 array([[1, 2, 3], 408 [1, 2, 3], 409 [1, 2, 3]]) 410 """ 411 return _broadcast_to(array, shape, subok=subok, readonly=True) 412 413 414def _broadcast_shape(*args): 415 """Returns the shape of the arrays that would result from broadcasting the 416 supplied arrays against each other. 417 """ 418 # use the old-iterator because np.nditer does not handle size 0 arrays 419 # consistently 420 b = np.broadcast(*args[:32]) 421 # unfortunately, it cannot handle 32 or more arguments directly 422 for pos in range(32, len(args), 31): 423 # ironically, np.broadcast does not properly handle np.broadcast 424 # objects (it treats them as scalars) 425 # use broadcasting to avoid allocating the full array 426 b = broadcast_to(0, b.shape) 427 b = np.broadcast(b, *args[pos:(pos + 31)]) 428 return b.shape 429 430 431@set_module('numpy') 432def broadcast_shapes(*args): 433 """ 434 Broadcast the input shapes into a single shape. 435 436 :ref:`Learn more about broadcasting here <basics.broadcasting>`. 437 438 .. versionadded:: 1.20.0 439 440 Parameters 441 ---------- 442 `*args` : tuples of ints, or ints 443 The shapes to be broadcast against each other. 444 445 Returns 446 ------- 447 tuple 448 Broadcasted shape. 449 450 Raises 451 ------ 452 ValueError 453 If the shapes are not compatible and cannot be broadcast according 454 to NumPy's broadcasting rules. 455 456 See Also 457 -------- 458 broadcast 459 broadcast_arrays 460 broadcast_to 461 462 Examples 463 -------- 464 >>> np.broadcast_shapes((1, 2), (3, 1), (3, 2)) 465 (3, 2) 466 467 >>> np.broadcast_shapes((6, 7), (5, 6, 1), (7,), (5, 1, 7)) 468 (5, 6, 7) 469 """ 470 arrays = [np.empty(x, dtype=[]) for x in args] 471 return _broadcast_shape(*arrays) 472 473 474def _broadcast_arrays_dispatcher(*args, subok=None): 475 return args 476 477 478@array_function_dispatch(_broadcast_arrays_dispatcher, module='numpy') 479def broadcast_arrays(*args, subok=False): 480 """ 481 Broadcast any number of arrays against each other. 482 483 Parameters 484 ---------- 485 `*args` : array_likes 486 The arrays to broadcast. 487 488 subok : bool, optional 489 If True, then sub-classes will be passed-through, otherwise 490 the returned arrays will be forced to be a base-class array (default). 491 492 Returns 493 ------- 494 broadcasted : list of arrays 495 These arrays are views on the original arrays. They are typically 496 not contiguous. Furthermore, more than one element of a 497 broadcasted array may refer to a single memory location. If you need 498 to write to the arrays, make copies first. While you can set the 499 ``writable`` flag True, writing to a single output value may end up 500 changing more than one location in the output array. 501 502 .. deprecated:: 1.17 503 The output is currently marked so that if written to, a deprecation 504 warning will be emitted. A future version will set the 505 ``writable`` flag False so writing to it will raise an error. 506 507 See Also 508 -------- 509 broadcast 510 broadcast_to 511 broadcast_shapes 512 513 Examples 514 -------- 515 >>> x = np.array([[1,2,3]]) 516 >>> y = np.array([[4],[5]]) 517 >>> np.broadcast_arrays(x, y) 518 [array([[1, 2, 3], 519 [1, 2, 3]]), array([[4, 4, 4], 520 [5, 5, 5]])] 521 522 Here is a useful idiom for getting contiguous copies instead of 523 non-contiguous views. 524 525 >>> [np.array(a) for a in np.broadcast_arrays(x, y)] 526 [array([[1, 2, 3], 527 [1, 2, 3]]), array([[4, 4, 4], 528 [5, 5, 5]])] 529 530 """ 531 # nditer is not used here to avoid the limit of 32 arrays. 532 # Otherwise, something like the following one-liner would suffice: 533 # return np.nditer(args, flags=['multi_index', 'zerosize_ok'], 534 # order='C').itviews 535 536 args = [np.array(_m, copy=False, subok=subok) for _m in args] 537 538 shape = _broadcast_shape(*args) 539 540 if all(array.shape == shape for array in args): 541 # Common case where nothing needs to be broadcasted. 542 return args 543 544 return [_broadcast_to(array, shape, subok=subok, readonly=False) 545 for array in args] 546