In [6]:
%pylab inline
from struct import *
In [87]:
INPUT_RESOLUTION=512
WEB_RESOLUTION=50
SMOOTH_SCALE=1.00
BOX_SIZE=100
PATH_TWEB="/home/forero/work/data/WebFinderComparison/Box%d/%d/Tweb/%d/"%(BOX_SIZE,INPUT_RESOLUTION,WEB_RESOLUTION)
PATH_FOF="/home/forero/work/data/WebFinderComparison/Box%d/%d/"%(BOX_SIZE,INPUT_RESOLUTION)
SNAP_NAME="snap%dMpc%d_z0"%(BOX_SIZE, INPUT_RESOLUTION)
FOF_NAME="fof_%dMpc_%d_z0_b0p2"%(BOX_SIZE, INPUT_RESOLUTION)
ENV_CUBE_NAME="CUBE_%d_%d_%s.txt"%(BOX_SIZE, INPUT_RESOLUTION, "ForeroRomero")
ENV_FOF_NAME="FOF_%d_%d_%s.txt"%(BOX_SIZE, INPUT_RESOLUTION, "ForeroRomero")
FILE_ENV_CUBE="%s%s"%(PATH_TWEB, ENV_CUBE_NAME)
FILE_ENV_FOF="%s%s"%(PATH_TWEB, ENV_FOF_NAME)
FILE_ENV_FOF_COLUMN="%s%s.column"%(PATH_TWEB, ENV_FOF_NAME)
FILE_FOF="%s%s"%(PATH_FOF, FOF_NAME)
FILE_DENSITY_TWEB="%s%s.CIC"%(PATH_TWEB, SNAP_NAME)
FILE_DELTA_TWEB_SMOOTH="%s%s.s%.2f.DELTA"%(PATH_TWEB, SNAP_NAME, SMOOTH_SCALE)
FILE_EIGENV1_TWEB="%s%s.eigen_1"%(PATH_TWEB,SNAP_NAME)
FILE_EIGENV2_TWEB="%s%s.eigen_2"%(PATH_TWEB,SNAP_NAME)
FILE_EIGENV3_TWEB="%s%s.eigen_3"%(PATH_TWEB,SNAP_NAME)
FILE_EIGENVEC1_TWEB="%s%s.eigenvec_1"%(PATH_TWEB, SNAP_NAME)
FILE_EIGENVEC2_TWEB="%s%s.eigenvec_2"%(PATH_TWEB, SNAP_NAME)
FILE_EIGENVEC3_TWEB="%s%s.eigenvec_3"%(PATH_TWEB, SNAP_NAME)
FILE_EIGENV1_TWEB_SMOOTH="%s%s.s%.2f.eigen_1"%(PATH_TWEB,SNAP_NAME,SMOOTH_SCALE)
FILE_EIGENV2_TWEB_SMOOTH="%s%s.s%.2f.eigen_2"%(PATH_TWEB,SNAP_NAME,SMOOTH_SCALE)
FILE_EIGENV3_TWEB_SMOOTH="%s%s.s%.2f.eigen_3"%(PATH_TWEB,SNAP_NAME,SMOOTH_SCALE)
FILE_EIGENVEC1_TWEB_SMOOTH="%s%s.s%.2f.eigenvec_1"%(PATH_TWEB, SNAP_NAME, SMOOTH_SCALE)
FILE_EIGENVEC2_TWEB_SMOOTH="%s%s.s%.2f.eigenvec_2"%(PATH_TWEB, SNAP_NAME, SMOOTH_SCALE)
FILE_EIGENVEC3_TWEB_SMOOTH="%s%s.s%.2f.eigenvec_3"%(PATH_TWEB, SNAP_NAME, SMOOTH_SCALE)
In [88]:
FOF_ID=0
FOF_PX=1
FOF_PY=2
FOF_PZ=3
FOF_VX=4
FOF_VY=5
FOF_VZ=6
FOF_NP=7
FOF_MASS=8
ENV_KNOT=3
ENV_FILAMENT=2
ENV_SHEET=1
ENV_VOID=0
In [89]:
def read_CIC_scalar(filename):
f = open(filename, "rb")
dumb = f.read(38)
dumb = f.read(4)
n_x = f.read(4)
n_y = f.read(4)
n_z = f.read(4)
nodes = f.read(4)
x0 = f.read(4)
y0 = f.read(4)
z0 = f.read(4)
dx = f.read(4)
dy = f.read(4)
dz = f.read(4)
dumb = f.read(4)
n_x = (unpack('i', n_x))[0]
n_y = (unpack('i', n_y))[0]
n_z = (unpack('i', n_z))[0]
nodes = (unpack('i', nodes))[0]
dx = (unpack('f', dx))[0]
dy = (unpack('f', dy))[0]
dz = (unpack('f', dz))[0]
x0 = (unpack('f', x0))[0]
y0 = (unpack('f', y0))[0]
z0 = (unpack('f', z0))[0]
print n_x, n_y, n_z, nodes, dx, dy, dz
total_nodes = n_x * n_y *n_z
dumb = f.read(4)
array_data = f.read(total_nodes*4)
dumb = f.read(4)
format_s = str(total_nodes)+'f'
array_data = unpack(format_s, array_data)
f.close()
array_data = np.array(array_data)
array_data.resize(n_z,n_y,n_x)
array_data = array_data.transpose()
return array_data
In [90]:
def get_environment(lambda_1, lambda_2, lambda_3,lambda_th=0.3):
env = int_(lambda_1.copy())
index = where(lambda_3>lambda_th)
env[index] = ENV_KNOT
index = where((lambda_2>lambda_th) & (lambda_3<lambda_th))
env[index] = ENV_FILAMENT
index = where((lambda_1>lambda_th) & (lambda_2<lambda_th) & (lambda_3<lambda_th))
env[index] = ENV_SHEET
index = where(lambda_1<lambda_th)
env[index] = ENV_VOID
return env
In [91]:
lambda1=read_CIC_scalar(FILE_EIGENV1_TWEB_SMOOTH)
lambda2=read_CIC_scalar(FILE_EIGENV2_TWEB_SMOOTH)
lambda3=read_CIC_scalar(FILE_EIGENV3_TWEB_SMOOTH)
fof_data=loadtxt(FILE_FOF, skiprows=3)
grid_env = get_environment(lambda1, lambda2, lambda3)
In [92]:
trace=lambda1 + lambda2 +lambda3
pos_x = fof_data[:,FOF_PX]
pos_y = fof_data[:,FOF_PY]
pos_z = fof_data[:,FOF_PZ]
mass = fof_data[:,FOF_MASS]
In [93]:
n_x = size(lambda1[0,0,:])
boxsize = BOX_SIZE * 1.0
delta_x = boxsize/n_x
slice_point = 0.7*boxsize
n_slice = int_(slice_point/delta_x)
print delta_x, n_x
In [94]:
min_z = n_slice*delta_x
max_z = (n_slice+1.0)*delta_x
slice_id = where((pos_z>min_z) & (pos_z<max_z))
slice_pos_x = pos_x[slice_id]
slice_pos_y = pos_y[slice_id]
slice_pos_z = pos_z[slice_id]
slice_mass = mass[slice_id]
In [95]:
imshow(log10(lambda1[::-1,:,n_slice]+1.0), cmap="bone", extent=[0,boxsize,0,boxsize])
scatter(slice_pos_y,slice_pos_x, color='red',s=5)
Out[95]:
In [96]:
n_halos = size(pos_x)
fof_trace = ones(n_halos)
fof_env = ones(n_halos)
fof_lambda1 = ones(n_halos)
fof_lambda2 = ones(n_halos)
fof_lambda3 = ones(n_halos)
for i in range(n_halos):
ix = pos_x[i]/delta_x
iy = pos_y[i]/delta_x
iz = pos_z[i]/delta_x
fof_trace[i] = trace[ix,iy,iz]
fof_lambda1[i] = lambda1[ix,iy,iz]
fof_lambda2[i] = lambda2[ix,iy,iz]
fof_lambda3[i] = lambda3[ix,iy,iz]
fof_env[i] = grid_env[ix,iy,iz]
#fof_env = get_environment(fof_lambda1, fof_lambda2, fof_lambda3,lambda_th=0.3)
In [97]:
scatter(log10(mass), log10(fof_trace+1.1),s=0.1)
Out[97]:
In [98]:
index = where(fof_env==ENV_KNOT)
scatter(log10(mass[index]), log10(fof_trace[index]+1.0),s=0.1)
Out[98]:
In [99]:
fof_out = open(FILE_ENV_FOF_COLUMN, "w")
for i in range(2):
fof_out.write("#\n")
fof_out.write("WebID\n")
for i in range(n_halos):
fof_out.write("%d\n"%fof_env[i])
fof_out.close()
print FILE_ENV_FOF
!paste $FILE_FOF $FILE_ENV_FOF_COLUMN > $FILE_ENV_FOF
!cp $FILE_ENV_FOF ../data/
In [100]:
cube_out = open(FILE_ENV_CUBE, "w")
cube_out.write("# x y z WebID\n")
for i in range(n_x):
for j in range(n_x):
for k in range(n_x):
cube_out.write("%d %d %d %d\n"%((i + 0.5)*delta_x,(j+0.5)*delta_x,(k+0.5)*delta_x,grid_env[i,j,k]))
cube_out.close()
!cp $FILE_ENV_CUBE ../data/
In [101]:
grid_env = int_(grid_env)
print size(where(grid_env==ENV_KNOT)[0])*1.0/size(grid_env)
print size(where(grid_env==ENV_FILAMENT)[0])*1.0/(1.0*size(grid_env))
print size(where(grid_env==ENV_SHEET)[0])*1.0/size(grid_env)
print size(where(grid_env==ENV_VOID)[0])*1.0/(1.0*size(grid_env))
In [102]:
print FILE_ENV_CUBE
In [ ]: