[Numpy-svn] r5341 - in trunk/numpy/lib: . tests
numpy-svn@scip...
numpy-svn@scip...
Thu Jul 3 01:42:30 CDT 2008
Author: rkern
Date: 2008-07-03 01:42:28 -0500 (Thu, 03 Jul 2008)
New Revision: 5341
Added:
trunk/numpy/lib/stride_tricks.py
trunk/numpy/lib/tests/test_stride_tricks.py
Modified:
trunk/numpy/lib/__init__.py
Log:
ENH: Add broadcast_arrays() function to expose broadcasting to pure Python functions that cannot be made to be ufuncs.
Modified: trunk/numpy/lib/__init__.py
===================================================================
--- trunk/numpy/lib/__init__.py 2008-07-03 06:23:15 UTC (rev 5340)
+++ trunk/numpy/lib/__init__.py 2008-07-03 06:42:28 UTC (rev 5341)
@@ -5,6 +5,7 @@
from index_tricks import *
from function_base import *
from shape_base import *
+from stride_tricks import *
from twodim_base import *
from ufunclike import *
@@ -24,6 +25,7 @@
__all__ += index_tricks.__all__
__all__ += function_base.__all__
__all__ += shape_base.__all__
+__all__ += stride_tricks.__all__
__all__ += twodim_base.__all__
__all__ += ufunclike.__all__
__all__ += polynomial.__all__
Added: trunk/numpy/lib/stride_tricks.py
===================================================================
--- trunk/numpy/lib/stride_tricks.py 2008-07-03 06:23:15 UTC (rev 5340)
+++ trunk/numpy/lib/stride_tricks.py 2008-07-03 06:42:28 UTC (rev 5341)
@@ -0,0 +1,109 @@
+""" Utilities that manipulate strides to achieve desirable effects.
+"""
+import numpy as np
+
+__all__ = ['broadcast_arrays']
+
+class DummyArray(object):
+ """ Dummy object that just exists to hang __array_interface__ dictionaries
+ and possibly keep alive a reference to a base array.
+ """
+ def __init__(self, interface, base=None):
+ self.__array_interface__ = interface
+ self.base = base
+
+def as_strided(x, shape=None, strides=None):
+ """ Make an ndarray from the given array with the given shape and strides.
+ """
+ interface = dict(x.__array_interface__)
+ if shape is not None:
+ interface['shape'] = tuple(shape)
+ if strides is not None:
+ interface['strides'] = tuple(strides)
+ return np.asarray(DummyArray(interface, base=x))
+
+def broadcast_arrays(*args):
+ """ Broadcast any number of arrays against each other.
+
+ Parameters
+ ----------
+ *args : arrays
+
+ Returns
+ -------
+ broadcasted : list of arrays
+ These arrays are views on the original arrays. They are typically not
+ contiguous. Furthermore, more than one element of a broadcasted array
+ may refer to a single memory location. If you need to write to the
+ arrays, make copies first.
+
+ Examples
+ --------
+ >>> x = np.array([[1,2,3]])
+ >>> y = np.array([[1],[2],[3]])
+ >>> np.broadcast_arrays(x, y)
+ [array([[1, 2, 3],
+ [1, 2, 3],
+ [1, 2, 3]]), array([[1, 1, 1],
+ [2, 2, 2],
+ [3, 3, 3]])]
+
+ Here is a useful idiom for getting contiguous copies instead of
+ non-contiguous views.
+
+ >>> map(np.array, np.broadcast_arrays(x, y))
+ [array([[1, 2, 3],
+ [1, 2, 3],
+ [1, 2, 3]]), array([[1, 1, 1],
+ [2, 2, 2],
+ [3, 3, 3]])]
+
+ """
+ args = map(np.asarray, args)
+ shapes = [x.shape for x in args]
+ if len(set(shapes)) == 1:
+ # Common case where nothing needs to be broadcasted.
+ return args
+ shapes = [list(s) for s in shapes]
+ strides = [list(x.strides) for x in args]
+ nds = [len(s) for s in shapes]
+ biggest = max(nds)
+ # Go through each array and prepend dimensions of length 1 to each of the
+ # shapes in order to make the number of dimensions equal.
+ for i in range(len(args)):
+ diff = biggest - nds[i]
+ if diff > 0:
+ shapes[i] = [1] * diff + shapes[i]
+ strides[i] = [0] * diff + strides[i]
+ # Chech each dimension for compatibility. A dimension length of 1 is
+ # accepted as compatible with any other length.
+ common_shape = []
+ for axis in range(biggest):
+ lengths = [s[axis] for s in shapes]
+ unique = set(lengths + [1])
+ if len(unique) > 2:
+ # There must be at least two non-1 lengths for this axis.
+ raise ValueError("shape mismatch: two or more arrays have "
+ "incompatible dimensions on axis %r." % (axis,))
+ elif len(unique) == 2:
+ # There is exactly one non-1 length. The common shape will take this
+ # value.
+ unique.remove(1)
+ new_length = unique.pop()
+ common_shape.append(new_length)
+ # For each array, if this axis is being broadcasted from a length of
+ # 1, then set its stride to 0 so that it repeats its data.
+ for i in range(len(args)):
+ if shapes[i][axis] == 1:
+ shapes[i][axis] = new_length
+ strides[i][axis] = 0
+ else:
+ # Every array has a length of 1 on this axis. Strides can be left
+ # alone as nothing is broadcasted.
+ common_shape.append(1)
+
+ # Construct the new arrays.
+ broadcasted = [as_strided(x, shape=sh, strides=st) for (x,sh,st) in
+ zip(args, shapes, strides)]
+ return broadcasted
+
Property changes on: trunk/numpy/lib/stride_tricks.py
___________________________________________________________________
Name: svn:eol-style
+ native
Added: trunk/numpy/lib/tests/test_stride_tricks.py
===================================================================
--- trunk/numpy/lib/tests/test_stride_tricks.py 2008-07-03 06:23:15 UTC (rev 5340)
+++ trunk/numpy/lib/tests/test_stride_tricks.py 2008-07-03 06:42:28 UTC (rev 5341)
@@ -0,0 +1,206 @@
+from nose.tools import assert_raises
+import numpy as np
+from numpy.testing import assert_array_equal
+
+from numpy.lib.stride_tricks import broadcast_arrays
+
+
+def assert_shapes_correct(input_shapes, expected_shape):
+ """ Broadcast a list of arrays with the given input shapes and check the
+ common output shape.
+ """
+ inarrays = [np.zeros(s) for s in input_shapes]
+ outarrays = broadcast_arrays(*inarrays)
+ outshapes = [a.shape for a in outarrays]
+ expected = [expected_shape] * len(inarrays)
+ assert outshapes == expected
+
+def assert_incompatible_shapes_raise(input_shapes):
+ """ Broadcast a list of arrays with the given (incompatible) input shapes
+ and check that they raise a ValueError.
+ """
+ inarrays = [np.zeros(s) for s in input_shapes]
+ assert_raises(ValueError, broadcast_arrays, *inarrays)
+
+def assert_same_as_ufunc(shape0, shape1, transposed=False, flipped=False):
+ """ Broadcast two shapes against each other and check that the data layout
+ is the same as if a ufunc did the broadcasting.
+ """
+ x0 = np.zeros(shape0, dtype=int)
+ # Note that multiply.reduce's identity element is 1.0, so when shape1==(),
+ # this gives the desired n==1.
+ n = int(np.multiply.reduce(shape1))
+ x1 = np.arange(n).reshape(shape1)
+ if transposed:
+ x0 = x0.T
+ x1 = x1.T
+ if flipped:
+ x0 = x0[::-1]
+ x1 = x1[::-1]
+ # Use the add ufunc to do the broadcasting. Since we're adding 0s to x1, the
+ # result should be exactly the same as the broadcasted view of x1.
+ y = x0 + x1
+ b0, b1 = broadcast_arrays(x0, x1)
+ assert_array_equal(y, b1)
+
+
+def test_same():
+ x = np.arange(10)
+ y = np.arange(10)
+ bx, by = broadcast_arrays(x, y)
+ assert_array_equal(x, bx)
+ assert_array_equal(y, by)
+
+def test_one_off():
+ x = np.array([[1,2,3]])
+ y = np.array([[1],[2],[3]])
+ bx, by = broadcast_arrays(x, y)
+ bx0 = np.array([[1,2,3],[1,2,3],[1,2,3]])
+ by0 = bx0.T
+ assert_array_equal(bx0, bx)
+ assert_array_equal(by0, by)
+
+def test_same_input_shapes():
+ """ Check that the final shape is just the input shape.
+ """
+ data = [
+ (),
+ (1,),
+ (3,),
+ (0,1),
+ (0,3),
+ (1,0),
+ (3,0),
+ (1,3),
+ (3,1),
+ (3,3),
+ ]
+ for shape in data:
+ input_shapes = [shape]
+ # Single input.
+ yield assert_shapes_correct, input_shapes, shape
+ # Double input.
+ input_shapes2 = [shape, shape]
+ yield assert_shapes_correct, input_shapes2, shape
+ # Triple input.
+ input_shapes3 = [shape, shape, shape]
+ yield assert_shapes_correct, input_shapes3, shape
+
+def test_two_compatible_by_ones_input_shapes():
+ """ Check that two different input shapes (of the same length but some have
+ 1s) broadcast to the correct shape.
+ """
+ data = [
+ [[(1,), (3,)], (3,)],
+ [[(1,3), (3,3)], (3,3)],
+ [[(3,1), (3,3)], (3,3)],
+ [[(1,3), (3,1)], (3,3)],
+ [[(1,1), (3,3)], (3,3)],
+ [[(1,1), (1,3)], (1,3)],
+ [[(1,1), (3,1)], (3,1)],
+ [[(1,0), (0,0)], (0,0)],
+ [[(0,1), (0,0)], (0,0)],
+ [[(1,0), (0,1)], (0,0)],
+ [[(1,1), (0,0)], (0,0)],
+ [[(1,1), (1,0)], (1,0)],
+ [[(1,1), (0,1)], (0,1)],
+ ]
+ for input_shapes, expected_shape in data:
+ yield assert_shapes_correct, input_shapes, expected_shape
+ # Reverse the input shapes since broadcasting should be symmetric.
+ yield assert_shapes_correct, input_shapes[::-1], expected_shape
+
+def test_two_compatible_by_prepending_ones_input_shapes():
+ """ Check that two different input shapes (of different lengths) broadcast
+ to the correct shape.
+ """
+ data = [
+ [[(), (3,)], (3,)],
+ [[(3,), (3,3)], (3,3)],
+ [[(3,), (3,1)], (3,3)],
+ [[(1,), (3,3)], (3,3)],
+ [[(), (3,3)], (3,3)],
+ [[(1,1), (3,)], (1,3)],
+ [[(1,), (3,1)], (3,1)],
+ [[(1,), (1,3)], (1,3)],
+ [[(), (1,3)], (1,3)],
+ [[(), (3,1)], (3,1)],
+ [[(), (0,)], (0,)],
+ [[(0,), (0,0)], (0,0)],
+ [[(0,), (0,1)], (0,0)],
+ [[(1,), (0,0)], (0,0)],
+ [[(), (0,0)], (0,0)],
+ [[(1,1), (0,)], (1,0)],
+ [[(1,), (0,1)], (0,1)],
+ [[(1,), (1,0)], (1,0)],
+ [[(), (1,0)], (1,0)],
+ [[(), (0,1)], (0,1)],
+ ]
+ for input_shapes, expected_shape in data:
+ yield assert_shapes_correct, input_shapes, expected_shape
+ # Reverse the input shapes since broadcasting should be symmetric.
+ yield assert_shapes_correct, input_shapes[::-1], expected_shape
+
+def test_incompatible_shapes_raise_valueerror():
+ """ Check that a ValueError is raised for incompatible shapes.
+ """
+ data = [
+ [(3,), (4,)],
+ [(2,3), (2,)],
+ [(3,), (3,), (4,)],
+ [(1,3,4), (2,3,3)],
+ ]
+ for input_shapes in data:
+ yield assert_incompatible_shapes_raise, input_shapes
+ # Reverse the input shapes since broadcasting should be symmetric.
+ yield assert_incompatible_shapes_raise, input_shapes[::-1]
+
+def test_same_as_ufunc():
+ """ Check that the data layout is the same as if a ufunc did the operation.
+ """
+ data = [
+ [[(1,), (3,)], (3,)],
+ [[(1,3), (3,3)], (3,3)],
+ [[(3,1), (3,3)], (3,3)],
+ [[(1,3), (3,1)], (3,3)],
+ [[(1,1), (3,3)], (3,3)],
+ [[(1,1), (1,3)], (1,3)],
+ [[(1,1), (3,1)], (3,1)],
+ [[(1,0), (0,0)], (0,0)],
+ [[(0,1), (0,0)], (0,0)],
+ [[(1,0), (0,1)], (0,0)],
+ [[(1,1), (0,0)], (0,0)],
+ [[(1,1), (1,0)], (1,0)],
+ [[(1,1), (0,1)], (0,1)],
+ [[(), (3,)], (3,)],
+ [[(3,), (3,3)], (3,3)],
+ [[(3,), (3,1)], (3,3)],
+ [[(1,), (3,3)], (3,3)],
+ [[(), (3,3)], (3,3)],
+ [[(1,1), (3,)], (1,3)],
+ [[(1,), (3,1)], (3,1)],
+ [[(1,), (1,3)], (1,3)],
+ [[(), (1,3)], (1,3)],
+ [[(), (3,1)], (3,1)],
+ [[(), (0,)], (0,)],
+ [[(0,), (0,0)], (0,0)],
+ [[(0,), (0,1)], (0,0)],
+ [[(1,), (0,0)], (0,0)],
+ [[(), (0,0)], (0,0)],
+ [[(1,1), (0,)], (1,0)],
+ [[(1,), (0,1)], (0,1)],
+ [[(1,), (1,0)], (1,0)],
+ [[(), (1,0)], (1,0)],
+ [[(), (0,1)], (0,1)],
+ ]
+ for input_shapes, expected_shape in data:
+ yield assert_same_as_ufunc, input_shapes[0], input_shapes[1]
+ # Reverse the input shapes since broadcasting should be symmetric.
+ yield assert_same_as_ufunc, input_shapes[1], input_shapes[0]
+ # Try them transposed, too.
+ yield assert_same_as_ufunc, input_shapes[0], input_shapes[1], True
+ # ... and flipped for non-rank-0 inputs in order to test negative
+ # strides.
+ if () not in input_shapes:
+ yield assert_same_as_ufunc, input_shapes[0], input_shapes[1], False, True
+ yield assert_same_as_ufunc, input_shapes[0], input_shapes[1], True, True
Property changes on: trunk/numpy/lib/tests/test_stride_tricks.py
___________________________________________________________________
Name: svn:eol-style
+ native
More information about the Numpy-svn
mailing list