Skip to content

Commit

Permalink
Merge pull request #88 from deepskies/dev
Browse files Browse the repository at this point in the history
functionality to specify which cls to output + other (param/option name) updates for clarity
  • Loading branch information
samueldmcdermott committed Nov 30, 2023
2 parents b819b82 + abc456d commit 60f1bb8
Show file tree
Hide file tree
Showing 3 changed files with 64 additions and 38 deletions.
66 changes: 44 additions & 22 deletions deepcmbsim/camb_power_spectrum.py
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ def get_noise(self):
provides noise for the TT power spectrum and the polarization power spectra;
shape is (2, max_l_use)
"""
if self.UserParams['noise_type'] == 'white':
if self.UserParams['noise_type'] == 'detector-white':
t_noise = noise.detector_white_noise(self.UserParams['noise_uKarcmin'], self.UserParams['beamfwhm_arcmin'], self.max_l_use,
TT=True)
eb_noise = noise.detector_white_noise(self.UserParams['noise_uKarcmin'], self.UserParams['beamfwhm_arcmin'], self.max_l_use,
Expand All @@ -72,7 +72,7 @@ def get_noise(self):
elif self.UserParams['noise_type'] is None:
return np.zeros((2, self.max_l_use))
else:
print("only detector white noise is currently implemented, via `noise_type = 'white'` in `user_config.yaml`")
print("only detector white noise is currently implemented, via `noise_type = 'detector-white'` in `user_config.yaml`")
return np.zeros((2, self.max_l_use))

def get_cls(self, save_to_dict=None, user_params=True):
Expand All @@ -98,29 +98,51 @@ def get_cls(self, save_to_dict=None, user_params=True):
# main calculation: https://camb.readthedocs.io/en/latest/camb.html#camb.get_results
results = camb.get_results(self.CAMBparams)

outdict = { 'l': range(self.max_l_use + 1) }
if self.UserParams['cls_to_output'] == 'all':
cls_needed = ['clTT', 'clEE', 'clBB', 'clTE', 'clPP', 'clPT', 'clPE']
else:
cls_needed = self.UserParams['cls_to_output']

# https://camb.readthedocs.io/en/latest/results.html#camb.results.CAMBdata.get_total_cls
tt, ee, bb, te = results.get_total_cls(raw_cl=self.normalize_cls, CMB_unit=self.TT_units)[:self.max_l_use + 1].T
if self.UserParams['noise_type'] is not None:
_noise = self.get_noise()
tt += _noise[0]
ee += _noise[1]
bb += _noise[1]
te += _noise[1]
if ('clTT' in cls_needed) or ('clEE' in cls_needed) or ('clEB' in cls_needed) or ('clTE' in cls_needed):
# need to run things to get one/some/all of tt, ee, bb, te
tt, ee, bb, te = results.get_total_cls(raw_cl=self.normalize_cls, CMB_unit=self.TT_units)[:self.max_l_use + 1].T
# add noise
if self.UserParams['noise_type'] is not None:
_noise = self.get_noise()
tt += _noise[0]
ee += _noise[1]
bb += _noise[1]
te += _noise[1]
# now add to outdict
for key in ['clTT', 'clEE', 'clBB', 'clTE']:
if key in cls_needed:
if key == 'clTT':
outdict[key] = tt
elif key == 'clEE':
outdict[key] = ee
elif key == 'clBB':
outdict[key] = bb
elif key == 'clTE':
outdict[key] = te
else:
raise ValueError('somethings wrong.')

# https://camb.readthedocs.io/en/latest/results.html#camb.results.CAMBdata.get_lens_potential_cls
pp, pt, pe = results.get_lens_potential_cls(raw_cl=self.normalize_cls)[:self.max_l_use + 1].T
lvals = range(self.max_l_use + 1)

outdict = {
'l': lvals,
'clTT': tt,
'clEE': ee,
'clBB': bb,
'clTE': te,
'clPP': pp,
'clPT': pt,
'clPE': pe
}
if ('clPP' in cls_needed) or ('clPT' in cls_needed) or ('clPE' in cls_needed):
pp, pt, pe = results.get_lens_potential_cls(raw_cl=self.normalize_cls)[:self.max_l_use + 1].T
# now add to outdict
for key in ['clPP', 'clPT', 'clPE']:
if key in cls_needed:
if key == 'clPP':
outdict[key] = pp
elif key == 'clPT':
outdict[key] = pt
elif key == 'clPE':
outdict[key] = pe
else:
raise ValueError('somethings wrong.')

if bool(self.UserParams["verbose"]):
time_end = dt.now()
Expand Down
26 changes: 13 additions & 13 deletions deepcmbsim/cl_plotting.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,8 @@ class flatmap:
number of pixels on each side of a flat map
degrees : float
number of degrees on each side of a flat map
seed : int, default -1
random seed to provide `namaster` (potential bug in `namaster`, does not always have the
namaster_seed : int, default -1
random namaster_seed to provide `namaster` (potential bug in `namaster`, does not always have the
behavior that is described at https://namaster.readthedocs.io/en/latest/pymaster.html
cl_dict : dict, optional
dictionary of power spectra. Necessary to use the primary method of this class.
Expand All @@ -40,7 +40,7 @@ class flatmap:
function for plotting a realization of a flat-sky map with the given number of pixels
representing the given size of sky patch based on the provided power spectra
"""
def __init__(self, pixels, degrees, seed = -1, cl_dict = None):
def __init__(self, pixels, degrees, namaster_seed = -1, cl_dict = None):
"""
initialize a map-making class
Parameters
Expand All @@ -49,32 +49,32 @@ def __init__(self, pixels, degrees, seed = -1, cl_dict = None):
number of pixels on each side of a flat map
degrees : float
number of degrees on each side of a flat map
seed : int, default -1
random seed to provide `namaster` (potential bug in `namaster`, does not always have the
namaster_seed : int, default -1
random namaster_seed to provide `namaster` (potential bug in `namaster`, does not always have the
behavior that is described at https://namaster.readthedocs.io/en/latest/pymaster.html
cl_dict : dict, optional
dictionary of power spectra. Necessary to use the primary method of this class.
If not specified, can still use the hidden mapping function.
"""
self.pixels, self.degrees = pixels, degrees # pixels on each side, degrees on each side
self.seed, self.cl_dict = seed, cl_dict # random seed if applicable, dictionary of CLs if you want to include them
self.namaster_seed, self.cl_dict = namaster_seed, cl_dict # random namaster_seed if applicable, dictionary of CLs if you want to include them
self.reso = self.degrees / self.pixels
self.nx, self.lx, self.lx_rad, self.ny, self.ly, self.ly_rad = [int(self.pixels), int(self.degrees), self.degrees*np.pi/180]*2
self.ticks, self.labels = np.linspace(0, self.pixels, self.degrees+1), np.arange(0, self.degrees+1, dtype=float)

