Как ускорить работу метода четырехъядерного интегрирования, применяемого к большому массиву numpy, который должен быть построен с помощью matplotlib?

У меня есть функция многих постоянных параметров, но только одна независимая переменная m. Эта независимая переменная изменяется от некоторого минимального значения, low, до некоторого максимального значения, up. Эти границы известны априори. Теперь я хотел бы определить отдельную функцию, которая вычисляет площадь под первой функцией в 10 ячейках по м. Фактически, я должен интегрировать первую функцию в любой из этих 10 интервалов, а затем сохранить результат в другом массиве, который позже нужно будет построить относительно независимой переменной m. Я поместил левую часть моих значений x в один массив, а правую часть моих значений x - в другой массив, так что эти два будут входными данными функции, которая будет построена, а именно нижняя и верхняя границы интеграла от первой функции. Вот мой код для этого на Python:

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import quad, dblquad
import mpmath as mp


low, up  = 5.630e5, 1.167e12
alpha, xo = 1.05   , 2.15e10     
beta = 274


def g(x, low, up, beta):
    return mp.gamma(-2/3) * (mp.gammainc(-2/3, beta*(x/low)**3) - mp.gammainc(-2/3, beta*(x/up)**3))
gv = np.vectorize(g)

def Integrand1(x, low, up, xo, alpha, beta):
    return pow(x/xo, alpha) * g(x, low, up, beta)

def Integrand2(x, low, up, xo, alpha, beta):
    return g(x, low, up, beta)


def PDF(x, low, up, xo, alpha, beta):
    Integral1 = quad(Integrand1, low, xo, args=(low, up, xo, alpha, beta))
    Integral2 = quad(Integrand2, xo, up, args=(low, up, xo, alpha, beta))
    A=(Integral1[0]+Integral2[0])**(-1)
    y = np.piecewise(x, 
                    [x < xo], [lambda x: A * pow(x/xo, alpha) * gv(x, low, up, beta), 
                                   lambda x: A * gv(x, low, up, beta)
                                  ]
                    )
    return y



x_array = np.array(np.logspace(8.2, 11.8, 10))
deltam = 0.4
lolims, uplims = np.array(np.logspace(8.2-deltam/2, 11.8-deltam/2, len(x_array))), np.array(np.logspace(8.2+deltam/2, 11.8+deltam/2, len(x_array)))


def myfunc(M1, M2, low, up, xo, alpha, beta):
    return np.array([quad(PDF, M1, M2, args=(low, up, xo, alpha, beta)) for M1 in lolims for M2 in uplims])


y_array = np.array([myfunc(M1, M2, low, up, xo, alpha, beta) for M1 in lolims for M2 in uplims]) 
plt.plot(x_array, y_array, color='green', linestyle='-')
plt.gca().autoscale(False)
plt.vlines([xo], linestyles='dashed', color='k', label='')

plt.legend(loc='upper left', ncol=1, fontsize=12)
plt.xlim([1e8, M_max])
plt.ylim([2e-6, 2e-1])
plt.ylabel('x', fontsize=12)
plt.xlabel('y', fontsize=12)
plt.xscale("log", nonposx='clip')
plt.yscale("log", nonposy='clip')
plt.show()

Проблема в том, что код работает вечно. Я не уверен, что он правильный, но не чистый и лаконичный, или что-то серьезное не так, что не дает ожидаемого сюжета.


arrays function matplotlib tuples integral
person Ash    schedule 21.11.2017    source источник
comment
Вы пробовали использовать профилировщик, такой как runnakerun, чтобы проверить, какие вызовы функций на самом деле являются дорогостоящими?   -  person alexblae    schedule 22.11.2017
comment
Я не знаю профилировщика. Но после печати некоторых результатов в разных местах кода я считаю, что последняя функция (myfunc), которая фактически вычисляет площадь под функцией PDF, занимает много времени, особенно когда две границы (low и up) расположены по разные стороны от x = xo (low ‹xo и вверх ›xo). Но это происходит только один раз, хотя на это уходит очень много времени. Частично задержка также связана с тем, что мне нужно как минимум 10 интервалов (чтобы увидеть, как выглядит общий тренд) или ~ 100 интервалов, чтобы определить точный паттерн (если в нем есть какая-то особенность). Спасибо еще раз,   -  person Ash    schedule 22.11.2017
comment
извините, я имел в виду M1 и M2, а не low и up (которые исправлены).   -  person Ash    schedule 22.11.2017

Ответы (1)


arrow_upward
0
arrow_downward

Так как выходом функции quad является кортеж, первый элемент которого является результатом, а второй элемент - ошибкой результата, необходимо сохранить только первые элементы и работать с ними как с массивом numpy. Однако узкое место проблемы заключается в способе определения myfunc. Предполагается, что эта функция берет свою переменную M1 из нижних пределов и свою переменную M2 из массивов верхних пределов. Однако неясно, как это написано, в каком порядке происходят эти присвоения, что означает, что они смешиваются между двумя массивами до такой степени, что вызывает путаницу для кода, чтобы знать, какие именно элементы он должен принимать для (M1 < / em>, M2) кортеж.

Частичный ответ:

Используйте функцию zip, чтобы поместить ваш диапазон M1 и M2 в упорядоченном виде, чтобы myfunc знал в каждом цикле, что следует принимать в качестве своего входы. Для этого нужно изменить определение myfunc и следующего y_array на следующие операторы:

y_tuples = [quad(PDF, M1, M2, args=(low, up, xo, alpha, beta)) for (M1, M2) in zip(lolims, uplims)]
y_array = np.array([i[0] for i in y_tuples])

Таким образом, проблема с кодом решена, но проблема того, что он медленный, все еще сохраняется и ожидает ответа.

person Ash    schedule 22.11.2017