创建和操作适合的文件

Question

我编写了一个程序，操作VLA观测到的数据，并创建FITS文件，用于创建高红移射电星系的光谱能量分布。我真的试着尽可能地提高效率。我列出了必须使用的变量和列表。我尽我所能地做了所有的事，但结果还是令人毛骨悚然。

#!/usr/bin/python

import math
from astropy.io import ascii
import astropy.cosmology as cosmo
import pyfits as fits
import matplotlib.pyplot as plt
from matplotlib import lines
import numpy as np
import os

# removing any duplicates from the current directory
try:
    os.system("rm -f sed_and_alpha.png")
except:
    pass
try:
    os.system("rm -f sed_and_alpha.fits")
except:
    pass
#constants
CBAND = 4.86*10**9
XBAND = 8.46*10**9
CF = (10**-26)

#functions
def get_luminosity_distance(redshift):
    '''finds luminosity distance'''
    lum_dist = cosmo.luminosity_distance(redshift)
    tmp = str(lum_dist)
    tmp = tmp.replace("Mpc","")
    tmp = float(tmp)
    return tmp*3.09*10**22
    print lum_dist

def get_alpha(s1,s2,v1,v2):
    '''gets alpha from two flux densities'''
    top = math.log10(s1/s2)
    bottom = math.log10(v2/v1)
    return top/bottom

def get_luminosity(D_l,z,a,nu1,nu2,Sv2):
    '''gets luminosity from alpha'''
    Sv2 = (Sv2/1000)*CF
    lum = (4*math.pi*((D_l**2)))*(1/((1+z)**(1+a)))*((nu1/nu2)**a)*(Sv2)
    return lum


#lists
available_redshifts = [] # For making sure there are no blank spaces
luminosity_distribution = [] # which calculated luminosities are available
alpha_list = [] # available alphas
new_luminosity_distribution = [] # luminosity distr. after str[-2:]
default_data = {}
list22 = []
list23 = []
list24 = []
list25 = []
list26 = []
list27 = []
list28 = []
exponent_list = [
                    list22,
                    list23,
                    list24,
                    list25,
                    list26,
                    list27,
                    list28,
                ]
exponent_list_names = [
                        "list22",
                        "list23",
                        "list24",
                        "list25",
                        "list26",
                        "list27",
                        "list28",
                    ]
range_column_list = [10.0**22, 10.0**23, 10.0**24,10.0**25,10.0**26,10.0**27,10.0**28]
ylist = []
num_galaxies_list = []
filename = "/home/vhx/Documents/code/python/fits/radio_template.fits/apj384116t3_mrt.txt"
table_data = ascii.read(filename) # getting the table data
# Determining which redshifts are available
for i in xrange(246):
    blank = True
    tmp = (table_data["z"][i])
    if tmp != 'NO' and table_data["SX"][i] != 'NO' and table_data["SC"][i] != 'NO':
        available_redshifts.append(i)
    else:
        pass
# get the spectral indices 
for i in range(len(available_redshifts)-1):
    alpha_list.append(get_alpha(table_data["SX"][available_redshifts[i]],table_data["SC"][available_redshifts[i]],XBAND,CBAND))
    print table_data["z"][available_redshifts[i]], table_data["SX"][available_redshifts[i]], table_data["SC"][available_redshifts[i]]


# solving for luminosity
for i in range(len(available_redshifts)-1):
    x = get_luminosity_distance(table_data["z"][available_redshifts[i]])
    luminosity_distribution.append(get_luminosity(x,table_data["z"][available_redshifts[i]],alpha_list[i],CBAND,XBAND,table_data["SX"][available_redshifts[i]]))
# removing the units from the luminosity distribution so it can be converted
for z in range(len(luminosity_distribution)-1):
    tmp = str(luminosity_distribution[z])
    new_tmp = tmp.replace("Mpc","")
    print new_tmp
    try:
        new_luminosity_distribution.append(float(new_tmp))
    except ValueError:
        print 'Found blank space!'
        pass
new_luminosity_distribution.pop(0)
# checking stuff
alpha_list[:] = [x for x in alpha_list if math.isnan(x) != True]
print alpha_list
k = 0
for item in alpha_list:
    k += item
print k/(len(alpha_list)-1)
print max(alpha_list),min(alpha_list)
print '\n'*5
for item in new_luminosity_distribution:
    print '%e' % item
print '%e, %e' % (max(new_luminosity_distribution), min(new_luminosity_distribution))
sum = 0
for item in new_luminosity_distribution:
    sum += item
print sum/len(new_luminosity_distribution)
new_luminosity_distribution_array = np.array(new_luminosity_distribution)
alpha_list_array = np.array(alpha_list)
numbers = np.random.normal(size = 1000)
width = 1
zed = 0
#making new_luminosity_dist. and alpha_list equal in length
print len(new_luminosity_distribution), len(alpha_list)
alpha_list.pop()
alpha_list.pop()

