Optimization of numpy mesh creation for efficient interpolation

Question

I am reading magnetic field data from a text file. My goal is to correctly and efficiently load the mesh points (in 3 dimensions) and the associated fields (for simplicity I will assume below that I have a scalar field).

I managed to make it work, however I feel that some steps might not be necessary. In particular, reading the numpy doc it might be that "broadcasting" would be able to work its magic to my advantage.

import numpy as np
from scipy import interpolate
# Loaded from a text file, here the sampling over each dimension is identical but it is not required
x = np.array([-1.0, -0.5, 0.0, 0.5, 1.0])
y = np.array([-1.0, -0.5, 0.0, 0.5, 1.0])
z = np.array([-1.0, -0.5, 0.0, 0.5, 1.0])
# Create a mesh explicitely
mx, my, mz = np.meshgrid(x, y, z, indexing='ij')  # I have to switch from 'xy' to 'ij'
# These 3 lines seem odd
mx = mx.reshape(np.prod(mx.shape))
my = my.reshape(np.prod(my.shape))
mz = mz.reshape(np.prod(mz.shape))
# Loaded from a text file
field = np.random.rand(len(mx))
# Put it all together
data = np.array([mx, my, mz, field]).T
# Interpolate
interpolation_points = np.array([[0, 0, 0]])
interpolate.griddata(data[:, 0:3], data[:, 3], interpolation_points, method='linear')

Is it really necessary to construct the mesh like this? Is it possible to make it more efficient?

Feel free to let me know if this should be moved to codereview.SE — Cedric H.
– Cedric H., Commented Mar 17, 2019 at 10:44
Your reshape lines could use mx = mx.ravel() or reshape(-1). And data = np.stack([mx,my,mz,field], axis=1) is an alternative to np.array. — hpaulj
– hpaulj, Commented Mar 17, 2019 at 16:05

Divakar · Accepted Answer · 2019-03-17 11:08:39Z

1

Here's one with broadcasted-assignment to generate data directly from x,y,z and hence avoid the memory overhead of creating all the mesh-grids and hopefully lead to better performance -

m,n,r = len(x),len(y),len(z)
out = np.empty((m,n,r,4))
out[...,0]  = x[:,None,None]
out[...,1]  = y[:,None]
out[...,2]  = z
out[...,3]  = np.random.rand(m,n,r)
data_out = out.reshape(-1,out.shape[-1])

answered Mar 17, 2019 at 11:08

Divakar

222k19 gold badges273 silver badges374 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

Cedric H. Over a year ago

Thanks for the suggestion! I tested it: it works, but on my machine and for the size of data that I have I could not see a significant gain. This is still more elegant I believe than what I had. I think it's a pity that meshgrid does not provide that kind of output directly.

Divakar Over a year ago

@CedricH. Smells like others are the bottleneck, I see interpolate.griddata being one of them.

Collectives™ on Stack Overflow

Optimization of numpy mesh creation for efficient interpolation

1 Answer 1

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related