2-D function and meshgrid
Hello, I have a function that receives a array of shape (2,) and returns a number (a function from R^2 -> R). It basically looks like this: def weirdDistance2(x): return dot(dot(weirdMatrix, x), x) (weirdMatrix is a "global" (2,2) array) I want to see its level sets in the box [0, 1] x [0, 1], hence I have to create a meshgrid and then compute it at each point of the mesh: x = linspace(0, 1, 200) y = x.copy() X, Y = meshgrid(x, y) My problem is how to actually compute the function at each point of the mesh. I have come out with two solutions. One very short and clear, but slow, and another longer and more convoluted (it has a loop, I hate loops in numpy code), but faster. Does anyone know a "no explicit-loops" and fast solution? Solution1: def myDistance(a, b): return weirdDistance(np.array((a, b))) vecDistance = np.vectorize(myDistance) return vecDistance(X, Y) Solution 2: nPoints = X.size result = np.zeros(nPoints) points = np.array( [X.ravel(), Y.ravel()] ).T for i in xrange(nPoints): result[i] = weirdDistance(points[i]) result = result.reshape(X.shape) Of course, the first one is slow because the myDistance function creates an array at each call. The second one, even with a loop, avoids the array creations. Best, Paulo
On Fri, Jan 9, 2009 at 11:45 AM, Paulo J. S. Silva <pjssilva@ime.usp.br>wrote:
Hello,
I have a function that receives a array of shape (2,) and returns a number (a function from R^2 -> R). It basically looks like this:
def weirdDistance2(x): return dot(dot(weirdMatrix, x), x)
(weirdMatrix is a "global" (2,2) array)
I want to see its level sets in the box [0, 1] x [0, 1], hence I have to create a meshgrid and then compute it at each point of the mesh:
x = linspace(0, 1, 200) y = x.copy() X, Y = meshgrid(x, y)
My problem is how to actually compute the function at each point of the mesh. I have come out with two solutions. One very short and clear, but slow, and another longer and more convoluted (it has a loop, I hate loops in numpy code), but faster. Does anyone know a "no explicit-loops" and fast solution?
Solution1:
def myDistance(a, b): return weirdDistance(np.array((a, b))) vecDistance = np.vectorize(myDistance) return vecDistance(X, Y)
Solution 2:
nPoints = X.size result = np.zeros(nPoints) points = np.array( [X.ravel(), Y.ravel()] ).T for i in xrange(nPoints): result[i] = weirdDistance(points[i]) result = result.reshape(X.shape)
Of course, the first one is slow because the myDistance function creates an array at each call. The second one, even with a loop, avoids the array creations.
Try IDLE 1.2.4
import numpy as np pts = np.random.rand(5,2) mat = np.random.rand(2,2) res = (np.dot(pts,mat)*pts).sum(axis=1) res
array([ 0.63018561, 0.30829864, 0.23173343, 1.79972127, 0.69498856])
for row in pts : np.dot(row,np.dot(mat,row))
0.63018560596590589 0.30829864146737423 0.23173343333294744 1.7997212735553192 0.69498855520540959 Chuck
Chuck, Thanks, your version is much faster. I would prefer a solution that doesn't force me to re-implement weirdDistance (as my two solutions were). But the function is so simple that it is easier just to re-write it for speed as you did. By the way, I came out with one more solution that looks more Pythonic and does not need to re-write weirdDistance (and hence can be used in more complicated cases). It is also a tad faster than the fastest solution from my first post: Solution 3 points = np.vstack( [x.ravel(), y.ravel()] ).T results = np.array([weirDistance(p) for p in points]) return results.reshape(x.shape) (This is basically solution 2 using list comprehensions to make to code clearer) best, Paulo
Try
IDLE 1.2.4
import numpy as np pts = np.random.rand(5,2) mat = np.random.rand(2,2) res = (np.dot(pts,mat)*pts).sum(axis=1) res
array([ 0.63018561, 0.30829864, 0.23173343, 1.79972127, 0.69498856])
for row in pts : np.dot(row,np.dot(mat,row))
0.63018560596590589 0.30829864146737423 0.23173343333294744 1.7997212735553192 0.69498855520540959
Chuck
participants (2)
-
Charles R Harris -
Paulo J. S. Silva