# checking the range of exponent values
for i in range(len(new_luminosity_distribution)-1):
    tmp = str(new_luminosity_distribution[i])
    tmp = tmp[-2:]
    if int(tmp) == 22:
        list22.append(i)
    elif int(tmp) == 23:
        list23.append(i)
    elif int(tmp) == 24:
        list24.append(i)
    elif int(tmp) == 25:
        list25.append(i)
    elif int(tmp) == 26:
        list26.append(i)
    elif int(tmp) == 27:
        list27.append(i)
    elif int(tmp) == 28:
        list28.append(i)
for lists in exponent_list:
    print lists

# grouping average alphas
for i in range(len(exponent_list)-1):
    sum = 0
    for item in exponent_list[i]:
        sum += alpha_list[item]
    try:
        sum /= len(exponent_list[i])
        ylist.append(sum)
    except ZeroDivisionError:
        sum = 0
    average_dictionary["average of "+str(exponent_list_names[i])] = sum
print average_dictionary
# making the bar graph
for zed in range(len(exponent_list)-1):
    current_list = exponent_list[zed]
    item = 0
    try:
        x = str(new_luminosity_distribution[current_list[item]])
        x = int(x[-2:])
    except:
        x = 0
    height = len(current_list)
    plt.bar(x,height,width,color = 'gray')

for item in exponent_list:
    num_galaxies_list.append(len(item))

#annotating the alphas
plt.annotate("alpha = \n"+str(average_dictionary['average of list22']),xy = (22,5),xytext = (22.3,10),rotation = 90)
plt.annotate("alpha = \n"+str(average_dictionary['average of list24']),xy = (24,5),xytext = (24.3,15),rotation = 90)
plt.annotate("alpha = \n"+str(average_dictionary['average of list25']),xy = (25,5),xytext = (25.3,20),rotation = 90)
plt.annotate("alpha = \n"+str(average_dictionary['average of list26']),xy = (22,5),xytext = (26.3,20),rotation = 90)
plt.annotate("alpha = \n"+str(average_dictionary['average of list22']),xy = (22,5),xytext = (27.3,20),rotation = 90)
plt.savefig("sed_and_alpha")

# Putting everything into a FITS file
hdu = fits.PrimaryHDU()
hdulist = fits.HDUList([hdu])
range_column = fits.Column(name = 'Range', format = 'FLOAT',
                                        array = range_column_list)
num_galaxies = fits.Column(name = 'Number of Galaxies', format = 'FLOAT',array = num_galaxies_list)
mean_alpha_column = fits.Column(name = 'mean_alpha', format = 'FLOAT', array = ylist)
mean_alpha = fits.ColDefs([range_column,mean_alpha_column,num_galaxies])
alpha_change = fits.new_table(mean_alpha)
alpha_change.header.update('EXTNAME','ALPHACHANGE','mapping the change of spectral index')
primary_header = fits.Header()
primary_header.append('LUMINOS')
for item in range_column_list:
    primary_header.append(str(item))
primary_hdu = fits.PrimaryHDU(header = primary_header)
thdulist = fits.HDUList([primary_hdu,alpha_change])
thdulist.writeto("sed_and_alpha.fits")
print '<====DONE====>'

问题：

• 有什么可以让它更有效率的吗？
• 列表理解可以用来缩短时间吗？
• 密码好吗？我是否已经从我过去的岗位上有所改善了？

注意:不能做任何事情来缩短plt.annotate块--这些值是任意的(但不是随机的)。

Answer 1

这里发生了很多事情，我真的不知道从哪里开始从…开始也许可以指出你在先前的问题中已经被告知的一些错误：

• 异常处理；
• 空白区；
• 清单-理解。

异常处理

不要使用裸露的except子句。在考虑使用try之前，您必须知道预期会出现什么样的异常。如果没有，那么就没有必要尝试从“错误”中恢复，因为您将不知道如何处理所说的“错误”。

即使你知道该期待哪种异常，如果你不知道如何处理它，尝试这样做也是没有意义的。如果它知道如何执行…，最好让异常冒泡，并让代码的其他部分来处理它。

这些显然是一般性的注释，但关键是，只有在您现在既期待什么，又如何处理它的情况下，您才应该使用try ... except子句。

白空间

把你的代号给你。午餐店代码。使用空格使读取更容易，并在单个(复杂)表达式上分隔逻辑部分。

Python的官方编码风格是PEP8。读一读，它有一些关于白空间的部分。

List-comprehensions

和关于迭代

的一些事情

在您的代码中有几次，您的构造看起来像

a_list = []
for i in range(len(some_other_list)-1):
    value = call_function(some_other_list[i])
    a_list.append(value)

