In [1]:
import localgroup
import iso_criteria
import sys
sys.path.append('/u/ki/yymao/pyscripts')
from helpers.SimulationAnalysis import readHlist
from fast3tree import fast3tree
In [98]:
import numpy as np
from sdfpy import load_sdf, SDFIndex, SDFRead
from itertools import product
fname = "/nfs/slac/g/ki/darksky/simulations/ds14_b/ds14_b_halos_1.0000"
halos_data = SDFRead(fname)
halos_midx = SDFRead(fname + '.midx7')
halos = SDFIndex(halos_data, halos_midx)
fields = 'x y z vx vy vz mvir rvir vmax id mvir_pid'.split()
pad = np.array([2.0, 2.0, 2.0]) #in Mpc/h
def get_data(level, cell_iarr, pad, fields, cut=None, cut_fields=None):
dtype = np.dtype([(f, int if f.endswith('id') else float) for f in fields])
count = 0
for h in halos.iter_padded_bbox_data(level, idx, pad, fields):
if hasattr(cut, '__call__'):
flag = cut(*(h[f] for f in cut_fields))
count_this = np.count_nonzero(flag)
else:
flag = slice(None)
count_this = len(h[fields[0]])
if count_this:
if count:
data.resize(count+count_this)
else:
data = np.empty(count_this, dtype)
for f in fields:
data[f][count:] = h[f][flag]
count += count_this
if not count:
data = np.empty(0, dtype)
return data
level = 3
for idx in product(xrange(2**level), repeat=3):
data = get_data(level, idx, pad, fields, lambda x: x > 40.0, ['vmax'])
break
In [100]:
data['mvir']
Out[100]:
In [2]:
level = 3
for idx in product(xrange(2**level), repeat=3):
print idx
In [ ]:
In [ ]:
In [2]:
#path = '/nfs/slac/g/ki/ki20/cosmo/behroozi/Consuelo/4001/hlists/hlist_1.00000.list'
#halos = readHlist(path)
In [3]:
#np.save('/lustre/ki/pfs/mwillia1/LG_project/Consuelo_Boxes/4001hlist.npy', halos)
In [9]:
halos = np.load('/lustre/ki/pfs/mwillia1/LG_project/Consuelo_Boxes/4001hlist.npy')
In [16]:
halos = halos[halos['vmax']>100]
In [17]:
halos.shape
Out[17]:
In [19]:
host_flag = np.where(halos['upid'] == -1)[0]
In [40]:
def find_pairs(halos, host_flag, D):
pairs = []
points = halos[host_flag][list('xyz')].view(float).reshape(-1, 3)
with fast3tree(points) as tree:
for i, p in enumerate(points):
idx = tree.query_radius(p, D, True, 'index')
idx = idx[idx > i]
pairs.extend(((i, j) if halos['mvir'][host_flag[i]] > halos['mvir'][host_flag[j]] else (j, i) for j in idx))
return pairs
def find_periodic_midpoint(p1, p2, box_size, periodic=True):
mid = (p1 + p2)*0.5
if not periodic:
return mid
half_box_size = box_size*0.5
mid[np.fabs(p1 - p2) > half_box_size] += half_box_size
mid = np.fmod(mid, box_size)
return mid
def isolated(pairs, halos, host_flag, D_iso, D_M33, box_size, periodic=True):
print 'before: ', len(pairs)
# iso_pairs = []
MW_M31_larger = []
M31_M31_larger = []
M33_M31_larger = []
M33_M33_larger = []
M31_M33_larger = []
MW_M33_larger = []
with fast3tree(halos[list('xyz')].view(float).reshape(-1, 3)) as tree:
for pair in pairs:
print pair
h1, h2 = halos[host_flag[list(pair)]] #h1 is the larger halo
p1 = np.fromiter((h1[ax] for ax in 'xyz'), float)
p2 = np.fromiter((h2[ax] for ax in 'xyz'), float)
mid = find_periodic_midpoint(p1, p2, box_size, periodic)
idx = tree.query_radius(mid, D_iso, periodic, 'index')
biggest = halos['mvir'][idx].argmax()
idx = np.delete(idx, biggest)
biggest = halos['mvir'][idx].argmax() #biggest is now second-biggest
# the second index has a smaller mass
if h2['mvir'] == halos['mvir'][idx[biggest]]:
# iso_pairs.append(pair) #satisfying iso
M31_M31_larger.append(pair[0])
MW_M31_larger.append(pair[1])
else:
continue #no need to search M33
#find M33, in two ways
idx = tree.query_radius(p1, D_M33, periodic, 'index')
# idx = idx[halos['vmax'][idx] < 80]
# M33_ind = idx[halos['mvir'][idx].argmax()] if len(idx) else -1
m31_delete_ind = halos['mvir'][idx].argmax()
idx = np.delete(idx, m31_delete_ind) #So we don't identify M31 as M33
if len(idx):
m33_candidate = halos['mvir'][idx].argmax()
if halos[idx[m33_candidate]]['id']==halos[host_flag][pair[1]]['id']: #If MW is identified as M33
idx = np.delete(idx, m33_candidate)
M33_ind = idx[halos['mvir'][idx].argmax()] if len(idx) else -1
if halos['vmax'][M33_ind] > 80: M33_ind = -1
M33_M31_larger.append(M33_ind)
idx = tree.query_radius(p2, D_M33, periodic, 'index')
# idx = idx[halos['vmax'][idx] < 80]
# LMC_ind = idx[halos['mvir'][idx].argmax()] if len(idx) else -1
m31_delete_ind = halos['mvir'][idx].argmax()
idx = np.delete(idx, m31_delete_ind) #So we don't identify M31 as M33
if len(idx):
m33_candidate = halos['mvir'][idx].argmax()
if halos[idx[m33_candidate]]['id']==halos[host_flag][pair[1]]['id']: #If MW is identified as M33
idx = np.delete(idx, m33_candidate)
LMC_ind = idx[halos['mvir'][idx].argmax()] if len(idx) else -1
if halos['vmax'][LMC_ind] > 80: LMC_ind = -1
M33_mass = halos['mvir'][M33_ind] if M33_ind != -1 else 0
LMC_mass = halos['mvir'][LMC_ind] if LMC_ind != -1 else 0
if LMC_mass > M33_mass:
M33_ind = LMC_ind
M31_M33_larger.append(pair[1]) #append the smaller mass index in the pair -> M31 since LMC > M33
MW_M33_larger.append(pair[0]) #append the larger mass index in the pair -> MW since LMC > M33
else:
M31_M33_larger.append(pair[0]) #append the smaller mass index in the pair -> M31 since LMC < M33
MW_M33_larger.append(pair[1]) #append the larger mass index in the pair -> MW since LMC < M33
M33_M33_larger.append(M33_ind)
print 'after: ', len(MW_M31_larger)
return MW_M31_larger, M31_M31_larger, M33_M31_larger, M33_M33_larger, M31_M33_larger, MW_M33_larger
In [ ]:
In [20]:
pairs = find_pairs(halos, host_flag, 1.0)
In [38]:
len(pairs)
Out[38]:
In [41]:
MW_M31_larger, M31_M31_larger, M33_M31_larger, M33_M33_larger, M31_M33_larger, MW_M33_larger = isolated(pairs, halos, host_flag, 2.0/0.7, 0.4/0.7, 420, periodic=True)
In [22]:
a=np.array(M33_M31_larger)
len(a[a[:] > 0])
Out[22]:
In [ ]:
In [88]:
#hosts_subs = halos[halos['vmax']>80]
#zall = np.vstack((hosts_subs['x'], hosts_subs['y'], hosts_subs['z'])).T
#zall = zall[1100000:1200000]
#masses_subs = hosts_subs['mvir'][1100000:1200000]
In [5]:
#hosts = halos[halos['upid']==-1]
#hosts = hosts[hosts['vmax']>80]
#hosts.shape
Out[5]:
In [6]:
#z = np.vstack((hosts['x'], hosts['y'], hosts['z'])).T
#z = z[1100000:1200000]
#masses = hosts['mvir'][1100000:1200000]
In [12]:
#for i in range(z.shape[1]):
# print z[:,i].max()
In [59]:
#pairs = iso_criteria.criteria.find_pairs(z, 1.0)
In [60]:
len(pairs)
Out[60]:
In [13]:
isopairs = iso_criteria.criteria.isolated(pairs, 2, z, masses, 420, True)
In [84]:
MW_inds, M31_inds = iso_criteria.criteria.label_MW_M31(isopairs, masses)
In [90]:
M33_inds = iso_criteria.criteria.find_M33(zall, 0.4, masses_subs, M31_inds, z, True)
In [94]:
M33_inds = np.array(M33_inds)
a=np.array(M33_inds)[M33_inds>0]
In [104]:
len(M33_inds),len(M31_inds),len(MW_inds)
Out[104]:
In [82]:
def get_data(halos, host_flag, MW_inds, M31_inds, M33_inds):
fn = ['MW_'+nm for nm in halos.dtype.names]
fn1 = ['M31_'+nm for nm in halos.dtype.names]
fn2 = ['M33_'+nm for nm in halos.dtype.names]
fields = np.hstack((np.array(fn), np.array(fn1), np.array(fn2)))
dtype = np.dtype([(f, int if f.endswith('id') else float) for f in fields])
num_sys = len(MW_inds)
data = np.empty(num_sys, dtype)
no_M33_data = np.copy(halos[0])
for f in no_M33_data.dtype.names: no_M33_data[f]=-1
mw_dat = halos[host_flag[MW_M31_larger]]
m31_dat = halos[host_flag[M31_M31_larger]]
m33_dat = np.array([halos[i] if i>0 else no_M33_data for i in M33_M31_larger])
for i in range(len(data)):
data[i] = tuple(mw_dat[i])+tuple(m31_dat[i])+tuple(m33_dat[i])
return data
In [37]:
mw_dat = halos[host_flag[MW_M31_larger]]
m31_dat = halos[host_flag[M31_M31_larger]]
m33_dat = halos[M33_M31_larger]
In [65]:
no_M33_data = np.copy(halos[0])
for f in no_M33_data.dtype.names: no_M33_data[f]=-1
no_M33_data
Out[65]:
In [77]:
m33_dat = np.array([halos[i] if i>0 else no_M33_data for i in M33_M31_larger])
In [39]:
fn = ['MW_'+nm for nm in halos.dtype.names]
fn1 = ['M31_'+nm for nm in halos.dtype.names]
fn2 = ['M33_'+nm for nm in halos.dtype.names]
fields = np.hstack((np.array(fn), np.array(fn1), np.array(fn2)))
dtype = np.dtype([(f, int if f.endswith('id') else float) for f in fields])
num_sys = len(MW_M31_larger)
data = np.empty(num_sys, dtype)
In [78]:
for i in range(len(data)):
data[i] = tuple(mw_dat[i])+tuple(m31_dat[i])+tuple(m33_dat[i])
In [81]:
data['M31_id']
Out[81]:
In [83]:
dat = get_data(halos, host_flag, MW_M31_larger, M31_M31_larger, M33_M31_larger)
In [96]:
t = np.copy(dat[0:3])
u = np.copy(dat[0:3])
t
Out[96]:
In [97]:
np.random.shuffle(u)
u
Out[97]:
In [199]:
M33_inds.index(977)
Out[199]:
In [202]:
hosts[4568]
Out[202]:
In [207]:
dat[0]['MWid']
Out[207]:
In [118]:
-1*np.ones(11)
Out[118]:
In [125]:
hosts[0]
Out[125]:
In [136]:
hosts.dtype
Out[136]:
In [141]:
hosts.dtype.names
Out[141]:
In [186]:
fn = ['MW'+nm for nm in hosts.dtype.names]
fn1 = ['M31'+nm for nm in hosts.dtype.names]
fn2 = ['M33'+nm for nm in hosts.dtype.names]
In [188]:
np.hstack((np.array(fn), np.array(fn1), np.array(fn2)))
Out[188]:
In [185]:
np.array(fn).flatten()
Out[185]:
In [162]:
t = np.copy(hosts[0])
In [169]:
for f in t.dtype.names: t[f]=-1
In [171]:
t
Out[171]:
In [208]:
def sq(x):
return x**2
In [209]:
sq(3.0/2)
Out[209]:
In [6]:
a = np.array([1,2,3,4,5,6])
In [8]:
a
Out[8]:
In [15]:
halos[host_flag][0]['id']
Out[15]:
In [ ]: