numpy - 在网格点上评估函数

Question

有什么好的方法可以生成一个包含在n维点网格上求值函数的numpy数组呢？

例如，假设我想要求解由以下函数定义：

def func(x, y):
    return <some function of x and y>

假设我想在一个二维点数组上进行评估，其中x值从0到4,分为十个步骤，y值从-1到1，分为20个步骤。在numpy中有什么好的方法可以做到这一点？

P.S. 这个问题以各种形式在StackOverflow上被问过很多次，但我找不到一个简洁明了的问题和答案。我发布这篇文章是为了提供一个简明的简单解决方案（如下）。

- DanielSank

---

你还希望我在我的回答中做出更多的改进吗？如果没有，您可以接受它作为答案吗？ - usethedeathstar

@usethedeathstar：我不喜欢你的方法的一个方面是，你必须明确地键入诸如x [:，None，：]之类的内容。这是不可扩展或脚本化的。另一方面，使用meshgrid，您只需放置尽可能多的参数，它就会执行正确的操作。如果您可以修改您的答案以处理该情况，我一定会接受它。 - DanielSank

修改了我的答案以使其按照你的要求工作。这种方法在numpy 1.8之前也可以使用，而meshgrid在numpy 1.8之前的版本中无法在超过2个维度上工作。 - usethedeathstar

Answer 1

更短、更快、更清晰的答案，避免使用meshgrid：

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是在一个1D数组上完成的（而不是2D数组），只有在执行"+"时才会增加维度。对于meshgrid，x^2将在一个2D数组上完成，其中每一行基本上都是相同的，导致时间大幅增加。

编辑： "x[:,None]"使x成为一个2D数组，但第二个维度为空。这个"None"相当于使用"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,:])

如果您希望，可以使用尽可能多的None或:来扩展到 n 个维度。每个:都会形成一个维度，而每个None则表示一个“空”维度。下一个示例展示了这些空维度的工作方式。正如您所看到的，在下一个示例中使用None会改变形状，但是只有在乘以实际在这些维度上具有内容的对象时，空维度才会被填充（听起来很复杂，但下一个示例会解释）。

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的第一个参数作为第一个完整维度，第二个参数作为第二个完整维度等等- 它与之前使用的“硬编码”None语法相同。

简短解释：执行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。

- usethedeathstar

---

你能解释一下这里的索引是怎么回事吗？func(xaxis[:,None], yaxis[None,:]) - rerx

@rerx在答案中添加了解释。 - usethedeathstar

@DanielSank 更新以在三维中显示更多内容。要转到4个或更多维度，则类似。 - usethedeathstar

np.newaxis 存在是为了使这更清晰明了。 - endolith

而且result = func(x[:,None], y[None,:]) 应该改为 result = func(xaxis[:, None], yaxis[None, :])，对吗？ - endolith

Answer 2

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)

- DanielSank

---

这不是一个很好的例子，因为 np.outer(ys, np.sin(xs)) 可能会更快。 - behzad.nouri

哦，你的意思是我应该使用一个与y不成线性关系的函数？ - DanielSank

如果 f(x,y) 不能分解为 g(x) 乘以 h(y)，则说明它不可分。 - behzad.nouri

Answer 3

我使用这个函数来获取X、Y、Z值以便绘图：

def npmap2d(fun, xs, ys, doPrint=False):
  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, np.linspace(0, 0.5, 21), np.linspace(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)

- grabantot

Answer 4

如果您的函数实际上接受一个元组，其中包含d个元素，即f((x1,x2,x3,...xd))（例如scipy.stats.multivariate_normal function），并且您想要在N个变量的N^d个组合/网格上评估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 计算的二维元组。

- Yibo Yang

---

如果函数f()期望参数是一个元组(x1,x2,x3,...,xd)，那么它是如何工作的呢？np.dstack()将返回所有的网格坐标，而不是一个单独的点。 - hertzsprung

为了使其适用于任意函数，我发现有必要先将函数向量化，例如 np.vectorize(f, signature='(d)->()')(np.dstack(np.meshgrid(x, y)))，假设 f() 返回一个标量。 - hertzsprung

Answer 5

我的看法：

    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)

- rainman

Answer 6

如果您的函数中有任何非向量操作，那么您很可能会遇到以下错误：

TypeError: 只有大小为1的数组可以转换为Python标量

为了解决这个问题，您需要使用numpy.vectorize。以下是一个完整的示例：

import numpy as np
import math

def f(x, y):
    return x - math.log(2+y)

X, Y = np.meshgrid(np.linspace(0, 4, 10), np.linspace(-1, 1, 20))

# f(X, Y)  # TypeError: only size-1 arrays can be converted to Python scalars
v = np.vectorize(f)
v(X, Y)

- feetwet