On Di, 2014-07-15 at 10:22 +0100, Neil Hodgson wrote: > Hi, > > We came across this bug while using np.cross on 3D arrays of 2D > vectors.
Hi, which numpy version are you using? Until recently, the cross product simply did *not* work in a broadcasting manner (3d arrays of 2d vectors), it did something, but usually not the right thing. This is fixed in recent versions (not sure if 1.8 or only now with 1.9) - Sebastian > The first example shows the problem and we looked at the source for > np.cross and believe we found the bug - an unnecessary swapaxes when > returning the output (comment inserted in the code). > > Thanks > Neil > > # Example > > shape = (3,5,7,2) > > > # These are effectively 3D arrays (3*5*7) of 2D vectors > data1 = np.random.randn(*shape) > data2 = np.random.randn(*shape) > > > # The cross product of data1 and data2 should produce a (3*5*7) array > of scalars > cross_product_longhand = > data1[:,:,:,0]*data2[:,:,:,1]-data1[:,:,:,1]*data2[:,:,:,0] > print 'longhand output shape:',cross_product_longhand.shape # and it > does > > > cross_product_numpy = np.cross(data1,data2) > print 'numpy output shape:',cross_product_numpy.shape # It seems to > have transposed the last 2 dimensions > > > if (cross_product_longhand == np.transpose(cross_product_numpy, > (0,2,1))).all(): > print 'Unexpected transposition in numpy.cross (numpy version %s)'% > np.__version__ > > > # np.cross L1464 > if axis is not None: > axisa, axisb, axisc=(axis,)*3 > a = asarray(a).swapaxes(axisa, 0) > b = asarray(b).swapaxes(axisb, 0) > msg = "incompatible dimensions for cross product\n"\ > "(dimension must be 2 or 3)" > if (a.shape[0] not in [2, 3]) or (b.shape[0] not in [2, 3]): > raise ValueError(msg) > if a.shape[0] == 2: > if (b.shape[0] == 2): > cp = a[0]*b[1] - a[1]*b[0] > if cp.ndim == 0: > return cp > else: > ## WE SHOULD NOT SWAPAXES HERE! > ## For 2D vectors the first axis has been > > ## collapsed during the cross product > return cp.swapaxes(0, axisc) > > > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
