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