Почему мой код намного медленнее в форме функции? - PullRequest
Новая
       29

Почему мой код намного медленнее в форме функции?

0 голосов
/ 04 мая 2018

Я недавно написал код, который обрабатывает спектры, используя данные из текстовых файлов и выполняя вычисления на них. Я начал с кода, который просто выполняет все построчно без каких-либо функций, и, несмотря на то, что он длинный, он завершается за 2,11 секунды (согласно %%timeit). Ниже приведен исходный код, помеченный как таковой.

Однако я хотел вместо этого поместить свой код в функции, чтобы в будущем его было легче читать и использовать с различными моделями. Несмотря на то, что я использую все те же шаги, что и раньше (но на этот раз в своих функциях), это намного медленнее. Этот код также ниже. Теперь мне нужно подождать минут 15-20, чтобы получить те же результаты. Почему это намного медленнее, и есть ли способ, которым я могу сделать это значительно быстрее, но все же использовать функции?

Оригинальный код: 

import re
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate
filename = 'bpass_spectra.txt'
extinctionfile = 'ExtinctionLawPoints.txt' # from R_V = 4.0
pointslist = []
datalist = []
speclist = []

# Constants
Msun = 1.98892e30 # solar mass [kg]
h = 4.1357e-15 # Planck's constant [eV s]
c = float(3e8) # speed of light [m/s]

# Read spectra file
f = open(filename, 'r')
rawspectra = f.readlines()
met = re.findall('Z\s=\s(\d*\.\d+)', rawspectra[0])
del rawspectra[0]
for i in range(len(rawspectra)):
    newlist = rawspectra[i].split(' ')
    datalist.append(newlist)

# Read extinction curve data file
rawpoints = open(extinctionfile, 'r').readlines()
for i in range(len(rawpoints)):
    newlst = re.split('(?!\S)\s(?=\S)|(?!\S)\s+(?=\S)', rawpoints[i])
    pointslist.append(newlst)
pointslist = pointslist[3:]
lambdalist = [float(item[0]) for item in pointslist]
k_abslist = [float(item[4]) for item in pointslist]
xvallist = [(c*h)/(lamb*1e-6) for lamb in lambdalist]
k_interp = scipy.interpolate.interp1d(xvallist, k_abslist)

# Create new lists
Elist = [float(item[0]) for item in datalist]
speclambdalist = [h*c*1e9/E for E in Elist]
z1list = [float(item[1]) for item in datalist]
speclist.extend(z1list)
met = met[0]
klist = [None]*len(speclist)
Loutlist = [None]*len(speclist)
Tlist = [None]*len(speclist)

# Define parameters
b = 2.0
R = 1.0
z = 1.0
Mgas = 1.0 # mass of gas, input
Mhalo = 2e41 # mass of dark matter halo, known
if float(met) > 0.0052:
    DGRlist = [50.0*np.exp(-2.21)*float(met)]*len(speclist)
elif float(met) <= 0.0052:
    DGRlist = [((50.0*float(met))**3.15)*np.exp(-0.96)]*len(speclist)
for i in range(len(speclist)):
    if float(Elist[i]) <= 4.1357e-3: # frequencies <= 10^12 Hz
        klist[i] = 0.1*(float(Elist[i])/(1000.0*h))**b # extinction law [cm^2/g]
    elif float(Elist[i]) > 4.1357e-3: # frequencies > 10^12 Hz
        klist[i] = k_interp(Elist[i]) # interpolated function's value at Elist[i]
Mdustlist = [Mgas*DGR for DGR in DGRlist] # dust mass
Rhalo = 0.784*(0.27**2.0)*(0.7**(-2.0/3.0))*float(10.0/(1.0+z))*((Mhalo/(1e8*Msun))**(1.0/3.0))
Rdust = 0.018*Rhalo # [kpc]

for i in range(len(speclist)):
    Tlist[i] = 3*Mdustlist[i]*klist[i]/(4*np.pi*Rdust)

Linlist = [float(spectra)*R for spectra in speclist]

# Outgoing luminosity as function of wavelength
for i in range(len(Linlist)):
    Loutlist[i] = Linlist[i]*np.exp(-Tlist[i])

# Test the calculation
print "LIN ELEMENTS 0 AND 1000:", Linlist[0], Linlist[1000]
print "LOUT ELEMENTS 0 AND 1000:", Loutlist[0], Loutlist[1000]

Новый «функциональный» код (намного медленнее): 

import re
import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate

