Skip to content

modified post-processing code to be able to choose axis #107

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
andreufont wants to merge 1 commit into
base: master
Choose a base branch
from
Open

modified post-processing code to be able to choose axis #107

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 35 additions & 27 deletions post_process/py/extract_skewers.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,12 +9,14 @@
import camb_cosmo
import thermal_model

def get_skewers_filename(num,n_skewers,width_Mpc,scale_T0=None,
def get_skewers_filename(num,n_skewers,width_Mpc,axis=None,scale_T0=None,
scale_gamma=None):
"""Filename storing skewers for a particular temperature model"""

filename='skewers_'+str(num)+'_Ns'+str(n_skewers)
filename+='_wM'+str(int(1000*width_Mpc)/1000)
if axis is not None:
filename+='_axis'+str(axis)
if scale_T0:
filename+='_sT'+str(int(1000*scale_T0)/1000)
if scale_gamma:
Expand All @@ -23,11 +25,13 @@ def get_skewers_filename(num,n_skewers,width_Mpc,scale_T0=None,
return filename


def get_snapshot_json_filename(num,n_skewers,width_Mpc):
def get_snapshot_json_filename(num,n_skewers,width_Mpc,axis=None):
"""Filename describing the set of skewers for a given snapshot"""

filename='snap_skewers_'+str(num)+'_Ns'+str(n_skewers)
filename+='_wM'+str(int(1000*width_Mpc)/1000)
if axis is not None:
filename+='_axis'+str(axis)
filename+='.json'
return filename

Expand Down Expand Up @@ -55,8 +59,9 @@ def thermal_broadening_Mpc(T_0,dkms_dMpc):
return sigma_T_Mpc


def rescale_write_skewers_z(simdir,num,skewers_dir=None,n_skewers=50,
width_Mpc=0.1,scales_T0=None,scales_gamma=None):
def rescale_write_skewers_z(simdir,num,skewers_dir=None,
n_skewers=50,width_Mpc=0.1,axis=None,
scales_T0=None,scales_gamma=None):
"""Extract skewers for a given snapshot, for different temperatures."""

# don't rescale unless asked to
Expand All @@ -82,6 +87,8 @@ def rescale_write_skewers_z(simdir,num,skewers_dir=None,n_skewers=50,
'width_Mpc':width_Mpc, 'width_kms':width_kms,
'T0_ini':T0_ini, 'gamma_ini':gamma_ini,
'scales_T0':scales_T0, 'scales_gamma':scales_gamma}
if axis is not None:
sim_info['axis']=axis

# will also store measured values
sim_T0=[]
Expand All @@ -96,19 +103,20 @@ def rescale_write_skewers_z(simdir,num,skewers_dir=None,n_skewers=50,
for scale_gamma in scales_gamma:
T0=T0_ini*scale_T0
gamma=gamma_ini*scale_gamma
sk_filename=get_skewers_filename(num,n_skewers,width_Mpc,
scale_T0,scale_gamma)
sk_filename=get_skewers_filename(num,n_skewers=n_skewers,
width_Mpc=width_Mpc,axis=axis,
scale_T0=scale_T0,scale_gamma=scale_gamma)

# avoid (if possible) to use set_T0, might break fake_spectra
if (scale_T0==1.0) and (scale_gamma==1.0):
skewers=get_skewers_snapshot(simdir,skewers_dir,num,
n_skewers=n_skewers,width_kms=width_kms,
skewers_filename=sk_filename)
axis=axis,skewers_filename=sk_filename)
else:
skewers=get_skewers_snapshot(simdir,skewers_dir,num,
n_skewers=n_skewers,width_kms=width_kms,
set_T0=T0,set_gamma=gamma,
skewers_filename=sk_filename)
axis=axis,skewers_filename=sk_filename)

# call mean flux, so that the skewers are really computed
mf=skewers.get_mean_flux()
Expand All @@ -131,7 +139,7 @@ def rescale_write_skewers_z(simdir,num,skewers_dir=None,n_skewers=50,
sim_info['sim_scale_gamma']=sim_scale_gamma
sim_info['sk_files']=sk_files

snapshot_filename=get_snapshot_json_filename(num,n_skewers,width_Mpc)
snapshot_filename=get_snapshot_json_filename(num,n_skewers,width_Mpc,axis)
sim_info['snapshot_filename']=snapshot_filename
json_file = open(skewers_dir+'/'+snapshot_filename,"w")
json.dump(sim_info,json_file)
Expand All @@ -143,33 +151,33 @@ def rescale_write_skewers_z(simdir,num,skewers_dir=None,n_skewers=50,


def get_skewers_snapshot(simdir,skewers_dir,snap_num,n_skewers=50,width_kms=10,
set_T0=None,set_gamma=None,skewers_filename=None):
axis=None,set_T0=None,set_gamma=None,
skewers_filename='test_skewers.hdf5'):
"""Extract skewers for a particular snapshot"""

if not skewers_filename:
skewers_filename="skewers_"+str(snap_num)+"_"+str(n_skewers)
skewers_filename+="_"+str(width_kms)
if set_T0:
skewers_filename+='_T0_'+str(set_T0)
if set_gamma:
skewers_filename+='_gamma_'+str(set_gamma)
skewers_filename+='.hdf5'

# check that spectra file does not exist yet (should crash)
if os.path.exists(skewers_dir+'/'+skewers_filename):
print(skewers_filename,'already exists in',skewers_dir)

# avoid (if possible) to use set_T0, not always present in fake_spectra
if (set_T0 is None) and (set_gamma is None):
# new functionality to specify simulation axis for skewers
if axis is not None:
# for now only working without temperature rescalings
assert (set_T0 is None) and (set_gamma is None), 'get_skewers_snapshot'
skewers = grid_spec.GriddedSpectra(snap_num,simdir+'/output/',
nspec=n_skewers,res=width_kms,savefile=skewers_filename,
nspec=n_skewers,res=width_kms,axis=axis,
savefile=skewers_filename,
savedir=skewers_dir,reload_file=True)
else:
skewers = grid_spec.GriddedSpectra(snap_num,simdir+'/output/',
nspec=n_skewers,res=width_kms,savefile=skewers_filename,
savedir=skewers_dir,reload_file=True,
set_T0=set_T0,set_gamma=set_gamma)

# avoid (if possible) to use set_T0, not always present in fake_spectra
if (set_T0 is None) and (set_gamma is None):
skewers = grid_spec.GriddedSpectra(snap_num,simdir+'/output/',
nspec=n_skewers,res=width_kms,savefile=skewers_filename,
savedir=skewers_dir,reload_file=True)
else:
skewers = grid_spec.GriddedSpectra(snap_num,simdir+'/output/',
nspec=n_skewers,res=width_kms,savefile=skewers_filename,
savedir=skewers_dir,reload_file=True,
set_T0=set_T0,set_gamma=set_gamma)

return skewers

23 changes: 14 additions & 9 deletions post_process/py/snapshot_admin.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,20 +34,20 @@ def __init__(self,snap_json,scales_tau=None,kF_Mpc=None):
self.scales_tau=[1.0]


def get_all_flux_power(self,simdir=None):
def get_all_flux_power(self,simdir=None,skewers_dir=None):
"""Loop over all skewers, and return flux power for each"""

if simdir is not None:
# get box size from GenIC file, to normalize power
if 'simdir' in self.data:
simdir = self.data['simdir']
elif 'basedir' in self.data:
simdir = self.data['basedir']
if simdir is None:
simdir = self.data['simdir']

if skewers_dir is None:
skewers_dir=self.data['skewers_dir']

# get box size from GenIC file, to normalize power
genic_file=simdir+'/paramfile.genic'
L_Mpc=read_genic.L_Mpc_from_paramfile(genic_file,verbose=True)

skewers_dir=simdir+"output/skewers/"
# get snapshot number from JSON data
snap_num=self.data['snap_num']

# will loop over all temperature models in snapshot
Expand Down Expand Up @@ -96,6 +96,10 @@ def get_p1d_json_filename(self,p1d_label):

filename=p1d_label+'_'+str(num)+'_Ns'+str(n_skewers)
filename+='_wM'+str(int(1000*width_Mpc)/1000)

if 'axis' in self.data:
filename+='_axis'+str(self.data['axis'])

filename+='.json'

return filename
Expand All @@ -107,7 +111,8 @@ def write_p1d_json(self,p1d_label='p1d'):
if p1d_label is None:
p1d_label='p1d'

filename=self.data['simdir']+'/'+self.get_p1d_json_filename(p1d_label)
filename=self.data['skewers_dir']+'/'+self.get_p1d_json_filename(p1d_label)

print('will print P1d to file',filename)

# make sure we have already computed P1D
Expand Down
16 changes: 5 additions & 11 deletions post_process/scripts/arxiv_flux_power.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,9 +14,10 @@
parser = argparse.ArgumentParser()
parser.add_argument('--simdir', type=str, help='Base simulation directory',required=True)
parser.add_argument('--snap_num', type=int, help='Snapshop number',required=True)
parser.add_argument('--skewers_dir', type=str, help='Store skewers in this folder',required=False)
parser.add_argument('--skewers_dir', type=str, help='Store skewers in this folder',required=True)
parser.add_argument('--n_skewers', type=int, default=10, help='Number of skewers per side',required=False)
parser.add_argument('--width_Mpc', type=float, default=0.1, help='Cell width (in Mpc)',required=False)
parser.add_argument('--axis', type=int, help='Axis to use to extract skewers (1,2,3)', required=False)
parser.add_argument('--scales_tau', type=str, default='1.0', help='Comma-separated list of optical depth scalings to use.',required=False)
parser.add_argument('--p1d_label', type=str, default=None, help='String identifying P1D measurement and / or tau scaling.',required=False)
parser.add_argument('--verbose', action='store_true', help='Print runtime information',required=False)
Expand All @@ -25,19 +26,12 @@
verbose=args.verbose
print('verbose =',verbose)

# main simulation folder
simdir=args.simdir
if args.skewers_dir:
skewers_dir=args.skewers_dir
else:
skewers_dir=simdir+'/output/skewers/'

scales_tau=[float(scale) for scale in args.scales_tau.split(',')]
if verbose:
print('will scale tau by',scales_tau)

# try to read information about filtering length in simulation
kF_json=simdir+'/filtering_length.json'
kF_json=args.simdir+'/filtering_length.json'
if os.path.isfile(kF_json):
# read json file with filtering data
with open(kF_json) as json_data:
Expand All @@ -48,9 +42,9 @@
kF_Mpc=None

# read file containing information of all temperature rescalings in snapshot
snap_filename=skewers_dir+'/'+extract_skewers.get_snapshot_json_filename(
snap_filename=args.skewers_dir+'/'+extract_skewers.get_snapshot_json_filename(
num=args.snap_num,n_skewers=args.n_skewers,
width_Mpc=args.width_Mpc)
width_Mpc=args.width_Mpc,axis=args.axis)
if verbose:
print('setup snapshot admin from file',snap_filename)

Expand Down
5 changes: 3 additions & 2 deletions post_process/scripts/extract_default_skewers.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,8 @@
parser = argparse.ArgumentParser()
parser.add_argument('--simdir', type=str, help='Base simulation directory',required=True)
parser.add_argument('--skewers_dir', type=str, help='Store skewers in this folder',required=True)
parser.add_argument('--snap_num', type=int, default=8, help='Snapshot number',required=False)
parser.add_argument('--snap_num', type=int, help='Snapshot number',required=True)
parser.add_argument('--axis', type=int, help='Axis to use to extract skewers (1,2,3)',required=False)
parser.add_argument('--n_skewers', type=int, default=10, help='Number of skewers per side',required=False)
parser.add_argument('--width_Mpc', type=float, default=0.1, help='Cell width (in Mpc)',required=False)
parser.add_argument('--verbose', action='store_true', help='Print runtime information',required=False)
Expand All @@ -24,7 +25,7 @@
# extract skewers for one snapshot, and return some info
info=extract_skewers.rescale_write_skewers_z(simdir=args.simdir,
num=args.snap_num, skewers_dir=args.skewers_dir,
n_skewers=args.n_skewers,width_Mpc=args.width_Mpc)
axis=args.axis,n_skewers=args.n_skewers,width_Mpc=args.width_Mpc)

if args.verbose:
print('print extra info about skewers')
Expand Down