def _flatmap(self, cl_array, spin_array = [0], seed = None):
seed = seed if seed is not None else self.seed
return nmt.synfast_flat(self.nx, self.ny, self.lx_rad, self.ly_rad, np.array([x for x in cl_array]), spin_array, seed = seed)
def _flatmap(self, cl_array, spin_array = [0], namaster_seed = None):
namaster_seed = namaster_seed if namaster_seed is not None else self.namaster_seed
return nmt.synfast_flat(self.nx, self.ny, self.lx_rad, self.ly_rad, np.array([x for x in cl_array]), spin_array, namaster_seed = namaster_seed)

def flatmap(self, maps_out, seed = None):
def flatmap(self, maps_out, namaster_seed = None):
"""
Parameters
----------
maps_out : ["T", "E", "B", "P", "TT", "EE", "BB", "TE", "PP", "PT", "PE", "clTT", "clEE", "clBB", "clTE", "clPP", "clPT", "clPE", "TEB", "TQU"]
map(s) that you would like to produce
seed : int, optional
random seed to provide `namaster` (potential bug in `namaster`, does not always have the
namaster_seed : int, optional
random namaster_seed to provide `namaster` (potential bug in `namaster`, does not always have the
behavior that is described at https://namaster.readthedocs.io/en/latest/pymaster.html
Returns
-------
Expand All @@ -98,6 +98,6 @@ def flatmap(self, maps_out, seed = None):
else:
print("not a valid map specification")
return None
return self._flatmap(cl_arr, spin_arr, seed=seed)
return self._flatmap(cl_arr, spin_arr, namaster_seed=namaster_seed)
else:
print("if you don't want to restrict to a `cl_dict` dictionary, use `self._flatmap` instead")
10 changes: 7 additions & 3 deletions deepcmbsim/settings/user_config.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,18 @@ FORCAMB: # these simply overwrite any values in BASECAMBPARAMS
ITERABLES: # these provide iterables that can overwrite their corresponding values in BASECAMBPARAMS in loop settings
InitPower.r : [0.01, 0.1, 3] # for nested CAMBparams attributes, use the dot structure
Alens : [.8, 1.2, 3]
seed: 0
namaster_seed: 0 # seed for the map realization using namaster
verbose: 1
normalize_cls: False #raw_cl – return Cl rather than l*(l+1)*Cl/2π (Cl alone is not conventional)
TT_units: muK #return uK**2 units for the TT, EE, BB, TE Cls
outfile_dir: "outfiles"
noise_type: "white"
noise_type: "detector-white" # only option for now; adds detector white noise to the output; see noise.detector_white_noise for more details.
noise_uKarcmin: 5 # noise level in uK*arcmin
beamfwhm_arcmin: 3 # size of beam in arcmin
extra_l: 300
max_l_use: 10000 # max_l_use will differ from max_l and max_l_tensor by "extra_l" because
# according to the CAMB documentation errors affect the last "100 or so" multipoles
# according to the CAMB documentation errors affect the last "100 or so" multipoles
# ---
# decide on what to output; either 'all' or a subset list of
# ['clTT', 'clEE', 'clBB', 'clTE', 'clPP', 'clPT', 'clPE']
cls_to_output: 'all'

0 comments on commit 60f1bb8

Please sign in to comment.