Python, оптимизирующий вычисления с помощью scipy.integrate.quad (занимает очень много времени) - PullRequest
0 голосов
/ 24 октября 2018

В настоящее время я пишу программу на python для расчета общей спектральной излучательной способности (инфракрасных волн) любого данного материала при различных температурах (200K - 500K) на основе данных измерений, полученных путем измерения направленно-полусферической излучательной способности материалана разных длинах волн с использованием ИК - спектроскопа.Расчет выполняется путем интегрирования измеренной интенсивности по всем длинам волн с использованием закона Планка в качестве функции взвешивания (все это на самом деле не имеет значения для самого моего вопроса; я просто хочу объяснить фон, чтобы код был проще для понимания),Это мой код:

from scipy import integrate
from scipy.interpolate import CubicSpline
import numpy as np
import math as m


def planck_blackbody(lambda_, T):  # wavelength, temperature
    h = float(6.6260755e-34)
    c = float(2.99792458e+8)
    k = float(1.380658e-23)
    try:
        a = 2.0 * h * (c ** 2)
        b = h * c / (lambda_ * k * T)
        intensity = a / ((lambda_ ** 5) * (m.exp(b) - 1.0))
        return float(intensity)
    except OverflowError:  # for lower temperatures
        pass


def spectral_emissivity(emifilename, t, lambda_1, lambda_2):
    results = []

    with open(emifilename, 'r') as emifile:
       emilines = emifile.readlines()
       try:
           w = [float(x.split('\t')[0].strip('\n')) * 1e-6 for x in emilines]
           e = [float(x.split('\t')[1].strip('\n')) for x in emilines]
       except ValueError:
           pass
    w = np.asarray(w)  # wavelength
    e = np.asarray(e)  # measured intensity

    def part_1(lambda_, T):
        E = interp1d(w, e, fill_value = 'extrapolate')(lambda_)
        return E * planck_blackbody(lambda_, T)

    def E_complete(T):
        E_complete_part_1 = integrate.quad(part_1, lambda_1, lambda_2, args=T, limit=50)
        E_complete_part_2 = integrate.quad(planck_blackbody, lambda_1, lambda_2, args=T, limit=50)
        return E_complete_part_1[0] / E_complete_part_2[0]

    for T in t:
        results.append([T, E_complete(T)])

    with open("{}.plk".format(emifilename[:-4]), 'w') as resultfile:
        for item in results:
            resultfile.write("{}\t{}\n".format(item[0], item[1]))

t = np.arange(200, 501, 1)
spectral_emissivity('C:\test.dat', t, 1.4e-6, 35e-6)

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

Когда я запускаю этот код, пока он дает правильные результаты, я все еще сталкиваюсь с 2 проблемами:

  1. Я получаю сообщение об ошибке от scipy.integrate.quad:

      IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
      If increasing the limit yields no improvement it is advised to analyze 
      the integrand in order to determine the difficulties.  If the position of a 
      local difficulty can be determined (singularity, discontinuity) one will 
      probably gain from splitting up the interval and calling the integrator 
      on the subranges.  Perhaps a special-purpose integrator should be used.
      warnings.warn(msg, IntegrationWarning)
    

Может кто-нибудь объяснить мне, что именно это означает?Я понимаю, что integrate.quad - это числовой метод итераций, и что моим функциям, похоже, требуется более 50 итераций, но есть ли способ обойти это?я попытался увеличить предел, но даже с 200 я все еще получаю это сообщение об ошибке ... это особенно странно, учитывая, что подынтегральные функции довольно простые функции ...

тесно связана с первой проблемой: этой программе требуется около 5 минут, чтобы закончить один файл, но мне нужно обрабатывать много файлов каждый час.cProfile показывает, что 98% этого времени тратится внутри функции интеграции.Программа MathCad, которая делает то же самое и выдает одинаковые результаты, занимает всего несколько секунд.Несмотря на то, что на прошлой неделе я потратил время на поиски решения, мне просто не удалось ускорить эту программу, и, похоже, ни у кого другого в stackoverflow и в других местах нет сравнимых временных проблем с integrate.quad.

Итак, наконец, мой вопрос: есть ли очевидный способ оптимизировать этот код для его ускорения (кроме компиляции в C + или чего-то подобного)?Я попытался уменьшить все поплавки до 6 цифр (я не могу опускаться ниже по точности), но это ничего не изменило.

Обновление: , глядя на это еще немного, я понялчто большую часть времени фактически не потребляла сама Интеграция, а операция CubicSpline, которую я использовал для интерполяции моих данных.Я опробовал разные методы, и CubicSpline, по какой-то причине, оказался единственным работающим (хотя мои данные монотонно увеличивались, я получал ошибки от каждого другого метода, который я пробовал, говоря, что некоторые значения были выше или ниже диапазона интерполяции).).То есть, пока я не узнал об экстраполяции with scipy.interpolate.interp1d и (fill_value = 'extrapolate'.Это помогло мне, позволив мне использовать гораздо менее трудоемкий метод interp1d и эффективно сократить время выполнения моей программы с 280 до 49 секунд (также добавлено понимание списка для w и e).Хотя это большое улучшение, я все еще удивляюсь, почему моей программе требуется около 1 минуты для вычисления некоторых интегралов ... и я все еще получаю вышеупомянутый IntegrationWarning.Так что любые советы очень ценятся!

(кстати, так как я довольно новичок в python, я рад любым советам или критике, которые я могу получить!)

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...