# Required files and lists
filename = 'bpass_spectra.txt' # number of columns = 4
extinctionfile = 'ExtinctionLawPoints.txt' # R_V = 4.0
datalist = []
if filename == 'bpass_spectra.txt':
    filetype = 4
else:
    filetype = 1
if extinctionfile == 'ExtinctionLawPoints.txt':
    R_V = 4.0
else:
    R_V = 1.0 #to be determined

# Constants
M_sun = 1.98892e30 # solar mass [kg]
h = 4.1357e-15 # Planck's constant [eV s]
c = float(3e8) # speed of light [m/s]

# Inputs
beta = 2.0
R = 1.0
z = 1.0
M_gas = 1.0
M_halo = 2e41

# Read spectra file
f = open(filename, 'r')
rawlines = f.readlines()
met = re.findall('Z\s=\s(\d*\.\d+)', rawlines[0])
del rawlines[0]
for i in range(len(rawlines)):
    newlist = rawlines[i].split(' ')
    datalist.append(newlist)

# Read extinction curve data file
rawpoints = open(extinctionfile, 'r').readlines()
def interpolate(R_V, rawpoints, Elist, j):
    pointslist = []
    if R_V == 4.0:
        for i in range(len(rawpoints)):
            newlst = re.split('(?!\S)\s(?=\S)|(?!\S)\s+(?=\S)', rawpoints[i])
            pointslist.append(newlst)
        pointslist = pointslist[3:]
        lambdalist = [float(item[0]) for item in pointslist]
        k_abslist = [float(item[4]) for item in pointslist]
        xvallist = [(c*h)/(lamb*1e-6) for lamb in lambdalist]
        k_interp = scipy.interpolate.interp1d(xvallist, k_abslist)
        return k_interp(Elist[j])

# Dust extinction function
def dust(interpolate, filetype, datalist, beta, R, z, M_gas, M_halo, met):
    speclist = []
    if filetype == 4:
        metallicity = float(met[0])
        Elist = [float(item[0]) for item in datalist]
        speclambdalist = [h*c*1e9/E for E in Elist]
        met1list = [float(item[1]) for item in datalist]
        speclist.extend(met1list)
        klist, Tlist = [None]*len(speclist), [None]*len(speclist)
        if metallicity > 0.0052:
            DGRlist = [50.0*np.exp(-2.21)*metallicity]*len(speclist) # dust to gas ratio
        elif metallicity <= 0.0052:
            DGRlist = [((50.0*metallicity)**3.15)*np.exp(-0.96)]*len(speclist)
        for i in range(len(speclist)):
            if Elist[i] <= 4.1357e-3: # frequencies <= 10^12 Hz
                klist[i] = 0.1*(float(Elist[i])/(1000.0*h))**beta # extinction law [cm^2/g]
            elif Elist[i] > 4.1357e-3: # frequencies > 10^12 Hz
                klist[i] = interpolate(R_V, rawpoints, Elist, i) # interpolated function's value at Elist[i]
        Mdustlist = [M_gas*DGR for DGR in DGRlist] # dust mass
        R_halo = 0.784*(0.27**2.0)*(0.7**(-2.0/3.0))*float(10/(1+z))*((M_halo/(1e8*M_sun))**(1.0/3.0))
        R_dust = 0.018*R_halo # [kpc]

        # Optical depth calculation
        Tlist = [3*Mdustlist[i]*klist[i]/(4*np.pi*R_dust) for i in range(len(speclist))]

        # Ingoing and outgoing luminosities as functions of wavelength
        Linlist = [float(spectra)*R for spectra in speclist]
        Loutlist = [Linlist[i]*np.exp(-Tlist[i]) for i in range(len(speclist))]

        return speclambdalist, Linlist, Loutlist

print dust(interpolate, filetype, datalist, beta, R, z, M_gas, M_halo, met)

Даже когда у меня есть только функция, возвращающая Loutlist вместо набора из 3 списков, она все еще очень медленная. Есть идеи, почему это так? Кроме того, я хочу вернуть кортеж и затем построить speclambdalist против Linlist, а также построить speclambdalist против Loutlist на том же графике. Но у меня сложилось впечатление, что каждый раз, когда я вызываю dust(interpolate, filetype, datalist, beta, R, z, M_gas, M_halo, met)[i], где i = 0, 1, or 2 (я буду делать это несколько раз), придется каждый раз запускать функцию снова. Есть ли способ обойти эти дополнительные прогоны для дальнейшего увеличения скорости? Спасибо!

...