0.这里的EAGLE是什么
EAGLE (Evolution and Assembly of GaLaxies and their Environments) 是一系列不同大小及模型的宇宙学模拟实验，旨在研究星系的形成与演化。EAGLE的详细介绍可以参考其主页(EAGLE主页)及相关paper（ Schaye et al.2015; Crain et al.2015）(citation如果打错了那就是不会)。在接下来的内容中，我们将从EAGLE官方给出的代码示例出发，讨论如何高效的读取particle中需要的星系数据(当然，与其说是讨论，不如说是个人的记录)。
1.EAGLE数据及官方文档读取particle数据
EAGLE公开的数据包括星系数据库与粒子数据，前者为EAGLE团队对各个snap数据中星系的统计数据，例如星系的半质量半径，不同波段光度，形态学的参数。星系数据库的数据通过其自带的SQL数据库进行调用，也提供相应的python函数库实现在python上访问。而后者则需要下载到本地之后再进行读取，也是此处想要讨论的数据。
EAGLE的particle数据由一系列的snap（一般为28个）组成，每个snap对应一个红移，每一个snap下的数据包括一系列的hdf文件组成，每个hdf文件内的格式相同，相当于每个hdf文件包括了模拟盒子中的一部分粒子的数据，例如对于50Mpc，中等分辨率的盒子，每个snap有着128个大小接近的hdf数据文件。而为了读取盒子里的一个星系，都需要遍历所有的hdf文件读取数据。hdf文件的具体结构为 PartType/Attribute，前者为粒子类型，后者为这一粒子的某一属性，例如PartType=‘0’，Attribute=‘Coordinates’为该文件下气体粒子的坐标。EAGLE由gadget-3（gadget-2的变种）完成，gadget提供了5种粒子类型，而EAGLE使用了0（气体），1（暗物质），4（恒星），5（黑洞）。值得注意的是，暗物质粒子的质量是一样的，为此不能也不用单独读取。EAGLE中的特定星系的粒子具有相同的GroupNumber 与SubGroupNumber，这两个属性同时也标识着这个星系所处的子晕。
在EAGLE的官方文档（arXiv:1706.09899）中，提供了对particle数据提取的python代码示例，主要由 read_dataset与Galaxy类下的read_galaxy 两个函数完成。同时也提供了一个单独的python库进行读写，然而作为win10用户，果不其然安装不上，我也不知道那个库有多万能，所以还是接着看看代码示例。
import numpy as np
import h5py
import astropy.units as u
from astropy. constants import G
def read_dataset(itype, att, nfiles=16):
""" Read a selected dataset, itype is the PartType and att is the attribute name. """
# Output array.
data = []
# Loop over each file and extract the data.
for i in range(nfiles):
f = h5py.File('./data/snap_028_z000p000.%i.hdf5'%i, 'r')
tmp = f['PartType%i/%s'%(itype, att)][...]
data.append(tmp)
# Get conversion factors.
cgs = f['PartType%i/%s'%(itype, att)].attrs.get('CGSConversionFactor')
aexp = f['PartType%i/%s'%(itype, att)].attrs.get('aexp-scale-exponent')
hexp = f['PartType%i/%s'%(itype, att)].attrs.get('h-scale-exponent')
# Get expansion factor and Hubble parameter from the header.
a = f['Header'].attrs.get('Time')
h = f['Header'].attrs.get('HubbleParam')
f.close()
# Combine to a single array.
if len(tmp.shape) > 1:
data = np.vstack(data)
else:
data = np.concatenate(data)
# Convert to physical.
if data.dtype != np.int32 and data.dtype != np.int64:
data = np.multiply(data, cgs * a**aexp * h**hexp, dtype='f8')
return data
def read_galaxy(self, itype, gn, sgn, centre):
""" For a given galaxy (defined by its GroupNumber and SubGroupNumber)
extract the coordinates and mass of all particles of a selected type.
Coordinates are then wrapped around the centre to account for periodicity. """
data = {}
# Load data, then mask to selected GroupNumber and SubGroupNumber.
gns = read_dataset(itype, 'GroupNumber')
sgns = read_dataset(itype, 'SubGroupNumber')
mask = np.logical_and(gns == gn, sgns == sgn)
if itype == 1:
data['mass'] = read_dataset_dm_mass()[mask] * u.g.to(u.Msun)
else:
data['mass'] = read_dataset(itype, 'Mass')[mask] * u.g.to(u.Msun)
data['coords'] = read_dataset(itype, 'Coordinates')[mask] * u.cm.to(u.Mpc)
# Periodic wrap coordinates around centre.
boxsize = self.boxsize/self.h
data['coords'] = np.mod(data['coords']-centre+0.5*boxsize,boxsize)+centre-0.5*boxsize
return data
read_dataset的作用是读取一整个snap数据中，粒子类型为itype，属性为att的所有数据，而在read_galaxy中通过读取GroupNumber(gn)与SubGroupNumber(sgn)，与指定的gn与sgn对比，确定星系粒子数据。最后对数据进行单位转化。直接读取的数据为无量纲数据，在read_dataset中完成了无量纲数据到高斯单位制的转化，而在read_galaxy中使用astropy转化为了更实用的单位，同时处理坐标的周期性边界。
这样的读取存在两个显著的不足，首先是读取某一粒子的某一属性时，需要读取整个盒子里该粒子的属性，对一个50Mpc中等分辨率的盒子而言，数据的总大小是50G，假设除黑洞外的三种粒子数目一致，每种粒子对应28个属性，每个属性全部读取对应的大小就是50/3/28=0.60G，而实际上，暗物质的粒子数会更多一些，坐标这样的属性是n*3的矩阵，占用的内存会更大，对于我自己的电脑，可用内存开机后有着接近10G，依然有读取时内存撑爆的风险，而换到更大的盒子或者更高分辨率的模拟内存根本撑不住；另一方面在于读取的效率低下，read_galaxy中反复调用read_dataset来完成一个星系的读取，而每调用一次，就相当于要全部访问一遍数据，同样的，对于50Mpc中等分辨率的盒子，数据存储在电脑的硬盘里，读取一个星系的平均时间在3分钟以上，如果数据存储在外接硬盘，同时usb接口速度一般的话，这个时间会继续增加。三分钟对单个星系不算长，但当需要获得的星系量超过一百个时，超过三百分钟的读取时间能直接给人戴上痛苦面具。
2.针对官方文档代码示例中存在问题的一些解决方案
对于第一个问题可以少量调整read_dataset函数来解决，基本思路就是每打开一个数据文件就读取这个文件中目标星系的部分，也就是加上gn与sgn的判断，这样对内存的占用就会减小到一个数据文件中粒子的属性加上读出的粒子属性，大幅减小了内存撑爆的风险。
def read_dataset (itype , att , gn , sgn , nfiles=128,add='F:/RefL0050N0752/snapshot_028_z000p000'):
data = []
for i in range(nfiles):
nex=add[25:]
f = h5py.File(add+'/snap'+nex+'.%i.hdf5'%i, 'r')
gns=f['PartType%i/%s'%(itype , 'GroupNumber' )][...]
sgns=f['PartType%i/%s'%(itype , 'SubGroupNumber' )][...]
#读取GroupNumber和SubGroupNumber
mask=np.logical_and(gns==gn , sgns==sgn)
tmp = f['PartType%i/%s'%(itype , att )][...][mask]
#读出当前文件下需要的粒子的属性，并同时限制到需要的gn和sgn
data.append(tmp)
cgs = f['PartType%i/%s'%(itype , att )].attrs.get('CGSConversionFactor')
aexp = f['PartType%i/%s'%(itype , att )].attrs.get('aexp-scale-exponent')
hexp = f['PartType%i/%s'%(itype , att )].attrs.get('h-scale-exponent')
a = f['Header'].attrs.get('Time')
h = f['Header'].attrs.get('HubbleParam')
f.close ()
if len(tmp.shape) > 1:
data = np.vstack(data)
else:
data = np. concatenate (data)
if data.dtype != np.int32 and data.dtype != np.int64:
data = np.multiply (data , cgs * a**aexp * h**hexp , dtype='f8')
return data具体的实现可以参照上述代码，它与原本的read_dataset仅仅在中间额外读取gn与sgn处存在差别。当然，修改read_dataset同时也需要修改read_galaxy，不过基本上都是修改一下传的参数。此外，当前的函数存在潜在风险，当tmp为空时，可能会无法正确识别vstack和concatenate。
然而此时读取时间过长，反复打开数据文件的问题依然没有解决，甚至在理论上，由于反复读取gn与sgn，读取速度应当小于原先版本。显然，为了减少读取时间，应当在打开一个文件时，尽可能多的读出需要的星系的数据。为此需要对read_dataset与read_galaxy的功能进行合并，在下方将先给出代码，之后对其逻辑与细节进行说明。
import numpy as np
import h5py
from astropy.constants import G
import astropy.units as u
dic={}
dic['0']=['Mass','Coordinates','Velocity']
dic['4']=['Mass','Coordinates','Velocity']
dic['5']=['Mass','Coordinates','Velocity']
'''
dic 的结构为 keys：particle type
keys 下的元素为attribute，不用单独包含 GroupNumber 和 SubGroupNumber
'''
def mapping(itype):
if itype=='0':
return 'gas'
elif itype =='4':
return 'stars'
elif itype=='5':
return 'bh'
def read_dataset_custom(dic,gn,sgn,nfiles =128,add='F:/RefL0050N0752/snapshot_028_z000p000'):
#dic 格式： itype：att
keys=dic.keys()
data={}
for i in keys:
name=mapping(i)
data[name]={}
data[name]['gns']=[]
data[name]['sgns']=[]
for i in keys:
for att in dic[i]:
name=mapping(i)
data[name][att]=[]
#存储数据的字典初始化
for i in range(nfiles):
nex=add[-13:]
f = h5py.File(add+'/snap'+nex+'.%i.hdf5'%i, 'r')
for j in keys:
itype=int(j)
gns=f['PartType%i/%s'%(itype , 'GroupNumber' )][...]
sgns=f['PartType%i/%s'%(itype , 'SubGroupNumber' )][...]
tag=0
for k1,k2 in zip(gn,sgn):
if tag==0:
tag=1
mask=np.logical_and(gns==k1,sgns==k2)
else:
mask=mask|np.logical_and(gns==k1 ,sgns==k2)
#对于一个粒子类型 j，获取所有需要的gn与sgn的mask
name=mapping(j)
data[name]['gns'].append(gns[mask])
data[name]['sgns'].append(sgns[mask])
#存储gn，sgn
for att in dic[j]:
tmp = f['PartType%i/%s'%(itype , att )][...][mask]
data[name][att].append(tmp)
#存储对应属性
f.close()
for j in keys:
name=mapping(j)
data[name]['gns']=np.concatenate(data[name]['gns'])
data[name]['sgns']=np.concatenate(data[name]['sgns'])
#合并gn与sgn的数组
nex=add[-13:]
f = h5py.File(add+'/snap'+nex+'.%i.hdf5'%0, 'r')
for j in keys:
itype=int(j)
name=mapping(j)
for att in dic[j]:
cgs = f['PartType%i/%s'%(itype , att )].attrs.get('CGSConversionFactor')
aexp = f['PartType%i/%s'%(itype , att )].attrs.get('aexp-scale-exponent')
hexp = f['PartType%i/%s'%(itype , att )].attrs.get('h-scale-exponent')
a = f['Header'].attrs.get('Time')
h = f['Header'].attrs.get('HubbleParam')
#转化为高斯单位制
for kk in data[name][att]:
if len(kk)!=0:
tmp=kk
if len(tmp.shape) > 1:
data[name][att] = np.vstack(data[name][att])
else:
data[name][att] = np. concatenate (data[name][att])
#合并属性数组
if data[name][att].dtype != np.int32 and data[name][att].dtype != np.int64:
data[name][att] = np.multiply (data[name][att] , cgs * a**aexp * h**hexp , dtype='f8')
if att=='Mass':
data[name][att]=data[name][att]*u.g.to(u.Msun)
elif att=='Coordinates':
data[name][att]=data[name][att]*u.cm.to(u.Mpc)
#转化为常见单位
f.close()
return data如果需要快速而简明的看懂这段代码，建议把name看成gas，att看作Mass。
我们先定义了一个字典用于标记需要读取的粒子类型与属性，之后用mapping函数将对应的itype转化为对应的名称，用于减小下方函数的复制粘贴量。在read_dataset_custom中，传入的gn与sgn是所有需要的gn与sgn——值得注意的是，对于0号子晕通常有着更多多的粒子，尤其是对于gn靠前的晕，大量读取这类晕，特别是包括暗物质读取时，需要小心内存爆炸。存储的data字典初始化之后的处理逻辑十分简单，无非就是打开一个文件，将所有需要的gn与sgn做一个mask（我也不知道这种布尔值的数组叫什么），然后对应的写入需要的属性。在整个读取结束之后，再进行单位的转化和数组的合并（这一步被留到最后，因为如果紧跟着读取就完成数组合并，数组没有.append可以用，当然，np.r_应该能解决这个问题）。基本上就是先前read_dataset与read_galaxy缝到一起，将模拟盒子里需要的部分单独放到自己的一个盒子里，为此也能方便的看到自己使用的数据的宏观分布（虽然用不上）。同时，在这里我们并没有处理周期性边界对于坐标的影响。
如果仍需单独分离出单个星系，因为gn与sgn的数据被一并读了出来，可以直接在读出的data上进行操作，也可以单独存成独立文件，具体方案取决于需求。
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
data=read_dataset_custom(dic,[2,2,2,2,2],[1,2,3,4,5])
coords=data['stars']['Coordinates']
fig=plt.figure(figsize=(10,10))
ax=Axes3D(fig)
ax.scatter3D(coords.T[0],coords.T[1],coords.T[2])在上方的调用中，我们将gn为2的5个子晕存到了data中，再将其中恒星的粒子画到3d图上，作为一个简单展示。
图中有三个星系可能是由于视角的原因重叠在了一起，同时颜色的深浅是不代表星系真实大小的，甚至实际上也不存在深浅之分。
3.总结
我们通过对EAGLE官方文档给出的代码示例进行改造，完成了一个用于批量提取所需星系数据的函数，大幅减小了IO的时间，即使对于需要全部检查的情况，分批进行批量读取在处理速度上也有着优势。
写完博客之后呢，觉得写的其实没什么难度，这个想法非常的自然，可能是EAGLE给出的代码示例多少给人一种思维惯性，以至于有一天发现自己在坐牢才开始思考怎么改。虽然但是把读取的时间缩短后血压显著下降，终于不用开着电脑，听着风扇坐牢了
reference
Schaye J., Crain R.~A., Bower R.~G., Furlong M., Schaller M., Theuns T., Dalla Vecchia C., et al., 2015, MNRAS, 446, 521.
Crain R.~A., Schaye J., Bower R.~G., Furlong M., Schaller M., Theuns T., Dalla Vecchia C., et al., 2015, MNRAS, 450, 1937
The EAGLE team, 2017, arXiv, arXiv:1706.09899
扫一扫
614
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
