On 21/11/2016 12:44, Peter Otten wrote:
After a look into the source this is no longer a big surprise (numpy 1.8.2):
if axis is not None:
axisa, axisb, axisc=(axis,)*3
a = asarray(a).swapaxes(axisa, 0)
b = asarray(b).swapaxes(axisb, 0)
The situation may be different when you process vectors in bulk, i. e.
instead of [cross(a, bb) for bb in b] just say cross(a, b).
I thought that numpy was supposed to be fast? Also that the critical
bits were not implemented in Python?
Anyway I tried your code (put into a function as shown below) in the
same test program, and it was *still* 3 times as fast as numpy!
(mycross() was 3 times as fast as np.cross().)
Explain that...
---------------------------------------
def mycross(a,b,axisa=-1, axisb=-1, axisc=-1, axis=None):
if axis is not None:
axisa, axisb, axisc=(axis,)*3
a = np.asarray(a).swapaxes(axisa, 0)
b = np.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:
return cp.swapaxes(0, axisc)
else:
x = a[1]*b[2]
y = -a[0]*b[2]
z = a[0]*b[1] - a[1]*b[0]
elif a.shape[0] == 3:
if (b.shape[0] == 3):
x = a[1]*b[2] - a[2]*b[1]
y = a[2]*b[0] - a[0]*b[2]
z = a[0]*b[1] - a[1]*b[0]
else:
x = -a[2]*b[1]
y = a[2]*b[0]
z = a[0]*b[1] - a[1]*b[0]
cp = np.array([x, y, z])
if cp.ndim == 1:
return cp
else:
return cp.swapaxes(0, axisc)
---------------------------------------
Tested as:
x=np.array([1,2,3])
y=np.array([4,5,6])
start=time.clock()
for i in range(loops):
z=mycross(x,y)
print "Calc, %s loops: %.2g seconds" %(loops,time.clock()-start)
--
Bartc
--
https://mail.python.org/mailman/listinfo/python-list
