Numpy - оценивать функцию на сетке точек

Что такое хороший способ создания массива numpy, содержащего значения функции, вычисляемой на n-мерной сетке точек?

Например, предположим, что я хочу оценить функцию, определенную

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

Предположим, что я хочу оценить его на двухмерном массиве точек с x значениями, идущими от 0 до 4 за десять шагов, а значения y идут от -1 до 1 за двадцать шагов. Какой хороший способ сделать это в numpy?

P.S. Это неоднократно задавалось в разных формах StackOverflow, но я не мог найти кратко сформулированный вопрос и ответ. Я разместил это, чтобы предоставить краткое простое решение (ниже).

Ответы

Ответ 1

короче, быстрее и яснее ответ, избегая сетки:

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 будет выполняться на двумерном массиве, в котором практически все строки одинаковы, что приводит к значительному увеличению времени.

Редактировать: "x [:, None]", превращает x в двумерный массив, но с пустым вторым измерением. Это "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,:])

Таким образом, вы можете легко расширить до 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 -slicing для создания дополнительных пустых измерений, делая первый аргумент, который вы даете ndm как первое полное измерение, а второй как второе полное измерение etc- делает то же самое, что и "жестко закодированный" Нетипизированный синтаксис, используемый ранее.

Краткое объяснение: выполнение x2, y2, z2 = ndm(xaxis, yaxis, zaxis) аналогично выполнению

x2 = xaxis[:,None,None]
y2 = yaxis[None,:,None]
z2 = zaxis[None,None,:]

но метод ndm также должен работать для большего количества измерений, без необходимости жестко кодировать None -slices в несколько строк, как только что показано. Это также будет работать в numpy версиях до 1.8, в то время как numpy.meshgrid работает только для более чем 2-х измерений, если у вас numpy 1.8 или выше.

Ответ 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)

Ответ 3

В случае, когда ваша функция фактически берет кортеж элементов d, то есть f((x1,x2,x3,...xd)) (например, функция scipy.stats.multivariate_normal), и вы хотите оценить f на N ^ d комбинациях/сетке из N переменных, вы также можете сделать следующее (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.

Ответ 4

Я использую эту функцию для получения значений 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)

Тот же результат может быть достигнут с помощью карты:

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)

Ответ 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)

Ответ 6

Начинающий вопрос: что я тут не так делаю? если f (x, y) = x + y, то, очевидно, это работает. что отличается в max (x, y)?

def f(x,y):
    return max(x,y)

x=np.arange(5)
y=np.arange(5)
X,Y=np.meshgrid(x,y)
Z=f(X,Y)

Traceback (последний вызов последний):

Файл "", строка 7, в Z = f (X, Y)

Файл "", строка 2, в f возвращает max (x, y)

ValueError: Значение истинности массива с более чем одним элементом неоднозначно. Используйте a.any() или a.all()