有三件事不对：

• range(x)生成整数，从0(包含)到x (独占)。我很难在最后看到-1的需要，因为它每次删除最后一项。
• Python的for循环意味着对序列的元素进行迭代(如果我们想要迭代，那么任何可迭代的，但是序列是可迭代的)。你不需要索引来一个接一个地得到元素。您甚至不希望索引，因为使用它们的索引访问元素要比让for循环提供给您: a_list = [] for元素在some_other_list: a_list.append(call_function(元素))中要慢，如果您认为需要索引一起迭代多个序列(…)再想一想，了解zip。如果确实需要索引和元素，请使用enumerate。
• 对列表的append比较慢，最好直接使用列表理解来构建列表: a_list = 打电话_函数(E)在某些_其他_列表

写函数

从目前来看，您的代码几乎只是一大堆指令，没有任何相关的意义。将事物组织成函数可能会导致您的顶层代码沿着

remove_old_files()
alpha, luminosity = process_redshifts("/home/vhx/Documents/code/python/fits/radio_template.fits/apj384116t3_mrt.txt")
exponent_list = separate(luminosity, alpha)
plot(exponent_list, "sed_and_alpha.png")
save(exponent_list, "sed_and_alpha.fits")

这是一个粗略的草图，我不确定这里计算的是什么，但是在顶层使用这种代码提供了更多关于正在发生的事情的信息。

此外，顶层代码最好放在if __name__ == '__main__':语句下。它避免在交互式会话(或任何其他测试环境)中import文件的事件中运行代码。

使用函数还具有赋予参数更多意义的优点，而不是让它们作为全局变量在代码中浮动。

总结一下我到目前为止说的话，我们可以用以下方式编写process_redshifts：

from itertools import izip

def process_redshifts(input_file_name):
    data = ascii.read(input_file_name)
    # Canot generate two lists at once with a list-comprehension
    alpha = []
    luminosity = []

    for z, sx, sc in izip(data['z'], data['SX'], data['SC']):
        if 'NO' not in (z, sx, sc):
            alpha_value = get_alpha(sx, sc, XBAND, CBAND)
            alpha.append(alpha_value)
            x = get_luminosity_distance(z)
            luminosity.append(get_luminosity(z, x, CBAND, XBAND, sx))

    return alpha, luminosity

我在这里使用izip而不是zip来处理三胞胎，而不是预先构建一个完整的三胞胎列表。

os operations

要从计算机中删除文件吗？那么，您应该使用os.remove，它正是为此目的构建的。

def remove_old_files():
    for filename in ('sed_and_alpha.png', 'sed_and_alpha.fits'):
        try:
            os.remove(filename)
        except FileNotFoundError:
            pass

使用不稳定的内置功能

我一点也不精通，但是我看了一下他们在github上的代码。事实证明，使用float转换为str和移除'Mpc'的对象都是astropy.units.Quantity对象。

这些对象通过在其单元前面形成__str__来定义它们的self.value方法。它们还提供了一个__float__方法。因此，您有两种内置的方法来获取要查找的浮点值：

1. the_value_i_want = astropy_object.value
2. the_value_i_want = float(astropy_object)

停止使用str进行此类转换。

打印

它们有很多，它们看起来都像调试消息。此时，您应该已经对代码进行了足够长的测试，以便对其所做的工作有信心，因此您能够删除它们。它们的数量增加了评审人员理解代码的难度。

您甚至碰巧在您的一个函数中，在print之后有一个return行。

在print '<====DONE====>'的末尾还有一个特殊的…这是完全无用的，因为脚本的结尾将通过用户的提示符再次显示出来。

Answer 2

首先，我要做的非常简单的事情是整理你使用科学符号的方式。与编写4.86*10**9不同，您可以编写4.86e9等等。它更容易阅读和减少错误的发生。

其次，您想要做的一些事情在库函数中是重复的。例如，您的get_luminosity_distance函数可以替换为cosmo.luminosity_distance(redshift).si.value。(至少在目前的版本中是如此。)我不知道你用的是哪个版本)。

我会把你的exponent_lists放在字典里，你已经在用exponent_list_names列表了。要检查指数值的范围，可以执行以下操作：

for idx, luminosity in enumerate(new_luminosity_distribution):
    exponent_dict[int(np.log10(luminosity))].append(idx)

它比

# checking the range of exponent values
for i in range(len(new_luminosity_distribution)-1):
    tmp = str(new_luminosity_distribution[i])
    tmp = tmp[-2:]
    if int(tmp) == 22:
        list22.append(i)
    elif int(tmp) == 23:
        list23.append(i)
    elif int(tmp) == 24:
        list24.append(i)
    elif int(tmp) == 25:
        list25.append(i)
    elif int(tmp) == 26:
        list26.append(i)
    elif int(tmp) == 27:
        list27.append(i)
    elif int(tmp) == 28:
        list28.append(i)

