numpy meshgrid操作问题


问题内容
Y, X = np.mgrid[-3:-3:10j, -3:3:10j]

我注意到,在像上面的网格网格上应用某些操作时,我收到一个错误,因为这些操作可能与numpy不兼容。有时可能有一个numpy函数替代sin,cos,但不是所有函数(例如scipy.integrate中的quad函数)。

我如何解决这个问题?我需要在整个网格上应用操作。


问题答案:

您的问题(带有后续注释)可以至少采用两种不同的方式:

  1. 您具有多个参数的函数,并且希望能够以与numpy本机支持的广播调用在语法上相似的方式调用该函数。性能不是问题,只是函数的调用语法。

  2. 您具有要在一系列numpy数组上求值的多个参数的函数,但是该函数的实现方式无法利用numpy数组的连续内存布局。性能是问题;您将很乐意遍历numpy数组,并以无聊的普通的for循环样式调用函数,但这样做太慢了。

对于第1项,numpy提供了一个便捷函数,vectorize该函数接受一个常规可调用对象并返回一个可调用对象,该可调用对象可以使用numpy数组作为参数来调用,并且将遵循numpy的广播规则。

考虑这个人为的例子:

def my_func(x, y):
    return x + 2*y

现在,假设我需要在二维网格中的任何地方评估此函数。这是普通的无聊方式:

Y, X  =  np.mgrid[0:10:1, 0:10:1]
Z = np.zeros_like(Y)

for i in range(Y.shape[0]):
    for j in range(Y.shape[1]):
        Z[i,j] = my_func(X[i,j], Y[i,j])

如果我们有几个类似的函数my_func,最好将这个过程概括为一个将给定函数“映射”到二维数组上的函数。

import itertools
def array_map(some_func, *arg_arrays):
    output = np.zeros_like(arg_arrays[0])
    coordinates = itertools.imap(range, output.shape)
    for coord in itertools.product(coordinates):
        args = [arg_array[coord] for arg_array in arg_arrays]
        output[coord] = some_func(*args)
    return output

现在我们可以看到它的array_map(my_func, X, Y)行为就像嵌套的for循环:

In [451]: array_map(my_func, X, Y)
Out[451]: 
array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9],
       [ 2,  3,  4,  5,  6,  7,  8,  9, 10, 11],
       [ 4,  5,  6,  7,  8,  9, 10, 11, 12, 13],
       [ 6,  7,  8,  9, 10, 11, 12, 13, 14, 15],
       [ 8,  9, 10, 11, 12, 13, 14, 15, 16, 17],
       [10, 11, 12, 13, 14, 15, 16, 17, 18, 19],
       [12, 13, 14, 15, 16, 17, 18, 19, 20, 21],
       [14, 15, 16, 17, 18, 19, 20, 21, 22, 23],
       [16, 17, 18, 19, 20, 21, 22, 23, 24, 25],
       [18, 19, 20, 21, 22, 23, 24, 25, 26, 27]])

现在,如果我们可以调用array_map(my_func)并省去多余的数组参数,那不是很好吗?取而代之的是取回一个新函数,该函数正等待执行所需的for循环。

我们可以用functools.partial-做到这一点,所以我们可以编写一个方便的矢量化器,如下所示:

import functools
def vectorizer(regular_function):
    awesome_function = functools.partial(array_map, regular_function)
    return awesome_function

并进行测试:

In [453]: my_awesome_func = vectorizer(my_func)

In [454]: my_awesome_func(X, Y)
Out[454]: 
array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9],
       [ 2,  3,  4,  5,  6,  7,  8,  9, 10, 11],
       [ 4,  5,  6,  7,  8,  9, 10, 11, 12, 13],
       [ 6,  7,  8,  9, 10, 11, 12, 13, 14, 15],
       [ 8,  9, 10, 11, 12, 13, 14, 15, 16, 17],
       [10, 11, 12, 13, 14, 15, 16, 17, 18, 19],
       [12, 13, 14, 15, 16, 17, 18, 19, 20, 21],
       [14, 15, 16, 17, 18, 19, 20, 21, 22, 23],
       [16, 17, 18, 19, 20, 21, 22, 23, 24, 25],
       [18, 19, 20, 21, 22, 23, 24, 25, 26, 27]])

现在的my_awesome_func行为就好像您可以直接在ndarrays上调用它一样!

在将这个玩具版本称为vectorizer…时,我忽略了许多额外的性能细节,边界检查等内容,但幸运的是,在numpyvectorize中已经有这样做了!

In [455]: my_vectorize_func = np.vectorize(my_func)

In [456]: my_vectorize_func(X, Y)
Out[456]: 
array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9],
       [ 2,  3,  4,  5,  6,  7,  8,  9, 10, 11],
       [ 4,  5,  6,  7,  8,  9, 10, 11, 12, 13],
       [ 6,  7,  8,  9, 10, 11, 12, 13, 14, 15],
       [ 8,  9, 10, 11, 12, 13, 14, 15, 16, 17],
       [10, 11, 12, 13, 14, 15, 16, 17, 18, 19],
       [12, 13, 14, 15, 16, 17, 18, 19, 20, 21],
       [14, 15, 16, 17, 18, 19, 20, 21, 22, 23],
       [16, 17, 18, 19, 20, 21, 22, 23, 24, 25],
       [18, 19, 20, 21, 22, 23, 24, 25, 26, 27]])

再次强调,正如我之前对OP的评论以及该文档中所强调的那样vectorize-这 不是
速度优化。实际上,在某些情况下,额外的函数调用开销比直接编写for循环要慢。但是,对于速度不是问题的情况,此方法可以使您的自定义函数遵循与numpy相同的调用约定-
这样可以提高库接口的统一性,并使代码更一致和更具可读性。

关于第2项,还有很多其他文章。如果您的问题是您需要优化功能以利用连续的内存块并绕过重复的动态类型检查(numpy数组添加到Python列表的主要功能)
),那么以下一些链接可能会有所帮助:

  1. < http://pandas.pydata.org/pandas-docs/stable/enhancingperf.html >
  2. < http://csl.name/C-functions-from-Python/ >
  3. < https://jakevdp.github.io/blog/2014/05/09/why-python-is-slow >
  4. <nbviewer.ipython.org/url/jakevdp.github.io/downloads/notebooks/NumbaCython.ipynb>