CODEWITHSUNDEEP
pythonnumpy
mrkwjc
mrkwjc

Reputation: 1219

Numpy "multi meshgrid"

Is there any obvious way in numpy to replace something like:

for x in X:
    xi, xj = meshgrid(x, x, indexing='ij')

with a single (and possibly more efficient) operation like:

Xi, Xj = multi_meshgrid(X, X, indexing='ij')

The example of X is the following:

X = np.array([[0,1,2,3,4,5], [5,6,7,8,9,10], [11,12,13,14,15], ...])

The main problem is that i can have tens and hundreds of thousands of entries in X and possibly operation is often repeated.

The problem arises from assembling global stiffness matrix K in finite element method. For each entry in X of length n i have a matrix "n x n" which i have to inscribe into this global matrix. This matrix is in scipy.sparse coordinate format.

Regards, Marek

Upvotes: 0

Views: 236

Answers (1)

user2379410
user2379410

Reputation:

I think this answers the question, though I'm not sure if this is the best for constructing the sparse matrix in the end.. Anyway the following code creates a "view" into X, so it's very efficient both computationally and memory-wise. Try it :)

from numpy.lib.stride_tricks import as_strided

m = 3
n = 4
X = np.arange(m*n).reshape((m,n))

sz = X.itemsize
Xi = as_strided(X, shape=(m,n,n), strides=(n*sz, sz, 0))
Xj = as_strided(X, shape=(m,n,n), strides=(n*sz, 0, sz))

This does, however, not work when X is not a regular matrix. E.g. in your example the third row has 5 elements whereas the others have 6.

Upvotes: 2

Related Questions