_stride_tricks_impl.py 19 KB

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