生成包含在n维点网格上评估的函数值的numpy数组有什么好方法?
例如,假设我想评估由其定义的函数
def func(x, y):
return <some function of x and y>
假设我想在二维点阵列上对它进行评估,其中x值从10到4分为10步,y值从-1到1分为20步。在numpy中这样做的好方法是什么?
附:这已经多次在StackOverflow上以各种形式提出,但我找不到简明扼要的问答。我发布了这个提供简洁的简单解决方案(如下)。
更短,更快,更清晰的答案,避免网格网格:
import numpy as np
def func(x, y):
return np.sin(y * x)
xaxis = np.linspace(0, 4, 10)
yaxis = np.linspace(-1, 1, 20)
result = func(xaxis[:,None], yaxis[None,:])
如果你得到像x ^ 2 + y这样的函数,那么内存会更快,因为x ^ 2是在一维数组(而不是2D数组)上完成的,并且尺寸的增加只发生在你执行“ +”。对于meshgrid,x ^ 2将在2D数组上完成,其中基本上每行都相同,导致大量时间增加。
编辑:“x [:,None]”,使x成为2D数组,但第二维空。这个“无”与使用“x [:,numpy.newaxis]”相同。使用Y完成同样的事情,但是制作一个空的第一维。
编辑:3维:
def func2(x, y, z):
return np.sin(y * x)+z
xaxis = np.linspace(0, 4, 10)
yaxis = np.linspace(-1, 1, 20)
zaxis = np.linspace(0, 1, 20)
result2 = func2(xaxis[:,None,None], yaxis[None,:,None],zaxis[None,None,:])
这样,您可以根据需要轻松扩展到n维,使用尽可能多的None
或:
。每个:
形成一个维度,每个None
形成一个“空”维度。下一个示例更多地展示了这些空白维度的工作原理。正如您所看到的,如果您使用None
,形状会发生变化,在下一个示例中显示它是一个3D对象,但只有当您与实际具有这些维度内容的对象相乘时,空尺寸才会被填满(听起来很复杂,但下一个例子显示了我的意思)
In [1]: import numpy
In [2]: a = numpy.linspace(-1,1,20)
In [3]: a.shape
Out[3]: (20,)
In [4]: a[None,:,None].shape
Out[4]: (1, 20, 1)
In [5]: b = a[None,:,None] # this is a 3D array, but with the first and third dimension being "empty"
In [6]: c = a[:,None,None] # same, but last two dimensions are "empty" here
In [7]: d=b*c
In [8]: d.shape # only the last dimension is "empty" here
Out[8]: (20, 20, 1)
编辑:无需自己键入None
def ndm(*args):
return [x[(None,)*i+(slice(None),)+(None,)*(len(args)-i-1)] for i, x in enumerate(args)]
x2,y2,z2 = ndm(xaxis,yaxis,zaxis)
result3 = func2(x2,y2,z2)
通过这种方式,您可以使None
切片创建额外的空白维度,方法是将第一个参数赋予ndm作为第一个完整维度,第二个参数作为第二个完整维度等 - 它与“硬编码”无关 - 之前使用的类型语法。
简短的解释:做x2, y2, z2 = ndm(xaxis, yaxis, zaxis)
和做的一样
x2 = xaxis[:,None,None]
y2 = yaxis[None,:,None]
z2 = zaxis[None,None,:]
但是ndm方法也应该适用于更多维度,而不需要像刚刚显示的那样在多行中对None
切片进行硬编码。这也适用于1.8之前的numpy版本,而numpy.meshgrid仅适用于高于2维的情况,如果你有1.8或更高的numpy。
import numpy as np
def func(x, y):
return np.sin(y * x)
xaxis = np.linspace(0, 4, 10)
yaxis = np.linspace(-1, 1, 20)
x, y = np.meshgrid(xaxis, yaxis)
result = func(x, y)
如果您的函数实际上采用了d
元素的元组,即f((x1,x2,x3,...xd))
(例如scipy.stats.multivariate_normal function),并且您想要在N ^ d组合/ N个变量的网格上评估f
,您还可以执行以下操作(2D情况):
x=np.arange(-1,1,0.2) # each variable is instantiated N=10 times
y=np.arange(-1,1,0.2)
Z=f(np.dstack(np.meshgrid(x,y))) # result is an NxN (10x10) matrix, whose entries are f((xi,yj))
在这里,np.dstack(np.meshgrid(x,y))
创建了一个10x10“矩阵”(技术上是一个10x10x2 numpy数组),其条目是由f
评估的二维元组。
我使用此函数来获取X,Y,Z值以供绘图:
def npmap2d(fun, x_spec, y_spec, doPrint=False):
xs = np.linspace(*x_spec)
ys = np.linspace(*y_spec)
Z = np.empty(len(xs) * len(ys))
i = 0
for y in ys:
for x in xs:
Z[i] = fun(x, y)
if doPrint: print([i, x, y, Z[i]])
i += 1
X, Y = np.meshgrid(xs, ys)
Z.shape = X.shape
return X, Y, Z
用法:
def f(x, y):
# ...some function that can't handle numpy arrays
X, Y, Z = npmap2d(f, (0, 0.5, 21), (0.6, 0.4, 41))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(X, Y, Z)
使用map可以获得相同的结果:
xs = np.linspace(0, 4, 10)
ys = np.linspace(-1, 1, 20)
X, Y = np.meshgrid(xs, ys)
Z = np.fromiter(map(f, X.ravel(), Y.ravel()), X.dtype).reshape(X.shape)
我的两分钱:
import numpy as np
x = np.linspace(0, 4, 10)
y = np.linspace(-1, 1, 20)
[X, Y] = np.meshgrid(x, y, indexing = 'ij', sparse = 'true')
def func(x, y):
return x*y/(x**2 + y**2 + 4)
# I have defined a function of x and y.
func(X, Y)