Распараллелить вычисление матричных элементов с помощью scipy.integrate.quad

Мне нужно оценить каждый элемент матрицы, используя функцию с числовым интегралом (scipy.integrate.quad). Элементами матрицы являются пиксели серого изображения размером 5202 x 3465.
У меня есть доступ к графическому процессору, и я хотел бы параллельно оценивать как можно больше элементов, потому что прямо сейчас, при линейном программировании, все вычисления занимают более 24 часа.
Вот пример кода:

for i in range(0, rows):
    for j in range(0, columns):
        img[i, j] = myFun(constant_args, i, j)

def myFunc(constant_args, i, j):
    new_pixel = quad(integrand, constant_args, i, j)
    ... other calculations ... 
    return new_pixel

Я пытался использовать многопроцессорность (как mp) следующим образом:

arows = list(range(0, rows))
acolumns = list(range(0, columns))
with mp.Pool() as pool:
    img = pool.map(myFunc, (constant_args, arows, acolumns))

или с img = pool.map(myFunc(constant_args), (стрелки, столбцы))
но это дает мне:
TypeError: myFunc() отсутствует 2 обязательных позиционных аргумента: 'j' и 'i'

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

Заранее спасибо за вашу помощь!


python scipy parallel-processing multiprocessing quad
person elena faillace    schedule 08.10.2019    source источник
comment
Обычно первым шагом является ускорение интеграции. Можете ли вы предоставить полный рабочий пример хотя бы интеграции? stackoverflow.com/q/45823212/4045774   -  person max9111    schedule 08.10.2019
comment
Я был бы рад помочь вам с этим (бесплатно, конечно), но мне нужен доступ к коду, который вы пытаетесь оптимизировать. Не стесняйтесь обращаться ко мне в Твиттере, если вы хотите узнать больше twitter.com/walkingrandomly   -  person WalkingRandomly    schedule 08.10.2019

Ответы (2)


arrow_upward
0
arrow_downward

Во-первых, ошибка в вашем вызове операции map-. Это должно быть:

arows = list(range(0, rows))
acolumns = list(range(0, columns))
with mp.Pool() as pool:
    img = pool.map(myFunc, constant_args, arows, acolumns)

Однако это может не привести к тому, что вы ищете, поскольку это просто проходит через 3 аргумента (которые должны быть списками). Он не проходит через их комбинации, особенно arows и acolumns. Например, если constant_args имеет 3 элемента, Pool.map остановится после 3 итераций, не проходя более длинные списки arows или acolumns.

Сначала вы должны сделать все комбинации индексов строк и столбцов

from itertools import product, repeat
comb = list(product(arows, acolumns))

Это производит что-то вроде (все возможные комбинации)

[(1, 1), (1, 2), (1, 3), (2, 1), (2, 2), (2, 3), (3, 1), (3, 2), (3, 3)]

Затем я бы заархивировал эти комбинации с вашим constant_args

constant_args = [10, 11]
arguments = list(zip(comb , repeat(constant_args)))

Он создает список кортежей, каждый из которых содержит два элемента. Первое — это ваша позиция в пикселях, а второе — ваше constant_args.

[((1, 1), [10, 11]),
 ((1, 2), [10, 11]),
 ((1, 3), [10, 11]),
 ((2, 1), [10, 11]),
 ((2, 2), [10, 11]),
 ((2, 3), [10, 11]),
 ((3, 1), [10, 11]),
 ((3, 2), [10, 11]),
 ((3, 3), [10, 11])]

Теперь нам нужно немного изменить ваш myFunc:

def myFunc(pix_id, constant_args):
    new_pixel = quad(integrand, constant_args, pix_id[0], pix_id[1])
    ... other calculations ... 
    return new_pixel

Наконец, мы используем Pool.starmap для волшебства (см. здесь: использование звездной карты):

with mp.Pool() as pool:
    img = pool.starmap(myFunc, arguments )

Что происходит, так это то, что starmap берет список кортежей и предоставляет их в качестве входных данных для функций. Однако starmap автоматически распаковывает список кортежей в отдельные аргументы для вашей функции. Первый аргумент — pix_id, состоящий из двух элементов, а второй аргумент — constant_args.

person RaJa    schedule 11.10.2019

arrow_upward
0
arrow_downward

Вы можете использовать quadpy (один из моих проектов). Он выполняет векторизованные вычисления, так что это будет работать очень быстро. Пример с выводом формы 2x2:

import quadpy


def f(x):
    return [[x ** 2, x**3], [x**4, x**5]]


val, err = quadpy.quad(f, 0, 1)
print(val)

[[0.33333333 0.25      ]
 [0.2        0.16666667]]

Выход f должен иметь форму (..., x.shape), а ... может быть любым кортежем, например, (5202, 3465).

person Nico Schlömer    schedule 30.10.2019