AHF Example
We mentioned that identifying halos and galaxies for a simulation box is a challenging task. Most state-of-the-art halo finders have no direct Python wrappers, yet pynbody allows to overcome this limitation by wrapping to most major halo finders, including the Amiga Halo Finder (AHF), Rockstar and SubFind.
For this reason we demonstrate in this example how one can use the s.halos() command of pynbody to identify halos in a simulation box. For pynbody’s AHF-invoking s.halos() command to work with a Gadget output, one needs to install AHF and place the executable into ~/bin (or extend the $PATH variable)
$ wget popia.ft.uam.es/AHF/files/ahf-v1.0-111.tgz
$ tar zxf ahf-v1.0-111.tgz; cd ahf-v1.0-111 # extract
$ make AHF # build
$ cp bin/./AHF-v1.0-111 ~/bin/ # copy into $PATH
$ echo 'PATH=$PATH:~/bin' >> ~/.bashrc # make $PATH aware of AHF
Then, modify _run_ahf(self, sim) in /path/to/pynbody/halo/ahf.py to:
def _run_ahf(self, sim):
typecode = '60' # '61'
import pynbody.units as units
# find AHFstep
groupfinder = config_parser.get('AHFCatalogue', 'Path')
if groupfinder == 'None':
for directory in os.environ["PATH"].split(os.pathsep):
ahfs = glob.glob(os.path.join(directory, "AHF*"))
for iahf, ahf in enumerate(ahfs):
# if there are more AHF*'s than 1, it's not the last one, and
# it's AHFstep, then continue, otherwise it's OK.
if ((len(ahfs) > 1) & (iahf != len(ahfs) - 1) &
(os.path.basename(ahf) == 'AHFstep')):
continue
else:
groupfinder = ahf
break
if not os.path.exists(groupfinder):
raise RuntimeError("Path to AHF (%s) is invalid" % groupfinder)
if (os.path.basename(groupfinder) == 'AHFstep'):
isAHFstep = True
else:
isAHFstep = False
# build units file
if isAHFstep:
f = open('tipsy.info', 'w')
f.write(str(sim.properties['omegaM0']) + "\n")
f.write(str(sim.properties['omegaL0']) + "\n")
f.write(str(sim['pos'].units.ratio(
units.kpc, a=1) / 1000.0 * sim.properties['h']) + "\n")
f.write(
str(sim['vel'].units.ratio(units.km / units.s, a=1)) + "\n")
f.write(str(sim['mass'].units.ratio(units.Msol)) + "\n")
f.close()
# make input file
f = open('AHF.in', 'w')
f.write(sim._filename + " " + str(typecode) + " 1\n")
f.write(sim._filename + "\n256\n5\n5\n0\n0\n0\n0\n")
f.close()
else:
# make input file
f = open('AHF.in', 'w')
lgmax = np.min([int(2 ** np.floor(np.log2(
1.0 / np.min(sim['eps'])))), 131072])
#lgmax = np.min([int(2 ** np.floor(np.log2(
# 1.0 / 0.19))), 131072])
# hardcoded maximum 131072 might not be necessary
print(config_parser.get('AHFCatalogue', 'Config', vars={
'filename': str(sim._filename),
'typecode': int(typecode),
'gridmax': int(lgmax)
}), file=f)
print(config_parser.get('AHFCatalogue', 'ConfigGadget', vars={
'omega0': sim.properties['omegaM0'],
'lambda0': sim.properties['omegaL0'],
'boxsize': sim['pos'].units.ratio('Mpc a h^-1', **sim.conversion_context()),
'vunit': sim['vel'].units.ratio('km s^-1 a', **sim.conversion_context()),
'munit': sim['mass'].units.ratio('Msol h^-1', **sim.conversion_context()),
'eunit': 0.03 # surely this can't be right?
}), file=f)
f.close()
if (not os.path.exists(sim._filename)):
os.system("gunzip " + sim._filename + ".gz")
# determine parallel possibilities
if os.path.exists(groupfinder):
# run it
os.system(groupfinder + " AHF.in")
return
and the [AHFCatalogue] section in /path/to/pynbody/config.ini to:
[AHFCatalogue]
# settings for the AHF Catalogue reader
AutoRun: True
# automatically attempt to run AHF if no catalogue can be found
# on disk
Path: None
# /path/to/AHF, or None to attempt to find it in your $PATH
AutoGrp: False
# set to true to automatically create a 'grp' array on load
# The grp array
AutoPid: False
# set to true to automatically create a 'pid' array on load
# the PID array is another way to get the particle IDs in the ancestor snapshot,
# but the framework provides h[n].get_index_list(f) for halo catalogue h and
# base snapshot f, so you probably don't need AutoPid
Config: [AHF]
ic_filename = %(filename)s
ic_filetype = %(typecode)s
outfile_prefix = %(filename)s
LgridDomain = 128
LgridMax = %(gridmax)s
NperDomCell = 5
NperRefCell = 5
VescTune = 1.5
NminPerHalo = 50
RhoVir = 0
Dvir = 200
MaxGatherRad = 10.0
ConfigGadget: [GADGET]
GADGET_MUNIT = 1.0e10
GADGET_LUNIT = 1.0e-3
In this example, we generate a mock universe using nbodykit, save the universe to a Gadget 2 file, load the Gadget 2 file with pynbody, identify halos with AHF and estimate shape profiles with CosmicProfiles.
If pynbody.plot.image(halos[2].d, width = '500 kpc', cmap=plt.cm.Greys, units = 'Msol kpc^-2') fails, modify the argument cen_size in the center() function of /path/to/pynbody/analysis/halo.py to something like cen_size="10 kpc".
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Mar 16 18:23:43 2021
"""
from mpi4py import MPI
import matplotlib.pyplot as plt
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()
import pynbody
import os
import sys
import h5py
import inspect
import subprocess
currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
sys.path.append(os.path.join(currentdir, '..')) # Only needed if cosmic_profiles is not installed
import numpy as np
from cosmic_profiles import DensShapeProfs, getMassDMParticle, updateInUnitSystem, updateOutUnitSystem
from nbodykit.lab import cosmology, LogNormalCatalog
def createLogNormUni(L_BOX, nbar, redshift, Nmesh, UNIT_MASS):
""" Create mock simulation box by Poisson-sampling a lognormal density distribution
The Poisson-sampled distribution is evolved according to the Zeldovich (1LPT) prescription
up until redshift ``redshift`` under the constraint of an 'EisensteinHu' power spectrum
:param L_BOX: size of to-be-obtained simulation box
:type L_BOX: float
:param nbar: number density of points (i.e. sampling density / resolution) in box, units: 1/(Mpc/h)**3
Note: ``nbar`` is assumed to be constant across the box
:type nbar: float
:param redshift: redshift of interest
:type redshift: float
:param Nmesh: the mesh size to use when generating the density and displacement fields,
which are Poisson-sampled to particles
:type Nmesh: int
:param UNIT_MASS: in units of solar masses / h. Returned masses will have units UNIT_MASS*(solar_mass)/h
:type UNIT_MASS: float
:return: total number of particles, xyz-coordinates of DM particles, xyz-values of DM particle velocities,
masses of the DM particles (all identical)
:rtype: int, (N,) floats, (N,) floats, (N,) floats, (N,) floats, (N,) floats, (N,) floats, (N,) floats"""
print('Starting createLogNormUni()')
if rank == 0:
# Generating LogNormal Catalog
redshift = redshift
cosmo = cosmology.Planck15
Plin = cosmology.LinearPower(cosmo, redshift, transfer='EisensteinHu')
cat = LogNormalCatalog(Plin=Plin, nbar=nbar, BoxSize=L_BOX, Nmesh=Nmesh, bias=2.0, seed=42)
x_vec = np.float32(np.array(cat['Position'][:,0])) # Mpc/h
y_vec = np.float32(np.array(cat['Position'][:,1]))
z_vec = np.float32(np.array(cat['Position'][:,2]))
x_vel = np.float32(np.array(cat['Velocity'][:,0]))
y_vel = np.float32(np.array(cat['Velocity'][:,1]))
z_vel = np.float32(np.array(cat['Velocity'][:,2]))
N = int(round(len(x_vec)**(1/3)))
N_tot = len(x_vec)
dm_mass = getMassDMParticle(N, L_BOX)/UNIT_MASS
return N_tot, x_vec, y_vec, z_vec, x_vel, y_vel, z_vel, np.ones((len(x_vec),),dtype = np.float32)*dm_mass
else:
return None, None, None, None, None, None, None, None
############## Parameters ####################################################################################
updateInUnitSystem(length_in_cm = 'Mpc/h', mass_in_g = 'Msun/h', velocity_in_cm_per_s = 1e5, little_h = 0.6774)
updateOutUnitSystem(length_in_cm = 'kpc/h', mass_in_g = 'Msun/h', velocity_in_cm_per_s = 1e5, little_h = 0.6774)
L_BOX = np.float32(10) # Mpc/h
nbar = 8e+3 # If too small, e.g. 5e+3: pynbody later yields OSError: Corrupt header record. If too large, need many GBs of RAM.
Nmesh = 256
redshift = 5.5
h = 0.6774
UNIT_MASS = 10**10 # in M_sun
SNAP = '025'
MIN_NUMBER_PTCS = 200
CENTER = 'mode'
CAT_DEST = "./cat"
VIZ_DEST = "./viz"
katz_config = {
'ROverR200': np.logspace(-1.5,0,70),
'IT_TOL': 1e-2,
'IT_WALL': 100,
'IT_MIN': 10,
'REDUCED': False,
'SHELL_BASED': False
}
# Create CAT_DEST and VIZ_DEST if not available
subprocess.call(['mkdir', '-p', '{}'.format(CAT_DEST)], cwd=os.path.join(currentdir))
subprocess.call(['mkdir', '-p', '{}'.format(VIZ_DEST)], cwd=os.path.join(currentdir))
#############################################################################################################
def AHFEx():
############## Generate mock universe #######################################################################
N_tot, dm_x, dm_y, dm_z, vel_x, vel_y, vel_z, mass_array = createLogNormUni(L_BOX, nbar, redshift, Nmesh, UNIT_MASS) # mass_array in UNIT_MASS * solar masses * h^-1
#############################################################################################################
############## Save mock universe to Gadget 2 file #########################################################
f = h5py.File('{}/snap_{}.hdf5'.format(CAT_DEST, SNAP), 'w')
config = f.create_group('Config')
header = f.create_group('Header')
param = f.create_group('Parameters')
ptype1 = f.create_group('PartType1')
# Populate Config
config.attrs['NTYPES'] = 6.0
# Populate Header
header.attrs['MassTable'] = np.array([0, 1, 0, 0, 0, 0])
header.attrs['Redshift'] = redshift
header.attrs['BoxSize'] = L_BOX # in Mpc/h
header.attrs['Omega0'] = 0.3089 # Planck 2015
header.attrs['OmegaBaryon'] = 0.048
header.attrs['OmegaLambda'] = 0.6911
header.attrs['UnitLength_in_cm'] = 3.085678e+24 # 1 Mpc
header.attrs['UnitVelocity_in_cm_per_s'] = 100000.0 # 1 km/s
header.attrs['UnitMass_in_g'] = 1.989e+43 # 10**10 M_sun
header.attrs['HubbleParam'] = h
header.attrs['Time'] = 1/(redshift+1)
header.attrs['NumFilesPerSnapshot'] = 1
header.attrs['NumPart_ThisFile'] = np.array([0, N_tot, 0, 0, 0, 0])
header.attrs['NumPart_Total'] = np.array([0, N_tot, 0, 0, 0, 0])
header.attrs['MassTable'] = np.array([0, mass_array[0], 0, 0, 0, 0])
# Populate Parameters
param.attrs['MassTable'] = np.array([0, mass_array[0], 0, 0, 0, 0])
param.attrs['Redshift'] = redshift
param.attrs['BoxSize'] = L_BOX # in Mpc/h
param.attrs['Omega0'] = 0.3089 # Planck 2015
param.attrs['OmegaBaryon'] = 0.048
param.attrs['OmegaLambda'] = 0.6911
param.attrs['UnitLength_in_cm'] = 3.085678e+24 # 1 Mpc
param.attrs['UnitVelocity_in_cm_per_s'] = 100000.0 # 1 km/s
param.attrs['UnitMass_in_g'] = 1.989e+43 # 10**10 M_sun
param.attrs['HubbleParam'] = h
param.attrs['NumFilesPerSnapshot'] = 1
param.attrs['SofteningComovingType1'] = 0.19/1000 # Mpc
param.attrs['SofteningMaxPhysType1'] = 0.19/1000
param.attrs['SofteningTypeOfPartType1'] = 1
# Populate PartType1
dm_xyz = np.hstack((np.reshape(dm_x, (dm_x.shape[0],1)), np.reshape(dm_y, (dm_y.shape[0],1)), np.reshape(dm_z, (dm_z.shape[0],1))))
vel_xyz = np.hstack((np.reshape(vel_x, (vel_x.shape[0],1)), np.reshape(vel_y, (vel_y.shape[0],1)), np.reshape(vel_z, (vel_z.shape[0],1))))
ptype1_coords = ptype1.create_dataset('Coordinates', dtype = np.float32, data = dm_xyz) # Mpc/h
ptype1_coords.attrs['a_scaling'] = 1.0
ptype1_coords.attrs['h_scaling'] = -1.0
ptype1_ids = ptype1.create_dataset('ParticleIDs', dtype = np.uint32, data = np.arange(N_tot, dtype = np.uint32))
ptype1_ids.attrs['a_scaling'] = 0.0
ptype1_ids.attrs['h_scaling'] = 0.0
ptype1_vels = ptype1.create_dataset('Velocities', dtype = np.float32, data = vel_xyz)
ptype1_vels.attrs['a_scaling'] = 0.5
ptype1_vels.attrs['h_scaling'] = 0.0
f.close()
s = pynbody.load('{}/snap_{}.hdf5'.format(CAT_DEST, SNAP))
s['pos']; s['mass']; s['vel']; s['iord']
s.write(fmt=pynbody.gadget.GadgetSnap,
filename='{}/snap_{}'.format(CAT_DEST, SNAP))
#############################################################################################################
############## Load Gadget 2 file with pynbody ###############################################################
s = pynbody.load('{}/snap_{}'.format(CAT_DEST, SNAP))
print("The snaphost properties read", s.properties)
#############################################################################################################
############## Find halos with e.g. AHF #####################################################################
s['eps'] = 0.01 # The smaller, the larger AHF's LgridMax: SpatialResolution ≃ BoxSize/LgridMax
halos = s.halos() # Install e.g. AHF and place executable in ~/bin (or extend $PATH variable)
# Modify /halo/ahf.py (Gadget, not Tipsy) and config.ini (Gadget, not Tipsy)
#############################################################################################################
############## Viz some halos and retrieve basic properties #################################################
print("The length of halo 1 and 2 are {} and {}, respectively.".format(len(halos[1]), len(halos[2])))
print("The second halo has mass {} 1e12 Msol".format(halos[2]['mass'].sum().in_units('1e12 Msol')))
print("The first 5 particles of halo 2 are located at (in kpc)", halos[2]['pos'].in_units('kpc')[:5])
pynbody.analysis.halo.center(halos[2]) # Modify cen_size="1 kpc" argument in analysis/halo.py's center() method to cen_size="10 kpc" if resolution is too low (i.e. nbar too low)
pynbody.plot.image(halos[2].d, width = '500 kpc', cmap=plt.cm.Greys, units = 'Msol kpc^-2')
plt.savefig('{}/RhoHalo2.pdf'.format(VIZ_DEST))
#############################################################################################################
############## Extract R_vir, halo indices and halo sizes for e.g. first 5 halos ############################
ahf_halos = np.loadtxt('{}/snap_{}.z{}.AHF_halos'.format(CAT_DEST, SNAP, format(redshift, '.3f')), unpack=True)
h_sizes = np.int32(ahf_halos[4])[:5] # 5th column
r_vir = np.float32(ahf_halos[11])[:5] # 12th column
ahf_ptcs = np.loadtxt('{}/snap_{}.z{}.AHF_particles'.format(CAT_DEST, SNAP, format(redshift, '.3f')), skiprows=2, unpack = True, dtype = int)[0]
h_indices = [[] for i in range(len(h_sizes))]
offset = 0
for h_idx in range(len(h_sizes)):
h_indices[h_idx].extend(ahf_ptcs[offset:offset+h_sizes[h_idx]]) # True indices, no plus 1 etc
offset += h_sizes[h_idx] + 1 # Skip line containing number of particles in halo
#############################################################################################################
############## Run cosmic_profiles: define DensShapeProfs object ############################################
cprofiles = DensShapeProfs(dm_xyz, mass_array, h_indices, r_vir, L_BOX, SNAP, VIZ_DEST, CAT_DEST, MIN_NUMBER_PTCS = MIN_NUMBER_PTCS, CENTER = CENTER)
if rank == 0:
h_idx_cat_len = len(cprofiles.getIdxCat()[1])
else:
h_idx_cat_len = None
h_idx_cat_len = comm.bcast(h_idx_cat_len, root = 0)
obj_numbers = np.arange(h_idx_cat_len//2)
############## Create local halo shape catalogue ############################################################
cprofiles.dumpShapeCatLocal(obj_numbers, katz_config)
#############################################################################################################
############## Viz first halo ###############################################################################
cprofiles.vizLocalShapes(obj_numbers, katz_config)
#############################################################################################################
AHFEx()