这实际上引出了另一个问题:当您正在遍历列表中的元素时，而不是执行

for i in range(len(thelist)):
    ... etc. ...

要做的更多的是琵琶：

for item in thelist:
    ... etc. ...

如果您只需要这些项，或者使用enumerate (正如我前面所做的)，如果您也需要索引的话。

Answer 3

我不是python开发人员，但我会尽力而为的。

首先，list22，list23，list24很奇怪，你应该做这样的事情

exponent_list[0]  #for list22
exponent_list[1]  #for list23
exponent_list[2]  #for list24

等等。在这种方式下，它更抽象和有用，可以像现在这样保存来执行3或4个数组。

这个密码，

range_column_list = [10.0**22, 10.0**23, 10.0**24,10.0**25,10.0**26,10.0**27,10.0**28]

可以这样做

range_column_list = []
for i in range(22,28)
    range_column_list.append(10.0 ** i)

或者更好，10.0是一个常量

for i in xrange(246):

在这一行中，我认为你想把文件中的所有行都取出来，但是246对我来说是一个神奇的数字，如果明天要改变，你需要修改程序。

如果是必需的，或者始终是这样的，那么放置一个常量，但是如果不是这里，则是一个用于计算文件https://stackoverflow.com/questions/845058/how-to-get-line-count-cheaply-in-python行数的链接。

首先，为文件中的所有项运行一个循环，如果是红移，则将其添加到数组中。

在那个数组上循环之后，为什么不一起做其他的事情呢？我想这是相同的代码

# Determining which redshifts are available
for i in xrange(246):
    blank = True
    tmp = (table_data["z"][i])
    if tmp != 'NO' and table_data["SX"][i] != 'NO' and table_data["SC"][i] != 'NO':
        available_redshifts.append(i)
        #now you know that is a redshift so... I do more things
        # get the spectral indices 
        alpha_list.append(get_alpha(table_data["SX"][available_redshifts[i]],table_data["SC"][available_redshifts[i]],XBAND,CBAND))
        print table_data["z"][available_redshifts[i]], table_data["SX"][available_redshifts[i]], table_data["SC"][available_redshifts[i]]
        # now you iterate again for each redshift but I already have one, so... I don't need a loop
        # solving for luminosity
        x = get_luminosity_distance(table_data["z"][available_redshifts[i]])
        luminosity_distribution.append(get_luminosity(x,table_data["z"][available_redshifts[i]],alpha_list[i],CBAND,XBAND,table_data["SX"][available_redshifts[i]]))
        # I was done the last loop here too but I realise that all spectral indexes will not together so I do in a different loop like you    

    else:
        pass

你真的需要或验证数据之前，我看到了很多假设，我知道数据总是以相同的格式，但你永远不知道。请检查所有数组索引。

这段代码有点混乱。

new_luminosity_distribution.pop(0)

如果只是附加数据，为什么要删除第一项？

有很多这样的

print k/(len(alpha_list)-1)

我不知道这是什么，现在用文字来澄清这个值

alpha_list.pop()
alpha_list.pop()

你移除两样东西，这对我来说太神奇了

这

for item in new_luminosity_distribution:
    print '%e' % item
print '%e, %e' % (max(new_luminosity_distribution), min(new_luminosity_distribution))
sum = 0
for item in new_luminosity_distribution:
    sum += item

如果你能重复一次，为什么要重复两次？

sum = 0
for item in new_luminosity_distribution:
    print '%e' % item
    sum += item
print '%e, %e' % (max(new_luminosity_distribution), min(new_luminosity_distribution))

在这里，同样地，假设执行exponent_list而不是list22，您将得到下一个

# checking the range of exponent values
for i in range(len(new_luminosity_distribution)-1):
    tmp = str(new_luminosity_distribution[i])
    tmp = tmp[-2:]
    if int(tmp) == 22:
        list22.append(i)
    elif int(tmp) == 23:
        list23.append(i)
    elif int(tmp) == 24:
        list24.append(i)
    elif int(tmp) == 25:
        list25.append(i)
    elif int(tmp) == 26:
        list26.append(i)
    elif int(tmp) == 27:
        list27.append(i)
    elif int(tmp) == 28:
        list28.append(i)

像这样

# checking the range of exponent values
for i in range(len(new_luminosity_distribution)-1):
    tmp = str(new_luminosity_distribution[i])
    tmp = tmp[-2:]
    exponent_list[tmp-22].append(i)
    # tmp-22 will give you the index
    # 0 for 22, 1 for 23, 2 for 24 etc like list22, list23, list24 etc

如果您有两个循环，并且是相同的循环，但是内容是不同的，那么尝试减少循环的数量并执行优化，如果同时有一种方法，则查找do。