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

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

0 голосов
/ 08 октября 2019

Мне нужно оценить каждый элемент матрицы, используя функцию с числовым интегралом (scipy.integrate.quad). Элементами матрицы являются пиксели серого изображения 5202x3465.
У меня есть доступ к графическому процессору, и я хотел бы оценить как можно больше элементов параллельно, потому что сейчас при линейном программировании все вычисления занимают более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), (arows, acolumns))
но это дает мне:
TypeError: myFunc () отсутствует 2 обязательных позиционных аргумента: 'j' и 'i'

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

Заранее благодарю за помощь!

Ответы [ 2 ]

0 голосов
/ 30 октября 2019

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

import quadpy


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


scheme = quadpy.line_segment.gauss_legendre(5)
val = scheme.integrate(f, [0, 1])
print(val)

[[0.33333333 0.25      ]
 [0.2        0.16666667]]

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

0 голосов
/ 11 октября 2019

Сначала ошибка в вашем вызове операции 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чтобы применить магию (см. здесь: использование starmap ): 

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

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

...