From e37a90c8bd1bf99f89b537ce47492095d6ae147c Mon Sep 17 00:00:00 2001 From: "Per A. Brodtkorb" Date: Tue, 13 Dec 2011 23:01:22 +0100 Subject: ENH: enhance meshgrid to generate 3D grids, sparse grids, matrix indexing. --- numpy/lib/function_base.py | 140 +++++++++++++++++++++++++--------- numpy/lib/tests/test_function_base.py | 18 +++++ 2 files changed, 124 insertions(+), 34 deletions(-) diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py index f254bbacf..86ae3ab08 100644 --- a/numpy/lib/function_base.py +++ b/numpy/lib/function_base.py @@ -6,7 +6,7 @@ __all__ = ['select', 'piecewise', 'trim_zeros', 'copy', 'iterable', 'histogram', 'histogramdd', 'bincount', 'digitize', 'cov', 'corrcoef', 'msort', 'median', 'sinc', 'hamming', 'hanning', 'bartlett', 'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc', 'add_docstring', - 'meshgrid', 'delete', 'insert', 'append', 'interp', 'add_newdoc_ufunc'] + 'meshgrid', 'ndgrid', 'delete', 'insert', 'append', 'interp', 'add_newdoc_ufunc'] import warnings import types @@ -3200,62 +3200,134 @@ def add_newdoc(place, obj, doc): pass -# From matplotlib -def meshgrid(x,y): +# Based on scitools meshgrid +def meshgrid(*xi, **kwargs): """ - Return coordinate matrices from two coordinate vectors. + Return coordinate matrices from two or more coordinate vectors. + + Make N-D coordinate arrays for vectorized evaluations of + N-D scalar/vector fields over N-D grids, given + one-dimensional coordinate arrays x1, x2,..., xn. Parameters ---------- - x, y : ndarray - Two 1-D arrays representing the x and y coordinates of a grid. + x1, x2,..., xn : array_like + 1-D arrays representing the coordinates of a grid. + indexing : 'xy' or 'ij' (optional) + cartesian ('xy', default) or matrix ('ij') indexing of output + sparse : True or False (default) (optional) + If True a sparse grid is returned in order to conserve memory. + copy : True (default) or False (optional) + If False a view into the original arrays are returned in order to + conserve memory. Please note that sparse=False, copy=False will likely + return non-contiguous arrays. 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. Returns ------- - X, Y : ndarray - For vectors `x`, `y` with lengths ``Nx=len(x)`` and ``Ny=len(y)``, - return `X`, `Y` where `X` and `Y` are ``(Ny, Nx)`` shaped arrays - with the elements of `x` and y repeated to fill the matrix along - the first dimension for `x`, the second for `y`. + X1, X2,..., XN : ndarray + For vectors `x1`, `x2`,..., 'xn' with lengths ``Ni=len(xi)`` , + return ``(N1, N2, N3,...Nn)`` shaped arrays if indexing='ij' + or ``(N2, N1, N3,...Nn)`` shaped arrays if indexing='xy' + with the elements of `xi` repeated to fill the matrix along + the first dimension for `x1`, the second for `x2` and so on. + + Notes + ----- + This function supports both indexing conventions through the indexing keyword + argument. Giving the string 'ij' returns a meshgrid with matrix indexing, + while 'xy' returns a meshgrid with Cartesian indexing. The difference is + illustrated by the following code snippet: + + xv, yv = meshgrid(x, y, sparse=False, indexing='ij') + for i in range(nx): + for j in range(ny): + # treat xv[i,j], yv[i,j] + + xv, yv = meshgrid(x, y, sparse=False, indexing='xy') + for i in range(nx): + for j in range(ny): + # treat xv[j,i], yv[j,i] See Also -------- index_tricks.mgrid : Construct a multi-dimensional "meshgrid" - using indexing notation. + using indexing notation. index_tricks.ogrid : Construct an open multi-dimensional "meshgrid" - using indexing notation. + using indexing notation. Examples -------- - >>> X, Y = np.meshgrid([1,2,3], [4,5,6,7]) - >>> X - array([[1, 2, 3], - [1, 2, 3], - [1, 2, 3], - [1, 2, 3]]) - >>> Y - array([[4, 4, 4], - [5, 5, 5], - [6, 6, 6], - [7, 7, 7]]) + >>> nx, ny = (3, 2) + >>> x = np.linspace(0, 1, nx) + >>> y = np.linspace(0, 1, ny) + >>> xv, yv = meshgrid(x, y) + >>> xv + array([[ 0. , 0.5, 1. ], + [ 0. , 0.5, 1. ]]) + >>> yv + array([[ 0., 0., 0.], + [ 1., 1., 1.]]) + >>> xv, yv = meshgrid(x, y, sparse=True) # make sparse output arrays + >>> xv + array([[ 0. , 0.5, 1. ]]) + >>> yv + array([[ 0.], + [ 1.]]) `meshgrid` is very useful to evaluate functions on a grid. >>> x = np.arange(-5, 5, 0.1) >>> y = np.arange(-5, 5, 0.1) - >>> xx, yy = np.meshgrid(x, y) + >>> xx, yy = meshgrid(x, y, sparse=True) >>> z = np.sin(xx**2+yy**2)/(xx**2+yy**2) - + + >>> import matplotlib.pyplot as plt + >>> h = plt.contourf(x,y,z) """ - x = asarray(x) - y = asarray(y) - numRows, numCols = len(y), len(x) # yes, reversed - x = x.reshape(1,numCols) - X = x.repeat(numRows, axis=0) + copy_ = kwargs.get('copy', True) + args = np.atleast_1d(*xi) + ndim = len(args) + + if not isinstance(args, list) or ndim<2: + raise TypeError('meshgrid() takes 2 or more arguments (%d given)' % int(ndim>0)) + + sparse = kwargs.get('sparse', False) + indexing = kwargs.get('indexing', 'xy') + + s0 = (1,) * ndim + output = [x.reshape(s0[:i] + (-1,) + s0[i + 1::]) for i, x in enumerate(args)] + + shape = [x.size for x in output] - y = y.reshape(numRows,1) - Y = y.repeat(numCols, axis=1) - return X, Y + if indexing == 'xy': + # switch first and second axis + output[0].shape = (1, -1) + (1,)*(ndim - 2) + output[1].shape = (-1, 1) + (1,)*(ndim - 2) + shape[0], shape[1] = shape[1], shape[0] + + if sparse: + if copy_: + return [x.copy() for x in output] + else: + return output + else: + # Return the full N-D matrix (not only the 1-D vector) + if copy_: + mult_fact = np.ones(shape, dtype=int) + return [x * mult_fact for x in output] + else: + return np.broadcast_arrays(*output) + + +def ndgrid(*args, **kwargs): + """ + Same as calling meshgrid with indexing='ij' (see meshgrid for + documentation). + """ + kwargs['indexing'] = 'ij' + return meshgrid(*args, **kwargs) def delete(arr, obj, axis=None): """ diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py index ba6b336ff..af85be326 100644 --- a/numpy/lib/tests/test_function_base.py +++ b/numpy/lib/tests/test_function_base.py @@ -1025,6 +1025,24 @@ class TestMeshgrid(TestCase): [5, 5, 5], [6, 6, 6], [7, 7, 7]]))) + def test_indexing(self): + [X, Y] = meshgrid([1, 2, 3], [4, 5, 6, 7], indexing='ij') + assert_(all(X == array([[1, 1, 1, 1], + [2, 2, 2, 2], + [3, 3, 3, 3]]))) + assert_(all(Y == array([[4, 5, 6, 7], + [4, 5, 6, 7], + [4, 5, 6, 7]]))) + + def test_sparse(self): + [X, Y] = meshgrid([1, 2, 3], [4, 5, 6, 7], sparse=True) + assert_(all(X == array([[1, 2, 3]]))) + assert_(all(Y == array([[4], + [5], + [6], + [7]]))) + + class TestPiecewise(TestCase): -- cgit v1.2.1 From 37d723cb629e3e8cb4d9ed7b85702fbd6350818d Mon Sep 17 00:00:00 2001 From: Ralf Gommers Date: Tue, 13 Dec 2011 23:24:47 +0100 Subject: MAINT: clean up docstring and some minor items in meshgrid. Remove ndgrid. --- numpy/lib/function_base.py | 78 ++++++++++++++++------------------- numpy/lib/tests/test_function_base.py | 2 - 2 files changed, 36 insertions(+), 44 deletions(-) diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py index 86ae3ab08..abb239ffc 100644 --- a/numpy/lib/function_base.py +++ b/numpy/lib/function_base.py @@ -6,7 +6,7 @@ __all__ = ['select', 'piecewise', 'trim_zeros', 'copy', 'iterable', 'histogram', 'histogramdd', 'bincount', 'digitize', 'cov', 'corrcoef', 'msort', 'median', 'sinc', 'hamming', 'hanning', 'bartlett', 'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc', 'add_docstring', - 'meshgrid', 'ndgrid', 'delete', 'insert', 'append', 'interp', 'add_newdoc_ufunc'] + 'meshgrid', 'delete', 'insert', 'append', 'interp', 'add_newdoc_ufunc'] import warnings import types @@ -3213,16 +3213,18 @@ def meshgrid(*xi, **kwargs): ---------- x1, x2,..., xn : array_like 1-D arrays representing the coordinates of a grid. - indexing : 'xy' or 'ij' (optional) - cartesian ('xy', default) or matrix ('ij') indexing of output - sparse : True or False (default) (optional) + indexing : {'xy', 'ij'}, optional + Cartesian ('xy', default) or matrix ('ij') indexing of output. + sparse : bool, optional If True a sparse grid is returned in order to conserve memory. - copy : True (default) or False (optional) - If False a view into the original arrays are returned in order to - conserve memory. Please note that sparse=False, copy=False will likely - return non-contiguous arrays. 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. + Default is False. + copy : bool, optional + If False, a view into the original arrays are returned in + order to conserve memory. Default is True. Please note that + ``sparse=False, copy=False`` will likely return non-contiguous arrays. + 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. Returns ------- @@ -3232,23 +3234,23 @@ def meshgrid(*xi, **kwargs): or ``(N2, N1, N3,...Nn)`` shaped arrays if indexing='xy' with the elements of `xi` repeated to fill the matrix along the first dimension for `x1`, the second for `x2` and so on. - + Notes ----- - This function supports both indexing conventions through the indexing keyword - argument. Giving the string 'ij' returns a meshgrid with matrix indexing, - while 'xy' returns a meshgrid with Cartesian indexing. The difference is - illustrated by the following code snippet: - - xv, yv = meshgrid(x, y, sparse=False, indexing='ij') - for i in range(nx): - for j in range(ny): - # treat xv[i,j], yv[i,j] - - xv, yv = meshgrid(x, y, sparse=False, indexing='xy') - for i in range(nx): - for j in range(ny): - # treat xv[j,i], yv[j,i] + This function supports both indexing conventions through the indexing keyword + argument. Giving the string 'ij' returns a meshgrid with matrix indexing, + while 'xy' returns a meshgrid with Cartesian indexing. The difference is + illustrated by the following code snippet:: + + xv, yv = meshgrid(x, y, sparse=False, indexing='ij') + for i in range(nx): + for j in range(ny): + # treat xv[i,j], yv[i,j] + + xv, yv = meshgrid(x, y, sparse=False, indexing='xy') + for i in range(nx): + for j in range(ny): + # treat xv[j,i], yv[j,i] See Also -------- @@ -3260,9 +3262,9 @@ def meshgrid(*xi, **kwargs): Examples -------- >>> nx, ny = (3, 2) - >>> x = np.linspace(0, 1, nx) - >>> y = np.linspace(0, 1, ny) - >>> xv, yv = meshgrid(x, y) + >>> x = np.linspace(0, 1, nx) + >>> y = np.linspace(0, 1, ny) + >>> xv, yv = meshgrid(x, y) >>> xv array([[ 0. , 0.5, 1. ], [ 0. , 0.5, 1. ]]) @@ -3281,17 +3283,17 @@ def meshgrid(*xi, **kwargs): >>> x = np.arange(-5, 5, 0.1) >>> y = np.arange(-5, 5, 0.1) >>> xx, yy = meshgrid(x, y, sparse=True) - >>> z = np.sin(xx**2+yy**2)/(xx**2+yy**2) - - >>> import matplotlib.pyplot as plt + >>> z = np.sin(xx**2 + yy**2) / (xx**2 + yy**2) >>> h = plt.contourf(x,y,z) + """ copy_ = kwargs.get('copy', True) args = np.atleast_1d(*xi) ndim = len(args) - - if not isinstance(args, list) or ndim<2: - raise TypeError('meshgrid() takes 2 or more arguments (%d given)' % int(ndim>0)) + + if ndim < 2: + msg = 'meshgrid() takes 2 or more arguments (%d given)' % int(ndim > 0) + raise TypeError(msg) sparse = kwargs.get('sparse', False) indexing = kwargs.get('indexing', 'xy') @@ -3321,14 +3323,6 @@ def meshgrid(*xi, **kwargs): return np.broadcast_arrays(*output) -def ndgrid(*args, **kwargs): - """ - Same as calling meshgrid with indexing='ij' (see meshgrid for - documentation). - """ - kwargs['indexing'] = 'ij' - return meshgrid(*args, **kwargs) - def delete(arr, obj, axis=None): """ Return a new array with sub-arrays along an axis deleted. diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py index af85be326..83532e70d 100644 --- a/numpy/lib/tests/test_function_base.py +++ b/numpy/lib/tests/test_function_base.py @@ -1043,8 +1043,6 @@ class TestMeshgrid(TestCase): [7]]))) - - class TestPiecewise(TestCase): def test_simple(self): # Condition is single bool list -- cgit v1.2.1 From 4e571172a210612bdeba1db0135b7d5fbc44ee0c Mon Sep 17 00:00:00 2001 From: Ralf Gommers Date: Wed, 28 Dec 2011 10:01:42 +0100 Subject: BUG: meshgrid: raise error on single input. --- numpy/lib/function_base.py | 9 +++++---- numpy/lib/tests/test_function_base.py | 4 ++++ 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py index abb239ffc..8e2fc375e 100644 --- a/numpy/lib/function_base.py +++ b/numpy/lib/function_base.py @@ -3288,13 +3288,14 @@ def meshgrid(*xi, **kwargs): """ copy_ = kwargs.get('copy', True) + + if len(xi) < 2: + msg = 'meshgrid() takes 2 or more arguments (%d given)' % int(len(xi) > 0) + raise ValueError(msg) + args = np.atleast_1d(*xi) ndim = len(args) - if ndim < 2: - msg = 'meshgrid() takes 2 or more arguments (%d given)' % int(ndim > 0) - raise TypeError(msg) - sparse = kwargs.get('sparse', False) indexing = kwargs.get('indexing', 'xy') diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py index 83532e70d..eda919df6 100644 --- a/numpy/lib/tests/test_function_base.py +++ b/numpy/lib/tests/test_function_base.py @@ -1025,6 +1025,10 @@ class TestMeshgrid(TestCase): [5, 5, 5], [6, 6, 6], [7, 7, 7]]))) + + def test_single_input(self): + assert_raises(ValueError, meshgrid, np.arange(5)) + def test_indexing(self): [X, Y] = meshgrid([1, 2, 3], [4, 5, 6, 7], indexing='ij') assert_(all(X == array([[1, 1, 1, 1], -- cgit v1.2.1 From 3fb541ef51ca60cf5a903ed83f69120fe3f5373a Mon Sep 17 00:00:00 2001 From: Ralf Gommers Date: Wed, 28 Dec 2011 10:30:08 +0100 Subject: TST: meshgrid: test expected shapes for Cartesian and matrix indexing. --- numpy/lib/function_base.py | 16 +++++++++++----- numpy/lib/tests/test_function_base.py | 13 ++++++++++++- 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py index 8e2fc375e..39d1d1c87 100644 --- a/numpy/lib/function_base.py +++ b/numpy/lib/function_base.py @@ -3215,6 +3215,7 @@ def meshgrid(*xi, **kwargs): 1-D arrays representing the coordinates of a grid. indexing : {'xy', 'ij'}, optional Cartesian ('xy', default) or matrix ('ij') indexing of output. + See Notes for more details. sparse : bool, optional If True a sparse grid is returned in order to conserve memory. Default is False. @@ -3238,9 +3239,13 @@ def meshgrid(*xi, **kwargs): Notes ----- This function supports both indexing conventions through the indexing keyword - argument. Giving the string 'ij' returns a meshgrid with matrix indexing, - while 'xy' returns a meshgrid with Cartesian indexing. The difference is - illustrated by the following code snippet:: + argument. Giving the string 'ij' returns a meshgrid with matrix indexing, + while 'xy' returns a meshgrid with Cartesian indexing. In the 2-D case + with inputs of length M and N, the outputs are of shape (N, M) for 'xy' + indexing and (M, N) for 'ij' indexing. In the 3-D case with inputs of + length M, N and P, outputs are of shape (N, M, P) for 'xy' indexing and (M, + N, P) for 'ij' indexing. The difference is illustrated by the following + code snippet:: xv, yv = meshgrid(x, y, sparse=False, indexing='ij') for i in range(nx): @@ -3287,8 +3292,6 @@ def meshgrid(*xi, **kwargs): >>> h = plt.contourf(x,y,z) """ - copy_ = kwargs.get('copy', True) - if len(xi) < 2: msg = 'meshgrid() takes 2 or more arguments (%d given)' % int(len(xi) > 0) raise ValueError(msg) @@ -3296,8 +3299,11 @@ def meshgrid(*xi, **kwargs): args = np.atleast_1d(*xi) ndim = len(args) + copy_ = kwargs.get('copy', True) sparse = kwargs.get('sparse', False) indexing = kwargs.get('indexing', 'xy') + if not indexing in ['xy', 'ij']: + raise ValueError("Valid values for `indexing` are 'xy' and 'ij'.") s0 = (1,) * ndim output = [x.reshape(s0[:i] + (-1,) + s0[i + 1::]) for i, x in enumerate(args)] diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py index eda919df6..57bfee61b 100644 --- a/numpy/lib/tests/test_function_base.py +++ b/numpy/lib/tests/test_function_base.py @@ -1030,7 +1030,9 @@ class TestMeshgrid(TestCase): assert_raises(ValueError, meshgrid, np.arange(5)) def test_indexing(self): - [X, Y] = meshgrid([1, 2, 3], [4, 5, 6, 7], indexing='ij') + x = [1, 2, 3] + y = [4, 5, 6, 7] + [X, Y] = meshgrid(x, y, indexing='ij') assert_(all(X == array([[1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3, 3]]))) @@ -1038,6 +1040,15 @@ class TestMeshgrid(TestCase): [4, 5, 6, 7], [4, 5, 6, 7]]))) + # Test expected shapes: + z = [8, 9] + assert_(meshgrid(x, y)[0].shape == (4, 3)) + assert_(meshgrid(x, y, indexing='ij')[0].shape == (3, 4)) + assert_(meshgrid(x, y, z)[0].shape == (4, 3, 2)) + assert_(meshgrid(x, y, z, indexing='ij')[0].shape == (3, 4, 2)) + + assert_raises(ValueError, meshgrid, x, y, indexing='notvalid') + def test_sparse(self): [X, Y] = meshgrid([1, 2, 3], [4, 5, 6, 7], sparse=True) assert_(all(X == array([[1, 2, 3]]))) -- cgit v1.2.1 From d48b756b232c99b6624d76db3188090052e0db60 Mon Sep 17 00:00:00 2001 From: Ralf Gommers Date: Sun, 5 Feb 2012 15:57:00 +0100 Subject: STY: meshgrid: some minor changes to address review comments. --- numpy/lib/function_base.py | 2 +- numpy/lib/tests/test_function_base.py | 5 +---- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py index 39d1d1c87..4f2f7b433 100644 --- a/numpy/lib/function_base.py +++ b/numpy/lib/function_base.py @@ -3223,7 +3223,7 @@ def meshgrid(*xi, **kwargs): If False, a view into the original arrays are returned in order to conserve memory. Default is True. Please note that ``sparse=False, copy=False`` will likely return non-contiguous arrays. - Furthermore, more than one element of a broadcasted array may refer to + Furthermore, more than one element of a broadcast array may refer to a single memory location. If you need to write to the arrays, make copies first. diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py index 57bfee61b..b38b39c94 100644 --- a/numpy/lib/tests/test_function_base.py +++ b/numpy/lib/tests/test_function_base.py @@ -1052,10 +1052,7 @@ class TestMeshgrid(TestCase): def test_sparse(self): [X, Y] = meshgrid([1, 2, 3], [4, 5, 6, 7], sparse=True) assert_(all(X == array([[1, 2, 3]]))) - assert_(all(Y == array([[4], - [5], - [6], - [7]]))) + assert_(all(Y == array([[4], [5], [6], [7]]))) class TestPiecewise(TestCase): -- cgit v1.2.1