Re: [Numpy-discussion] which one is best?
Hi Arnar, Your two commands below aren't doing the same thing - one is doing a[i]*b[i] and the other is doing a[i]*b[j] for all i and j. As the second is harder, it takes longer. Cheers, David On Fri, 2008-09-19 at 09:08 -0500, [EMAIL PROTECTED] wrote: I think [x*y for x in a for y in b] feels pythonic, however it has a surprisingly lousy performance. In [30]: %timeit [ x * y for x,y in zip(a,b) ] 10 loops, best of 3: 3.96 ?s per loop In [31]: %timeit [ i*j for i in a for j in b ] 10 loops, best of 3: 6.53 ?s per loop In [32]: a = range(100) In [33]: b = range(100) In [34]: %timeit [ x * y for x,y in zip(a,b) ] 1 loops, best of 3: 51.9 ?s per loop In [35]: %timeit [ i*j for i in a for j in b ] 100 loops, best of 3: 2.78 ms per loop Arnar -- ** David M. Kaplan Charge de Recherche 1 Institut de Recherche pour le Developpement Centre de Recherche Halieutique Mediterraneenne et Tropicale av. Jean Monnet B.P. 171 34203 Sete cedex France Phone: +33 (0)4 99 57 32 27 Fax: +33 (0)4 99 57 32 95 http://www.ur097.ird.fr/team/dkaplan/index.html ** ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
[Numpy-discussion] New patch for new mgrid / ogrid functionality
Hi all, Attached is a newer version of my patch that adds new mgrid / ogrid functionality for working with arrays in addition to slices. In fact, I have attached two versions of the patch: index_tricks.patch, that is just the last version of the patch I sent, and index_tricks.new.patch, that has been modified so that it is backward compatible. In the last version, mgrid calls where all arguments are slices will return an array, otherwise it returns a list as ogrid does. This is the only reasonable way to have the new functionality and maintain backwards compatibility. My 2 cents - I personally think the version that always returns a list will ultimately be more transparent and cause fewer problems than the newer version. In either case, the plan should be to eventually have it always return a list as that is the only fully consistent option, the question is just when that switch should be made and by who. If it is done at the next major release, someone else will have to remember to ax the additional code and correct the documentation Other changes that would be nice: add a __call__ method, create an instance called ndgrid for matlab compatibility, and have meshgrid be reimplimented using an nd_grid instance. Cheers, David -- ** David M. Kaplan Charge de Recherche 1 Institut de Recherche pour le Developpement Centre de Recherche Halieutique Mediterraneenne et Tropicale av. Jean Monnet B.P. 171 34203 Sete cedex France Phone: +33 (0)4 99 57 32 27 Fax: +33 (0)4 99 57 32 95 http://www.ur097.ird.fr/team/dkaplan/index.html ** Index: numpy/lib/tests/test_index_tricks.py === --- numpy/lib/tests/test_index_tricks.py (revision 5834) +++ numpy/lib/tests/test_index_tricks.py (working copy) @@ -24,15 +24,21 @@ def test_nd(self): c = mgrid[-1:1:10j,-2:2:10j] d = mgrid[-1:1:0.1,-2:2:0.2] -assert(c.shape == (2,10,10)) -assert(d.shape == (2,20,20)) +assert(array(c).shape == (2,10,10)) +assert(array(d).shape == (2,20,20)) assert_array_equal(c[0][0,:],-ones(10,'d')) assert_array_equal(c[1][:,0],-2*ones(10,'d')) assert_array_almost_equal(c[0][-1,:],ones(10,'d'),11) assert_array_almost_equal(c[1][:,-1],2*ones(10,'d'),11) -assert_array_almost_equal(d[0,1,:]-d[0,0,:], 0.1*ones(20,'d'),11) -assert_array_almost_equal(d[1,:,1]-d[1,:,0], 0.2*ones(20,'d'),11) +assert_array_almost_equal(d[0][1,:]-d[0][0,:], 0.1*ones(20,'d'),11) +assert_array_almost_equal(d[1][:,1]-d[1][:,0], 0.2*ones(20,'d'),11) +def test_listargs(self): +e = mgrid[ :2, ['a', 'b', 'c'], [1,5,50,500] ] +assert( array(e).shape == (3,2,3,4) ) +assert_array_equal( e[0][:,1,1].ravel(), r_[:2] ) +assert_array_equal( e[1][1,:,1].ravel(), array(['a','b','c']) ) +assert_array_equal( e[2][1,1,:].ravel(), array([1,5,50,500]) ) class TestConcatenator(TestCase): def test_1d(self): Index: numpy/lib/index_tricks.py === --- numpy/lib/index_tricks.py (revision 5834) +++ numpy/lib/index_tricks.py (working copy) @@ -11,7 +11,7 @@ from numpy.core.numerictypes import find_common_type import math -import function_base +import function_base, shape_base import numpy.core.defmatrix as matrix makemat = matrix.matrix @@ -118,14 +118,28 @@ number of points to create between the start and stop values, where the stop value **is inclusive**. +One can also use lists or arrays as indexing arguments, in which case +these will be meshed out themselves instead of generating matrices from +the slice arguments. See examples below. + If instantiated with an argument of sparse=True, the mesh-grid is open (or not fleshed out) so that only one-dimension of each returned argument is greater than 1 +***IMPORTANT NOTE*** Indexing an nd_grid instance with +sparse=False will currently return an array N+1 axis array if all +arguments are slices (i.e., something like -4:5:20j or :20:0.5) +and there are N arguments. However, if any of the arguments is +not a slice (e.g., is an array or list), then the return is a list +of arrays. This is to maintain backwards compatibility. However, +this functionality will disappear during the next major release +(after today's date: 2008-09-19) so that returns will always be +lists in the future. + Examples mgrid = np.lib.index_tricks.nd_grid() - mgrid[0:5,0:5] + mgrid[0:5,0:5] # NOTE currently returns array, but will become a list array([[[0, 0, 0, 0, 0], [1, 1, 1, 1, 1], [2, 2, 2, 2, 2], @@ -139,6 +153,27 @@ [0, 1, 2, 3, 4]]]) mgrid[-1:1:5j] array([-1. , -0.5, 0. , 0.5, 1. ]) + mgrid[:2,[1,5,50],['a','b']] # Example
[Numpy-discussion] patch for new mgrid / ogrid functionality
Hi, A while back, I sent some changes to index_tricks.py that would allow mgrid and ogrid to mesh things other than slices. For example: mgrid[['a','b'],[float,int],:3] [array([[['a', 'a', 'a'], ['a', 'a', 'a']], [['b', 'b', 'b'], ['b', 'b', 'b']]], dtype='|S1'), array([[[type 'float', type 'float', type 'float'], [type 'int', type 'int', type 'int']], [[type 'float', type 'float', type 'float'], [type 'int', type 'int', type 'int']]], dtype=object), array([[[0, 1, 2], [0, 1, 2]], [[0, 1, 2], [0, 1, 2]]])] At the time, there wasn't much follow-up, but I am hoping that there is still interest in this functionality, as I have gone ahead and finished the patch including documentation changes and updates to test_index_tricks.py. Attached is a patch set to the latest subversion of the numpy trunk. I don't think I am allowed to commit the changes myself - correct me if I am wrong. This functionality seems like a nice addition to me as it allows one to mesh things that are not uniformly spaced and potentially not even numbers. The changes don't affect functionality that existed previously except for one minor change - instead of returning a numpy array of arrays, mgrid/ogrid now return a list of arrays. However, this is unlikely to be a problem as the majority of users generally unpack the results of mgrid/ogrid so that each matrix can be used individually. Comments welcome. Cheers, David -- ** David M. Kaplan Charge de Recherche 1 Institut de Recherche pour le Developpement Centre de Recherche Halieutique Mediterraneenne et Tropicale av. Jean Monnet B.P. 171 34203 Sete cedex France Phone: +33 (0)4 99 57 32 27 Fax: +33 (0)4 99 57 32 95 http://www.ur097.ird.fr/team/dkaplan/index.html ** Index: numpy/lib/tests/test_index_tricks.py === --- numpy/lib/tests/test_index_tricks.py (revision 5654) +++ numpy/lib/tests/test_index_tricks.py (working copy) @@ -24,15 +24,21 @@ def test_nd(self): c = mgrid[-1:1:10j,-2:2:10j] d = mgrid[-1:1:0.1,-2:2:0.2] -assert(c.shape == (2,10,10)) -assert(d.shape == (2,20,20)) +assert(array(c).shape == (2,10,10)) +assert(array(d).shape == (2,20,20)) assert_array_equal(c[0][0,:],-ones(10,'d')) assert_array_equal(c[1][:,0],-2*ones(10,'d')) assert_array_almost_equal(c[0][-1,:],ones(10,'d'),11) assert_array_almost_equal(c[1][:,-1],2*ones(10,'d'),11) -assert_array_almost_equal(d[0,1,:]-d[0,0,:], 0.1*ones(20,'d'),11) -assert_array_almost_equal(d[1,:,1]-d[1,:,0], 0.2*ones(20,'d'),11) +assert_array_almost_equal(d[0][1,:]-d[0][0,:], 0.1*ones(20,'d'),11) +assert_array_almost_equal(d[1][:,1]-d[1][:,0], 0.2*ones(20,'d'),11) +def test_listargs(self): +e = mgrid[ :2, ['a', 'b', 'c'], [1,5,50,500] ] +assert( array(e).shape == (3,2,3,4) ) +assert_array_equal( e[0][:,1,1].ravel(), r_[:2] ) +assert_array_equal( e[1][1,:,1].ravel(), array(['a','b','c']) ) +assert_array_equal( e[2][1,1,:].ravel(), array([1,5,50,500]) ) class TestConcatenator(TestCase): def test_1d(self): Index: numpy/lib/index_tricks.py === --- numpy/lib/index_tricks.py (revision 5654) +++ numpy/lib/index_tricks.py (working copy) @@ -11,7 +11,7 @@ from numpy.core.numerictypes import find_common_type import math -import function_base +import function_base, shape_base import numpy.core.defmatrix as matrix makemat = matrix.matrix @@ -118,6 +118,10 @@ number of points to create between the start and stop values, where the stop value **is inclusive**. +One can also use lists or arrays as indexing arguments, in which case +these will be meshed out themselves instead of generating matrices from +the slice arguments. See examples below. + If instantiated with an argument of sparse=True, the mesh-grid is open (or not fleshed out) so that only one-dimension of each returned argument is greater than 1 @@ -126,19 +130,38 @@ mgrid = np.lib.index_tricks.nd_grid() mgrid[0:5,0:5] -array([[[0, 0, 0, 0, 0], -[1, 1, 1, 1, 1], -[2, 2, 2, 2, 2], -[3, 3, 3, 3, 3], -[4, 4, 4, 4, 4]], -BLANKLINE - [[0, 1, 2, 3, 4], -[0, 1, 2, 3, 4], -[0, 1, 2, 3, 4], -[0, 1, 2, 3, 4], -[0, 1, 2, 3, 4]]]) +[array([[0, 0, 0, 0, 0], + [1, 1, 1, 1, 1], + [2, 2, 2, 2, 2], + [3, 3, 3, 3, 3], + [4, 4, 4, 4, 4]]), array([[0, 1, 2, 3, 4], + [0, 1, 2, 3, 4], + [0, 1, 2, 3, 4], + [0, 1, 2, 3, 4], + [0, 1, 2, 3, 4]])] mgrid[-1:1:5j] array([-1. , -0.5, 0. , 0.5
[Numpy-discussion] [SciPy-user] unique, sort, sortrows
Hi, Well, as usual there are compromises in everything and the mgrid/ogrid functionality is the way it currently is for some good reasons. The first reason is that python appears to be fairly sloppy about how it passes indexing arguments to the __getitem__ method. It passes a tuple containing the arguments in all cases except when it has one argument, in which case it just passes that argument. This means that it is hard to tell a tuple argument from several non-tuple arguments. For example, the following two produce exactly the same call to __getitem__ : mgrid[1,2,3] mgrid[(1,2,3)] (__getitem__ receives a single tuple (1,2,3)), but different from: mgrid[[1,2,3]] (__getitem__ receives a single list = [1,2,3]). This seems like a bug to me, but is probably considered a feature by somebody. In any case, this is workable, but a bit annoying in that tuple arguments just aren't going to work well. The second problem is that the current implementation is fairly efficient because it forces all arguments to the same type so as to avoid some unnecessary copies (I think). Once you allow non-slice arguments, this is hard to maintain. That being said, attached is a replacement for index_tricks.py that should implement a reasonable set of features, while only very slightly altering performance. I have only touched nd_grid. I haven't fixed the documentation string yet, nor have I added tests to test_index_tricks.py, but will do that if the changes will be accepted into numpy. With the new version, old stuff should work as usual, except that mgrid now returns a list of arrays instead of an array of arrays (note that this change will cause current test_index_tricks.py to fail). With the new changes, you can now do: mgrid[-2:5:10j,[4.5,6,7.1],12,['abc','def']] The following will work as expected: mgrid[:5,(1,2,3)] But this will not: mgrid[(1,2,3)] # same as mgrid[1,2,3], but different from mgrid[[1,2,3]] Given these limitations, this seems like a fairly useful addition. If this looks usable, I will clean up and add tests if desired. If not, I recommend adding a ndgrid function to numpy that does the equivalent of matlab [X,Y,Z,...] = ndgrid(x,y,z,...) and then making the current meshgrid just call that changing the order of the first two arguments. Cheers, David -- ** David M. Kaplan Charge de Recherche 1 Institut de Recherche pour le Developpement Centre de Recherche Halieutique Mediterraneenne et Tropicale av. Jean Monnet B.P. 171 34203 Sete cedex France Phone: +33 (0)4 99 57 32 27 Fax: +33 (0)4 99 57 32 95 http://www.ur097.ird.fr/team/dkaplan/index.html ** __all__ = ['unravel_index', 'mgrid', 'ogrid', 'r_', 'c_', 's_', 'index_exp', 'ix_', 'ndenumerate','ndindex'] import sys import numpy.core.numeric as _nx from numpy.core.numeric import asarray, ScalarType, array, dtype from numpy.core.numerictypes import find_common_type import math import function_base import numpy.core.defmatrix as matrix makemat = matrix.matrix # contributed by Stefan van der Walt def unravel_index(x,dims): Convert a flat index into an index tuple for an array of given shape. e.g. for a 2x2 array, unravel_index(2,(2,2)) returns (1,0). Example usage: p = x.argmax() idx = unravel_index(p,x.shape) x[idx] == x.max() Note: x.flat[p] == x.max() Thus, it may be easier to use flattened indexing than to re-map the index to a tuple. if x _nx.prod(dims)-1 or x 0: raise ValueError(Invalid index, must be 0 = x = number of elements.) idx = _nx.empty_like(dims) # Take dimensions # [a,b,c,d] # Reverse and drop first element # [d,c,b] # Prepend [1] # [1,d,c,b] # Calculate cumulative product # [1,d,dc,dcb] # Reverse # [dcb,dc,d,1] dim_prod = _nx.cumprod([1] + list(dims)[:0:-1])[::-1] # Indices become [x/dcb % a, x/dc % b, x/d % c, x/1 % d] return tuple(x/dim_prod % dims) def ix_(*args): Construct an open mesh from multiple sequences. This function takes n 1-d sequences and returns n outputs with n dimensions each such that the shape is 1 in all but one dimension and the dimension with the non-unit shape value cycles through all n dimensions. Using ix_() one can quickly construct index arrays that will index the cross product. a[ix_([1,3,7],[2,5,8])] returns the array a[1,2] a[1,5] a[1,8] a[3,2] a[3,5] a[3,8] a[7,2] a[7,5] a[7,8] out = [] nd = len(args) baseshape = [1]*nd for k in range(nd): new = _nx.asarray(args[k]) if (new.ndim != 1): raise ValueError, Cross index must be 1 dimensional if issubclass(new.dtype.type, _nx.bool_): new = new.nonzero()[0] baseshape[k] = len(new) new = new.reshape(tuple(baseshape)) out.append(new
