CODEWITHSUNDEEP
pythonnumpyscipy
Cedric H.
Cedric H.

Reputation: 8298

Optimization of numpy mesh creation for efficient interpolation

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?

Upvotes: 1

Views: 112

Answers (1)

Divakar
Divakar

Reputation: 221564

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])

Upvotes: 1

Related Questions