Оптимизация кода Python с использованием векторизованного итерационного раздела

Question

Я новичок в python и numpy, есть этот код PDE второго порядка, который я хочу векторизовать для запуска за меньшее время, но поскольку он использует функцию для заполнения каждого значения в сетке, я застрял.

def func(x, y):
    return (-2 * np.pi ** 2) * np.sin(np.pi * x) * np.sin(np.pi * y)


def f1t(x, y):
    return (np.sin(np.pi * x) * np.sin(np.pi * y))


def five_point(grid, i, j, h, grid_x, grid_y):
    return ((grid[i + 1, j] + grid[i - 1, j] + grid[i, j + 1] + grid[i, j - 1]) / 4
            - ((h ** 2) / 4) * func(grid_x[i, 0], grid_y[0, j]))


def five_point_fin_int(X=np.ones(1), Y=np.ones(1), n_x=32, K_max=1000,
                       tol=0.0001, tol_type="grid"):
    import time;
    t0 = time.clock()

    h = 1 / n_x
    X_max = int(X / h)
    Y_max = int(Y / h)

    grid_x, grid_y = np.mgrid[0:X + h:h, 0:Y + h:h]
    grid_true = np.zeros((X_max + 1, Y_max + 1))

    for i in range(1, X_max):
        for j in range(1, Y_max):
            grid_true[i, j] = f1t(grid_x[i, 0], grid_y[0, j])

    grid = np.zeros((X_max + 1, Y_max + 1))
    grid_err = np.zeros((X_max + 1, Y_max + 1))

    count = 0
    tol_max = False
    while ((count < K_max) & (tol_max == False)):
        count += 1

        for i in range(1, X_max):
            for j in range(1, Y_max):
                grid[i, j] = five_point(grid, i, j, h, grid_x, grid_y)
                grid_err[i, j] = (grid[i, j] - grid_true[i, j])

                if (tol_type.lower() == "grid" ):
                    if (np.amax(abs(grid_err)) < tol):
                        tol_max = True
                else:
                    if (abs(np.linalg.norm(grid) - np.linalg.norm(grid_true)) < tol):
                        tol_max = True

    cpu_time = time.clock() - t0

В конце я печатаю время вычислений, так как сейчас оно вложено в циклы, затраченное на это время - около 9 секунд, и я хотел бы импровизировать на этом.

Answer 1

numpy позволяет заменять циклы на векторные вызовы.Вы определенно можете сделать следующее:

grid_true = np.zeros((X_max + 1, Y_max + 1))
grid_true[1:X_max,1:Y_max]=f1t(*np.meshgrid(grid_x[1:X_max,0], grid_y[0,1:Y_max]))

И вы также можете попробовать следующее:

grid = np.zeros((X_max + 1, Y_max + 1))
grid[1:-1, 1:-1] = five_point(grid, *np.meshgrid(np.arange(1,X_max), np.arange(1, Y_max)),
h, grid_x, grid_y)

Однако это не чистая "восходящая" интеграция, как та, которую вы делаете, так какВы, по сути, рассчитываете всю сетку вместе на каждом шаге (ваш вызов!).

Вероятно, процедура минимизации могла бы быть лучше.Нет большой разницы в производительности между numpy и pure python для коротких циклов или небольших векторов.