Interpolate in multidimensional array in python with matrix as output-CodePudding

In my problem I have a function that varies over one or more parameters and outputs a matrix. This function is now sampled. So suppose the function F is dependent on x and y, then for some x-y grid, I have for every grid point a matrix.

In MATLAB I solve this using a 4D matrix, where the first two dimensions are the matrix itself, and the remaining are the coordinates of the x-y grid. See the example below:

n = 2;
% x-y grid 
xg = [-1, 1, 7];
sz_xg = size(xg,2);
yg = [3, 4, 9];
sz_yg = size(yg,2);

% some function:
fcn = @(x,y) [3*x sin(y), cos(x); exp(-x), 5 y]; % some n x n matrix function

% Define grid
A = zeros(n,n,sz_xg,sz_yg);
for i = 1:sz_xg
    for j = 1:sz_yg
        A(1:n, 1:n, i, j) = fcn(xg(i), yg(j));
    end
end

xpoint = 5;
ypoint = 6;

Ainterp = interpn(1:n, 1:n, xg, yg, A, 1:n, 1:n, xpoint, ypoint)

I tried to fix this as well in Python, but no luck. When I try to implement the exact same example I get a 2x1 vector...

my attempt:

import numpy as np
from scipy.interpolate import interpn

n = 2
xg = np.array([-1, 1, 7])
yg = np.array([3, 4, 9])

# some function
def fcn(x,y):
    return np.array([[3*x np.sin(y), np.cos(x)],[np.exp(-x), 5 y]])

A = np.zeros((n,n,xg.size,yg.size))
for i in range(xg.size):
    for j in range(yg.size):
        A[:,:,i,j]=fcn(xg[i],yg[j])

xpoint = 5
ypoint = 6

At = interpn((np.arange(n), np.arange(n), xg, yg),
             A,(np.arange(0,n),np.arange(0,n),xpoint,ypoint))

Could you help me out here?

CodePudding user response：

the implementation of interpn is slightly different between Matlab and scipy.

Matlab returns the permutations on the input vectors , so when you specify [1,2] for the first and second dimension, it returns values at

[(1,1) (1,2);
 (2,1) (2,2)]

however this is not the case for scipy interpn, scipy uses parallel arrays to determine the location of interpolation, which means the first output is constructed from the first input in both arrays, and second output is constructed from the second input in both arrays and so on.

in order to get the behavior you want, you have to pass the entire matrix listed above, using a more efficient method described below:

grid_points = np.meshgrid(np.arange(0,n),np.arange(0,n),sparse=True)

At = interpn((np.arange(n), np.arange(n), xg, yg),
             A,(grid_points[0],grid_points[1],xpoint,ypoint))

basically, meshgrid(sparse=True) will produce the matrix i wrote above in an efficient and memory compressed manner for you.

in my opinion, python method is more predictable than matlab, as it's easier to determine the dimensions of the output array, unlike in matlab, where the dimension of the output is based on how many of the inputs were not a single value, which makes it hard to predict the order of the dimensions in the output.