Очень медленный код для моделирования спектров оптического излучения (комбинационного рассеяния)

Question

Может ли кто-нибудь помочь мне оптимизировать мой код? В настоящее время для выполнения вычислений требуется до 15 минут, и я знаю, что это станет возможным менее чем за 10 секунд, поскольку это то, что ранее было достигнуто в Mathematica для чего-то очень похожего. Мне нужно моделировать много спектров, поэтому мне нужно намного более быстрое выполнение.

Я написал простой код для моделирования типа спектра оптического излучения (комбинационного рассеяния) для разных молекул. Эти спектры комбинационного рассеяния характеризуются длиной волны излучения и интенсивностью сигнала.

Я имею дело со сложной молекулой со многими спектральными линиями (3,5 миллиона), что может сделать код, описанный ниже, очень медленным.

Моя проблема возникает при преобразовании спектра палки (то есть пар длин волн и интенсивностей) в моделируемый спектр, который учитывает спектральное разрешение спектрометра. Применяя разрешение спектрометра при моделировании, каждая дискретная спектральная линия на определенной длине волны становится распределением интенсивности света в диапазоне длин волн. Форма линии каждой линии является нормальным (гауссовым) распределением. Как только форма линии рассчитана, сумма интенсивностей каждой спектральной линии на соответствующих длинах волн рассчитывается для определения общего спектра.

Этот набор уравнений для преобразования из линии в спектр палки приведен ниже:

vlow  = 2800.0 #cm-1
vhigh = 3050.0 #cm-1
dv    =    1 #cm-1
resolution = 2. #cm-1
wavelength = np.arange(vlow,vhigh,dv)

def gauss(wavelength,v0,resolution):
    I = np.exp(-((v0-wavelength)/(2*resolution))**2)
    I /= np.sum(I*dv)
    return I

I = np.zeros(len(wavelength))

def simulate(linelist):
    if simulate_spec == True:
        for line in linelist:
            wn_line,sig_line = line
            a = ((wn_line-vlow)  % dv)/dv
            b = int((wn_line-vlow) // dv)
            if b>=0 and b<=len(I)-2:    
                I[b  ] += sig_line * (1-a)
                I[b+1] += sig_line *    a

            halfwayvelength = 0.5*(wavelength[0]+wavelength[-1])
            Iconv = np.convolve(I,gauss(wavelength, halfwayvelength,
            resolution)*dv,'same')

    return(wavelength,Iconv)

Затем я просто применяю функцию имитации к списку данных, и спектры генерируются.

Пример данных приведен здесь:

[2988.051201, 2.1528594547391786],
[2965.874737, 9.832493405930505],
[2803.2158010000003, 0.00041319356649735474],
[2825.751201, 2.695035803386361e-06],
[2999.060832, 4.306173915906166],
[2826.208982, 0.00034324344563385156],
[2849.130431, 4.14575931249117e-06],
[2916.502513, 8.227955023939035e-07],
[3019.955186, 4.1128004082464935],
[2864.1083719999997, 4.7168887436524236e-07],
[2967.5610460000003, 2.3584443718262116],
[3011.59334, 8.484855475132457e-06],
[3013.1311739999996, 0.0017841866237857066],
[2810.959775, 1.6278821649562382e-05],
[2824.372448, 8.1122568270338e-06],
[2997.864286, 2.152180604795494],
[3044.539569, 0.000536291455511021],
[2966.429826, 3.932154056795204],
[3013.105109, 0.000979839907319629],
[2855.130571, 0.001800236014966169],
[2946.900868, 5.792331072221127],
[2992.882515, 0.016941583295928393],
[2830.159341, 4.1885065256659556e-05],
[2850.825762, 5.1388969852984256e-05],
[2916.542136, 5.781541962764192e-06],
[3018.591801, 3.0706663183839553],
[2843.208214, 6.239525546259207e-06],
[2945.257879, 3.3138244296836072],
[2991.205655, 0.0003465527450543083],
[2993.135932, 0.004361954391985438],
[2841.004186, 0.0008298581247320486],
[2849.245868, 0.00018996099882164256],
[2870.293698, 1.4363019628341073e-05],
[2875.746432, 2.100659304108956e-05],
[3039.611398, 1.0057174209518613],
[3041.328008, 5.840788673645993],
[2817.899557, 8.412353568584454e-05]