analysis/substructure.py at master · caterpillarproject/analysis · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
import numpy as np
import haloutils
from caterpillaranalysis import *
from scipy import interpolate
import matplotlib.pyplot as plt
#import grifflib as glib


# to Read in values
#substruct = SubstructurePlugin()
#nsubs_vmax, nsubs_vmax_ever,subhalomass,submass_frac,subs_inner_frac = substruct.read(hpath)


def distance(posA, posB,boxsize=100*1000.):
    dist = abs(posA-posB)
    tmp = dist > boxsize/2.0
    dist[tmp] = boxsize-dist[tmp]
    if dist.shape == (3,):
        return np.sqrt(np.sum(dist**2))
    else:
        return np.sqrt(np.sum(dist**2,axis=1))

class SubstructurePlugin(PluginBase):
    def __init__(self):
        super(SubstructurePlugin,self).__init__()
        self.filename='substructure_data.dat'
        self.xmin=0;     self.xmax=256
        self.ymin=10**8; self.ymax=5*10**11
        self.xlog= False; self.ylog = False
        self.xlabel='' ; self.ylabel=''

    def _analyze(self,hpath):
        if not haloutils.check_last_rockstar_exists(hpath):
            raise IOError("No rockstar")

        # what to quantify substructure
        snap_z0 = haloutils.get_numsnaps(hpath)-1
        cat = haloutils.load_rscat(hpath,snap_z0,rmaxcut=False)
        hostID = haloutils.load_zoomid(hpath)
        hosthalo = cat.ix[hostID]
        hostpos = np.array(hosthalo[['posX','posY','posZ']])
        subs = cat.get_subhalos_within_halo(hostID) # just 1 level deep
        subs_all = cat.get_all_subhalos_within_halo(hostID)

        # num halos above some vmax thresh currently
        maskv = subs_all['vmax']> 25
        nsubs_vmax = np.sum(maskv)

        # num halos that ever entered above some vmax thresh
        import MTanalysis2 as MT
        AE = MT.AllExtantData()
        AD = MT.AllDestroyedData()
        dataE = AE.read(hpath)
        dataD = AD.read(hpath)
        maskE = dataE['infall_vmax'] > 25
        maskD = dataD['infall_vmax'] > 25
        nsubs_vmax_ever = np.sum(maskE)+np.sum(maskD)

        # SHMF normalization (divided by halo mass)
        # need to extract from something else
        shmf = SHMFPlugin()
        x,y,sx,sy = shmf.read(hpath)
        # x in Msun, y in dN/dM

        # total subhalo mass?
        subhalomass = np.sum(subs['mgrav']/cat.h0)

        # subhalo mass fraction
        submass_frac = subhalomass/(hosthalo['mgrav']/cat.h0)

        # num subhalos within 50% of rvir vs outside of it
        subposns = np.array(subs_all[['posX','posY','posZ']])
        dist = distance(hostpos,subposns)*cat.scale/cat.h0*1000
        subs_inner_frac = np.sum(dist<.5*hosthalo['rvir']/cat.h0)/float(len(subs_all))

        nsubs = len(subs_all)

        g = open(hpath+'/'+self.OUTPUTFOLDERNAME+'/'+self.filename,'wb')
        np.array([nsubs,nsubs_vmax, nsubs_vmax_ever,subhalomass,submass_frac,subs_inner_frac]).tofile(g)
        g.close()

    def _read(self,hpath):
        data = np.fromfile(hpath+'/'+self.OUTPUTFOLDERNAME+'/'+self.filename)
        return data

def substructure_tabfn(hpath):
    if hpath==None: return None
    plug = SubstructurePlugin()
    data = plug.read(hpath)
    nsubs,nsubs_vmax, nsubs_vmax_ever,subhalomass,submass_frac,subs_inner_frac = data
    names = ['nsubs','nsubs_vmax','nsubs_vmax_ever','subhalomass','submass_frac','subs_inner_frac']
    formats = [int, int, int, np.float, np.float, np.float]
    return data,names,formats

def substructure_formation_tabfn(hpath):
    import formationtime
    if hpath==None: return None
    out1 = substructure_tabfn(hpath)
    if out1==None: return None
    out2 = formationtime.formation_tabfn(hpath)
    if out2==None: return None
    d1,n1,f1 = out1
    d2,n2,f2 = out2
    data = np.concatenate([d1,d2])
    names = np.concatenate([n1,n2])
    formats = np.concatenate([f1,f2])
    return data, names, formats

def plot_substructure_formation(filename=None,exclude_hids=[94687,388476,1599988]):
    df = haloutils.tabulate(substructure_formation_tabfn,lx=14,exclude_hids=exclude_hids)
    import seaborn.apionly as sns
    g = sns.PairGrid(df)
    g.map_diag(plt.hist)
    g.map_offdiag(plt.scatter)
    if filename != None:
        g.savefig(filename)
    return g

#hpath = '/bigbang/data/AnnaGroup/caterpillar/halos/H1599988/H1599988_EX_Z127_P7_LN7_LX14_O4_NV4'
#substruct = SubstructurePlugin()
#nsubs_vmax, nsubs_vmax_ever,subhalomass,submass_frac,subs_inner_frac = substruct.read(hpath)
#lx = 14
#figdir = '/bigbang/data/AnnaGroup/GregFigs/lx'+str(lx)+'/'
#halo_paths = haloutils.find_halo_paths(levellist=[lx],require_mergertree=True,require_subfind=False,verbose=False)

"""
halo_paths = ['/bigbang/data/AnnaGroup/caterpillar/halos/H95289/H95289_BB_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H1195448/H1195448_EC_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H1354437/H1354437_EA_Z127_P7_LN7_LX14_O4_NV5',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H447649/H447649_EB_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H5320/H5320_EB_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H581141/H581141_EB_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H1631506/H1631506_EA_Z127_P7_LN7_LX14_O4_NV5',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H1130025/H1130025_EB_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H1725139/H1725139_EB_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H581180/H581180_EB_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H1422331/H1422331_EX_Z127_P7_LN7_LX14_O4_NV5',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H1232164/H1232164_EX_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H1387186/H1387186_EB_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H1725272/H1725272_EB_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H1725372/H1725372_EB_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H264569/H264569_EB_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H1292085/H1292085_EX_Z127_P7_LN7_LX14_O4_NV5',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H94687/H94687_EB_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H1268839/H1268839_EB_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H196589/H196589_EX_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H1599988/H1599988_EX_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H1079897/H1079897_EX_Z127_P7_LN7_LX14_O4_NV4',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H94638/H94638_EX_Z127_P7_LN7_LX14_O4_NV5',
 '/bigbang/data/AnnaGroup/caterpillar/halos/H796175/H796175_EB_Z127_P7_LN7_LX14_O4_NV4']


for hpath in halo_paths:
    nsubs, nsubs_vmax, nsubs_vmax_ever,subhalomass,submass_frac,subs_inner_frac = substruct.read(hpath)
    print subs_inner_frac, 'inner fraction'
    print 'done with', hpath
"""