numpy 向量与矩阵的Map

y3bcpkx1  于 4个月前  发布在  其他
关注(0)|答案(2)|浏览(72)

我需要根据坐标计算一些函数的值,比如说

f(x_i,y_i)

字符串
结果应该看起来像一个矩阵:

f(x1,y1),f(x2,y1),...,f(xn,y1)
f(x1,y2),f(x2,y2),...,f(xn,y2)
..............................
f(x1,yn),f(x2,yn),...,f(xn,yn)


现在,由于f(x,y)的计算需要每个x_i和y_i,我认为最好的方法是:

X=(x1,x2,x3,...,xn)
Y=(y1,y2,y3,...,yn)


我可以计算f(x,y)而不需要循环。但是,我需要进一步的计算,这需要每个f(Xi,yi)的位置。具体来说,我需要计算一个2*2子矩阵的每个行列式:

d_{ij}=f(x_i,y_j)+f(x_{i+1},y_{j+1})-f(x_i,y_{j+1})-f(x_{i+1},y_j)


并将其存储在一个新的矩阵中。据我所知,解决上述问题的方法是:
1.得到X和Y
1.计算f(X,Y)
1.将f(X,Y)重新Map回矩阵结构
1.计算每个行列式
我知道如何做第一步:如何生成一个三维坐标矩阵numpy.meshgrid下面是一个最小的代码:

import numpy as np
import scipy as scp
def GG(r):
    return np.sqrt((y-y0)**2 + r**2 * np.sin(r))
def beta1(r):
    return ((y-y0) - r * np.cos(r)) / GG(r)
def f(x,y):
    return scp.integrate.quad_vec(lambda r: 1/(np.sinh(self.beta1(r))))
x = np.linspace(2,5,10)
y = np.linspace(-5,5,10)
y0 = 1
x0 = 1
a = np.stack(np.meshgrid(x, y, indexing="ij"), axis=-1)
xx = a[:, 0, 0]
yy = a[0, :, 1]
fxy = f(xx,yy)

wnrlj8wa

wnrlj8wa1#

正如你已经提出的:要获取X和Y值,你可以简单地使用np.meshgrid对于第二个问题,使用numpy你可以计算f(X,Y),根据你的逻辑,你可以简单地将@vectorize Package 器放在函数上,然后用你拥有的X和Y调用它。对于第三部分,这应该已经在步骤2中完成了。最后,使用卷积将是一种方法来实现你正在寻找的第4部分。
如果你正在寻找一种有限方法,我建议你研究一下numpy vectorize及其操作!

q3qa4bjr

q3qa4bjr2#

前面的答案展示了将两个坐标数组组合成一个数组的各种方法。但我不确定这些方法是否有助于您理解问题。
使用meshgrid可以得到一个数组的元组,这些数组相互“广播”:

In [156]: x, y = np.meshgrid(np.arange(5), np.arange(4),indexing='ij', sparse=True)

In [157]: x,y
Out[157]: 
(array([[0],
        [1],
        [2],
        [3],
        [4]]),
 array([[0, 1, 2, 3]]))

In [158]: x.shape, y.shape
Out[158]: ((5, 1), (1, 4))

字符串
如果没有sparse,我们得到的二维数组具有(5,4)形状。功能上它们是相同的:

In [159]: X,Y = np.meshgrid(np.arange(5), np.arange(4),indexing='ij'); X.shape
Out[159]: (5, 4)


这两种形式都可以相加(或相乘,相减等),以产生(5,4)结果:

In [160]: x+y
Out[160]: 
array([[0, 1, 2, 3],
       [1, 2, 3, 4],
       [2, 3, 4, 5],
       [3, 4, 5, 6],
       [4, 5, 6, 7]])

In [161]: X+Y
Out[161]: 
array([[0, 1, 2, 3],
       [1, 2, 3, 4],
       [2, 3, 4, 5],
       [3, 4, 5, 6],
       [4, 5, 6, 7]])


执行类似np.stack((X,Y),axis=-1)的操作会产生一个(5,4,2)数组,但对f(x,y)形式的函数没有帮助。
mgridogridindex是单独或组合两个数组的其他方法。甚至np.ix_
如果你的f不能处理广播(和x,y形状一样),但是可以处理2d数组,你可以给予X,Y

Z = f(X.ravel(), Y,ravel())


如果ZX.ravel()的形状相同,则可以将其重新整形为与X相同的形状。

相关问题