I am trying to take transpose of a matrix using the interpolation function in scipy. The goal is NOT to take the transpose, but this is a benchmark for something. The idea is, I define my new coordinates (primed, xGp etc below), which are x'=z, y'=x, z'=y. This should be the same as a 2,0,1 transpose. However, these results do not match and I see a large discrepancy. Also, it seems to match a 1,2,0 transpose for reasons not clear to me. I encounter the same problem with map_coordinates function of scipy. If anybody can point out what I am missing that would be great, thanks! from scipy.interpolate import interpn #N = 64 Nx=8 Ny=8 Nz=8 Lx=2*np.pi dx=Lx/Nx dy=Lx/Ny dz=Lx/Nz x = np.linspace(0,Lx-dx,Nx) y = np.linspace(0,Lx-dy,Ny) z = np.linspace(0,Lx-dz,Nz) xG,yG,zG = np.meshgrid(x,y,z, indexing='ij') #generate some random data data = np.cos(xG*4) + 1e-3*np.sin(yG+xG) + zG #These are the new primed coordinates xGp = zG/dz yGp = xG/dx zGp = yG/dy #new coordinates are (z,x,y) #this should effectively be the same as taking a tranpose #in order 2,0,1 xGp = xGp.flatten() yGp = yGp.flatten() zGp = zGp.flatten() #this is the list of coordinates #at which I will interpolate coords = np.zeros((Nx*Ny*Nz,3)) for i in range(Nx*Ny*Nz): coords[i,0] = xGp coords[i,1] = yGp coords[i,2] = zGp coords = coords.astype(int) out = interpn((x/dx,y/dy,z/dz), data, coords) out = out.reshape((Nx,Ny,Nz)) print(np.max(np.abs(out-np.transpose(data,(2,0,1))))) #the above is NOT 0 print(out.shape) Continue reading...
