CODEWITHSUNDEEP
pythonperformancenumpyvectorization
user211337
user211337

Reputation: 504

Cross product all points in 3D space

I have a vector, a, which I wish to cross with every point in a defined 3D space.

import numpy as np

# Grid 
x = np.arange(-4,4,0.1)
y = np.arange(-4,4,0.1)
z = np.arange(-4,4,0.1)

a = [1,0,0]

result = [[] for i in range(3)]

for j in range(len(x)):  # loop on x coords
    for k in range(len(y)):  # loop on y coords
        for l in range(len(z)):  # loop on z coords
            r = [x[j] , y[k], z[l]]
            result[0].append(np.cross(a, r)[0])
            result[1].append(np.cross(a, r)[1])
            result[2].append(np.cross(a, r)[2])

This produces an array which has taken the cross product of a with every point in space. However, the process takes far too long, due to the nested loops. Is there anyway to exploit vectors (meshgrid perhaps?) to make this process faster?

Upvotes: 2

Views: 456

Answers (2)

Divakar
Divakar

Reputation: 221574

Here's one vectorized approach -

np.cross(a, np.array(np.meshgrid(x,y,z)).transpose(2,1,3,0)).reshape(-1,3).T

Sample run -

In [403]: x = np.random.rand(4)
     ...: y = np.random.rand(5)
     ...: z = np.random.rand(6)
     ...: 

In [404]: result = original_app(x,y,z,a)

In [405]: out = np.cross(a, np.array(np.meshgrid(x,y,z)).\
                    transpose(2,1,3,0)).reshape(-1,3).T

In [406]: np.allclose(result[0], out[0])
Out[406]: True

In [407]: np.allclose(result[1], out[1])
Out[407]: True

In [408]: np.allclose(result[2], out[2])
Out[408]: True

Runtime test -

# Original setup used in the question
In [393]: # Grid 
     ...: x = np.arange(-4,4,0.1)
     ...: y = np.arange(-4,4,0.1)
     ...: z = np.arange(-4,4,0.1)
     ...: 

# Original approach
In [397]: %timeit original_app(x,y,z,a)
1 loops, best of 3: 21.5 s per loop

# @Denziloe's soln
In [395]: %timeit [np.cross(a, r) for r in product(x, y, z)]
1 loops, best of 3: 7.34 s per loop

# Proposed in this post
In [396]: %timeit np.cross(a, np.array(np.meshgrid(x,y,z)).\
                        transpose(2,1,3,0)).reshape(-1,3).T
100 loops, best of 3: 16 ms per loop

More than 1000x speedup over the original one and more than 450x over the loopy approach from other post.

Upvotes: 1

Denziloe
Denziloe

Reputation: 8131

This takes a couple of seconds to run on my machine:

from itertools import product
result = [np.cross(a, r) for r in product(x, y, z)]

I don't know if that's fast enough for you, but there are a lot of calculations involved. It's certainly cleaner, and there is at least some reduction of redundancy (e.g. calculating np.cross(a, r) three times). It also gives the result in a slightly different format, but this is the natural way to store the result and is hopefully fine for your purposes.

Upvotes: 1

Related Questions