From ef6da39dd08eba4056cb42bb68f1bdb0009cafea Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Wed, 1 Sep 2021 04:45:24 +0200 Subject: [PATCH 01/29] Step 1: transfer Anne's works to new branch --- xfaster/xfaster_class.py | 1 + 1 file changed, 1 insertion(+) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 58a42c53..7ab9e0f6 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -572,6 +572,7 @@ def _get_files( mask_type="hitsmask_tailored", signal_type="r0p03", signal_transfer_type=None, + signal_transfer_type_dust=None, signal_subset="*", noise_type="stationary", noise_subset="*", From fcabcd6369ace4e6676ad5dbd0faa31f2bec71f3 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Fri, 3 Sep 2021 22:24:14 +0200 Subject: [PATCH 02/29] finish adding Anne's diffs. Need to review code for bugs --- xfaster/xfaster_class.py | 279 ++++++++++++++++++++++++++++++++++++--- xfaster/xfaster_exec.py | 65 ++++++++- 2 files changed, 323 insertions(+), 21 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 7ab9e0f6..fad24b68 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -574,6 +574,7 @@ def _get_files( signal_transfer_type=None, signal_transfer_type_dust=None, signal_subset="*", + signal_subset_dust="None", noise_type="stationary", noise_subset="*", signal_type_sim=None, @@ -687,6 +688,19 @@ def find_sim_files(name, root=None, subset="*"): "signal_transfer", "signal_{}".format(signal_transfer_type), signal_subset ) + # transfer_dust + # do option to separate_dust_transfer + #if separate_dust_transfer: + if signal_transfer_type_dust is None: ++ signal_transfer_root_dust = None ++ signal_transfer_files_dust = None + else: + if signal_subset_dust is None: + signal_subset_dust = signal_subset + find_sim_files( + "signal_transfer_dust", "signal_{}".format(signal_transfer_type_dust), signal_subset_dust + ) + # apply suffix out = dict() for k, v in fs.items(): @@ -706,6 +720,8 @@ def get_files( signal_type="synfast", signal_type_sim=None, signal_transfer_type=None, + signal_transfer_type_dust=None, # transfer_dust + signal_subset_dust=None, data_root2=None, data_subset2=None, foreground_type_sim=None, @@ -900,6 +916,8 @@ def get_files( mask_type=mask_type, signal_type=signal_type, signal_transfer_type=signal_transfer_type, + signal_transfer_type_dust=signal_transfer_type_dust, + signal_subset_dust=signal_subset_dust, signal_subset=signal_subset, noise_subset=noise_subset, ) @@ -2043,7 +2061,7 @@ def process_index(idx): gmat[xname] = OrderedDict() for spec in self.specs: si, sj = spec_inds[spec] - f = (fsky_eff[xname][si, sj] + fsky_eff[xname][sj, si]) / 2.0 + f = (ky_eff[xname][si, sj] + fsky_eff[xname][sj, si]) / 2.0 gmat[xname][spec] = f if np.any(np.asarray([f for f in fsky.values()]) > 0.1): @@ -2741,6 +2759,8 @@ def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): signal_files = self.signal_transfer_files signal_files2 = self.signal_transfer_files2 if null_run else None num_signal = self.num_signal_transfer + signal_files_dust = self.signal_transfer_files_dust # how about null_run + num_signal_dust = self.num_signal_transfer_dust # SET UP THESE else: signal_files = self.signal_files signal_files2 = self.signal_files2 if null_run else None @@ -2755,6 +2775,7 @@ def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): save_attrs = [ "cls_signal", + "cls_signal_dust", "cls_noise", "cls_sim", "cls_med", @@ -2772,6 +2793,8 @@ def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): if transfer: save_name = "sims_xcorr_{}".format(self.signal_transfer_type) cp = "sims_transfer" + if self.signal_transfer_type_dust is not None: + save_name += "_" + self.signal_transfer_type_dust else: save_name = "sims_xcorr_{}".format(self.signal_type) cp = "sims" @@ -2784,6 +2807,8 @@ def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): shape=data_shape, shape_ref="cls_signal", ) + # load cls_signal_dust + if ret is not None: if do_noise: if self.cls_noise is not None: @@ -2803,6 +2828,7 @@ def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): # process signal, noise, and S+N cls_sig = OrderedDict() + cls_sig_dust = OrderedDict() cls_null_sig = OrderedDict() if null_run else None cls_noise = OrderedDict() if do_noise else None cls_null_noise = OrderedDict() if null_run and do_noise else None @@ -2815,16 +2841,16 @@ def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): cls_res = OrderedDict() if do_noise else None cls_null_res = OrderedDict() if null_run and do_noise else None if do_noise: - for k in ["nxn0", "nxn1", "sxn0", "sxn1", "nxs0", "nxs1"]: + for k in ["nxn0", "nxn1", "sxn0", "sxn1", "nxs0", "nxs1", "cls_signal_dust"]: cls_res[k] = OrderedDict() if null_run: cls_null_res[k] = OrderedDict() if num_noise != 0: - nsim_min = min([num_signal, num_noise]) + nsim_min = min([num_signal, num_signal_dust, num_noise]) else: nsim_min = num_signal - nsim_max = max([num_signal, num_noise]) + nsim_max = max([num_signal, num_signal_dust, num_noise]) cls_all = np.zeros( [nsim_max, len(map_pairs.items()), len(self.specs), self.lmax + 1] ) @@ -2892,10 +2918,12 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): return cache[idx] sig_cache = dict() + sig_dust_cache = dict() noise_cache = dict() for isim in range(nsim_max): sig_cache.clear() + sig_dust_cache.clear() noise_cache.clear() for xind, (xname, (idx, jdx)) in enumerate(map_pairs.items()): self.log( @@ -2921,6 +2949,15 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): if null_run: cls_null1t = np.copy(cls_null1_sig) + if isim < num_signal_dust: + dimap_alms, _ = process_index( + signal_files_dust, None, idx, isim, sig_dust_cache + ) + djmap_alms, djnull_alms = process_index( + signal_files_dust, None, jdx, isim, sig_dust_cache + ) + cls1_sig_dust = self.alm2cl(dimap_alms, djmap_alms) + if do_noise and isim < num_noise: cls1_res = OrderedDict() if null_run: @@ -2976,6 +3013,8 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): quants += [[cls_sig, cls1_sig]] if null_run: quants += [[cls_null_sig, cls_null1_sig]] + if isim < num_signal_dust: + quants += [[cls_sig_dust, cls1_sig_dust]] if do_noise and isim < num_noise: quants += [[cls_noise, cls1_noise]] @@ -3020,6 +3059,7 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): cls_null_med[spec][xname] = cls_null_med_arr[xind][s] self.cls_signal = cls_sig + self.cls_signal_dust = cls_sig_dust self.cls_signal_null = cls_null_sig self.cls_noise = cls_noise self.cls_noise_null = cls_null_noise @@ -3267,7 +3307,10 @@ def get_signal_shape( flat=None, signal_mask=None, transfer=False, + transfer_dust=False, save=True, + template_alpha90=None, + template_alpha150=None, ): """ Load a shape spectrum for input to the Fisher iteration algorithm. @@ -3285,12 +3328,13 @@ def get_signal_shape( Arguments --------- - filename : string + filename : string, or list of strings Filename for a spectrum on disk. If None, and ``r`` is None and ``flat`` is False, this will search for a spectrum stored in ``signal_/spec_signal_.dat``. Otherwise, if the filename is a relative path and not found, the config directory is searched. + If a list, first element is CMB and second element is dust. r : float If supplied and ``flat`` is False, a spectrum is computed using CAMB for the given ``r`` value. Overrides ``filename``. @@ -3309,6 +3353,10 @@ def get_signal_shape( and ``r`` is None and ``flat`` is False, will search for a spectrum stored in ``signal_/spec_signal_.dat``. + transfer_dust : bool + If True, this is a transfer function run for dust and dust+cmb. + If `filename` is None, will search for a spectrum stored in + `signal_/spec_signal_.dat`. save : bool If True, save signal shape dict to disk. @@ -3318,6 +3366,9 @@ def get_signal_shape( Dictionary keyed by spectrum (cmb_tt, cmb_ee, ... , fg), each entry containing a vector of length 2 * lmax + 1 """ + if isinstance(filename, list): + filename_dust = filename[1] + filename = filename[0] lmax_kern = 2 * self.lmax @@ -3361,6 +3412,27 @@ def get_signal_shape( ellfac = ell * (ell + 1) / 2.0 / np.pi cls_shape = OrderedDict() + def fix_camb_specs(arr, ltmp, filename): + # camb starts at l=2, so set first 2 ells to be 0 + # make sure order is TT, EE, BB, TE, no ell row. + new_arr = np.append([0, 0], arr[1, : ltmp - 2]) + if self.pol: + if np.any(arr[2, : ltmp - 2] < 0): + # this is true if TE is the third index instead of EE + # # (ell is 0th index for CAMB) + self.log( + "Old CAMB format in model file {}. Re-indexing.".format( + filename + ), + "detail", + ) + pol_specs = [3, 4, 2] + else: + pol_specs = [2, 3, 4] + for d in arr[pol_specs]: + new_arr = np.vstack([new_arr, np.append([0, 0], d[: ltmp - 2])]) + return new_arr + if flat is not None and flat is not False: if flat is True: flat = 2e-5 @@ -3403,6 +3475,36 @@ def get_signal_shape( if not os.path.exists(filename): raise OSError("Missing model file {}".format(filename)) + +# chech this + if transfer_dust: + if filename_dust is None: + signal_root = self.signal_transfer_root_dust + filename_dust = "spec_{}.dat".format(os.path.basename(signal_root)) + filename_dust = os.path.join(signal_root, filename_dust) + if not os.path.exists(filename_dust) and not os.path.isabs( + filename_dust + ): + filename_dust = os.path.abspath( + os.path.join( + os.path.dirname(__file__), "../data", filename_dust + ) + ) + if not os.path.exists(filename_dust): + raise OSError("Missing model file {}".format(filename_dust)) + tmpd = np.loadtxt(filename_dust, unpack=True) + + ltmpd = tmpd.shape[-1] + if lmax_kern + 1 < ltmpd: + ltmpd = lmax_kern + 1 + else: + raise ValueError( + "Require at least lmax={} in model file, found {}".format( + lmax_kern, ltmpd + ) + ) + tmpd = fix_camb_specs(tmpd, ltmpd, filename_dust) + tmp = np.loadtxt(filename, unpack=True) ltmp = tmp.shape[-1] @@ -3415,6 +3517,7 @@ def get_signal_shape( ) ) + tmp = fix_camb_specs(tmp, ltmp, filename) # camb starts at l=2, so set first 2 ells to be 0 cls_shape["cmb_tt"] = np.append([0, 0], tmp[1, : ltmp - 2]) if self.pol: @@ -3427,10 +3530,30 @@ def get_signal_shape( ) ) pol_specs = [3, 4, 2] + spec_order = {"tt": 0, "ee": 1, "bb": 2, "te": 3} + for spec in specs: + parts = spec.split("_") + comp = parts[0] + spec0 = parts[1] + self.log(spec) + if len(parts) == 3: + nom_freq = parts[2] + if comp == "cmb": + arr = tmp + elif comp == "dust": + arr = tmpd else: pol_specs = [2, 3, 4] for spec, d in zip(specs[1:4], tmp[pol_specs]): cls_shape[spec] = np.append([0, 0], d[: ltmp - 2]) + # Combine CMB and dust scaled by alpha**2 for frequency + if nom_freq == "90": + alpha = template_alpha90 + elif nom_freq == "150": + alpha = template_alpha150 + arr = tmp + alpha ** 2 * tmpd + cls_shape[spec] = arr[spec_order[spec0]] + # EB and TB flat l^2 * C_l if self.pol: @@ -3452,6 +3575,25 @@ def get_signal_shape( "Added foreground to cls shape {}".format(list(cls_shape)), "debug" ) + if transfer_dust: + cls_shape["dust_tt"] = 0 + #nspecs += 1 + # Need separate fields per freq for dust+CMB to get relative + # amplitudes right + if self.pol: + cls_shape["dust_ee"] = 0 # add values here + cls_shape["dust_bb"] = 0 + cls_shape["dust_te"] = 0 + # nspecs += 3 + for freq in np.unique(list(self.nom_freqs.values())): + cls_shape["cmbdust_tt_{}".format(freq)] = 0 + # nspecs += 1 + if self.pol: + for s in ["ee", "bb", "te"]: + cls_shape["cmbdust_{}_{}".format(s, freq)] = 0 + nspecs += 3 + + # divide out l^2/2pi for spec in specs: cls_shape[spec][2:] /= ellfac[2:] @@ -3881,6 +4023,8 @@ def bin_cl_template( transfer_run=False, beam_error=False, use_precalc=True, + dust=False, + cmbdust=False, ): """ Compute the Cbl matrix from the input shape spectrum. @@ -3964,6 +4108,10 @@ def bin_cl_template( comps += ["cmb"] if "fg" in cls_shape and not transfer_run: comps += ["fg"] + if transfer_run and dust: + comps += ["dust"] + if transfer_run and cmbdust: + comps += ["dustcmb"] if self.nbins_res > 0 and not transfer_run: comps += ["res"] cls_noise = self.cls_noise_null if self.null_run else self.cls_noise @@ -4095,7 +4243,9 @@ def bin_things(comp, d, md): # single foreground spectrum s_arr[si, xi] = cls_shape["fg"][lk] else: - s_arr[si, xi] = cls_shape["cmb_{}".format(spec)][lk] + # for comp in [cmb, dust, cmbdust] + s_arr[si, xi] = cls = cls_shape["{}_{}".format(comp, spec)][lk] + # cls_shape["cmb_{}".format(spec)][lk] # get cross spectrum kernel terms k_arr[si, xi] = mll[spec][xname][ls, lk] @@ -4340,7 +4490,18 @@ def get_model_spectra( return cls - def get_data_spectra(self, map_tag=None, transfer_run=False, do_noise=True): + #def get_data_spectra(self, map_tag=None, transfer_run=False, do_noise=True): + def get_data_spectra( + self, + map_tag=None, + transfer_run=False, + do_noise=True, + transfer_dust=False, + transfer_cmbdust=False, + template_alpha90=None, + template_alpha150=None, + ): + """ Return data and noise spectra for the given map tag(s). Data spectra and signal/noise sim spectra must have been precomputed or loaded from @@ -4384,6 +4545,23 @@ def get_data_spectra(self, map_tag=None, transfer_run=False, do_noise=True): # obs depends on what you're computing if transfer_run: obs_quant = self.cls_signal + if transfer_dust: + obs_quant = self.cls_signal_dust + elif transfer_cmbdust: + # need to add cmb and dust scaled by alphas + obs_quant = copy.deepcopy(self.cls_signal_dust) + adict = {"90": template_alpha90, "150": template_alpha150} + for k2 in self.cls_signal_dust[k1].keys(): + mtag1, mtag2 = k2.split(":") + freq1 = self.nom_freqs[mtag1] + freq2 = self.nom_freqs[mtag2] + obs_quant[k1][k2] = ( + self.cls_signal + + adict[freq1] * adict[freq2] * self.cls_signal_dust + ) + else: + obs_quant = self.cls_signal + elif self.null_run: if self.reference_subtracted: obs_quant = self.cls_data_sub_null @@ -5114,6 +5292,10 @@ def fisher_iterate( like_profile_sigma=3.0, like_profile_points=100, file_tag=None, + transfer_dust=dust, + transfer_cmbdust=cmbdust, + template_alpha90=tempalate_alpha90, + template_alpha150=template_alpha150, ): """ Iterate over the Fisher calculation to compute bandpower estimates @@ -5240,7 +5422,13 @@ def fisher_iterate( qb = qb_start obs, nell, debias = self.get_data_spectra( - map_tag=map_tag, transfer_run=transfer_run + #map_tag=map_tag, transfer_run=transfer_run, + map_tag=map_tag, + transfer_run=transfer_run, + transfer_dust=transfer_dust, + transfer_cmbdust=cmbdust, + template_alpha90=tempalate_alpha90, + template_alpha150=tempalate_alpha150, ) bin_index = pt.dict_to_index(self.bin_def) @@ -5625,6 +5813,10 @@ def get_transfer( iter_max=200, save_iters=False, fix_bb_transfer=False, + dust=False, + cmbdust=False, + template_alpha90=None, + template_alpha150=None, ): """ Compute the transfer function from signal simulations created using @@ -5683,6 +5875,11 @@ def get_transfer( ) save_name = "transfer_all" + if dust: + save_name = "{}_dust".format(save_name) + elif cmbdust: + save_name = "{}_cmbdust".format(save_name) + if self.weighted_bins: save_name = "{}_wbins".format(save_name) @@ -5717,7 +5914,15 @@ def expand_transfer(qb_transfer, bin_def): return transfer if ret is not None: - self.qb_transfer = ret["qb_transfer"] + #self.qb_transfer = ret["qb_transfer"] + if dust: + self.qb_transfer_dust = ret["qb_transfer"] + elif cmbdust: + self.qb_transfer_cmbdust = ret["qb_transfer"] + else: + self.qb_transfer = ret["qb_transfer"] + # return ret["qb_transfer"] + if "transfer" in ret: self.transfer = ret["transfer"] else: @@ -5725,16 +5930,31 @@ def expand_transfer(qb_transfer, bin_def): return self.transfer self.qb_transfer = OrderedDict() - for spec in self.specs: - self.qb_transfer["cmb_" + spec] = OrderedDict() + qb_transfer = OrderedDict() + if dust: + comp = "dust" + elif cmbdust: + comp = "cmbdust" + else: + comp = "cmb" + for spec in self.specs: + #self.qb_transfer["cmb_" + spec] = OrderedDict() + if dust: + qb_transfer["{}_{}".format(comp, spec)] = OrderedDict() + + # success = False success = True msg = "" for im0, m0 in enumerate(self.map_tags): if not self.fit_transfer[self.map_tags_orig[im0]]: - for spec in self.specs: - self.qb_transfer["cmb_{}".format(spec)][m0] = np.ones( + #for spec in self.specs: + # self.qb_transfer["cmb_{}".format(spec)][m0] = np.ones( + # self.nbins_cmb // len(self.specs) + # ) + for stag in self.qb_transfer.keys(): + self.qb_transfer[stag][m0] = np.ones( self.nbins_cmb // len(self.specs) ) self.log("Setting map {} transfer to unity".format(m0), "info") @@ -5747,7 +5967,10 @@ def expand_transfer(qb_transfer, bin_def): "info", ) self.clear_precalc() - cbl = self.bin_cl_template(map_tag=m0, transfer_run=True) + #cbl = self.bin_cl_template(map_tag=m0, transfer_run=True) + cbl = self.bin_cl_template( + cls_shape, m0, transfer_run=True, dust=dust, cmbdust=cmbdust + ) ret = self.fisher_iterate( cbl, m0, @@ -5755,6 +5978,10 @@ def expand_transfer(qb_transfer, bin_def): iter_max=iter_max, converge_criteria=converge_criteria, save_iters=save_iters, + transfer_dust=dust, + transfer_cmbdust=cmbdust, + template_alpha90=tempalate_alpha90, + template_alpha150=template_alpha150, ) qb = ret["qb"] @@ -5788,11 +6015,14 @@ def expand_transfer(qb_transfer, bin_def): # Set EB/TB qb transfers to geometric means of components if len(self.specs) > 4: - qb["cmb_eb"] = np.sqrt(np.abs(qb["cmb_ee"] * qb["cmb_bb"])) - qb["cmb_tb"] = np.sqrt(np.abs(qb["cmb_tt"] * qb["cmb_bb"])) + #qb["cmb_eb"] = np.sqrt(np.abs(qb["cmb_ee"] * qb["cmb_bb"])) + #qb["cmb_tb"] = np.sqrt(np.abs(qb["cmb_tt"] * qb["cmb_bb"])) + qb[comp + "_eb"] = np.sqrt(np.abs(qb[comp + "_ee"] * qb[comp + "_bb"])) + qb[comp + "_tb"] = np.sqrt(np.abs(qb[comp + "_tt"] * qb[comp + "_bb"])) for stag, qbdat in qb.items(): - self.qb_transfer[stag][m0] = qbdat + #self.qb_transfer[stag][m0] = qbdat + qb_transfer[stag][m0] = qbdat if success: self.transfer = expand_transfer(self.qb_transfer, self.bin_def) @@ -5809,8 +6039,16 @@ def expand_transfer(qb_transfer, bin_def): if not success: raise RuntimeError("Error computing transfer function: {}".format(msg)) + + if dust: + self.qb_transfer_dust = qb_transfer + elif cmbdust: + self.qb_transfer_cmbdust = qb_transfer + else: + self.qb_transfer = qb_transfer - return self.transfer + return qb_transfer + #return self.transfer # equivalent to self.qb_transfer in Anne's old code. def get_bandpowers( self, @@ -5876,6 +6114,11 @@ def get_bandpowers( Keep first CMB bandpowers fixed to input shape (qb=1). return_cls : bool If True, return C_ls rather than D_ls + dust : bool + If True, compute transfer function using dust sims + cmbdust : bool + If True, compute transfer function using CMB+dust sims, with + dust scaled by template_alpha90, template_alpha150 like_profiles : bool If True, compute profile likelihoods for each qb, leaving all others fixed at their maximum likelihood values. Profiles are diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index 0f5cfbaf..2b1cbb35 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -38,6 +38,8 @@ def xfaster_run( signal_type="synfast", signal_subset="*", signal_transfer_type=None, + signal_transfer_type_dust=None, + signal_subset_dust=None, signal_type_sim=None, noise_type="stationary", noise_subset="*", @@ -97,6 +99,7 @@ def xfaster_run( qb_file_sim=None, signal_spec=None, signal_transfer_spec=None, + signal_transfer_spec_dust=None, model_r=None, ref_freq=359.7, beta_ref=1.54, @@ -171,6 +174,20 @@ def xfaster_run( signal_transfer_type : str The variant of signal sims to use for computing the transfer function. If not set, defaults to ``signal_type``. + signal_transfer_type_dust : string + 1. The variant of signal simulation to use for transfer function + calculation for dust, typically identified by the input spectrum model + used to generate it, e.g 'power_law'. These are assumed to be at 353 + GHz and thus scaled by the map frequency's alpha. + This directory may also contain + a copy of the input spectrum, to make sure that the correct + spectrum is used to compute the transfer function. + 2. The variant of signal sims to use for transfer function for dust, + or, to add to signal_transfer type for raw (uncleaned) spectra, + multiplied by alpha90**2 or alpha150**2 + signal_transfer_dust: string + If not provided, use signal_subset values + Add more. signal_type_sim : str The variant of signal sims to use for sim_index data maps. This enables having a different noise sim ensemble to use for @@ -354,6 +371,11 @@ def xfaster_run( The spectrum data file used to generate signal sims. If not supplied, will search for ``spec_signal_.dat`` in the transfer signal sim directory. Used for computing transfer functions. + signal_transfer_spec_dust : string + The spectrum data file used to generate signal sims for dust transfer + function. If not supplied, will search for + `spec_signal_.dat` in the + dust transfer signal sim directory. Used for computing transfer functions. model_r : float The ``r`` value to use to compute a spectrum for estimating bandpowers. Overrides ``signal_spec``. @@ -499,6 +521,8 @@ def xfaster_run( signal_type=signal_type, signal_type_sim=signal_type_sim if sim_data_r is None else "r", signal_transfer_type=signal_transfer_type, + signal_transfer_type_dust=signal_transfer_type_dust, + signal_subset_dust=signal_subset_dust, data_root2=data_root2, data_subset2=data_subset2, foreground_type_sim=foreground_type_sim, @@ -579,6 +603,7 @@ def xfaster_run( fix_bb_transfer=fix_bb_transfer, signal_spec=signal_spec, signal_transfer_spec=signal_transfer_spec, + signal_transfer_spec_dust=signal_transfer_spec_dust, model_r=model_r, ref_freq=ref_freq, beta_ref=beta_ref, @@ -599,6 +624,7 @@ def xfaster_run( spec_opts.pop("multi_map") spec_opts.pop("signal_spec") spec_opts.pop("signal_transfer_spec") + spec_opts.pop("signal_transfer_spec_dust") spec_opts.pop("model_r") spec_opts.pop("qb_file") bandpwr_opts = spec_opts.copy() @@ -671,10 +697,32 @@ def xfaster_run( X.get_beams(pixwin=pixwin) X.log("Loading spectrum shape for transfer function...", "notice") - X.get_signal_shape(filename=signal_transfer_spec, transfer=True) + X.get_signal_shape( + #filename=signal_transfer_spec, transfer=True, + filename=[signal_transfer_spec, signal_transfer_spec_dust], + tbeb=False, + transfer=True, + transfer_dust=signal_transfer_type_dust is not None, + template_alpha90=template_alpha90, + template_alpha150=template_alpha150, + ) - X.log("Computing transfer functions...", "notice") + #X.log("Computing transfer functions...", "notice") + X.log("Computing transfer functions for CMB...", "task") X.get_transfer(**transfer_opts) + if signal_transfer_type_dust is not None: + X.log("Computing transfer functions for Dust...", "task") + X.get_transfer(cls_shape, fix_bb_xfer=fix_bb_xfer, dust=True, **fisher_opts) + X.log("Computing transfer functions for Dust+CMB...", "task") + # this is old code, may break + X.get_transfer( + cls_shape, + fix_bb_xfer=fix_bb_xfer, + cmbdust=True, + template_alpha90=template_alpha90, + template_alpha150=template_alpha150, + **fisher_opts + ) X.log("Computing sim ensemble averages...", "notice") X.get_masked_sims(qb_file=qb_file_sim) @@ -974,7 +1022,13 @@ def add_arg(P, name, argtype=None, default=None, short=None, help=None, **kwargs add_arg(G, "mask_type") add_arg(G, "signal_type") add_arg(G, "signal_subset") - add_arg(G, "signal_transfer_type") + #add_arg(G, "signal_transfer_type") + add_arg(G, "signal_transfer_type", + help="Signal sim variant for CMB transfer functions", + ) + add_arg(G, "signal_transfer_type_dust", + help="Signal sim variant for dust transfer functions", + ) add_arg(G, "signal_type_sim") add_arg(G, "noise_type") add_arg(G, "noise_subset") @@ -1057,6 +1111,11 @@ def add_arg(P, name, argtype=None, default=None, short=None, help=None, **kwargs add_arg(E, "signal_spec") add_arg(E, "model_r") add_arg(G, "signal_transfer_spec") + add_arg(G, "signal_transfer_spec_dust", + help="Power spectrum used to create transfer signal simulations for dust. " + "Defaults to the spec_signal_{signal_transfer_type_dust}.dat file found " + "in the signal directory.", + ) add_arg(G, "ref_freq") add_arg(G, "beta_ref", argtype=float) add_arg(G, "delta_beta_prior", argtype=float) From 698c9a7299020b726aa652ec44c40c172e1dcfac Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Mon, 20 Sep 2021 05:19:10 +0200 Subject: [PATCH 03/29] add dust transfer function options. Todo: add dust component into cls_model and dSdqb and make option for transfer function: 1/ from cmb signal, 2/ seperately from dust and cmb signal, 3/from (alpha*dust+cmb) signal --- xfaster/xfaster_class.py | 364 ++++++++++++++++++++++----------------- xfaster/xfaster_exec.py | 83 +++++---- 2 files changed, 250 insertions(+), 197 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index fad24b68..8398134a 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -572,9 +572,9 @@ def _get_files( mask_type="hitsmask_tailored", signal_type="r0p03", signal_transfer_type=None, - signal_transfer_type_dust=None, + signal_transfer_type_dust=None, signal_subset="*", - signal_subset_dust="None", + signal_subset_dust="None", noise_type="stationary", noise_subset="*", signal_type_sim=None, @@ -685,21 +685,21 @@ def find_sim_files(name, root=None, subset="*"): if signal_transfer_type is None: signal_transfer_type = signal_type find_sim_files( - "signal_transfer", "signal_{}".format(signal_transfer_type), signal_subset + "signal_transfer", "signal_{}".format(signal_transfer_type), + signal_subset ) - # transfer_dust - # do option to separate_dust_transfer - #if separate_dust_transfer: if signal_transfer_type_dust is None: -+ signal_transfer_root_dust = None -+ signal_transfer_files_dust = None - else: + find_sim_files("signal_transfer_dust", None) + else: if signal_subset_dust is None: - signal_subset_dust = signal_subset - find_sim_files( - "signal_transfer_dust", "signal_{}".format(signal_transfer_type_dust), signal_subset_dust - ) + signal_subset_dust = signal_subset # better choice ?? + # find reobserved scaled dust + # TODO: implement options for unscaled dust, start with signal_ + find_sim_files( + "signal_transfer_dust", "foreground_{}".format(signal_transfer_type_dust), + signal_subset_dust + ) # apply suffix out = dict() @@ -720,8 +720,8 @@ def get_files( signal_type="synfast", signal_type_sim=None, signal_transfer_type=None, - signal_transfer_type_dust=None, # transfer_dust - signal_subset_dust=None, + signal_transfer_type_dust=None, + signal_subset_dust=None, data_root2=None, data_subset2=None, foreground_type_sim=None, @@ -826,6 +826,13 @@ def get_files( to generate it, e.g 'synfast'. This directory may also contain a copy of the input spectrum, to make sure that the correct spectrum is used to compute the transfer function. + signal_transfer_type_dust : string + The variant of signal simulation to use for transfer function + calculation for dust, typically identified by the input spectrum model + used to generate it, e.g 'power_law'. These are assumed to be at 353 + GHz and thus scaled by the map frequency's alpha. + This directory may also contain a copy of the input spectrum, to make + sure that the correct spectrum is used to compute the transfer function. data_root2, data_subset2 : string The root and subset for a second set of data. If either of these is keywords is supplied, then the two data sets are treated as two @@ -915,11 +922,11 @@ def get_files( noise_type=noise_type, mask_type=mask_type, signal_type=signal_type, - signal_transfer_type=signal_transfer_type, - signal_transfer_type_dust=signal_transfer_type_dust, - signal_subset_dust=signal_subset_dust, + signal_transfer_type=signal_transfer_type, + signal_transfer_type_dust=signal_transfer_type_dust, signal_subset=signal_subset, noise_subset=noise_subset, + signal_subset_dust=signal_subset_dust, ) ref_opts = dict(data_subset=data_subset, **opts) if null_run: @@ -1858,7 +1865,12 @@ def alm2cl(self, m1, m2=None, lmin=2, lmax=None, symmetric=True): cls[..., :lmin] = 0 return np.atleast_2d(cls) - def get_mask_weights(self, apply_gcorr=False, reload_gcorr=False, gcorr_file=None): + def get_mask_weights( + self, + apply_gcorr=False, + reload_gcorr=False, + gcorr_file=None + ): """ Compute cross spectra of the masks for each data map. @@ -2061,7 +2073,7 @@ def process_index(idx): gmat[xname] = OrderedDict() for spec in self.specs: si, sj = spec_inds[spec] - f = (ky_eff[xname][si, sj] + fsky_eff[xname][sj, si]) / 2.0 + f = (fsky_eff[xname][si, sj] + fsky_eff[xname][sj, si]) / 2.0 gmat[xname][spec] = f if np.any(np.asarray([f for f in fsky.values()]) > 0.1): @@ -2699,7 +2711,12 @@ def process_index(idx): save_name, from_attrs=save_attrs, extra_tag="xcorr", data_opts=True ) - def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): + def get_masked_sims( + self, + transfer=False, + do_noise=True, + qb_file=None, + ): """ Compute average signal and noise spectra for a given ensemble of sim maps. The same procedure that is used for computing data cross spectra @@ -2731,6 +2748,8 @@ def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): cls_signal, cls_signal_null: Mean signal spectra + cls_signal_dust: + Mean dust signal spectra cls_noise, cls_noise_null: Mean noise spectra cls_sim, cls_sim_null: @@ -2759,12 +2778,14 @@ def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): signal_files = self.signal_transfer_files signal_files2 = self.signal_transfer_files2 if null_run else None num_signal = self.num_signal_transfer - signal_files_dust = self.signal_transfer_files_dust # how about null_run - num_signal_dust = self.num_signal_transfer_dust # SET UP THESE + signal_files_dust = self.signal_transfer_dust_files + num_signal_dust = self.num_signal_transfer_dust else: signal_files = self.signal_files signal_files2 = self.signal_files2 if null_run else None - num_signal = self.num_signal + num_signal = self.num_signal + num_signal_dust = 0 + noise_files = self.noise_files noise_files2 = self.noise_files2 if null_run else None @@ -2792,9 +2813,9 @@ def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): if transfer: save_name = "sims_xcorr_{}".format(self.signal_transfer_type) + if self.signal_transfer_type_dust is not None: + save_name += "_" + self.signal_transfer_type_dust cp = "sims_transfer" - if self.signal_transfer_type_dust is not None: - save_name += "_" + self.signal_transfer_type_dust else: save_name = "sims_xcorr_{}".format(self.signal_type) cp = "sims" @@ -2807,7 +2828,6 @@ def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): shape=data_shape, shape_ref="cls_signal", ) - # load cls_signal_dust if ret is not None: if do_noise: @@ -2826,9 +2846,9 @@ def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): if transfer and self.signal_transfer_type == self.signal_type: self.force_rerun["sims"] = False - # process signal, noise, and S+N + # process signal, signal dust, noise, and S+N cls_sig = OrderedDict() - cls_sig_dust = OrderedDict() + cls_sig_dust = OrderedDict() #if transfer_dust else None cls_null_sig = OrderedDict() if null_run else None cls_noise = OrderedDict() if do_noise else None cls_null_noise = OrderedDict() if null_run and do_noise else None @@ -2841,16 +2861,17 @@ def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): cls_res = OrderedDict() if do_noise else None cls_null_res = OrderedDict() if null_run and do_noise else None if do_noise: - for k in ["nxn0", "nxn1", "sxn0", "sxn1", "nxs0", "nxs1", "cls_signal_dust"]: + #for k in ["nxn0", "nxn1", "sxn0", "sxn1", "nxs0", "nxs1", "cls_signal_dust"]: + for k in ["nxn0", "nxn1", "sxn0", "sxn1", "nxs0", "nxs1"]: # cls_signal_dust seperate cls_res[k] = OrderedDict() if null_run: cls_null_res[k] = OrderedDict() if num_noise != 0: - nsim_min = min([num_signal, num_signal_dust, num_noise]) + nsim_min = min([num_signal, num_signal_dust, num_noise]) else: - nsim_min = num_signal - nsim_max = max([num_signal, num_signal_dust, num_noise]) + nsim_min = min([num_signal, num_signal_dust]) # not sure why we need if case here? + nsim_max = max([num_signal, num_signal_dust, num_noise]) cls_all = np.zeros( [nsim_max, len(map_pairs.items()), len(self.specs), self.lmax + 1] ) @@ -2918,7 +2939,7 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): return cache[idx] sig_cache = dict() - sig_dust_cache = dict() + sig_dust_cache = dict() noise_cache = dict() for isim in range(nsim_max): @@ -2949,14 +2970,14 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): if null_run: cls_null1t = np.copy(cls_null1_sig) - if isim < num_signal_dust: - dimap_alms, _ = process_index( - signal_files_dust, None, idx, isim, sig_dust_cache - ) - djmap_alms, djnull_alms = process_index( - signal_files_dust, None, jdx, isim, sig_dust_cache - ) - cls1_sig_dust = self.alm2cl(dimap_alms, djmap_alms) + if isim < num_signal_dust: + dimap_alms, _ = process_index( + signal_files_dust, None, idx, isim, sig_dust_cache + ) + djmap_alms, _ = process_index( + signal_files_dust, None, jdx, isim, sig_dust_cache + ) + cls1_sig_dust = self.alm2cl(dimap_alms, djmap_alms) if do_noise and isim < num_noise: cls1_res = OrderedDict() @@ -2982,10 +3003,10 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): # need non-symmetric since will potentially modify these # with different residuals for T, E, B - for k, cx, cy, cnx, cny in [ + for k, cx, cy, cnx, cny in [ ("nxn", nimap_alms, njmap_alms, ninull_alms, njnull_alms), ("sxn", simap_alms, njmap_alms, sinull_alms, njnull_alms), - ("nxs", nimap_alms, sjmap_alms, ninull_alms, sjnull_alms), + ("nxs", nimap_alms, sjmap_alms, ninull_alms, sjnull_alms), ]: k0, k1 = ["{}{}".format(k, i) for i in range(2)] cls1_res[k0] = self.alm2cl(cx, cy, symmetric=False) @@ -2994,7 +3015,7 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): cls_null1_res[k0] = self.alm2cl(cnx, cny, symmetric=False) cls_null1_res[k1] = self.alm2cl(cny, cnx, symmetric=False) - # construct noise model + # construct noise model cls1_noise = (cls1_res["nxn0"] + cls1_res["nxn1"]) / 2.0 if null_run: cls_null1_noise = ( @@ -3007,14 +3028,23 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): cls1t += cls1_res[k] / 2.0 if null_run: cls_null1t += cls_null1_res[k] / 2.0 + + # add dust to cls1t + if isim < num_signal_dust: + cls1t += cls1_sig_dust + # cls1_res["cls_signal_dust"] = cls1_sig_dust # not need, see quants + if null_run: + cls_null1t += cls1_sig_dust + #cls_null1_res["cls_signal_dust"] = cls1_sig_dust + quants = [] if isim < num_signal: quants += [[cls_sig, cls1_sig]] if null_run: quants += [[cls_null_sig, cls_null1_sig]] - if isim < num_signal_dust: - quants += [[cls_sig_dust, cls1_sig_dust]] + if isim < num_signal_dust: + quants += [[cls_sig_dust, cls1_sig_dust]] if do_noise and isim < num_noise: quants += [[cls_noise, cls1_noise]] @@ -3059,7 +3089,8 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): cls_null_med[spec][xname] = cls_null_med_arr[xind][s] self.cls_signal = cls_sig - self.cls_signal_dust = cls_sig_dust + if transfer: + self.cls_signal_dust = cls_sig_dust self.cls_signal_null = cls_null_sig self.cls_noise = cls_noise self.cls_noise_null = cls_null_noise @@ -3307,10 +3338,10 @@ def get_signal_shape( flat=None, signal_mask=None, transfer=False, - transfer_dust=False, + transfer_dust=False, save=True, - template_alpha90=None, - template_alpha150=None, + template_alpha90=None, + template_alpha150=None, ): """ Load a shape spectrum for input to the Fisher iteration algorithm. @@ -3366,13 +3397,14 @@ def get_signal_shape( Dictionary keyed by spectrum (cmb_tt, cmb_ee, ... , fg), each entry containing a vector of length 2 * lmax + 1 """ - if isinstance(filename, list): - filename_dust = filename[1] + if isinstance(filename, list): + filename_dust = filename[1] filename = filename[0] lmax_kern = 2 * self.lmax specs = list(self.specs) + if transfer: if "eb" in specs: specs.remove("eb") @@ -3387,6 +3419,19 @@ def get_signal_shape( if component == "fg": specs = ["fg"] + if transfer_dust: + specs += ["dust_tt"] + if self.pol: + for s in ["ee", "bb", "te"]: + specs += ["dust_{}".format(s)] + # Need separate fields per freq for dust+CMB to get relative + # amplitudes right + for freq in np.unique(list(self.nom_freqs.values())): + specs += ["cmbdust_tt_{}".format(freq)] + if self.pol: + for s in ["ee", "bb", "te"]: + specs += ["cmbdust_{}_{}".format(s, freq)] + nspecs = len(specs) if save: @@ -3412,7 +3457,7 @@ def get_signal_shape( ellfac = ell * (ell + 1) / 2.0 / np.pi cls_shape = OrderedDict() - def fix_camb_specs(arr, ltmp, filename): + def fix_camb_specs(arr, ltmp, filename): # camb starts at l=2, so set first 2 ells to be 0 # make sure order is TT, EE, BB, TE, no ell row. new_arr = np.append([0, 0], arr[1, : ltmp - 2]) @@ -3474,12 +3519,24 @@ def fix_camb_specs(arr, ltmp, filename): filename = os.path.join(self.config_root, filename) if not os.path.exists(filename): raise OSError("Missing model file {}".format(filename)) + + tmp = np.loadtxt(filename, unpack=True) + + ltmp = tmp.shape[-1] + if lmax_kern + 1 < ltmp: + ltmp = lmax_kern + 1 + else: + raise ValueError( + "Require at least lmax={} in model file, found {}".format( + lmax_kern, ltmp + ) + ) + tmp = fix_camb_specs(tmp, ltmp, filename) -# chech this - if transfer_dust: + if transfer_dust: if filename_dust is None: - signal_root = self.signal_transfer_root_dust + signal_root = self.signal_transfer_dust_root #self.signal_transfer_root_dust filename_dust = "spec_{}.dat".format(os.path.basename(signal_root)) filename_dust = os.path.join(signal_root, filename_dust) if not os.path.exists(filename_dust) and not os.path.isabs( @@ -3503,34 +3560,22 @@ def fix_camb_specs(arr, ltmp, filename): lmax_kern, ltmpd ) ) - tmpd = fix_camb_specs(tmpd, ltmpd, filename_dust) - - tmp = np.loadtxt(filename, unpack=True) - - ltmp = tmp.shape[-1] - if lmax_kern + 1 < ltmp: - ltmp = lmax_kern + 1 - else: - raise ValueError( - "Require at least lmax={} in model file, found {}".format( - lmax_kern, ltmp - ) - ) + tmpd = fix_camb_specs(tmpd, ltmpd, filename_dust) - tmp = fix_camb_specs(tmp, ltmp, filename) + # camb starts at l=2, so set first 2 ells to be 0 - cls_shape["cmb_tt"] = np.append([0, 0], tmp[1, : ltmp - 2]) - if self.pol: - if np.any(tmp[2, : ltmp - 2] < 0): + #cls_shape["cmb_tt"] = np.append([0, 0], tmp[1, : ltmp - 2]) + #if self.pol: + # if np.any(tmp[2, : ltmp - 2] < 0): # this is true if TE is the third index instead of EE # (ell is 0th index for CAMB) - self.warn( - "Old CAMB format in model file {}. Re-indexing.".format( - filename - ) - ) - pol_specs = [3, 4, 2] - spec_order = {"tt": 0, "ee": 1, "bb": 2, "te": 3} + # self.warn( + # "Old CAMB format in model file {}. Re-indexing.".format( + # filename + # ) + # ) + # pol_specs = [3, 4, 2] + spec_order = {"tt": 0, "ee": 1, "bb": 2, "te": 3} for spec in specs: parts = spec.split("_") comp = parts[0] @@ -3541,19 +3586,18 @@ def fix_camb_specs(arr, ltmp, filename): if comp == "cmb": arr = tmp elif comp == "dust": - arr = tmpd + arr = tmpd else: - pol_specs = [2, 3, 4] - for spec, d in zip(specs[1:4], tmp[pol_specs]): - cls_shape[spec] = np.append([0, 0], d[: ltmp - 2]) - # Combine CMB and dust scaled by alpha**2 for frequency + # pol_specs = [2, 3, 4] + #for spec, d in zip(specs[1:4], tmp[pol_specs]): + # cls_shape[spec] = np.append([0, 0], d[: ltmp - 2]) + # Combine CMB and dust scaled by alpha**2 for frequency if nom_freq == "90": alpha = template_alpha90 elif nom_freq == "150": alpha = template_alpha150 - arr = tmp + alpha ** 2 * tmpd - cls_shape[spec] = arr[spec_order[spec0]] - + arr = tmp + alpha ** 2 * tmpd + cls_shape[spec] = arr[spec_order[spec0]] # # EB and TB flat l^2 * C_l if self.pol: @@ -3575,25 +3619,6 @@ def fix_camb_specs(arr, ltmp, filename): "Added foreground to cls shape {}".format(list(cls_shape)), "debug" ) - if transfer_dust: - cls_shape["dust_tt"] = 0 - #nspecs += 1 - # Need separate fields per freq for dust+CMB to get relative - # amplitudes right - if self.pol: - cls_shape["dust_ee"] = 0 # add values here - cls_shape["dust_bb"] = 0 - cls_shape["dust_te"] = 0 - # nspecs += 3 - for freq in np.unique(list(self.nom_freqs.values())): - cls_shape["cmbdust_tt_{}".format(freq)] = 0 - # nspecs += 1 - if self.pol: - for s in ["ee", "bb", "te"]: - cls_shape["cmbdust_{}_{}".format(s, freq)] = 0 - nspecs += 3 - - # divide out l^2/2pi for spec in specs: cls_shape[spec][2:] /= ellfac[2:] @@ -4023,8 +4048,8 @@ def bin_cl_template( transfer_run=False, beam_error=False, use_precalc=True, - dust=False, - cmbdust=False, + dust=False, + cmbdust=False, ): """ Compute the Cbl matrix from the input shape spectrum. @@ -4052,6 +4077,11 @@ def bin_cl_template( If True, load pre-calculated terms stored from a previous iteration, and store for a future iteration. Otherwise, all calculations are repeated. + dust : bool + If True, compute transfer function using dust sims + cmbdust : bool + If True, compute transfer function using CMB+dust sims, with + dust scaled by template_alpha90, template_alpha150 Returns ------- @@ -4108,10 +4138,10 @@ def bin_cl_template( comps += ["cmb"] if "fg" in cls_shape and not transfer_run: comps += ["fg"] - if transfer_run and dust: - comps += ["dust"] - if transfer_run and cmbdust: - comps += ["dustcmb"] + if transfer_run and dust: + comps += ["dust"] + if transfer_run and cmbdust: + comps += ["dustcmb"] if self.nbins_res > 0 and not transfer_run: comps += ["res"] cls_noise = self.cls_noise_null if self.null_run else self.cls_noise @@ -4242,10 +4272,10 @@ def bin_things(comp, d, md): if comp == "fg": # single foreground spectrum s_arr[si, xi] = cls_shape["fg"][lk] - else: - # for comp in [cmb, dust, cmbdust] - s_arr[si, xi] = cls = cls_shape["{}_{}".format(comp, spec)][lk] - # cls_shape["cmb_{}".format(spec)][lk] + else: + # cls = cls_shape["cmb_{}".format(spec)][lk] + # for comp in [cmb, dust, cmbdust] + s_arr[si, xi] = cls_shape["{}_{}".format(comp, spec)][lk] # get cross spectrum kernel terms k_arr[si, xi] = mll[spec][xname][ls, lk] @@ -4496,9 +4526,9 @@ def get_data_spectra( map_tag=None, transfer_run=False, do_noise=True, - transfer_dust=False, - transfer_cmbdust=False, - template_alpha90=None, + transfer_dust=False, + transfer_cmbdust=False, + template_alpha90=None, template_alpha150=None, ): @@ -4544,8 +4574,7 @@ def get_data_spectra( # obs depends on what you're computing if transfer_run: - obs_quant = self.cls_signal - if transfer_dust: + if transfer_dust: obs_quant = self.cls_signal_dust elif transfer_cmbdust: # need to add cmb and dust scaled by alphas @@ -4560,7 +4589,7 @@ def get_data_spectra( + adict[freq1] * adict[freq2] * self.cls_signal_dust ) else: - obs_quant = self.cls_signal + obs_quant = self.cls_signal elif self.null_run: if self.reference_subtracted: @@ -5002,6 +5031,7 @@ def fisher_calc( Dmat1 = OrderedDict() if cls_model is None: + # TO DO: if dust, Eq 8 paper: model = dust + cmb cls_model = self.get_model_spectra( qb, cbl, delta=True, cls_noise=cls_noise, cond_noise=cond_noise ) @@ -5292,10 +5322,10 @@ def fisher_iterate( like_profile_sigma=3.0, like_profile_points=100, file_tag=None, - transfer_dust=dust, - transfer_cmbdust=cmbdust, - template_alpha90=tempalate_alpha90, - template_alpha150=template_alpha150, + transfer_dust=False, + transfer_cmbdust=False, + template_alpha90=None, + template_alpha150=None, ): """ Iterate over the Fisher calculation to compute bandpower estimates @@ -5421,12 +5451,16 @@ def fisher_iterate( else: qb = qb_start + # TO DO: GET TEMPLATE INTERNALLY INSTEAD OF AS INPUTS + #template_alpha90=self.template_alpha["90"] + #template_alpha150=self.template_alpha["150"] + + # reminder: obs = cmb + a*dust if transfer_cmbdust obs, nell, debias = self.get_data_spectra( - #map_tag=map_tag, transfer_run=transfer_run, map_tag=map_tag, transfer_run=transfer_run, transfer_dust=transfer_dust, - transfer_cmbdust=cmbdust, + transfer_cmbdust=transfer_cmbdust, template_alpha90=tempalate_alpha90, template_alpha150=tempalate_alpha150, ) @@ -5813,10 +5847,10 @@ def get_transfer( iter_max=200, save_iters=False, fix_bb_transfer=False, - dust=False, - cmbdust=False, - template_alpha90=None, - template_alpha150=None, + dust=False, + cmbdust=False, + template_alpha90=None, + template_alpha150=None, ): """ Compute the transfer function from signal simulations created using @@ -5839,6 +5873,11 @@ def get_transfer( fix_bb_transfer : bool If True, after transfer functions have been calculated, impose the BB xfer is exactly equal to the EE transfer. + dust : bool + If True, compute transfer function using dust sims + cmbdust : bool + If True, compute transfer function using CMB+dust sims, with + dust scaled by template_alpha90, template_alpha150 Returns ------- @@ -5875,10 +5914,10 @@ def get_transfer( ) save_name = "transfer_all" - if dust: - save_name = "{}_dust".format(save_name) - elif cmbdust: - save_name = "{}_cmbdust".format(save_name) + if dust: + save_name = "{}_dust".format(save_name) + elif cmbdust: + save_name = "{}_cmbdust".format(save_name) if self.weighted_bins: save_name = "{}_wbins".format(save_name) @@ -5915,7 +5954,7 @@ def expand_transfer(qb_transfer, bin_def): if ret is not None: #self.qb_transfer = ret["qb_transfer"] - if dust: + if dust: self.qb_transfer_dust = ret["qb_transfer"] elif cmbdust: self.qb_transfer_cmbdust = ret["qb_transfer"] @@ -5929,8 +5968,8 @@ def expand_transfer(qb_transfer, bin_def): self.transfer = expand_transfer(ret["qb_transfer"], ret["bin_def"]) return self.transfer - self.qb_transfer = OrderedDict() - qb_transfer = OrderedDict() + #self.qb_transfer = OrderedDict() + qb_transfer = OrderedDict() if dust: comp = "dust" elif cmbdust: @@ -5940,21 +5979,21 @@ def expand_transfer(qb_transfer, bin_def): for spec in self.specs: #self.qb_transfer["cmb_" + spec] = OrderedDict() - if dust: - qb_transfer["{}_{}".format(comp, spec)] = OrderedDict() + #if dust: This is not necessary ??? + qb_transfer["{}_{}".format(comp, spec)] = OrderedDict() - # success = False - success = True + # success = False + success = True msg = "" for im0, m0 in enumerate(self.map_tags): if not self.fit_transfer[self.map_tags_orig[im0]]: - #for spec in self.specs: + #for spec in self.specs: # self.qb_transfer["cmb_{}".format(spec)][m0] = np.ones( # self.nbins_cmb // len(self.specs) # ) - for stag in self.qb_transfer.keys(): - self.qb_transfer[stag][m0] = np.ones( + for stag in qb_transfer.keys(): + qb_transfer[stag][m0] = np.ones( self.nbins_cmb // len(self.specs) ) self.log("Setting map {} transfer to unity".format(m0), "info") @@ -5968,9 +6007,10 @@ def expand_transfer(qb_transfer, bin_def): ) self.clear_precalc() #cbl = self.bin_cl_template(map_tag=m0, transfer_run=True) - cbl = self.bin_cl_template( - cls_shape, m0, transfer_run=True, dust=dust, cmbdust=cmbdust + cbl = self.bin_cl_template( + m0, transfer_run=True, dust=dust, cmbdust=cmbdust ) + ret = self.fisher_iterate( cbl, m0, @@ -5978,10 +6018,10 @@ def expand_transfer(qb_transfer, bin_def): iter_max=iter_max, converge_criteria=converge_criteria, save_iters=save_iters, - transfer_dust=dust, + transfer_dust=dust, transfer_cmbdust=cmbdust, template_alpha90=tempalate_alpha90, - template_alpha150=template_alpha150, + template_alpha150=template_alpha150, ) qb = ret["qb"] @@ -6017,12 +6057,12 @@ def expand_transfer(qb_transfer, bin_def): if len(self.specs) > 4: #qb["cmb_eb"] = np.sqrt(np.abs(qb["cmb_ee"] * qb["cmb_bb"])) #qb["cmb_tb"] = np.sqrt(np.abs(qb["cmb_tt"] * qb["cmb_bb"])) - qb[comp + "_eb"] = np.sqrt(np.abs(qb[comp + "_ee"] * qb[comp + "_bb"])) - qb[comp + "_tb"] = np.sqrt(np.abs(qb[comp + "_tt"] * qb[comp + "_bb"])) + qb[comp + "_eb"] = np.sqrt(np.abs(qb[comp + "_ee"] * qb[comp + "_bb"])) + qb[comp + "_tb"] = np.sqrt(np.abs(qb[comp + "_tt"] * qb[comp + "_bb"])) for stag, qbdat in qb.items(): #self.qb_transfer[stag][m0] = qbdat - qb_transfer[stag][m0] = qbdat + qb_transfer[stag][m0] = qbdat if success: self.transfer = expand_transfer(self.qb_transfer, self.bin_def) @@ -6040,15 +6080,15 @@ def expand_transfer(qb_transfer, bin_def): if not success: raise RuntimeError("Error computing transfer function: {}".format(msg)) - if dust: + if dust: self.qb_transfer_dust = qb_transfer elif cmbdust: self.qb_transfer_cmbdust = qb_transfer else: - self.qb_transfer = qb_transfer + self.qb_transfer = qb_transfer - return qb_transfer - #return self.transfer # equivalent to self.qb_transfer in Anne's old code. + return qb_transfer + #return self.transfer def get_bandpowers( self, @@ -6068,6 +6108,8 @@ def get_bandpowers( like_profile_sigma=3.0, like_profile_points=100, file_tag=None, + dust=False, + cmbdust=False, ): """ Compute the maximum likelihood bandpowers of the data, assuming @@ -6114,9 +6156,9 @@ def get_bandpowers( Keep first CMB bandpowers fixed to input shape (qb=1). return_cls : bool If True, return C_ls rather than D_ls - dust : bool + dust : bool If True, compute transfer function using dust sims - cmbdust : bool + cmbdust : bool If True, compute transfer function using CMB+dust sims, with dust scaled by template_alpha90, template_alpha150 like_profiles : bool diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index 2b1cbb35..4d9e3b9a 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -38,7 +38,7 @@ def xfaster_run( signal_type="synfast", signal_subset="*", signal_transfer_type=None, - signal_transfer_type_dust=None, + signal_transfer_type_dust=None, signal_subset_dust=None, signal_type_sim=None, noise_type="stationary", @@ -99,7 +99,7 @@ def xfaster_run( qb_file_sim=None, signal_spec=None, signal_transfer_spec=None, - signal_transfer_spec_dust=None, + signal_transfer_spec_dust=None, model_r=None, ref_freq=359.7, beta_ref=1.54, @@ -174,20 +174,14 @@ def xfaster_run( signal_transfer_type : str The variant of signal sims to use for computing the transfer function. If not set, defaults to ``signal_type``. + signal_subset_dust: str + If not provided, use signal_subset values + Add more. signal_transfer_type_dust : string - 1. The variant of signal simulation to use for transfer function - calculation for dust, typically identified by the input spectrum model - used to generate it, e.g 'power_law'. These are assumed to be at 353 - GHz and thus scaled by the map frequency's alpha. - This directory may also contain - a copy of the input spectrum, to make sure that the correct - spectrum is used to compute the transfer function. - 2. The variant of signal sims to use for transfer function for dust, + The variant of signal sims to use for transfer function for dust, or, to add to signal_transfer type for raw (uncleaned) spectra, multiplied by alpha90**2 or alpha150**2 - signal_transfer_dust: string - If not provided, use signal_subset values - Add more. + If provided, will compute separate transfer functions for dust and cmb. signal_type_sim : str The variant of signal sims to use for sim_index data maps. This enables having a different noise sim ensemble to use for @@ -521,7 +515,7 @@ def xfaster_run( signal_type=signal_type, signal_type_sim=signal_type_sim if sim_data_r is None else "r", signal_transfer_type=signal_transfer_type, - signal_transfer_type_dust=signal_transfer_type_dust, + signal_transfer_type_dust=signal_transfer_type_dust, signal_subset_dust=signal_subset_dust, data_root2=data_root2, data_subset2=data_subset2, @@ -624,7 +618,7 @@ def xfaster_run( spec_opts.pop("multi_map") spec_opts.pop("signal_spec") spec_opts.pop("signal_transfer_spec") - spec_opts.pop("signal_transfer_spec_dust") + spec_opts.pop("signal_transfer_spec_dust") spec_opts.pop("model_r") spec_opts.pop("qb_file") bandpwr_opts = spec_opts.copy() @@ -688,17 +682,27 @@ def xfaster_run( X.log("Computing kernels...", "notice") X.get_kernels(window_lmax=window_lmax) - X.log("Computing sim ensemble averages for transfer function...", "notice") + X.log("Computing sim ensemble averages for transfer function...", + "notice") # Do all the sims at once to also get the S+N sim ensemble average do_noise = signal_transfer_type in [signal_type, None] - X.get_masked_sims(transfer=True, do_noise=do_noise, qb_file=qb_file_sim) + + X.get_masked_sims( + transfer=True, + do_noise=do_noise, + qb_file=qb_file_sim) X.log("Computing beam window functions...", "notice") X.get_beams(pixwin=pixwin) + # todo: load template_alpha in get_signal_shape + template_alpha90 = template_alpha["90"] + template_alpha150 = template_alpha["150a"] + X.log("Loading spectrum shape for transfer function...", "notice") - X.get_signal_shape( - #filename=signal_transfer_spec, transfer=True, + # load signal shapes for dust, cmb, and cmbdust + X.get_signal_shape( + #filename=signal_transfer_spec, transfer=True, filename=[signal_transfer_spec, signal_transfer_spec_dust], tbeb=False, transfer=True, @@ -707,22 +711,29 @@ def xfaster_run( template_alpha150=template_alpha150, ) - #X.log("Computing transfer functions...", "notice") - X.log("Computing transfer functions for CMB...", "task") - X.get_transfer(**transfer_opts) - if signal_transfer_type_dust is not None: - X.log("Computing transfer functions for Dust...", "task") - X.get_transfer(cls_shape, fix_bb_xfer=fix_bb_xfer, dust=True, **fisher_opts) - X.log("Computing transfer functions for Dust+CMB...", "task") - # this is old code, may break - X.get_transfer( - cls_shape, - fix_bb_xfer=fix_bb_xfer, - cmbdust=True, - template_alpha90=template_alpha90, - template_alpha150=template_alpha150, - **fisher_opts - ) + #X.log("Computing transfer functions for CMB...", "task") + #X.get_transfer( + # fix_bb_transfer=fix_bb_transfer, + # **transfer_opts) + if (signal_transfer_spec_dust is not None): + if separate_dust_transfer_function: + X.log("Computing transfer functions for Dust...", "task") + X.get_transfer( + fix_bb_transfer=fix_bb_transfer, + dust=True, + **transfer_opts) + else: + X.log("Computing transfer functions for Dust+CMB...", "task") + X.get_transfer( + cmbdust=True, + fix_bb_transfer=fix_bb_transfer, + template_alpha90=template_alpha90, + template_alpha150=template_alpha150, + **transfer_opts, + ) + else: + X.log("Computing same transfer function for dust and CMB...", "task") + X.get_transfer(**transfer_opts) X.log("Computing sim ensemble averages...", "notice") X.get_masked_sims(qb_file=qb_file_sim) @@ -1026,7 +1037,7 @@ def add_arg(P, name, argtype=None, default=None, short=None, help=None, **kwargs add_arg(G, "signal_transfer_type", help="Signal sim variant for CMB transfer functions", ) - add_arg(G, "signal_transfer_type_dust", + add_arg(G, "signal_transfer_type_dust", help="Signal sim variant for dust transfer functions", ) add_arg(G, "signal_type_sim") From 6fed6f0d8a36974147f22904dbe2b211cef50eb7 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Fri, 24 Sep 2021 04:48:17 +0200 Subject: [PATCH 04/29] adding seperate dust Fl and Fl for (cmb + alpha*dust) options --- xfaster/xfaster_class.py | 84 +++++++++++++++++++++++++++++++++------- xfaster/xfaster_exec.py | 38 +++++++++++++----- 2 files changed, 97 insertions(+), 25 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 8398134a..28a49d2f 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -3339,6 +3339,7 @@ def get_signal_shape( signal_mask=None, transfer=False, transfer_dust=False, + #separate_dust_Fl=False, save=True, template_alpha90=None, template_alpha150=None, @@ -3424,6 +3425,7 @@ def get_signal_shape( if self.pol: for s in ["ee", "bb", "te"]: specs += ["dust_{}".format(s)] + # CHECK: DO I NEED THIS FOR SEPARATE DUST Fl ??? # Need separate fields per freq for dust+CMB to get relative # amplitudes right for freq in np.unique(list(self.nom_freqs.values())): @@ -3592,6 +3594,7 @@ def fix_camb_specs(arr, ltmp, filename): #for spec, d in zip(specs[1:4], tmp[pol_specs]): # cls_shape[spec] = np.append([0, 0], d[: ltmp - 2]) # Combine CMB and dust scaled by alpha**2 for frequency + # TODO: GENERALIZE THIS FOR FREQs if nom_freq == "90": alpha = template_alpha90 elif nom_freq == "150": @@ -3962,7 +3965,7 @@ def get_bin_def( return self.bin_def - def kernel_precalc(self, map_tag=None, transfer_run=False): + def kernel_precalc(self, dust=False, map_tag=None, transfer_run=False): """ Compute the mixing kernels M_ll' = K_ll' * F_l' * B_l'^2. Called by ``bin_cl_template`` to pre-compute kernel terms. @@ -3973,6 +3976,8 @@ def kernel_precalc(self, map_tag=None, transfer_run=False): If supplied, the kernels are computed only for the given map tag (or cross if map_tag is map_tag1:map_tag2). Otherwise, it is computed for all maps and crosses. + dust: bool + If true, get mll with dust Fl transfer_run : bool If True, set transfer function to 1 to solve for transfer function qbs. @@ -3982,6 +3987,7 @@ def kernel_precalc(self, map_tag=None, transfer_run=False): mll : OrderedDict Dictionary of M_ll' matrices, keyed by spec and xname. """ + # add one dim to add map_pairs = None if map_tag is not None: @@ -4004,7 +4010,10 @@ def kernel_precalc(self, map_tag=None, transfer_run=False): for spec in specs: transfer[spec] = OrderedDict() for tag in map_tags: - transfer[spec][tag] = self.transfer[spec][tag] + if dust: + transfer[spec][tag] = self.transfer_dust[spec][tag] + else: + transfer[spec][tag] = self.transfer[spec][tag] lk = slice(0, lmax + 1) mll = OrderedDict() @@ -4125,6 +4134,14 @@ def bin_cl_template( else: mll = self.mll + if not hasattr(self, "transfer_dust") and not use_precalc: + mll_dust = self.kernel_precalc(map_tag=map_tag, dust=True, + transfer_run=transfer_run) + if use_precalc: + self.mll_dust = mll_dust + else: + mll_dust = self.mll_dust + if beam_error: beam_error = self.get_beam_errors() beam_keys = ["b1", "b2", "b3"] @@ -4141,7 +4158,8 @@ def bin_cl_template( if transfer_run and dust: comps += ["dust"] if transfer_run and cmbdust: - comps += ["dustcmb"] + #comps += ["dustcmb"] + comps += ["cmbdust"] if self.nbins_res > 0 and not transfer_run: comps += ["res"] cls_noise = self.cls_noise_null if self.null_run else self.cls_noise @@ -4278,10 +4296,17 @@ def bin_things(comp, d, md): s_arr[si, xi] = cls_shape["{}_{}".format(comp, spec)][lk] # get cross spectrum kernel terms - k_arr[si, xi] = mll[spec][xname][ls, lk] + if comp == "dust" + k_arr[si, xi] = mll_dust[spec][xname][ls, lk] + else: + k_arr[si, xi] = mll[spec][xname][ls, lk] + if spec in ["ee", "bb"]: mspec = spec + "_mix" - mk_arr[si - 1, xi] = mll[mspec][xname][ls, lk] + if comp == "dust": + mk_arr[si - 1, xi] = mll_dust[mspec][xname][ls, lk] + else: + mk_arr[si - 1, xi] = mll[mspec][xname][ls, lk] md = None if s_arr.ndim == 1: @@ -4301,7 +4326,9 @@ def bin_things(comp, d, md): return cbl def get_model_spectra( - self, qb, cbl, delta=True, res=True, cls_noise=None, cond_noise=None + self, qb, cbl, delta=True, res=True, cls_noise=None, cond_noise=None, + transfer_dust=False, + transfer_cmbdust=False, ): """ Compute unbinned model spectra from qb amplitudes and a Cbl matrix. @@ -4343,11 +4370,15 @@ def get_model_spectra( """ comps = [] - if any([k.startswith("cmb_") for k in qb]): + if transfer_dust: + comps = ["dust"] + elif transfer_cmbdust: + comps = ["cmbdust"] + elif any([k.startswith("cmb_") for k in qb]): comps = ["cmb"] delta_beta = 0.0 - if "delta_beta" in qb: + if "delta_beta" in qb and not transfer_dust: # Evaluate fg at spectral index pivot for derivative # in Fisher matrix, unless delta is True if delta: @@ -4372,6 +4403,7 @@ def get_model_spectra( specs = [] for spec in self.specs: + # cmb, dust, and cmbdust should have same specs (Fl cases) if "cmb_{}".format(spec) in cbl: # Don't add entries that won't be filled in later cls["total_{}".format(spec)] = OrderedDict() @@ -4400,7 +4432,7 @@ def get_model_spectra( tag1, tag2 = self.map_pairs[xname] # extract qb's for the component spectrum - if comp == "cmb": + if comp in ["cmb", "dust", "cmbdust"]: qbs = qb[stag] if spec in ["ee", "bb"]: qbm = qb[mstag] @@ -4474,7 +4506,7 @@ def get_model_spectra( cl1 /= 2.0 # compute model spectra - if comp in ["cmb", "fg"]: + if comp in ["cmb", "fg", "dust", "cmb_dust"]: if xname not in cbl[stag]: continue cbl1 = cbl[stag][xname] @@ -4749,6 +4781,8 @@ def fisher_precalc( cls_debias=None, likelihood=False, windows=False, + transfer_dust=False, + transfer_cmbdust=False, ): """ Pre-compute the D_ell and signal derivative matrices necessary for @@ -4792,7 +4826,13 @@ def fisher_precalc( pol_dim = 3 if self.pol else 1 dim1 = pol_dim * num_maps - comps = ["cmb"] + if transfer_dust: + comps = ["dust"] + elif transfer_cmbdust: + comps = ["cmbdust"] + else: + comps = ["cmb"] + if "fg_tt" in cbl: comps += ["fg"] if "res_tt" in cbl or "res_ee" in cbl: @@ -4873,7 +4913,8 @@ def fisher_precalc( mix_cbl = cbl[stag + "_mix"][xname] dSdqb[comp][xname][spec][mspec] = mix_cbl - if comp == "fg": + if comp == "fg" and not transfer_dust: + # currently don't solve both beta_d and dust power. # add delta beta bin for spectral index dSdqb.setdefault("delta_beta", OrderedDict()).setdefault( xname, OrderedDict() @@ -4911,6 +4952,8 @@ def fisher_calc( use_precalc=True, windows=False, inv_fish=None, + transfer_dust=False, + transfer_cmbdust=False, ): """ Re-compute the Fisher matrix and qb amplitudes based on @@ -4997,6 +5040,8 @@ def fisher_calc( cls_debias=cls_debias, likelihood=likelihood, windows=windows, + transfer_dust=transfer_dust, + transfer_cmbdust=transfer_cmbdust, ) if use_precalc: self.Dmat_obs = Dmat_obs @@ -5033,7 +5078,9 @@ def fisher_calc( if cls_model is None: # TO DO: if dust, Eq 8 paper: model = dust + cmb cls_model = self.get_model_spectra( - qb, cbl, delta=True, cls_noise=cls_noise, cond_noise=cond_noise + qb, cbl, delta=True, cls_noise=cls_noise, cond_noise=cond_noise, + transfer_dust=transfer_dust, + transfer_cmbdust=transfer_cmbdust, ) mkeys = list(cls_model) @@ -5092,7 +5139,9 @@ def fisher_calc( cond_iter = 0 while not well_cond: cls_model = self.get_model_spectra( - qb, cbl, delta=True, cls_noise=cls_noise, cond_noise=cond_noise + qb, cbl, delta=True, cls_noise=cls_noise, cond_noise=cond_noise, + transfer_dust=transfer_dust, + transfer_cmbdust=transfer_cmbdust, ) for xname, (m0, m1) in self.map_pairs.items(): @@ -5483,6 +5532,8 @@ def fisher_iterate( delta_beta_prior=delta_beta_prior, cond_criteria=cond_criteria, null_first_cmb=null_first_cmb, + transfer_dust=transfer_dust, + transfer_cmbdust=transfer_cmbdust, ) qb_arr = pt.dict_to_arr(qb, flatten=True) @@ -6065,7 +6116,10 @@ def expand_transfer(qb_transfer, bin_def): qb_transfer[stag][m0] = qbdat if success: - self.transfer = expand_transfer(self.qb_transfer, self.bin_def) + if dust: + self.transfer_dust = expand_transfer(self.qb_transfer, self.bin_def) + else: + self.transfer = expand_transfer(self.qb_transfer, self.bin_def) else: self.transfer = None diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index 4d9e3b9a..5a9fab2b 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -100,6 +100,7 @@ def xfaster_run( signal_spec=None, signal_transfer_spec=None, signal_transfer_spec_dust=None, + separate_dust_Fl=False, model_r=None, ref_freq=359.7, beta_ref=1.54, @@ -695,7 +696,7 @@ def xfaster_run( X.log("Computing beam window functions...", "notice") X.get_beams(pixwin=pixwin) - # todo: load template_alpha in get_signal_shape + # TODO: load template_alpha in get_signal_shape template_alpha90 = template_alpha["90"] template_alpha150 = template_alpha["150a"] @@ -704,7 +705,7 @@ def xfaster_run( X.get_signal_shape( #filename=signal_transfer_spec, transfer=True, filename=[signal_transfer_spec, signal_transfer_spec_dust], - tbeb=False, + #tbeb=False, transfer=True, transfer_dust=signal_transfer_type_dust is not None, template_alpha90=template_alpha90, @@ -715,25 +716,40 @@ def xfaster_run( #X.get_transfer( # fix_bb_transfer=fix_bb_transfer, # **transfer_opts) + dust = False + cmbdust = False if (signal_transfer_spec_dust is not None): - if separate_dust_transfer_function: + if separate_dust_Fl: + dust = True + X.log("Computing transfer functions for CMB...", "task") + X.get_transfer( + fix_bb_transfer=fix_bb_transfer, + **transfer_opts, + ) X.log("Computing transfer functions for Dust...", "task") X.get_transfer( fix_bb_transfer=fix_bb_transfer, - dust=True, - **transfer_opts) + dust=dust, + template_alpha90=template_alpha90, + template_alpha150=template_alpha150, + **transfer_opts, + ) else: + cmbdust=True X.log("Computing transfer functions for Dust+CMB...", "task") X.get_transfer( - cmbdust=True, + cmbdust=cmbdust, fix_bb_transfer=fix_bb_transfer, template_alpha90=template_alpha90, template_alpha150=template_alpha150, **transfer_opts, - ) + ) else: - X.log("Computing same transfer function for dust and CMB...", "task") - X.get_transfer(**transfer_opts) + X.log("Computing transfer function...", "task") + X.get_transfer( + fix_bb_transfer=fix_bb_transfer, + **transfer_opts + ) X.log("Computing sim ensemble averages...", "notice") X.get_masked_sims(qb_file=qb_file_sim) @@ -757,7 +773,9 @@ def xfaster_run( if multi_map: X.log("Computing multi-map bandpowers...", "notice") - qb, inv_fish = X.get_bandpowers(return_qb=True, **bandpwr_opts) + qb, inv_fish = X.get_bandpowers(return_qb=True, + dust=dust, cmbdust=cmbdust, + **bandpwr_opts) if likelihood: X.log("Computing multi-map likelihood...", "notice") From 5619024a33c10a710b67fa7c103677117e32801e Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Tue, 12 Oct 2021 22:50:23 +1300 Subject: [PATCH 05/29] cleanup foreground model functions, add general load_camb_cl function for reading spectra from disk --- xfaster/spec_tools.py | 253 +++++++++++++++++++++++++++++++----------- 1 file changed, 190 insertions(+), 63 deletions(-) diff --git a/xfaster/spec_tools.py b/xfaster/spec_tools.py index 79063a6a..9051a283 100644 --- a/xfaster/spec_tools.py +++ b/xfaster/spec_tools.py @@ -2,117 +2,168 @@ from __future__ import absolute_import from __future__ import division import numpy as np +import warnings __all__ = [ "wigner3j", "get_camb_cl", + "load_camb_cl", "scale_dust", + "dust_model", ] +k = 1.38064852e-23 # Boltzmann constant +h = 6.626070040e-34 # Planck constant +hk = h / k * 1.0e9 # conversion from GHz / K to unitless value +T_dust = 19.6 # dust temperature in K +T_cmb = 2.72548 # CMB blackbody temperature in K -def blackbody(nu, ref_freq=353.0): + +def blackbody(nu, nu0=None, T=T_cmb): """ - The ratio of the blackbody function for dust at frequency nu - over the value for reference frequency ref_freq + The blackbody function at temperature ``T`` and frequency ``nu``, optionally + relative to the value at reference frequency ``nu0``. Arguments --------- nu : float Frequency in GHz. - ref_freq : float - Reference frequency in GHz. + nu0 : float + Reference frequency in GHz. If not supplied, return unreferenced + blackbody. + T : float + Blackbody temperature. Returns ------- - blackbody_ratio : float - B(nu, T_dust) / B(nu_ref, T_dust) + blackbody : float + B(nu, T) / B(nu0, T) """ - k = 1.38064852e-23 # Boltzmann constant - h = 6.626070040e-34 # Planck constant - T = 19.6 - nu_ref = ref_freq * 1.0e9 - nu *= 1.0e9 # GHz -> Hz - x = h * nu / k / T - x_ref = h * nu_ref / k / T - return x ** 3 / x_ref ** 3 * (np.exp(x_ref) - 1) / (np.exp(x) - 1) + x = hk * nu / T + if nu0 is None: + return x ** 3 / (np.exp(x) - 1) + x0 = hk * nu0 / T + return (x / x0) ** 3 * (np.exp(x0) - 1) / (np.exp(x) - 1) -def rj2cmb(nu_in): +def rj2bb(nu, nu0=None, T=T_cmb): """ - Conversion from Rayleigh-Jeans units to CMB temperature units + Conversion from Rayleigh-Jeans units to Blackbody temperature units Arguments --------- - nu_in : float + nu : float Frequency in GHz. + nu0 : float + Reference frequency in GHz. Ignored if not supplied. + T : float + Blakbody temperature Returns ------- cal_fac : float - Number by which to multiply a RJ temperature to get a CMB temp + Value by which to multiply a RJ temperature to get a CMB temperature. + If ``nu0`` is given, return the value relative to that at the reference + frequency. """ - k = 1.38064852e-23 # Boltzmann constant - h = 6.626070040e-34 # Planck constant - T = 2.72548 # Cmb BB temp in K - nu = nu_in * 1.0e9 # GHz -> Hz - x = h * nu / k / T - return (np.exp(x) - 1.0) ** 2 / (x ** 2 * np.exp(x)) + x = hk * nu / T + if nu0 is None: + return ((np.exp(x) - 1.0) / x) ** 2 / np.exp(x) + x0 = hk * nu0 / T + return ((np.exp(x) - 1.0) / (np.exp(x0) - 1.0) * x0 / x) ** 2 * np.exp(x0 - x) -def scale_dust(freq0, freq1, ref_freq, beta, delta_beta=None, deriv=False): +def scale_dust(nu1, nu2=None, nu0=359.7, beta=1.54, delta=False): """ Get the factor by which you must multiply the cross spectrum from maps of - frequencies freq0 and freq1 to match the dust power at ref_freq given - spectra index beta. - - If deriv is True, return the frequency scaling at the reference beta, - and the first derivative w.r.t. beta. - - Otherwise if delta_beta is given, return the scale factor adjusted - for a linearized offset delta_beta from the reference beta. + frequencies ``nu1`` and ``nu2`` to match the dust power at ``nu0`` given + spectral index ``beta``. If ``nu2`` is not given, compute the map-domain + factor for ``nu1`` alone. Optionally include terms to construct the + derivative with respect to ``beta`` if ``delta=True``. Arguments --------- - freq0 : float + nu1 : float Frequency of map0 in GHz. - freq1 : float - Frequency of map1 in GHz. - ref_freq : float + nu2 : float + Frequency of map1 in GHz. If not supplied, return the map-domain + multiplicative factor for ``nu1`` only. + nu0 : float Reference frequency from which to compute relative scaling in GHz. beta : float Dust spectral index. - delta_beta : float - Difference from beta-- scaling computed as a first order Taylor - expansion from original beta-scaling. - deriv : bool - If true, return the frequency scaling at the reference beta, along with - the first derivative w.r.t. beta at the reference beta. + delta : bool + If True, also return the multiplicative beta scaling and its first + derivative. Returns ------- freq_scale : float - The relative scaling factor for the dust cross spectrum-- multiply by - this number to get the dust spectrum at the reference frequency - -- or -- - freq_scale, deriv : floats - The relative scaling factor and its derivative + The relative scaling factor for the dust map or spectrum-- multiply by + this number to get the dust map or spectrum at the reference frequency + beta_scale : float + The multiplicative scaling factor for beta-- multiply the + frequency-scaled map or spectrum by ``beta_scale ** delta_beta`` to + obtain the dust map or spectrum at the adjusted spectral index ``beta + + delta_beta``. + log_beta_scale : float + The first derivative of the ``beta`` scaling-- multiply the + frequency-scaled map or spectrum by this number to obtain the derivative + of the frequency-scaled map or spectrum with respect to ``beta``. """ - freq_scale = ( - rj2cmb(freq0) - * rj2cmb(freq1) - / rj2cmb(ref_freq) ** 2.0 - * blackbody(freq0, ref_freq=ref_freq) - * blackbody(freq1, ref_freq=ref_freq) - * (freq0 * freq1 / ref_freq ** 2) ** (beta - 2.0) + bs = nu1 / nu0 + fs = ( + rj2bb(nu1, nu0=nu0, T=T_cmb) + * blackbody(nu1, nu0=nu0, T=T_dust) + * bs ** (beta - 2.0) ) - if deriv or delta_beta is not None: - delta = np.log(freq0 * freq1 / ref_freq ** 2) - if deriv: - return (freq_scale, freq_scale * delta) - return freq_scale * (1 + delta * delta_beta) + if nu2 is None: + if delta: + return (fs, bs, np.log(bs)) + return fs + + if delta: + fs2, bs2, log_bs2 = scale_dust(nu2, nu0=nu0, beta=beta, delta=True) + return (fs * fs2, bs * bs2, np.log(bs * bs2)) + + return fs * scale_dust(nu2, nu0=nu0, beta=beta, delta=False) + + +def dust_model(ell, pivot=80, amp=34.0, index=-2.28, lfac=True): + """ + Construct a power-law dust power spectrum. Default parameter values are for + the Planck LIV best-fit EE dust spectrum. + + The functional form is:: + + model = amp * (ell / pivot) ** index - return freq_scale + Arguments + --------- + ell : array_like + Multipoles over which to compute the model + pivot : scalar + Pivot ell at which to refence the amplitude + amp : scalar + Model amplitude + index : scalar + Model spectral index + lfac: bool + If True, multiply Cls by ell*(ell+1)/2/pi + + Returns + ------- + model : array_like + The dust model computed for each input ell value. + """ + ell = np.asarray(ell) + model = np.zeros(len(ell), dtype=float) + model[ell > 1] = amp * (ell[ell > 1] / float(pivot)) ** (index + 2.0) + if not lfac: + ellfac = ell * (ell + 1) / 2.0 / np.pi + model[ell > 1] /= ellfac[ell > 1] + return model def wigner3j(l2, m2, l3, m3): @@ -179,7 +230,7 @@ def get_camb_cl(r, lmax, nt=None, spec="total", lfac=True): Returns ------- cls : array_like - Array of spectra of shape (lmax + 1, nspec). + Array of spectra of shape (nspec, lmax + 1). Diagonal ordering (TT, EE, BB, TE). """ # Set up a new set of parameters for CAMB @@ -217,3 +268,79 @@ def get_camb_cl(r, lmax, nt=None, spec="total", lfac=True): totCL = powers[spec][: lmax + 1, :4].T return totCL + + +def load_camb_cl(filename, lmax=None, pol=None, lfac=True): + """ + Load a CAMB spectrum from a text file. Expects a file with at + least two columns, where the first column contains ell values, and + the rest are spectrum components (TT, EE, BB, TE). Old-style CAMB + files, where the spectra are ordered as (TT, TE, EE, BB), are + re-indexed into the new-style ordering. + + Arguments + --------- + filename : str + Path to spectrum file on disk. + lmax : int + If supplied, return spectrum up to this ell. Otherwise, all ell's are + included in the output. + pol : bool + If True, include polarization spectra in the output, if present in the + file. If False, include only the TT spectrum. Otherwise, include all + available spectra. + lfac : bool + If True, multiply Cls by ell*(ell+1)/2/pi + + Returns + ------- + cls : array_like + Array of spectra of shape (nspec, lmax + 1). + Diagonal ordering (TT, EE, BB, TE, EB, TB). + """ + + data = np.loadtxt(filename, unpack=True) + ell, data = data[0].astype(int), data[1:] + + # limit to lmax + if lmax is not None: + if lmax > ell.max(): + raise ValueError( + "Require at least lmax={} in CAMB file, found {}".format( + lmax, ell.max() + ) + ) + idx = ell <= lmax + ell = ell[idx] + data = data[..., idx] + + # include monopole and dipole + if ell.min() > 0: + n = ell.min() + ell = np.append(np.arange(n), ell) + data = np.append(np.zeros((len(data), n)), data, axis=1) + + # check polarization + if pol is None: + pol = len(data) > 1 + + # select columns + if pol: + if len(data) == 1: + raise ValueError( + "Missing polarization spectra in CAMB file {}".format(filename) + ) + if np.any(data[1, : ell.max() - 2] < 0): + warnings.warn("Old CAMB format file {}. Re-indexing.".format(filename)) + order = [0, 2, 3, 1] if len(data) == 4 else [0, 3, 5, 1, 4, 2] + data = data[order, :] + else: + data = data[[0], :] + + # check scaling + if not lfac: + f = ell * (ell + 1) / 2.0 / np.pi + f[ell == 0] = 1 + data /= f[None, :] + + return data From b64f07c9d21ffcf2f416319e45283dc40f9b4b1f Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Tue, 12 Oct 2021 22:51:20 +1300 Subject: [PATCH 06/29] * add foreground sims and data to example map ensemble * use separate seeds for each sim and component * clean up code to not rebuild existing files --- example/make_example_maps.py | 200 +++++++++++++++++++++++------------ 1 file changed, 134 insertions(+), 66 deletions(-) diff --git a/example/make_example_maps.py b/example/make_example_maps.py index 0d802a73..982ee3c3 100644 --- a/example/make_example_maps.py +++ b/example/make_example_maps.py @@ -8,26 +8,32 @@ xfaster_example.py script, including fake data maps, signal sims, noise sims, and masks. -Script should take <10 minutes and write about 240 M of maps to disk +Script should take <10 minutes and write about 350 M of maps to disk """ # set up options nsim = 100 nside = 256 -seed0 = 0 -np.random.seed(seed0) +do_fg = True + +betad = 1.54 +ref_freq = 359.7 tags = [95, 150] +freqs = [94.7, 151.0] +scales = [xf.spec_tools.scale_dust(f, nu0=ref_freq, beta=betad) for f in freqs] fwhms = [41.0, 29.0] # arcmin gnoise = [4.5, 3.5] # Temp std, uK_CMB data_noise_scale = 0.85 # make data noise less than sims by 15% mask_path = "maps_example/masks_rectangle" data_path = "maps_example/data_raw/full" +data_fg_path = "maps_example/data_cmbfg/full" sig_path = "maps_example/signal_synfast/full" noise_path = "maps_example/noise_gaussian/full" +fg_path = "maps_example/foreground_gaussian/full" -for p0 in [mask_path, data_path, sig_path, noise_path]: +for p0 in [mask_path, data_path, data_fg_path, sig_path, noise_path, fg_path]: if not os.path.exists(p0): os.makedirs(p0) @@ -39,83 +45,145 @@ # spectrum for signal sims-- synfast expects Cls spec_file = os.path.join(*sig_path.split("/")[:-1], "spec_signal_synfast.dat") if os.path.exists(spec_file): - dls = np.loadtxt(spec_file, unpack=True)[1:] - cls = np.hstack([np.zeros((len(dls), 2)), dls / lfac[2:]]) + cls = xf.spec_tools.load_camb_cl(spec_file, lfac=False) else: cls = xf.get_camb_cl(r=1.0, lmax=2000, lfac=False) # write to disk for transfer function dls = np.vstack([ell[2:], (lfac * cls)[:, 2:]]) np.savetxt(spec_file, dls.T) +# spectrum for foreground sims +fg_spec_file = os.path.join(*fg_path.split("/")[:-1], "spec_foreground_gaussian.dat") +if os.path.exists(spec_file): + cls_fg = xf.spec_tools.load_camb_cl(spec_file, lfac=False) +else: + cls_fg = xf.spec_tools.dust_model(ell, lfac=False) + # write to disk for transfer function + dls_fg = np.vstack([ell[2:], (lfac * cls_fg)[:, 2:]]) + np.savetxt(fg_spec_file, dls_fg.T) + # load a filter transfer function to smooth by fl = np.loadtxt("maps_example/transfer_example.txt") +mask = None + + +def read_map(filename, fill=None): + fields = 0 if "mask" in filename else (0, 1, 2) + data = np.asarray(hp.read_map(filename, field=fields)) + if fill is not None: + data[hp.mask_bad(data)] = fill + return data + + +def write_map(filename, data, mask=None): + if mask is not None: + data[..., ~mask] = hp.UNSEEN + hp.write_map(filename, data, partial=True, overwrite=True) + + +mask = None +mask_write = None # make celestial rectangular coordinate mask -latrange = [-55, -15] -lonrange = [18, 80] -theta, phi = hp.pix2ang(nside, np.arange(hp.nside2npix(nside))) -lat = 90 - np.rad2deg(theta) -lon = np.rad2deg(phi) -lon[lon > 180] -= 360 -mask = np.logical_and( - np.logical_and(lon > lonrange[0], lon < lonrange[1]), - np.logical_and(lat > latrange[0], lat < latrange[1]), -) -mask_write = mask.copy().astype(float) -mask_write[~mask] = hp.UNSEEN # to reduce size of file on disk -# Because using polarization, need I/Q/U mask -mask_write = np.tile(mask_write, [3, 1]) for tag in tags: - hp.write_map( - os.path.join(mask_path, "mask_map_{}.fits".format(tag)), - mask_write, - partial=True, - overwrite=True, - ) + mask_file = os.path.join(mask_path, "mask_map_{}.fits".format(tag)) + if os.path.exists(mask_file): + mask = read_map(mask_file, fill=0).astype(bool) + break + + if mask is None: + latrange = [-55, -15] + lonrange = [18, 80] + theta, phi = hp.pix2ang(nside, np.arange(hp.nside2npix(nside))) + lat = 90 - np.rad2deg(theta) + lon = np.rad2deg(phi) + lon[lon > 180] -= 360 + mask = np.logical_and( + np.logical_and(lon > lonrange[0], lon < lonrange[1]), + np.logical_and(lat > latrange[0], lat < latrange[1]), + ) + mask_write = mask.copy().astype(float) + mask_write[~mask] = hp.UNSEEN # to reduce size of file on disk + # Because using polarization, need I/Q/U mask + mask_write = np.tile(mask_write, [3, 1]) + + write_map(mask_file, mask_write) print("Wrote mask map {}".format(tag)) -# make data, signal, noise sims -for i in range(nsim + 1): - sig = hp.synfast(cls, nside=nside, new=True, pixwin=True) - noise = np.zeros_like(sig) +sig_cache = {} +fg_cache = {} + - # smooth both by filter transfer function - sig = hp.smoothing(sig, beam_window=np.sqrt(fl)) +def sim_signal(isim, itag): + if isim not in sig_cache: + np.random.seed(isim) + sig = hp.synfast(cls, nside=nside, new=True, pixwin=True) + sig = hp.smoothing(sig, beam_window=np.sqrt(fl)) + sig_cache[isim] = np.asarray(sig) + return hp.smoothing(sig_cache[isim], np.deg2rad(fwhms[itag] / 60.0), verbose=False) + + +def sim_noise(isim, itag): + np.random.seed(nsim * (1 + itag) + isim) + noise = np.tile((3, 1), np.zeros(mask.size, dtype=float)) + noise[0][mask] = np.random.normal(scale=gnoise[itag], size=np.sum(mask)) + noise[1][mask] = np.random.normal( + scale=gnoise[itag] * np.sqrt(2), size=np.sum(mask) # pol + ) + noise[2][mask] = np.random.normal( + scale=gnoise[itag] * np.sqrt(2), size=np.sum(mask) # pol + ) + return hp.smoothing(noise, beam_window=np.sqrt(fl)) + + +def sim_fg(isim, itag): + if isim not in fg_cache: + np.random.seed(nsim * (1 + len(tags)) + isim) + fg = hp.synfast(cls_fg, nside=nside, new=True, pixwin=True) + fg = hp.smoothing(fg, beam_window=np.sqrt(fl)) + fg_cache[isim] = np.asarray(fg) + return scales[itag] * hp.smoothing( + fg_cache[isim], np.deg2rad(fwhms[itag] / 60.0), verbose=False + ) + + +# make data, signal, noise, fg sims +for i in range(nsim + 1): for ti, tag in enumerate(tags): - sig_smooth = hp.smoothing(sig, np.deg2rad(fwhms[ti] / 60.0), verbose=False) - noise[0][mask] = np.random.normal(scale=gnoise[ti], size=np.sum(mask)) - noise[1][mask] = np.random.normal( - scale=gnoise[ti] * np.sqrt(2), size=np.sum(mask) # pol - ) - noise[2][mask] = np.random.normal( - scale=gnoise[ti] * np.sqrt(2), size=np.sum(mask) # pol - ) - noise_smooth = hp.smoothing(noise, beam_window=np.sqrt(fl)) if i == 0: - # call this data - dat = sig_smooth + data_noise_scale * noise_smooth - dat[:, ~mask] = hp.UNSEEN - hp.write_map( - os.path.join(data_path, "map_{}.fits".format(tag)), - dat, - partial=True, - overwrite=True, - ) - print("Wrote data map {}".format(tag)) + data_file = os.path.join(data_path, "map_{}.fits".format(tag)) + data_fg_file = os.path.join(data_fg_path, "map_{}.fits".format(tag)) + data = None + if not os.path.exists(data_file): + data = sim_sig(i, ti) + data_noise_scale * sim_noise(i, ti) + write_map(data_file, data, mask) + print("Wrote data signal+noise map {}".format(tag)) + if do_fg and not os.path.exists(data_fg_file): + if data is None: + data = read_map(data_file) + write_map(data_fg_file, data + sim_fg(i, ti), mask) + print("Wrote data signal+noise+fg map {}".format(tag)) else: - # signal and noise - sig_smooth[:, ~mask] = hp.UNSEEN - noise_smooth[:, ~mask] = hp.UNSEEN - hp.write_map( - os.path.join(sig_path, "map_{}_{:04}.fits".format(tag, i - 1)), - sig_smooth, - partial=True, - overwrite=True, - ) - hp.write_map( - os.path.join(noise_path, "map_{}_{:04}.fits".format(tag, i - 1)), - noise_smooth, - partial=True, - overwrite=True, + sig_file = os.path.join(sig_path, "map_{}_{:04}.fits".format(tag, i - 1)) + noise_file = os.path.join( + noise_path, "map_{}_{:04}.fits".format(tag, i - 1) ) - print("Wrote signal and noise map {} {} / {}".format(tag, i - 1, nsim)) + fg_file = os.path.join(fg_path, "map_{}_{:04}.fits".format(tag, i - 1)) + + comps = [] + if not os.path.exists(sig_file): + write_map(sig_file, sim_sig(i, ti), mask) + comps += ["sig"] + if not os.path.exists(noise_file): + write_map(noise_file, sim_noise(i, ti), mask) + comps += ["noise"] + if do_fg and not os.path.exists(fg_file): + write_map(fg_file, sim_fg(i, ti), mask) + comps += ["fg"] + if len(comps): + print( + "Wrote {} map {} {} / {}".format(", ".join(comps), tag, i - 1, nsim) + ) + + sig_cache.pop(i, None) + fg_cache.pop(i, None) From 5ffc0c020e9ec0e534d32eefadc55d519579f6b3 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Tue, 12 Oct 2021 22:53:03 +1300 Subject: [PATCH 07/29] overhaul foreground fitting code * consistent foreground_type and foreground_type_transfer options for handling foreground components in covariance model and transfer function calculation * fix_fg_transfer option fixes transfer function to be the same as CMB, interpolated onto the foreground binning * beta_fit option enables or disables fitting for delta_beta in addition to foreground amplitudes * correct calculation of window functions for foreground bins --- xfaster/parse_tools.py | 24 +- xfaster/xfaster_class.py | 1471 +++++++++++++++++--------------------- xfaster/xfaster_exec.py | 227 +++--- 3 files changed, 809 insertions(+), 913 deletions(-) diff --git a/xfaster/parse_tools.py b/xfaster/parse_tools.py index 8717146b..a6ffe08b 100644 --- a/xfaster/parse_tools.py +++ b/xfaster/parse_tools.py @@ -299,6 +299,9 @@ def load_and_parse(filename, check_version=True): data.pop("raw_root") data.pop("raw_files") + # In version 1, these parameters referred to ensembles used for + # sim_index runs. In version 2+, these parameters refer to ensembles + # used for the foreground component of the covariance model if "foreground_type" in data: data["foreground_type_sim"] = data.pop("foreground_type") if "foreground_root" in data: @@ -391,6 +394,22 @@ def load_and_parse(filename, check_version=True): if "fix_bb_xfer" in data: data["fix_bb_transfer"] = data.pop("fix_bb_xfer") + if version >= 1: + + if "ref_freq" in data: + data["freq_ref"] = data.pop("ref_freq") + + if "foreground_type_sim" in data and "foreground_type" not in data: + for k in ["type", "root", "files", "root2", "files2"]: + data["foreground_" + k] = None + data["foreground_transfer_" + k] = None + data["num_foreground"] = 0 + data["num_foreground_transfer"] = 0 + data["foreground_subset"] = "*" + + if "cls_sim" in data and "cls_fg" not in data: + data["cls_fg"] = None + # update data version in memory data["data_version"] = version = dv @@ -655,13 +674,14 @@ def dict_to_dsdqb_mat(dsdqb_dict, bin_def): in the Fisher iteration. """ # get the unique map tags in order from the keys map1:map2 - mtags = [x.split(":")[0] for x in dsdqb_dict["cmb"]] + mkeys = list(list(dsdqb_dict.values())[0].keys()) + mtags = [x.split(":")[0] for x in mkeys] _, uind = np.unique(mtags, return_index=True) map_tags = np.asarray(mtags)[sorted(uind)] map_pairs = tag_pairs(map_tags, index=True) nmaps = len(map_tags) - pol_dim = 3 if "cmb_ee" in bin_def else 1 + pol_dim = 3 if any(["ee" in x.split("_")[1] for x in bin_def]) else 1 inds = spec_index() bin_index = dict_to_index(bin_def) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 63a7f366..f9e2c140 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -572,11 +572,12 @@ def _get_files( mask_type="hitsmask_tailored", signal_type="r0p03", signal_transfer_type=None, - signal_transfer_type_dust=None, signal_subset="*", - signal_subset_dust="None", noise_type="stationary", noise_subset="*", + foreground_type=None, + foreground_transfer_type=None, + foreground_subset="*", signal_type_sim=None, noise_type_sim=None, foreground_type_sim=None, @@ -624,10 +625,7 @@ def _get_files( map_pairs_orig = pt.tag_pairs(map_tags, index=map_tags_orig) # make a dictionary of map freqs for each unique map tag - map_freqs_dict = {} - for im0, m0 in enumerate(map_tags): - map_freqs_dict[m0] = map_freqs[im0] - map_freqs = map_freqs_dict + map_freqs = dict(zip(map_tags, map_freqs)) # find all corresponding masks if mask_type is None: @@ -680,26 +678,29 @@ def find_sim_files(name, root=None, subset="*"): find_sim_files("noise", "noise_{}".format(noise_type), noise_subset) else: find_sim_files("noise", None) + if foreground_type is not None: + find_sim_files( + "foreground", "foreground_{}".format(foreground_type), foreground_subset + ) + else: + find_sim_files("foreground", None) # signal ensembles for transfer function covariance model if signal_transfer_type is None: signal_transfer_type = signal_type find_sim_files( - "signal_transfer", "signal_{}".format(signal_transfer_type), - signal_subset + "signal_transfer", "signal_{}".format(signal_transfer_type), signal_subset ) - - if signal_transfer_type_dust is None: - find_sim_files("signal_transfer_dust", None) - else: - if signal_subset_dust is None: - signal_subset_dust = signal_subset # better choice ?? - # find reobserved scaled dust - # TODO: implement options for unscaled dust, start with signal_ - find_sim_files( - "signal_transfer_dust", "foreground_{}".format(signal_transfer_type_dust), - signal_subset_dust + if foreground_transfer_type is None: + foreground_transfer_type = foreground_type + if foreground_transfer_type is not None: + find_sim_files( + "foreground_transfer", + "foreground_{}".format(foreground_transfer_type), + foreground_subset, ) + else: + find_sim_files("foreground_transfer", None) # apply suffix out = dict() @@ -713,6 +714,7 @@ def get_files( data_subset="full/*0", signal_subset="*", noise_subset="*", + foreground_subset="*", data_type="raw", noise_type="gaussian", noise_type_sim=None, @@ -720,11 +722,11 @@ def get_files( signal_type="synfast", signal_type_sim=None, signal_transfer_type=None, - signal_transfer_type_dust=None, - signal_subset_dust=None, + foreground_type=None, + foreground_transfer_type=None, + foreground_type_sim=None, data_root2=None, data_subset2=None, - foreground_type_sim=None, template_type=None, template_noise_type=None, template_type_sim=None, @@ -795,6 +797,12 @@ def get_files( that is added onto the data_subset path to indicate which sims to use. For example, for all, use ``'*'``. For the first 300 sims, use ``'0[0-2]*'``. + foreground_subset : string + Subset of map tags to use for spectrum estimation for foreground + sims. This should be a string that is parseable using ``glob`` + that is added onto the data_subset path to indicate which sims + to use. For example, for all, use ``'*'``. For the first 300 sims, + use ``'0[0-2]*'``. data_type : string The type of data to use, default: "raw" noise_type: string @@ -826,23 +834,26 @@ def get_files( to generate it, e.g 'synfast'. This directory may also contain a copy of the input spectrum, to make sure that the correct spectrum is used to compute the transfer function. - signal_transfer_type_dust : string - The variant of signal simulation to use for transfer function - calculation for dust, typically identified by the input spectrum model - used to generate it, e.g 'power_law'. These are assumed to be at 353 - GHz and thus scaled by the map frequency's alpha. - This directory may also contain a copy of the input spectrum, to make - sure that the correct spectrum is used to compute the transfer function. + foreground_type : string + The variant of foreground simulation to use, typically identified by + the input spectrum model used to generate it, e.g 'power_law'. + foreground_transfer_type : string + The variant of foreground simulation to use for transfer function + calculation, typically identified by the input spectrum model used + to generate it, e.g 'power_law'. These are assumed to be at 353 GHz. + This directory may also contain a copy of the input spectrum, to + make sure that the correct spectrum is used to compute the transfer + function. + foreground_type_sim : string + Tag for directory (foreground_) where + foreground sims are that should be added to the signal and noise + sims when running in sim_index mode. data_root2, data_subset2 : string The root and subset for a second set of data. If either of these is keywords is supplied, then the two data sets are treated as two halves of a null test. In this case, XFaster computes the sum and difference spectra for each map tag in order to estimate a null spectrum. - foreground_type_sim : string - Tag for directory (foreground_) where - foreground sims are that should be added to the signal and noise - sims when running in sim_index mode. template_type : string Tag for directory (templates_) containing templates (e.g. a foreground model) to be scaled by a scalar value per map tag @@ -881,8 +892,14 @@ def get_files( if signal_transfer_type is None: signal_transfer_type = signal_type + if foreground_transfer_type is None: + foreground_transfer_type = foreground_type + # make sure sims get rerun correctly - if signal_transfer_type == signal_type: + if ( + signal_transfer_type == signal_type + and foreground_transfer_type == foreground_type + ): # if signal types match, then sims are run before computing the # transfer function, so need to set the correct checkpoint to rerun if self.checkpoint == "sims": @@ -894,6 +911,9 @@ def get_files( if noise_type_sim is None: noise_type_sim = noise_type + if foreground_type_sim is None: + foreground_type_sim = foreground_type + if template_noise_type is None: template_noise_type = template_type @@ -922,11 +942,12 @@ def get_files( noise_type=noise_type, mask_type=mask_type, signal_type=signal_type, - signal_transfer_type=signal_transfer_type, - signal_transfer_type_dust=signal_transfer_type_dust, + signal_transfer_type=signal_transfer_type, + foreground_type=foreground_type, + foreground_transfer_type=foreground_transfer_type, signal_subset=signal_subset, noise_subset=noise_subset, - signal_subset_dust=signal_subset_dust, + foreground_subset=foreground_subset, ) ref_opts = dict(data_subset=data_subset, **opts) if null_run: @@ -1866,12 +1887,7 @@ def alm2cl(self, m1, m2=None, lmin=2, lmax=None, symmetric=True): cls[..., :lmin] = 0 return np.atleast_2d(cls) - def get_mask_weights( - self, - apply_gcorr=False, - reload_gcorr=False, - gcorr_file=None - ): + def get_mask_weights(self, apply_gcorr=False, reload_gcorr=False, gcorr_file=None): """ Compute cross spectra of the masks for each data map. @@ -2074,7 +2090,7 @@ def process_index(idx): gmat[xname] = OrderedDict() for spec in self.specs: si, sj = spec_inds[spec] - f = (fsky_eff[xname][si, sj] + fsky_eff[xname][sj, si]) / 2.0 + f = (fsky_eff[xname][si, sj] + fsky_eff[xname][sj, si]) / 2.0 gmat[xname][spec] = f if np.any(np.asarray([f for f in fsky.values()]) > 0.1): @@ -2712,12 +2728,7 @@ def process_index(idx): save_name, from_attrs=save_attrs, extra_tag="xcorr", data_opts=True ) - def get_masked_sims( - self, - transfer=False, - do_noise=True, - qb_file=None, - ): + def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): """ Compute average signal and noise spectra for a given ensemble of sim maps. The same procedure that is used for computing data cross spectra @@ -2749,8 +2760,8 @@ def get_masked_sims( cls_signal, cls_signal_null: Mean signal spectra - cls_signal_dust: - Mean dust signal spectra + cls_fg: + Mean foreground spectra cls_noise, cls_noise_null: Mean noise spectra cls_sim, cls_sim_null: @@ -2779,25 +2790,26 @@ def get_masked_sims( signal_files = self.signal_transfer_files signal_files2 = self.signal_transfer_files2 if null_run else None num_signal = self.num_signal_transfer - signal_files_dust = self.signal_transfer_dust_files - num_signal_dust = self.num_signal_transfer_dust + foreground_files = self.foreground_transfer_files + num_foreground = self.num_foreground_transfer else: signal_files = self.signal_files signal_files2 = self.signal_files2 if null_run else None - num_signal = self.num_signal - num_signal_dust = 0 - + num_signal = self.num_signal + foreground_files = self.foreground_files + num_foreground = self.num_foreground noise_files = self.noise_files noise_files2 = self.noise_files2 if null_run else None num_noise = self.num_noise + do_fg = num_foreground > 0 and not null_run if do_noise: do_noise = noise_files is not None save_attrs = [ "cls_signal", - "cls_signal_dust", + "cls_fg", "cls_noise", "cls_sim", "cls_med", @@ -2814,11 +2826,13 @@ def get_masked_sims( if transfer: save_name = "sims_xcorr_{}".format(self.signal_transfer_type) - if self.signal_transfer_type_dust is not None: - save_name += "_" + self.signal_transfer_type_dust + if self.foreground_transfer_type: + save_name += "_" + self.foreground_transfer_type cp = "sims_transfer" else: save_name = "sims_xcorr_{}".format(self.signal_type) + if self.foreground_type: + save_name += "_" + self.foreground_type cp = "sims" ret = self.load_data( @@ -2847,9 +2861,9 @@ def get_masked_sims( if transfer and self.signal_transfer_type == self.signal_type: self.force_rerun["sims"] = False - # process signal, signal dust, noise, and S+N + # process signal, foreground, noise, and S+N cls_sig = OrderedDict() - cls_sig_dust = OrderedDict() #if transfer_dust else None + cls_fg = OrderedDict() if do_fg else None cls_null_sig = OrderedDict() if null_run else None cls_noise = OrderedDict() if do_noise else None cls_null_noise = OrderedDict() if null_run and do_noise else None @@ -2862,17 +2876,16 @@ def get_masked_sims( cls_res = OrderedDict() if do_noise else None cls_null_res = OrderedDict() if null_run and do_noise else None if do_noise: - #for k in ["nxn0", "nxn1", "sxn0", "sxn1", "nxs0", "nxs1", "cls_signal_dust"]: - for k in ["nxn0", "nxn1", "sxn0", "sxn1", "nxs0", "nxs1"]: # cls_signal_dust seperate + for k in ["nxn0", "nxn1", "sxn0", "sxn1", "nxs0", "nxs1"]: cls_res[k] = OrderedDict() if null_run: cls_null_res[k] = OrderedDict() if num_noise != 0: - nsim_min = min([num_signal, num_signal_dust, num_noise]) + nsim_min = min([num_signal, num_noise]) else: - nsim_min = min([num_signal, num_signal_dust]) # not sure why we need if case here? - nsim_max = max([num_signal, num_signal_dust, num_noise]) + nsim_min = num_signal + nsim_max = max([num_signal, num_foreground, num_noise]) cls_all = np.zeros( [nsim_max, len(map_pairs.items()), len(self.specs), self.lmax + 1] ) @@ -2927,7 +2940,7 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): m_alms = get_map_alms(filename, mask, rls=rls) mn_alms = None - if null_run: + if null_run and files2 is not None: # second null half filename2 = files2[idx] if idx2 is not None: @@ -2940,12 +2953,12 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): return cache[idx] sig_cache = dict() - sig_dust_cache = dict() + fg_cache = dict() noise_cache = dict() for isim in range(nsim_max): sig_cache.clear() - sig_dust_cache.clear() + fg_cache.clear() noise_cache.clear() for xind, (xname, (idx, jdx)) in enumerate(map_pairs.items()): self.log( @@ -2971,14 +2984,14 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): if null_run: cls_null1t = np.copy(cls_null1_sig) - if isim < num_signal_dust: - dimap_alms, _ = process_index( - signal_files_dust, None, idx, isim, sig_dust_cache - ) - djmap_alms, _ = process_index( - signal_files_dust, None, jdx, isim, sig_dust_cache - ) - cls1_sig_dust = self.alm2cl(dimap_alms, djmap_alms) + if do_fg and isim < num_foreground: + fimap_alms, _ = process_index( + foreground_files, None, idx, isim, fg_cache + ) + fjmap_alms, _ = process_index( + foreground_files, None, jdx, isim, fg_cache + ) + cls1_fg = self.alm2cl(fimap_alms, fjmap_alms) if do_noise and isim < num_noise: cls1_res = OrderedDict() @@ -3004,10 +3017,10 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): # need non-symmetric since will potentially modify these # with different residuals for T, E, B - for k, cx, cy, cnx, cny in [ + for k, cx, cy, cnx, cny in [ ("nxn", nimap_alms, njmap_alms, ninull_alms, njnull_alms), ("sxn", simap_alms, njmap_alms, sinull_alms, njnull_alms), - ("nxs", nimap_alms, sjmap_alms, ninull_alms, sjnull_alms), + ("nxs", nimap_alms, sjmap_alms, ninull_alms, sjnull_alms), ]: k0, k1 = ["{}{}".format(k, i) for i in range(2)] cls1_res[k0] = self.alm2cl(cx, cy, symmetric=False) @@ -3016,7 +3029,7 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): cls_null1_res[k0] = self.alm2cl(cnx, cny, symmetric=False) cls_null1_res[k1] = self.alm2cl(cny, cnx, symmetric=False) - # construct noise model + # construct noise model cls1_noise = (cls1_res["nxn0"] + cls1_res["nxn1"]) / 2.0 if null_run: cls_null1_noise = ( @@ -3029,23 +3042,14 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): cls1t += cls1_res[k] / 2.0 if null_run: cls_null1t += cls_null1_res[k] / 2.0 - - # add dust to cls1t - if isim < num_signal_dust: - cls1t += cls1_sig_dust - # cls1_res["cls_signal_dust"] = cls1_sig_dust # not need, see quants - if null_run: - cls_null1t += cls1_sig_dust - #cls_null1_res["cls_signal_dust"] = cls1_sig_dust - quants = [] if isim < num_signal: quants += [[cls_sig, cls1_sig]] if null_run: quants += [[cls_null_sig, cls_null1_sig]] - if isim < num_signal_dust: - quants += [[cls_sig_dust, cls1_sig_dust]] + if do_fg and isim < num_foreground: + quants += [[cls_fg, cls1_fg]] if do_noise and isim < num_noise: quants += [[cls_noise, cls1_noise]] @@ -3090,8 +3094,7 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): cls_null_med[spec][xname] = cls_null_med_arr[xind][s] self.cls_signal = cls_sig - if transfer: - self.cls_signal_dust = cls_sig_dust + self.cls_fg = cls_fg self.cls_signal_null = cls_null_sig self.cls_noise = cls_noise self.cls_noise_null = cls_null_noise @@ -3337,13 +3340,12 @@ def get_signal_shape( r=None, component=None, flat=None, + filename_fg=None, + freq_ref=359.7, + beta_ref=1.54, signal_mask=None, transfer=False, - transfer_dust=False, - #separate_dust_Fl=False, save=True, - template_alpha90=None, - template_alpha150=None, ): """ Load a shape spectrum for input to the Fisher iteration algorithm. @@ -3361,98 +3363,95 @@ def get_signal_shape( Arguments --------- - filename : string, or list of strings - Filename for a spectrum on disk. If None, and ``r`` is None and - ``flat`` is False, this will search for a spectrum stored in - ``signal_/spec_signal_.dat``. - Otherwise, if the filename is a relative path and not found, - the config directory is searched. - If a list, first element is CMB and second element is dust. + filename : string + Filename for a signal spectrum on disk. If None, and ``r`` is None + and ``flat`` is False, this will search for a spectrum stored in + ``signal_/spec_signal_.dat``. Otherwise, + if the filename is a relative path and not found, the config + directory is searched. r : float If supplied and ``flat`` is False, a spectrum is computed using CAMB for the given ``r`` value. Overrides ``filename``. component : 'scalar', 'tensor', 'fg' If 'scalar', and ``r`` is not None, return just the r=0 scalar terms in the signal model. If 'tensor', return just the tensor component - scaled by the input ``r`` value. If 'fg', return just fg term + scaled by the input ``r`` value. If 'fg', return just fg terms flat : float If given, a spectrum that is flat in ell^2 Cl is returned, with amplitude given by the supplied value. Overrides all other options. + filename_fg : string + Filename for a foreground spectrum on disk. If None, this will + search for a spectrum stored in + ``foreground_/spec_foreground_.dat``. + Otherwise, if the filename is a relative path and not found, the + config directory is searched. + freq_ref : float + In GHz, reference frequency for dust model. Dust bandpowers output + will be at this reference frequency. + beta_ref : float + The spectral index of the dust model. This is a fixed value, with + an additive deviation from this value fit for in foreground fitting + mode. signal_mask: str array Include only these spectra, others set to zero. Options: TT, EE, BB, TE, EB, TB transfer : bool If True, this is a transfer function run. If ``filename`` is None - and ``r`` is None and ``flat`` is False, will search for a spectrum - stored in + and ``r`` is None and ``flat`` is False, will search for a signal + spectrum stored in ``signal_/spec_signal_.dat``. - transfer_dust : bool - If True, this is a transfer function run for dust and dust+cmb. - If `filename` is None, will search for a spectrum stored in - `signal_/spec_signal_.dat`. + Also, if ``filename_fg`` is None, will search for a foreground + spectrum stored in + ``foreground_/spec_foreground_.dat``. save : bool If True, save signal shape dict to disk. Returns ------- cls : OrderedDict - Dictionary keyed by spectrum (cmb_tt, cmb_ee, ... , fg), each - entry containing a vector of length 2 * lmax + 1 + Dictionary keyed by spectrum (cmb_tt, cmb_ee, ... , fg_tt, ...), + each entry containing a vector of length 2 * lmax + 1. """ - if isinstance(filename, list): - filename_dust = filename[1] - filename = filename[0] - lmax_kern = 2 * self.lmax specs = list(self.specs) + comps = [] + if "cmb_tt" in self.bin_def or "cmb_ee" in self.bin_def: + comps += ["cmb"] + if not self.null_run and "fg_tt" in self.bin_def: + comps += ["fg"] + if component in ["cmb", "fg"]: + comps = [component] + if transfer: if "eb" in specs: specs.remove("eb") if "tb" in specs: specs.remove("tb") tbeb = "eb" in specs and "tb" in specs and r is None - specs = ["cmb_{}".format(spec) for spec in specs] - - if not self.null_run and "fg_tt" in self.bin_def: - specs.append("fg") - - if component == "fg": - specs = ["fg"] - - if transfer_dust: - specs += ["dust_tt"] - if self.pol: - for s in ["ee", "bb", "te"]: - specs += ["dust_{}".format(s)] - # CHECK: DO I NEED THIS FOR SEPARATE DUST Fl ??? - # Need separate fields per freq for dust+CMB to get relative - # amplitudes right - for freq in np.unique(list(self.nom_freqs.values())): - specs += ["cmbdust_tt_{}".format(freq)] - if self.pol: - for s in ["ee", "bb", "te"]: - specs += ["cmbdust_{}_{}".format(s, freq)] - - nspecs = len(specs) if save: - shape = (nspecs, lmax_kern + 1) + shape = (len(specs) * len(comps), lmax_kern + 1) save_name = "shape_transfer" if transfer else "shape" - opts = dict( - filename=filename, - r=r, - flat=flat, - signal_mask=signal_mask, - ) + opts = dict(signal_mask=signal_mask) + if "cmb" in comps: + opts.update(filename=filename, r=r, flat=flat) + if "fg" in comps: + opts.update( + filename_fg=filename_fg, freq_ref=freq_ref, beta_ref=beta_ref + ) ret = self.load_data( save_name, save_name, shape_ref="cls_shape", shape=shape, value_ref=opts ) if ret is not None: - if r is not None: + if "cmb" in comps and r is not None: self.r_model = ret["r_model"] + if "fg" in comps: + self.freq_ref = ret["freq_ref"] + self.beta_ref = ret["beta_ref"] + self.fg_scales = ret["fg_scales"] self.cls_shape = ret["cls_shape"] return ret["cls_shape"] @@ -3460,173 +3459,121 @@ def get_signal_shape( ellfac = ell * (ell + 1) / 2.0 / np.pi cls_shape = OrderedDict() - def fix_camb_specs(arr, ltmp, filename): - # camb starts at l=2, so set first 2 ells to be 0 - # make sure order is TT, EE, BB, TE, no ell row. - new_arr = np.append([0, 0], arr[1, : ltmp - 2]) - if self.pol: - if np.any(arr[2, : ltmp - 2] < 0): - # this is true if TE is the third index instead of EE - # # (ell is 0th index for CAMB) - self.log( - "Old CAMB format in model file {}. Re-indexing.".format( - filename - ), - "detail", - ) - pol_specs = [3, 4, 2] + if "cmb" in comps: + if flat is not None and flat is not False: + if flat is True: + flat = 2e-5 + # flat spectrum for null tests + for spec in specs: + cls_shape["cmb_" + spec] = flat * np.ones_like(ell) + + elif r is not None: + # cache model components + if not hasattr(self, "r_model") or self.r_model is None: + # scalar CAMB spectrum + scal = st.get_camb_cl(r=0, lmax=lmax_kern) + # tensor CAMB spectrum for r=1, scales linearly with r + tens = st.get_camb_cl(r=1, lmax=lmax_kern, nt=0, spec="tensor") + self.r_model = {"scalar": scal, "tensor": tens} + if save: + opts["r_model"] = self.r_model + # CAMB spectrum for given r value + component = str(component).lower() + if component == "scalar": + cls_camb = self.r_model["scalar"] + elif component == "tensor": + cls_camb = r * self.r_model["tensor"] else: - pol_specs = [2, 3, 4] - for d in arr[pol_specs]: - new_arr = np.vstack([new_arr, np.append([0, 0], d[: ltmp - 2])]) - return new_arr - - if flat is not None and flat is not False: - if flat is True: - flat = 2e-5 - # flat spectrum for null tests - for spec in specs: - cls_shape[spec] = flat * np.ones_like(ell) - - elif r is not None: - # cache model components - if not hasattr(self, "r_model") or self.r_model is None: - # scalar CAMB spectrum - scal = st.get_camb_cl(r=0, lmax=lmax_kern) - # tensor CAMB spectrum for r=1, scales linearly with r - tens = st.get_camb_cl(r=1, lmax=lmax_kern, nt=0, spec="tensor") - self.r_model = {"scalar": scal, "tensor": tens} - if save: - opts["r_model"] = self.r_model - # CAMB spectrum for given r value - component = str(component).lower() - if component == "scalar": - cls_camb = self.r_model["scalar"] - elif component == "tensor": - cls_camb = r * self.r_model["tensor"] - else: - cls_camb = self.r_model["scalar"] + r * self.r_model["tensor"] - ns, _ = cls_camb.shape - for s, spec in enumerate(specs[:ns]): - cls_shape[spec] = cls_camb[s] - - elif any([x.startswith("cmb") for x in specs]): - # signal sim model or custom filename - if filename is None: - signal_root = ( - self.signal_transfer_root if transfer else self.signal_root - ) - filename = "spec_{}.dat".format(os.path.basename(signal_root)) - filename = os.path.join(signal_root, filename) - elif not os.path.exists(filename): - filename = os.path.join(self.config_root, filename) - if not os.path.exists(filename): - raise OSError("Missing model file {}".format(filename)) - - tmp = np.loadtxt(filename, unpack=True) - - ltmp = tmp.shape[-1] - if lmax_kern + 1 < ltmp: - ltmp = lmax_kern + 1 + cls_camb = self.r_model["scalar"] + r * self.r_model["tensor"] + cls_shape.update({"cmb_" + s: cls for s, cls in zip(specs, cls_camb)}) + else: - raise ValueError( - "Require at least lmax={} in model file, found {}".format( - lmax_kern, ltmp + # signal sim model or custom filename + if filename is None: + signal_root = ( + self.signal_transfer_root if transfer else self.signal_root ) - ) + filename = "spec_{}.dat".format(os.path.basename(signal_root)) + filename = os.path.join(signal_root, filename) + elif not os.path.exists(filename): + filename = os.path.join(self.config_root, filename) + if not os.path.exists(filename): + raise OSError("Missing model file {}".format(filename)) - tmp = fix_camb_specs(tmp, ltmp, filename) - - if transfer_dust: - if filename_dust is None: - signal_root = self.signal_transfer_dust_root #self.signal_transfer_root_dust - filename_dust = "spec_{}.dat".format(os.path.basename(signal_root)) - filename_dust = os.path.join(signal_root, filename_dust) - if not os.path.exists(filename_dust) and not os.path.isabs( - filename_dust - ): - filename_dust = os.path.abspath( - os.path.join( - os.path.dirname(__file__), "../data", filename_dust - ) - ) - if not os.path.exists(filename_dust): - raise OSError("Missing model file {}".format(filename_dust)) - tmpd = np.loadtxt(filename_dust, unpack=True) + cls_cmb = st.load_camb_cl(filename, lmax=lmax_kern, pol=self.pol) + cls_shape.update({"cmb_" + s: cls for s, cls in zip(specs, cls_cmb)}) - ltmpd = tmpd.shape[-1] - if lmax_kern + 1 < ltmpd: - ltmpd = lmax_kern + 1 - else: - raise ValueError( - "Require at least lmax={} in model file, found {}".format( - lmax_kern, ltmpd - ) + # EB and TB flat l^2 * C_l + if self.pol: + if tbeb and (flat is None or flat is False): + tbeb_flat = np.abs(cls_shape["cmb_bb"][100]) * ellfac[100] * 1e-4 + tbeb_flat = np.ones_like(cls_shape["cmb_bb"]) * tbeb_flat + tbeb_flat[:2] = 0 + cls_shape["cmb_eb"] = np.copy(tbeb_flat) + cls_shape["cmb_tb"] = np.copy(tbeb_flat) + elif not tbeb and not transfer: + cls_shape["cmb_eb"] = np.zeros_like(ell, dtype=float) + cls_shape["cmb_tb"] = np.zeros_like(ell, dtype=float) + + if "fg" in comps: + # dust signal sim model or custom filename + # XXX optionally load freq_ref and beta_ref from file + if filename_fg is None: + fg_root = ( + self.foreground_transfer_root if transfer else self.foreground_root + ) + if fg_root is not None: + filename_fg = "spec_{}.dat".format(os.path.basename(fg_root)) + filename_fg = os.path.join(fg_root, filename_fg) + elif not os.path.exists(filename_fg): + filename_fg = os.path.join(self.config_root, filename_fg) + + if filename_fg is None: + # From Planck LIV EE dust + cls_dust = st.dust_model(ell, lfac=True) + for spec in specs: + cls_shape["fg_" + spec] = cls_dust + self.log( + "Added simple foregrounds to cls shape {}".format(list(cls_shape)), + "debug", + ) + else: + if not os.path.exists(filename_fg): + raise OSError( + "Missing foreground model file {}".format(filename_fg) ) - tmpd = fix_camb_specs(tmpd, ltmpd, filename_dust) - - - # camb starts at l=2, so set first 2 ells to be 0 - #cls_shape["cmb_tt"] = np.append([0, 0], tmp[1, : ltmp - 2]) - #if self.pol: - # if np.any(tmp[2, : ltmp - 2] < 0): - # this is true if TE is the third index instead of EE - # (ell is 0th index for CAMB) - # self.warn( - # "Old CAMB format in model file {}. Re-indexing.".format( - # filename - # ) - # ) - # pol_specs = [3, 4, 2] - spec_order = {"tt": 0, "ee": 1, "bb": 2, "te": 3} - for spec in specs: - parts = spec.split("_") - comp = parts[0] - spec0 = parts[1] - self.log(spec) - if len(parts) == 3: - nom_freq = parts[2] - if comp == "cmb": - arr = tmp - elif comp == "dust": - arr = tmpd - else: - # pol_specs = [2, 3, 4] - #for spec, d in zip(specs[1:4], tmp[pol_specs]): - # cls_shape[spec] = np.append([0, 0], d[: ltmp - 2]) - # Combine CMB and dust scaled by alpha**2 for frequency - # TODO: GENERALIZE THIS FOR FREQs - if nom_freq == "90": - alpha = template_alpha90 - elif nom_freq == "150": - alpha = template_alpha150 - arr = tmp + alpha ** 2 * tmpd - cls_shape[spec] = arr[spec_order[spec0]] # - - # EB and TB flat l^2 * C_l - if self.pol: - if tbeb and (flat is None or flat is False): - tbeb_flat = np.abs(cls_shape["cmb_bb"][100]) * ellfac[100] * 1e-4 - tbeb_flat = np.ones_like(cls_shape["cmb_bb"]) * tbeb_flat - tbeb_flat[:2] = 0 - cls_shape["cmb_eb"] = np.copy(tbeb_flat) - cls_shape["cmb_tb"] = np.copy(tbeb_flat) - elif not tbeb and not transfer: - cls_shape["cmb_eb"] = np.zeros_like(ell, dtype=float) - cls_shape["cmb_tb"] = np.zeros_like(ell, dtype=float) - - if "fg" in specs: - # From Planck LIV EE dust - cls_dust = 34.0 * (ell[2:] / 80.0) ** (-2.28 + 2.0) - cls_shape["fg"] = np.append([0, 0], cls_dust) - self.log( - "Added foreground to cls shape {}".format(list(cls_shape)), "debug" - ) + + cls_fg = st.load_camb_cl(filename_fg, lmax=lmax_kern, pol=None) + if len(cls_fg) == 1: + cls_fg = np.tile(cls_fg, (len(specs), 1)) + cls_shape.update({"fg_" + s: cls for s, cls in zip(specs, cls_fg)}) + + if self.pol and "fg_eb" not in cls_shape: + if tbeb: + tbeb_flat = np.abs(cls_shape["fg_ee"][100]) * 1e-4 + tbeb_flat = np.ones_like(cls_shape["fg_ee"]) * tbeb_flat + cls_shape["fg_eb"] = np.copy(tbeb_flat) + cls_shape["fg_tb"] = np.copy(tbeb_flat) + elif not transfer: + cls_shape["fb_eb"] = np.zeros_like(ell, dtype=float) + cls_shape["fb_tb"] = np.zeros_like(ell, dtype=float) + + # frequency scaling for cross spectra + fg_scales = OrderedDict() + for xname, (m0, m1) in self.map_pairs.items(): + f0 = self.map_freqs[m0] + f1 = self.map_freqs[m1] + fg_scales[xname] = st.scale_dust(f0, f1, freq_ref, beta_ref, delta=True) + self.fg_scales = fg_scales + self.freq_ref = freq_ref + self.beta_ref = beta_ref + if save: + opts["fg_scales"] = self.fg_scales # divide out l^2/2pi - for spec in specs: - cls_shape[spec][2:] /= ellfac[2:] - cls_shape[spec][:2] = 0.0 + for cl in cls_shape.values(): + cl[2:] /= ellfac[2:] + cl[:2] = 0.0 if signal_mask is not None: self.log("Masking {} spectra".format(signal_mask), "debug") @@ -3801,6 +3748,7 @@ def get_bin_def( res_specs=None, bin_width_res=25, foreground_fit=False, + beta_fit=False, bin_width_fg=25, ): """ @@ -3846,6 +3794,9 @@ def get_bin_def( foreground_fit : bool If True, construct bin definitions for foreground components as well. + beta_fit : bool + If True, include ``delta_beta`` in the foreground fitting parameters, + along with the foreground amplitudes. bin_width_fg : int or list of ints Width of each foreground spectrum bin. If a scalar, the same width is applied to all cross spectra. Otherwise, must be a list of up to @@ -3909,9 +3860,11 @@ def get_bin_def( bins = np.append(bins, self.lmax + 1) bin_def["fg_{}".format(spec)] = np.column_stack((bins[:-1], bins[1:])) nbins_fg += len(bins) - 1 - bin_def["delta_beta"] = np.array([[0, 0]]) + if beta_fit: + bin_def["delta_beta"] = np.array([[0, 0]]) self.log( - "Added {} foreground bins to bin_def".format(nbins_fg + 1), "debug" + "Added {} foreground bins to bin_def".format(nbins_fg + int(beta_fit)), + "debug", ) # Do the same for residual bins @@ -3966,7 +3919,7 @@ def get_bin_def( return self.bin_def - def kernel_precalc(self, dust=False, map_tag=None, transfer_run=False): + def kernel_precalc(self, map_tag=None, component=None, transfer=False): """ Compute the mixing kernels M_ll' = K_ll' * F_l' * B_l'^2. Called by ``bin_cl_template`` to pre-compute kernel terms. @@ -3977,18 +3930,18 @@ def kernel_precalc(self, dust=False, map_tag=None, transfer_run=False): If supplied, the kernels are computed only for the given map tag (or cross if map_tag is map_tag1:map_tag2). Otherwise, it is computed for all maps and crosses. - dust: bool - If true, get mll with dust Fl - transfer_run : bool - If True, set transfer function to 1 to solve for transfer function - qbs. + component : str + If supplied, the matrix is computed only for the given component(s). + Required if ``transfer=True``. + transfer : bool + If True, assumes a unity transfer function for all bins for the + given component(s). Returns ------- mll : OrderedDict - Dictionary of M_ll' matrices, keyed by spec and xname. + Dictionary of M_ll' matrices, keyed by component spec and xname. """ - # add one dim to add map_pairs = None if map_tag is not None: @@ -4005,33 +3958,41 @@ def kernel_precalc(self, dust=False, map_tag=None, transfer_run=False): specs = list(self.specs) lmax = self.lmax # 2 * lmax - if not transfer_run: - # collect transfer function terms - transfer = OrderedDict() - for spec in specs: - transfer[spec] = OrderedDict() - for tag in map_tags: - if dust: - transfer[spec][tag] = self.transfer_dust[spec][tag] - else: - transfer[spec][tag] = self.transfer[spec][tag] - lk = slice(0, lmax + 1) mll = OrderedDict() - for spec in specs: - mll[spec] = OrderedDict() + all_comps = ["cmb", "fg"] + if component is None: + component = all_comps + elif component not in all_comps: + raise ValueError( + "Invalid component, must be one of {}".format(", ".join(all_comps)) + ) + + for stag in self.bin_def: + comp, spec = stag.split("_", 1) + if comp not in ["cmb", "fg"]: + continue + + if transfer and comp not in component: + continue + + mll[stag] = OrderedDict() if spec in ["ee", "bb"]: - mspec = "{}_mix".format(spec) - mll[mspec] = OrderedDict() + mstag = "{}_mix".format(stag) + mll[mstag] = OrderedDict() + + beam = self.beam_windows[spec] + if not transfer: + xfer = self.transfer[stag] for xname, (m0, m1) in map_pairs.items(): # beams - fb2 = self.beam_windows[spec][m0][lk] * self.beam_windows[spec][m1][lk] + fb2 = beam[m0][lk] * beam[m1][lk] # transfer function - if not transfer_run: - fb2 *= np.sqrt(transfer[spec][m0][lk] * transfer[spec][m1][lk]) + if not transfer: + fb2 *= np.sqrt(xfer[m0][lk] * xfer[m1][lk]) # kernels if spec == "tt": @@ -4045,9 +4006,9 @@ def kernel_precalc(self, dust=False, map_tag=None, transfer_run=False): k = self.pkern[xname][:, lk] - self.mkern[xname][:, lk] # store final product - mll[spec][xname] = k * fb2 + mll[stag][xname] = k * fb2 if spec in ["ee", "bb"]: - mll[mspec][xname] = mk * fb2 + mll[mstag][xname] = mk * fb2 return mll @@ -4055,17 +4016,16 @@ def bin_cl_template( self, cls_shape=None, map_tag=None, - transfer_run=False, + component=None, + transfer=False, beam_error=False, use_precalc=True, - dust=False, - cmbdust=False, ): """ Compute the Cbl matrix from the input shape spectrum. This method requires beam windows, kernels and transfer functions - (if ``transfer_run`` is False) to have been precomputed. + (if ``transfer=False``) to have been precomputed. Arguments --------- @@ -4076,10 +4036,14 @@ def bin_cl_template( If supplied, the Cbl is computed only for the given map tag (or cross if map_tag is map_tag1:map_tag2). Otherwise, it is computed for all maps and crosses. - transfer_run : bool - If True, this assumes a unity transfer function for all bins, and - the output Cbl is used to compute the transfer functions that are - then loaded when this method is called with ``transfer_run = False``. + component : str + If supplied, the Cbl is computed only for the given component(s). + Required if ``transfer=True``. + transfer : bool + If True, assumes a unity transfer function for all bins for the + given component, and the output Cbl is used to compute the transfer + functions for that component, which are then loaded when this method + is called with ``transfer=False``. beam_error : bool If True, use beam error envelope instead of beam to get cbls that are 1 sigma beam error envelope offset of signal terms. @@ -4087,11 +4051,6 @@ def bin_cl_template( If True, load pre-calculated terms stored from a previous iteration, and store for a future iteration. Otherwise, all calculations are repeated. - dust : bool - If True, compute transfer function using dust sims - cmbdust : bool - If True, compute transfer function using CMB+dust sims, with - dust scaled by template_alpha90, template_alpha150 Returns ------- @@ -4119,30 +4078,30 @@ def bin_cl_template( map_pairs = pt.tag_pairs(map_tags) specs = list(self.specs) - if transfer_run: + if transfer: + assert component in [ + "cmb", + "fg", + ], "Argument `component` required if `transfer=True`" if "eb" in specs: specs.remove("eb") if "tb" in specs: specs.remove("tb") + if component == "fg" and "te" in specs: + specs.remove("te") lmax = self.lmax lmax_kern = lmax # 2 * self.lmax if getattr(self, "mll", None) is None or not use_precalc: - mll = self.kernel_precalc(map_tag=map_tag, transfer_run=transfer_run) + mll = self.kernel_precalc( + map_tag=map_tag, component=component, transfer=transfer + ) if use_precalc: self.mll = mll else: mll = self.mll - if not hasattr(self, "transfer_dust") and not use_precalc: - mll_dust = self.kernel_precalc(map_tag=map_tag, dust=True, - transfer_run=transfer_run) - if use_precalc: - self.mll_dust = mll_dust - else: - mll_dust = self.mll_dust - if beam_error: beam_error = self.get_beam_errors() beam_keys = ["b1", "b2", "b3"] @@ -4152,19 +4111,17 @@ def bin_cl_template( cbl = OrderedDict() comps = [] - if "cmb_tt" in cls_shape or "cmb_ee" in cls_shape: - comps += ["cmb"] - if "fg" in cls_shape and not transfer_run: - comps += ["fg"] - if transfer_run and dust: - comps += ["dust"] - if transfer_run and cmbdust: - #comps += ["dustcmb"] - comps += ["cmbdust"] - if self.nbins_res > 0 and not transfer_run: - comps += ["res"] - cls_noise = self.cls_noise_null if self.null_run else self.cls_noise - cls_res = self.cls_res_null if self.null_run else self.cls_res + if component: + comps = component if isinstance(component, list) else [component] + else: + if "cmb_tt" in cls_shape or "cmb_ee" in cls_shape: + comps += ["cmb"] + if "fg_tt" in cls_shape or "fg_ee" in cls_shape: + comps += ["fg"] + if self.nbins_res > 0: + comps += ["res"] + cls_noise = self.cls_noise_null if self.null_run else self.cls_noise + cls_res = self.cls_res_null if self.null_run else self.cls_res ell = np.arange(lmax_kern + 1) @@ -4237,11 +4194,12 @@ def binup2(d, bd, bw): continue # use correct shape spectrum + stag = "{}_{}".format(comp, spec) if comp == "fg": - # single foreground spectrum - s_arr[xi] = cls_shape["fg"][lk] + freq_scale = self.fg_scales[xname][0] + s_arr[xi] = freq_scale * cls_shape[stag][lk] else: - s_arr[xi] = cls_shape["{}_{}".format(comp, spec)][lk] + s_arr[xi] = cls_shape[stag][lk] # use correct beam error shape if beam_error: @@ -4249,11 +4207,10 @@ def binup2(d, bd, bw): b_arr["b2"][xi] = beam_error[spec][tag2][lk] # apply cross spectrum kernel terms - mll1 = mll_dust if comp == "dust" else mll - d_arr[xi, ls] = mll1[spec][xname][ls, lk] + d_arr[xi, ls] = mll[stag][xname][ls, lk] if spec in ["ee", "bb"]: - mspec = spec + "_mix" - md_arr[xi, ls] = mll1[mspec][xname][ls, lk] + mstag = stag + "_mix" + md_arr[xi, ls] = mll[mstag][xname][ls, lk] if "res" in comp: continue @@ -4295,11 +4252,7 @@ def binup2(d, bd, bw): return cbl - def get_model_spectra( - self, qb, cbl, delta=True, res=True, cls_noise=None, cond_noise=None, - transfer_dust=False, - transfer_cmbdust=False, - ): + def get_model_spectra(self, qb, cbl, cls_noise=None, cond_noise=None): """ Compute unbinned model spectra from qb amplitudes and a Cbl matrix. Requires pre-loaded bin definitions using ``get_bin_def`` or @@ -4313,11 +4266,6 @@ def get_model_spectra( Array of bandpowers for every spectrum bin. cbl : dict Cbl dict as computed by ``bin_cl_template``. - delta : bool - If True, evaluate the foreground model at the spectral - index offset by qb['delta_beta'] - res : bool - If True, include the residual noise model terms. cls_noise : OrderedDict If supplied, the noise spectrum is applied to the model spectrum. cond_noise : float @@ -4340,22 +4288,18 @@ def get_model_spectra( """ comps = [] - if transfer_dust: - comps = ["dust"] - elif transfer_cmbdust: - comps = ["cmbdust"] - elif any([k.startswith("cmb_") for k in qb]): - comps = ["cmb"] + if any([k.startswith("cmb_") for k in qb]): + comps += ["cmb"] + if any([k.startswith("fg_") for k in qb]): + comps += ["fg"] delta_beta = 0.0 - if "delta_beta" in qb and not transfer_dust: + if "delta_beta" in qb: # Evaluate fg at spectral index pivot for derivative # in Fisher matrix, unless delta is True - if delta: - delta_beta = qb["delta_beta"][0] - comps += ["fg"] + delta_beta = qb["delta_beta"][0] - if res and any([k.startswith("res_") for k in qb]): + if any([k.startswith("res_") for k in qb]): comps += ["res"] if cls_noise is not None: @@ -4373,7 +4317,6 @@ def get_model_spectra( specs = [] for spec in self.specs: - # cmb, dust, and cmbdust should have same specs (Fl cases) if "cmb_{}".format(spec) in cbl: # Don't add entries that won't be filled in later cls["total_{}".format(spec)] = OrderedDict() @@ -4402,25 +4345,17 @@ def get_model_spectra( tag1, tag2 = self.map_pairs[xname] # extract qb's for the component spectrum - if comp in ["cmb", "dust", "cmbdust"]: + if comp == "cmb" or (comp == "fg" and delta_beta == 0.0): qbs = qb[stag] if spec in ["ee", "bb"]: qbm = qb[mstag] elif comp == "fg": - # frequency scaling for foreground model - # I don't remember why delta beta was done this way. - # For likelihood, it makes sense to just use beta_ref+db - freq_scale = st.scale_dust( - self.map_freqs[tag1], - self.map_freqs[tag2], - ref_freq=self.ref_freq, - beta=self.beta_ref, - delta_beta=delta_beta, - ) - qbs = freq_scale * qb[stag] + # beta scaling for foreground model + beta_scale = self.fg_scales[xname][1] ** delta_beta + qbs = beta_scale * qb[stag] if spec in ["ee", "bb"]: - qbm = freq_scale * qb[mstag] + qbm = beta_scale * qb[mstag] elif comp == "res": # modify model by previously fit res, including @@ -4476,7 +4411,7 @@ def get_model_spectra( cl1 /= 2.0 # compute model spectra - if comp in ["cmb", "fg", "dust", "cmb_dust"]: + if comp in ["cmb", "fg"]: if xname not in cbl[stag]: continue cbl1 = cbl[stag][xname] @@ -4522,18 +4457,7 @@ def get_model_spectra( return cls - #def get_data_spectra(self, map_tag=None, transfer_run=False, do_noise=True): - def get_data_spectra( - self, - map_tag=None, - transfer_run=False, - do_noise=True, - transfer_dust=False, - transfer_cmbdust=False, - template_alpha90=None, - template_alpha150=None, - ): - + def get_data_spectra(self, map_tag=None, transfer_comp=False, do_noise=True): """ Return data and noise spectra for the given map tag(s). Data spectra and signal/noise sim spectra must have been precomputed or loaded from @@ -4544,21 +4468,23 @@ def get_data_spectra( map_tag : str If None, all map-map cross-spectra are included in the outputs. Otherwise, only the autospectra of the given map are included. - transfer_run : bool - If True, the data cls are the average of the signal simulations, and - noise cls are ignored. If False, the data cls are either - ``cls_data_null`` (for null tests) or ``cls_data``. See + transfer_comp : str + If supplied, the data cls are the average of the simulations for the + given component, and noise cls are ignored. If False, the data cls + are either ``cls_data_null`` (for null tests) or ``cls_data``. See ``get_masked_data`` for how these are computed. The input noise is similarly either ``cls_noise_null`` or ``cls_noise``. do_noise : bool - If True, return noise spectra along with data. + If True, return noise spectra and debiased spectra along with data. Returns ------- obs : OrderedDict Dictionary of data cross spectra nell : OrderedDict - Dictionary of noise cross spectra, or None if transfer_run is True. + Dictionary of noise cross spectra, or None if transfer_comp is True. + debias : OrderedDict + Dictionary of debiased data cross spectra """ # select map pairs if map_tag is not None: @@ -4568,31 +4494,26 @@ def get_data_spectra( map_pairs = pt.tag_pairs(map_tags) # select spectra - tbeb = "cmb_tb" in self.bin_def - if transfer_run or not tbeb: - specs = self.specs[:4] - else: - specs = self.specs + tbeb = "cmb_tb" in self.bin_def or "fg_tb" in self.bin_def + specs = list(self.specs) + if transfer_comp or not tbeb: + if "eb" in specs: + specs.remove("eb") + if "tb" in specs: + specs.remove("tb") + if transfer_comp == "fg" and "te" in specs: + specs.remove("te") # obs depends on what you're computing - if transfer_run: - if transfer_dust: - obs_quant = self.cls_signal_dust - elif transfer_cmbdust: - # need to add cmb and dust scaled by alphas - obs_quant = copy.deepcopy(self.cls_signal_dust) - adict = {"90": template_alpha90, "150": template_alpha150} - for k2 in self.cls_signal_dust[k1].keys(): - mtag1, mtag2 = k2.split(":") - freq1 = self.nom_freqs[mtag1] - freq2 = self.nom_freqs[mtag2] - obs_quant[k1][k2] = ( - self.cls_signal - + adict[freq1] * adict[freq2] * self.cls_signal_dust - ) + if transfer_comp: + if transfer_comp == "cmb": + obs_quant = self.cls_signal + elif transfer_comp == "fg": + obs_quant = self.cls_fg else: - obs_quant = self.cls_signal - + raise ValueError( + "Unknown transfer function component {}".format(transfer_comp) + ) elif self.null_run: if self.reference_subtracted: obs_quant = self.cls_data_sub_null @@ -4613,12 +4534,14 @@ def get_data_spectra( if not do_noise: return obs - nell = None - debias = None + elif transfer_comp: + return obs, None, None + + nell = OrderedDict() + debias = OrderedDict() + # Nulls are debiased by average of S+N sims - if self.null_run and not transfer_run and self.cls_sim_null is not None: - nell = OrderedDict() - debias = OrderedDict() + if self.null_run and self.cls_sim_null is not None: for spec in specs: nell[spec] = OrderedDict() debias[spec] = OrderedDict() @@ -4631,9 +4554,7 @@ def get_data_spectra( debias[spec][xname] = np.copy(self.cls_sim_null[spec][xname]) # Non-nulls are debiased by average of N sims - elif not transfer_run and self.cls_noise is not None: - nell = OrderedDict() - debias = OrderedDict() + elif self.cls_noise is not None: for spec in specs: nell[spec] = OrderedDict() debias[spec] = OrderedDict() @@ -4707,7 +4628,7 @@ def do_qb2cb(self, qb, inv_fish, wbl): for stag, wbl1 in wbl.items(): # compute conversion factors qb2cb[stag] = np.zeros((len(wbl1), len(wbl1))) - cls_shape = self.cls_shape["fg" if "fg" in stag else stag][: len(ell)] + cls_shape = self.cls_shape[stag][: len(ell)] v = norm * wbl1 * cls_shape bd = self.bin_def[stag] bw = self.bin_weights[stag] @@ -4751,8 +4672,6 @@ def fisher_precalc( cls_debias=None, likelihood=False, windows=False, - transfer_dust=False, - transfer_cmbdust=False, ): """ Pre-compute the D_ell and signal derivative matrices necessary for @@ -4796,13 +4715,9 @@ def fisher_precalc( pol_dim = 3 if self.pol else 1 dim1 = pol_dim * num_maps - if transfer_dust: - comps = ["dust"] - elif transfer_cmbdust: - comps = ["cmbdust"] - else: - comps = ["cmb"] - + comps = [] + if "cmb_tt" in cbl or "cmb_ee" in cbl: + comps += ["cmb"] if "fg_tt" in cbl: comps += ["fg"] if "res_tt" in cbl or "res_ee" in cbl: @@ -4819,8 +4734,11 @@ def fisher_precalc( if windows: Dmat_obs = None Mmat = OrderedDict() - for spec in specs: - Mmat[spec] = OrderedDict() + for comp in comps: + if comp not in ["cmb", "fg"]: + continue + for spec in specs: + Mmat["{}_{}".format(comp, spec)] = OrderedDict() mll = getattr(self, "mll", None) if mll is None: mll = self.kernel_precalc() @@ -4838,8 +4756,11 @@ def fisher_precalc( if likelihood: Dmat_obs_b[xname] = OrderedDict() elif windows: - for spec in specs: - Mmat[spec][xname] = OrderedDict() + for comp in comps: + if comp not in ["cmb", "fg"]: + continue + for spec in specs: + Mmat["{}_{}".format(comp, spec)][xname] = OrderedDict() else: Dmat_obs[xname] = OrderedDict() @@ -4848,10 +4769,14 @@ def fisher_precalc( # without bias subtraction for likelihood Dmat_obs_b[xname][spec] = cls_input[spec][xname] elif windows: - Mmat[spec][xname][spec] = mll[spec][xname] - if spec in ["ee", "bb"]: - mspec = "bb" if spec == "ee" else "ee" - Mmat[spec][xname][mspec] = mll["{}_mix".format(spec)][xname] + for comp in comps: + if comp not in ["cmb", "fg"]: + continue + stag = "{}_{}".format(comp, spec) + Mmat[stag][xname][spec] = mll[stag][xname] + if spec in ["ee", "bb"]: + mspec = "bb" if spec == "ee" else "ee" + Mmat[stag][xname][mspec] = mll["{}_mix".format(stag)][xname] else: if cls_debias is not None: Dmat_obs[xname][spec] = ( @@ -4883,16 +4808,6 @@ def fisher_precalc( mix_cbl = cbl[stag + "_mix"][xname] dSdqb[comp][xname][spec][mspec] = mix_cbl - if comp == "fg" and not transfer_dust: - # currently don't solve both beta_d and dust power. - # add delta beta bin for spectral index - dSdqb.setdefault("delta_beta", OrderedDict()).setdefault( - xname, OrderedDict() - ) - for spec in specs: - # this will be filled in in fisher_calc - dSdqb["delta_beta"][xname][spec] = OrderedDict() - return Dmat_obs, Dmat_obs_b, dSdqb, Mmat def clear_precalc(self): @@ -4922,8 +4837,6 @@ def fisher_calc( use_precalc=True, windows=False, inv_fish=None, - transfer_dust=False, - transfer_cmbdust=False, ): """ Re-compute the Fisher matrix and qb amplitudes based on @@ -4997,7 +4910,6 @@ def fisher_calc( well_cond = False pol_dim = 3 if self.pol else 1 - do_fg = "fg_tt" in cbl dkey = "Dmat_obs_b" if likelihood else "Mmat" if windows else "Dmat_obs" @@ -5010,8 +4922,6 @@ def fisher_calc( cls_debias=cls_debias, likelihood=likelihood, windows=windows, - transfer_dust=transfer_dust, - transfer_cmbdust=transfer_cmbdust, ) if use_precalc: self.Dmat_obs = Dmat_obs @@ -5035,9 +4945,6 @@ def fisher_calc( if "delta_beta" in qb: delta_beta = qb["delta_beta"][0] - if not likelihood: - dSdqb_mat1_freq = copy.deepcopy(dSdqb_mat1) - if likelihood or windows or not cond_noise: well_cond = True cond_noise = None @@ -5046,11 +4953,8 @@ def fisher_calc( Dmat1 = OrderedDict() if cls_model is None: - # TO DO: if dust, Eq 8 paper: model = dust + cmb cls_model = self.get_model_spectra( - qb, cbl, delta=True, cls_noise=cls_noise, cond_noise=cond_noise, - transfer_dust=transfer_dust, - transfer_cmbdust=transfer_cmbdust, + qb, cbl, cls_noise=cls_noise, cond_noise=cond_noise ) mkeys = list(cls_model) @@ -5071,47 +4975,45 @@ def fisher_calc( if well_cond: Dmat1_mat = pt.dict_to_dmat(Dmat1) - # Set up dSdqb frequency dependence - for xname, (m0, m1) in self.map_pairs.items(): - # transfer function does not have crosses - if xname not in cls_model[mkeys[0]]: - continue + # Set up dSdqb spectral index dependence + if "delta_beta" in qb and not likelihood: + # don't make changes to cached matrix + dSdqb_mat1 = copy.deepcopy(dSdqb_mat1) + dSdqb_mat1["delta_beta"] = OrderedDict() + + for xname in self.map_pairs: + # transfer function does not have crosses + if xname not in cls_model[mkeys[0]]: + continue + + dSdqb_mat1["delta_beta"][xname] = OrderedDict() - if do_fg and not likelihood: # Add spectral index dependence # dSdqb now depends on qb (spec index) because # model is non-linear so cannot be precomputed. # get foreground at pivot point spectral index - # and first derivative - freq_scale0, freq_scale_deriv = st.scale_dust( - self.map_freqs[m0], - self.map_freqs[m1], - ref_freq=self.ref_freq, - beta=self.beta_ref, - deriv=True, - ) - freq_scale = freq_scale0 + delta_beta * freq_scale_deriv - freq_scale_ratio = freq_scale_deriv / freq_scale + _, beta_scale, log_beta_scale = self.fg_scales[xname] + # separable beta correction + beta_scale = beta_scale ** delta_beta # non-linear + # beta_scale = (1 + log_beta_scale * delta_beta) # linear # scale foreground model by frequency scaling adjusted for beta - for s1, sdat in dSdqb_mat1_freq["fg"][xname].items(): + for s1, sdat in dSdqb_mat1["fg"][xname].items(): for s2, sdat2 in sdat.items(): - dSdqb_mat1_freq["fg"][xname][s1][s2] *= freq_scale + sdat2 *= beta_scale - # build delta_beta term from frequency scaled model, - # divide out frequeny scaling and apply derivative term - for s1, sdat in dSdqb_mat1_freq["delta_beta"][xname].items(): - sdat[s1] = cls_model["fg_{}".format(s1)][xname] * freq_scale_ratio + # with linearized derivative term + dSdqb_mat1["delta_beta"][xname][s1] = OrderedDict( + {s1: cls_model["fg_{}".format(s1)][xname] * log_beta_scale} + ) # Set up Dmat -- if it's not well conditioned, add noise to the # diagonal until it is. cond_iter = 0 while not well_cond: cls_model = self.get_model_spectra( - qb, cbl, delta=True, cls_noise=cls_noise, cond_noise=cond_noise, - transfer_dust=transfer_dust, - transfer_cmbdust=transfer_cmbdust, + qb, cbl, cls_noise=cls_noise, cond_noise=cond_noise ) for xname, (m0, m1) in self.map_pairs.items(): @@ -5162,7 +5064,11 @@ def fisher_calc( else: if not windows: Dmat_obs = pt.dict_to_dmat(Dmat_obs) - dSdqb_mat1_freq = pt.dict_to_dsdqb_mat(dSdqb_mat1_freq, self.bin_def) + # select only derivative terms for bins that are iterated + bin_def = OrderedDict( + [(k, v) for k, v in self.bin_def.items() if k in qb] + ) + dSdqb_mat1 = pt.dict_to_dsdqb_mat(dSdqb_mat1, bin_def) # apply ell limits if likelihood: @@ -5178,7 +5084,7 @@ def fisher_calc( else: if not windows: Dmat_obs = Dmat_obs[..., ell] - dSdqb_mat1_freq = dSdqb_mat1_freq[..., ell] + dSdqb_mat1 = dSdqb_mat1[..., ell] gmat = gmat[..., ell] lam, R = np.linalg.eigh(Dmat1.swapaxes(0, -1)) @@ -5213,7 +5119,7 @@ def fisher_calc( del lam, R eye = np.eye(len(gmat)) mat1 = np.einsum("ij...,jk...->ik...", eye, Dinv) - mat2 = np.einsum("klm...,ln...->knm...", dSdqb_mat1_freq, Dinv) + mat2 = np.einsum("klm...,ln...->knm...", dSdqb_mat1, Dinv) del Dinv mat = np.einsum("ik...,knm...->inm...", mat1, mat2) del mat1, mat2 @@ -5242,8 +5148,8 @@ def fisher_calc( # construct matrices for the qb and fisher terms, # and take the trace and sum over ell simultaneously if not windows or (windows and inv_fish is None): - fisher = np.einsum("iil,ijkl,jiml->km", gmat, mat, dSdqb_mat1_freq) / 2 - del dSdqb_mat1_freq + fisher = np.einsum("iil,ijkl,jiml->km", gmat, mat, dSdqb_mat1) / 2 + del dSdqb_mat1 if not windows: qb_vec = np.einsum("iil,ijkl,jil->k", gmat, mat, Dmat_obs) / 2.0 del gmat, mat, Dmat_obs @@ -5265,15 +5171,15 @@ def fisher_calc( # compute window functions for each spectrum for k, (left, right) in bin_index.items(): - comp, spec = k.split("_", 1) - if comp not in ["cmb", "fg"]: + if k not in Mmat: continue + comp, spec = k.split("_", 1) # select bins for corresponding spectrum sarg = arg[:, :, left:right] # construct Mll' matrix - marg = pt.dict_to_dmat(Mmat[spec], pol=self.pol)[..., ell, :] + marg = pt.dict_to_dmat(Mmat[k], pol=self.pol)[..., ell, :] # qb window function wbl1 = np.einsum("iil,ijkl,jilm->km", gmat, sarg, marg) / 2.0 / norm @@ -5285,7 +5191,7 @@ def fisher_calc( chi_bl[l:r] += w # check normalization - cls_shape = self.cls_shape["fg" if comp == "fg" else k][: len(norm)] + cls_shape = self.cls_shape[k][: len(norm)] self.log( "{} qb window function normalization: {}".format( k, np.sum(wbl1 * norm * chi_bl * cls_shape, axis=-1) @@ -5331,7 +5237,7 @@ def fisher_iterate( iter_max=200, converge_criteria=0.005, qb_start=None, - transfer_run=False, + transfer_comp=None, save_iters=False, null_first_cmb=False, delta_beta_prior=None, @@ -5341,10 +5247,6 @@ def fisher_iterate( like_profile_sigma=3.0, like_profile_points=100, file_tag=None, - transfer_dust=False, - transfer_cmbdust=False, - template_alpha90=None, - template_alpha150=None, ): """ Iterate over the Fisher calculation to compute bandpower estimates @@ -5369,9 +5271,9 @@ def fisher_iterate( qb_start : OrderedDict Initial guess at ``qb`` bandpower amplitudes. If None, unity is assumed for all bins. - transfer_run : bool - If True, the input Cls passed to ``fisher_calc`` are the average - of the signal simulations, and noise cls are ignored. + transfer_comp : str + If set, the input Cls passed to ``fisher_calc`` are the average + of the component simulations, and noise cls are ignored. If False, the input Cls are either ``cls_data_null`` (for null tests) or ``cls_data``. See ``get_masked_data`` for how these are computed. The input noise is similarly either @@ -5424,7 +5326,7 @@ def fisher_iterate( cls_shape : shape spectrum iters : number of iterations completed - If ``transfer_run`` is False, this dictionary also contains:: + If ``transfer_comp`` is not set, this dictionary also contains:: qb_transfer : transfer function amplitudes transfer : ell-by-ell transfer function @@ -5434,15 +5336,16 @@ def fisher_iterate( Notes ----- This method stores outputs to files with name 'transfer' for transfer - function runs (if ``transfer_run = True``), otherwise with name + function runs (if ``transfer_comp`` is set), otherwise with name 'bandpowers'. Outputs from each individual iteration, containing only the quantities that change with each step, are stored in separate files with the same name (and different index). """ - if transfer_run: + if transfer_comp: null_first_cmb = False - save_name = "transfer_wbins" if self.weighted_bins else "transfer" + save_name = "transfer_{}wbins" if self.weighted_bins else "transfer_{}" + save_name = save_name.format(transfer_comp) else: save_name = "bandpowers" @@ -5453,15 +5356,22 @@ def fisher_iterate( if qb_start is None: qb = OrderedDict() for k, v in self.bin_def.items(): - if transfer_run: - if "cmb" not in k or "eb" in k or "tb" in k: + # transfer function for signal component + if transfer_comp: + if "eb" in k or "tb" in k: + continue + if transfer_comp == "fg" and "te" in k: continue + if not k.startswith(transfer_comp): + continue + qb[k] = np.ones(len(v)) + continue if k == "delta_beta": # qb_delta beta is a coefficient on the change from beta, # so expect that it should be small if beta_ref is close # (zeroes cause singular matrix problems) qb[k] = [self.delta_beta_fix] - elif k.startswith("res_") or k.startswith("fg_"): + elif k.startswith("res_"): # res qb=0 means noise model is 100% accurate. qb[k] = 1e-5 * np.ones(len(v)) else: @@ -5470,18 +5380,8 @@ def fisher_iterate( else: qb = qb_start - # TO DO: GET TEMPLATE INTERNALLY INSTEAD OF AS INPUTS - #template_alpha90=self.template_alpha["90"] - #template_alpha150=self.template_alpha["150"] - - # reminder: obs = cmb + a*dust if transfer_cmbdust obs, nell, debias = self.get_data_spectra( - map_tag=map_tag, - transfer_run=transfer_run, - transfer_dust=transfer_dust, - transfer_cmbdust=transfer_cmbdust, - template_alpha90=tempalate_alpha90, - template_alpha150=tempalate_alpha150, + map_tag=map_tag, transfer_comp=transfer_comp ) bin_index = pt.dict_to_index(self.bin_def) @@ -5502,8 +5402,6 @@ def fisher_iterate( delta_beta_prior=delta_beta_prior, cond_criteria=cond_criteria, null_first_cmb=null_first_cmb, - transfer_dust=transfer_dust, - transfer_cmbdust=transfer_cmbdust, ) qb_arr = pt.dict_to_arr(qb, flatten=True) @@ -5533,9 +5431,7 @@ def fisher_iterate( # put qb_new in original dict qb = copy.deepcopy(qb_new) - cls_model = self.get_model_spectra( - qb, cbl, delta=True, cls_noise=nell, cond_noise=None - ) + cls_model = self.get_model_spectra(qb, cbl, cls_noise=nell, cond_noise=None) if "delta_beta" in qb: # get beta fit and beta error @@ -5561,7 +5457,6 @@ def fisher_iterate( inv_fish=inv_fish, cls_model=cls_model, cbl=cbl, - map_freqs=self.map_freqs, cls_signal=self.cls_signal, cls_noise=self.cls_noise, Dmat_obs=self.Dmat_obs, @@ -5571,13 +5466,17 @@ def fisher_iterate( if "fg_tt" in self.bin_def: out.update( - beta_fit=beta_fit, - beta_err=beta_err, - ref_freq=self.ref_freq, + map_freqs=self.map_freqs, + freq_ref=self.freq_ref, beta_ref=self.beta_ref, ) + if "delta_beta" in qb: + out.update( + beta_fit=beta_fit, + beta_err=beta_err, + ) - self.save_data(save_name, bp_opts=not transfer_run, **out) + self.save_data(save_name, bp_opts=not transfer_comp, **out) (nans,) = np.where(np.isnan(qb_new_arr)) if len(nans): @@ -5600,7 +5499,7 @@ def fisher_iterate( ) if np.nanmax(np.abs(fqb)) < converge_criteria: - if not transfer_run: + if not transfer_comp: # Calculate final fisher matrix without conditioning self.log("Calculating final Fisher matrix.", "info") _, inv_fish = self.fisher_calc( @@ -5629,13 +5528,15 @@ def fisher_iterate( else: msg = "{} {} did not converge in {} iterations".format( "Multi-map" if map_tag is None else "Map {}".format(map_tag), - "transfer function" if transfer_run else "spectrum", + "{} transfer function".format(transfer_comp) + if transfer_comp + else "spectrum", iter_max, ) # Check the slope of the last ten fqb_maxpoints. # If there's not a downward trend, adjust conditioning # criteria to help convergence. - if len(prev_fqb) <= 10 or transfer_run: + if len(prev_fqb) <= 10 or transfer_comp: continue m, b = np.polyfit(np.arange(10), prev_fqb[-10:], 1) if m > 0: # Not converging @@ -5677,7 +5578,6 @@ def fisher_iterate( iters=iter_idx, success=success, map_tags=self.map_tags, - map_freqs=self.map_freqs, converge_criteria=converge_criteria, cond_noise=cond_noise, cond_criteria=cond_criteria, @@ -5688,12 +5588,16 @@ def fisher_iterate( if "fg_tt" in self.bin_def: out.update( - delta_beta_prior=delta_beta_prior, - beta_fit=beta_fit, - beta_err=beta_err, - ref_freq=self.ref_freq, + map_freqs=self.map_freqs, + freq_ref=self.freq_ref, beta_ref=self.beta_ref, ) + if "delta_beta" in qb: + out.update( + delta_beta_prior=delta_beta_prior, + beta_fit=beta_fit, + beta_err=beta_err, + ) if self.debug: out.update( @@ -5707,7 +5611,7 @@ def fisher_iterate( Dmat_obs=self.Dmat_obs, ) - if not transfer_run: + if not transfer_comp: out.update( qb_transfer=self.qb_transfer, transfer=self.transfer, @@ -5715,7 +5619,7 @@ def fisher_iterate( if self.template_cleaned: out.update(template_alpha=self.template_alpha) - if success and not transfer_run: + if success and not transfer_comp: # do one more fisher calc that doesn't include sample variance # set qb=very close to 0. 0 causes singular matrix problems. # don't do this for noise residual bins @@ -5748,8 +5652,8 @@ def fisher_iterate( invfish_nosampvar=inv_fish_ns, ) - # compute window functions for CMB bins - self.log("Calculating window functions for CMB bins", "info") + # compute window functions for signal bins + self.log("Calculating signal window functions", "info") wbl_qb = self.fisher_calc( qb, cbl, @@ -5857,7 +5761,7 @@ def fisher_iterate( return self.save_data( save_name, map_tag=map_tag, - bp_opts=not transfer_run, + bp_opts=not transfer_comp, extra_tag=file_tag, **out ) @@ -5868,37 +5772,33 @@ def get_transfer( iter_max=200, save_iters=False, fix_bb_transfer=False, - dust=False, - cmbdust=False, - template_alpha90=None, - template_alpha150=None, + fix_fg_transfer=False, ): """ - Compute the transfer function from signal simulations created using - the same spectrum as the input shape. + Compute the transfer function from signal simulations created using the + same spectrum as the input shape. - This raises a ValueError if a negative transfer function amplitude - is found. + This raises a ValueError if a negative transfer function amplitude is + found. Arguments --------- converge_criteria : float - Maximum fractional change in qb that indicates convergence and - stops iteration. + Maximum fractional change in qb that indicates convergence and stops + iteration. iter_max : int - Maximum number of iterations to perform. if this limit is - reached, a warning is issued. + Maximum number of iterations to perform. if this limit is reached, + a warning is issued. save_iters : bool - If True, the output data from each Fisher iteration are stored - in an individual npz file. + If True, the output data from each Fisher iteration are stored in an + individual npz file. fix_bb_transfer : bool - If True, after transfer functions have been calculated, impose - the BB xfer is exactly equal to the EE transfer. - dust : bool - If True, compute transfer function using dust sims - cmbdust : bool - If True, compute transfer function using CMB+dust sims, with - dust scaled by template_alpha90, template_alpha150 + If True, after transfer functions have been calculated, impose the + BB xfer is exactly equal to the EE transfer. + fix_fg_transfer : bool + If True, the transfer function for the foreground component is fixed + to an interpolation of the CMB transfer function onto the foreground + binning. Returns ------- @@ -5907,8 +5807,8 @@ def get_transfer( Notes ----- - This method is called at the 'transfer' checkpoint, and loads or saves - a data dictionary named 'transfer_all' with the following entries: + This method is called at the 'transfer' checkpoint, and loads or saves a + data dictionary named 'transfer_all' with the following entries: nbins: number of bins @@ -5919,13 +5819,14 @@ def get_transfer( transfer: ell-by-ell transfer function for each map and spectrum component - Additionally the final output of ``fisher_iterate`` is stored - in a dictionary called ``transfer_map`` for each map. + Additionally the final output of ``fisher_iterate`` is stored in a + dictionary called ``transfer_map`` for each map. """ - transfer_shape = ( - self.num_maps * len(self.specs), - self.nbins_cmb / len(self.specs), - ) + comps = [] + if "cmb_tt" in self.bin_def or "cmb_ee" in self.bin_def: + comps += ["cmb"] + if "fg_tt" in self.bin_def or "fg_ee" in self.bin_def: + comps += ["fg"] opts = dict( converge_criteria=converge_criteria, @@ -5934,12 +5835,12 @@ def get_transfer( weighted_bins=self.weighted_bins, ) - save_name = "transfer_all" - if dust: - save_name = "{}_dust".format(save_name) - elif cmbdust: - save_name = "{}_cmbdust".format(save_name) + if "fg" in comps: + if self.cls_fg is None: + fix_fg_transfer = True + opts.update(fix_fg_transfer=fix_fg_transfer) + save_name = "transfer_all" if self.weighted_bins: save_name = "{}_wbins".format(save_name) @@ -5947,149 +5848,155 @@ def get_transfer( save_name, "transfer", to_attrs=False, - shape_ref="qb_transfer", - shape=transfer_shape, value_ref=opts, ) # function for converting qb's to ell-by-ell transfer function - def expand_transfer(qb_transfer, bin_def): - transfer = OrderedDict() - for spec in self.specs: - stag = "cmb_{}".format(spec) - transfer[spec] = OrderedDict() - bd = bin_def[stag] - for m0 in self.map_tags: - if spec == "bb" and fix_bb_transfer: - transfer[spec][m0] = transfer["ee"][m0] - elif spec in ["tb", "eb"]: - speca, specb = [s + s for s in spec] - transfer[spec][m0] = np.sqrt( - np.abs(transfer[speca][m0] * transfer[specb][m0]) - ) - else: - transfer[spec][m0] = pt.expand_qb( - qb_transfer[stag][m0], bd, self.lmax + 1 - ) - return transfer + def expand_transfer(qb_transfer, bin_def, check_only=False): + if not check_only: + transfer = OrderedDict() + if not check_only and fix_fg_transfer and "fg" in comps: + for stag in [k for k in qb_transfer if k.startswith("cmb_")]: + ftag = stag.replace("cmb_", "fg_") + for m0, v in qb_transfer[stag].items(): + cb = np.mean(bin_def[stag], axis=1) + fb = np.mean(bin_def[ftag], axis=1) + vint = np.interp(fb, np.append([0], cb), np.append([0], v)) + qb_transfer[ftag][m0] = vint + for comp in comps: + for spec in self.specs: + stag = "{}_{}".format(comp, spec) + if stag not in bin_def or stag not in qb_transfer: + raise KeyError(stag) + if not check_only: + transfer[stag] = OrderedDict() + bd = bin_def[stag] + for m0 in self.map_tags: + if m0 not in qb_transfer[stag]: + raise KeyError(m0) + lq = len(qb_transfer[stag][m0]) + lb = len(bin_def[stag]) + if lq != lb: + raise ValueError( + "Found {} transfer bins for component {} map {}, expected {}".format( + lq, stag, m0, lb, + ) + ) + if check_only: + continue + if spec == "bb" and fix_bb_transfer: + v = transfer["{}_ee".format(comp)][m0] + elif spec in ["te", "tb", "eb"] and spec not in qb_transfer: + staga, stagb = ["{}_{}{}".format(comp, s, s) for s in spec] + v = np.sqrt( + np.abs(transfer[staga][m0] * transfer[stagb][m0]) + ) + else: + v = pt.expand_qb(qb_transfer[stag][m0], bd, self.lmax) + transfer[stag][m0] = v + if not check_only: + return transfer if ret is not None: - #self.qb_transfer = ret["qb_transfer"] - if dust: - self.qb_transfer_dust = ret["qb_transfer"] - elif cmbdust: - self.qb_transfer_cmbdust = ret["qb_transfer"] + try: + check_only = ret.get("transfer", None) is not None + xfer = expand_transfer(ret["qb_transfer"], ret["bin_def"], check_only) + except: + ret = None else: self.qb_transfer = ret["qb_transfer"] - # return ret["qb_transfer"] - - if "transfer" in ret: - self.transfer = ret["transfer"] - else: - self.transfer = expand_transfer(ret["qb_transfer"], ret["bin_def"]) - return self.transfer - - #self.qb_transfer = OrderedDict() - qb_transfer = OrderedDict() - if dust: - comp = "dust" - elif cmbdust: - comp = "cmbdust" - else: - comp = "cmb" + self.transfer = xfer or ret["transfer"] + return self.transfer - for spec in self.specs: - #self.qb_transfer["cmb_" + spec] = OrderedDict() - #if dust: This is not necessary ??? - qb_transfer["{}_{}".format(comp, spec)] = OrderedDict() - - # success = False - success = True + self.qb_transfer = OrderedDict() + success = True msg = "" - for im0, m0 in enumerate(self.map_tags): - if not self.fit_transfer[self.map_tags_orig[im0]]: - #for spec in self.specs: - # self.qb_transfer["cmb_{}".format(spec)][m0] = np.ones( - # self.nbins_cmb // len(self.specs) - # ) - for stag in qb_transfer.keys(): - qb_transfer[stag][m0] = np.ones( - self.nbins_cmb // len(self.specs) + for comp in comps: + for spec in self.specs: + self.qb_transfer["{}_{}".format(comp, spec)] = OrderedDict() + + for im0, m0 in enumerate(self.map_tags): + if not self.fit_transfer[self.map_tags_orig[im0]]: + for spec in self.specs: + stag = "{}_{}".format(comp, spec) + nbins = self.nbins_cmb if comp == "cmb" else self.nbins_fg + self.qb_transfer[stag][m0] = np.ones(nbins // len(self.specs)) + self.log("Setting map {} transfer to unity".format(m0), "info") + continue + + if comp == "fg" and fix_fg_transfer: + # if no input spectrum given, set foreground transfer + # function equal to CMB transfer function, interpolated to + # the foreground bins + self.log( + "Setting map {} fg transfer to cmb transfer".format(m0), "info" ) - self.log("Setting map {} transfer to unity".format(m0), "info") - continue + continue - self.log( - "Computing transfer function for map {}/{}".format( - im0 + 1, self.num_maps - ), - "info", - ) - self.clear_precalc() - #cbl = self.bin_cl_template(map_tag=m0, transfer_run=True) - cbl = self.bin_cl_template( - m0, transfer_run=True, dust=dust, cmbdust=cmbdust + self.log( + "Computing {} transfer function for map {}/{}".format( + comp, im0 + 1, self.num_maps + ), + "info", ) - - ret = self.fisher_iterate( - cbl, - m0, - transfer_run=True, - iter_max=iter_max, - converge_criteria=converge_criteria, - save_iters=save_iters, - transfer_dust=dust, - transfer_cmbdust=cmbdust, - template_alpha90=tempalate_alpha90, - template_alpha150=template_alpha150, - ) - qb = ret["qb"] - - success = success and ret["success"] - if not ret["success"]: - msg = "Error in fisher_iterate for map {}".format(m0) - - # fix negative amplitude bins - for k, v in qb.items(): - if np.any(v < 0): - (negbin,) = np.where(v < 0) - self.warn( - "Transfer function amplitude {} < 0" - "for {} bin {} of map {}".format(v, k, negbin, m0) - ) - # XXX cludge - # This happens in first bin - # try linear interp between zero and next value - try: - qb[k][negbin] = qb[k][negbin + 1] / 2.0 + self.clear_precalc() + cbl = self.bin_cl_template(map_tag=m0, component=comp, transfer=True) + ret = self.fisher_iterate( + cbl, + m0, + transfer_comp=comp, + iter_max=iter_max, + converge_criteria=converge_criteria, + save_iters=save_iters, + ) + qb = ret["qb"] + + success = success and ret["success"] + if not ret["success"]: + msg = "Error in fisher_iterate for map {}".format(m0) + + # fix negative amplitude bins + for k, v in qb.items(): + if np.any(v < 0): + (negbin,) = np.where(v < 0) self.warn( - "Setting Transfer function in negative bin to small " - "positive. This is probably due to choice of bins or " - "insufficient number of signal sims" + "Transfer function amplitude {} < 0" + "for {} bin {} of map {}".format(v, k, negbin, m0) ) - except Exception as e: - msg = "Unable to adjust negative bins for map {}: {}".format( - m0, str(e) - ) - success = False + # XXX cludge + # This happens in first bin + # try linear interp between zero and next value + try: + qb[k][negbin] = qb[k][negbin + 1] / 2.0 + self.warn( + "Setting Transfer function in negative bin to small " + "positive. This is probably due to choice of bins or " + "insufficient number of signal sims" + ) + except Exception as e: + msg = ( + "Unable to adjust negative bins for map {}: {}".format( + m0, str(e) + ) + ) + success = False - # Set EB/TB qb transfers to geometric means of components - if len(self.specs) > 4: - #qb["cmb_eb"] = np.sqrt(np.abs(qb["cmb_ee"] * qb["cmb_bb"])) - #qb["cmb_tb"] = np.sqrt(np.abs(qb["cmb_tt"] * qb["cmb_bb"])) - qb[comp + "_eb"] = np.sqrt(np.abs(qb[comp + "_ee"] * qb[comp + "_bb"])) - qb[comp + "_tb"] = np.sqrt(np.abs(qb[comp + "_tt"] * qb[comp + "_bb"])) + # Set TE/EB/TB qb transfers to geometric means of components, if necessary + for spec in ["te", "eb", "tb"]: + if spec not in self.specs: + continue + stag = "{}_{}".format(comp, spec) + if stag in qb: + continue + staga, stagb = ["{}_{}{}".format(comp, s, s) for s in spec] + qb[stag] = np.sqrt(np.abs(qb[staga] * qb[stagb])) - for stag, qbdat in qb.items(): - #self.qb_transfer[stag][m0] = qbdat - qb_transfer[stag][m0] = qbdat + for stag, qbdat in qb.items(): + self.qb_transfer[stag][m0] = qbdat if success: - if dust: - self.transfer_dust = expand_transfer(self.qb_transfer, self.bin_def) - else: - self.transfer = expand_transfer(self.qb_transfer, self.bin_def) + self.transfer = expand_transfer(self.qb_transfer, self.bin_def) else: self.transfer = None @@ -6103,16 +6010,8 @@ def expand_transfer(qb_transfer, bin_def): if not success: raise RuntimeError("Error computing transfer function: {}".format(msg)) - - if dust: - self.qb_transfer_dust = qb_transfer - elif cmbdust: - self.qb_transfer_cmbdust = qb_transfer - else: - self.qb_transfer = qb_transfer - return qb_transfer - #return self.transfer + return self.transfer def get_bandpowers( self, @@ -6121,8 +6020,6 @@ def get_bandpowers( iter_max=200, return_qb=False, save_iters=False, - ref_freq=359.7, - beta_ref=1.54, delta_beta_prior=None, cond_noise=None, cond_criteria=None, @@ -6132,8 +6029,6 @@ def get_bandpowers( like_profile_sigma=3.0, like_profile_points=100, file_tag=None, - dust=False, - cmbdust=False, ): """ Compute the maximum likelihood bandpowers of the data, assuming @@ -6160,13 +6055,6 @@ def get_bandpowers( save_iters : bool If True, the output data from each Fisher iteration are stored in an individual npz file. - ref_freq : float - In GHz, reference frequency for dust model. Dust bandpowers output - will be at this reference frequency. - beta_ref : float - The spectral index of the dust model. This is a fixed value, with - an additive deviation from this value fit for in foreground fitting - mode. delta_beta_prior : float The width of the prior on the additive change from beta_ref. If you don't want the code to fit for a spectral index different @@ -6180,11 +6068,6 @@ def get_bandpowers( Keep first CMB bandpowers fixed to input shape (qb=1). return_cls : bool If True, return C_ls rather than D_ls - dust : bool - If True, compute transfer function using dust sims - cmbdust : bool - If True, compute transfer function using CMB+dust sims, with - dust scaled by template_alpha90, template_alpha150 like_profiles : bool If True, compute profile likelihoods for each qb, leaving all others fixed at their maximum likelihood values. Profiles are @@ -6229,16 +6112,14 @@ def get_bandpowers( # foreground fitting if "fg_tt" in self.bin_def: - # reference frequency and spectral index - self.ref_freq = ref_freq - self.beta_ref = beta_ref - # priors on frequency spectral index - self.delta_beta_fix = 1.0e-8 opts.update( - delta_beta_prior=delta_beta_prior, - ref_freq=self.ref_freq, + freq_ref=self.freq_ref, beta_ref=self.beta_ref, ) + if "delta_beta" in self.bin_def: + # priors on frequency spectral index + self.delta_beta_fix = 1.0e-8 + opts.update(delta_beta_prior=delta_beta_prior) if self.template_cleaned: opts.update(template_alpha=self.template_alpha) @@ -6262,12 +6143,11 @@ def get_bandpowers( self.clear_precalc() - cbl = self.bin_cl_template(map_tag=map_tag, transfer_run=False) + cbl = self.bin_cl_template(map_tag=map_tag) ret = self.fisher_iterate( cbl, map_tag, - transfer_run=False, iter_max=iter_max, converge_criteria=converge_criteria, save_iters=save_iters, @@ -6529,11 +6409,6 @@ def get_likelihood( self.bin_def[k] = bd[good] v = qb.pop(k)[good] num_res += len(v) - - # use average qb res in good range per map - # self.bin_def[k] = np.array([[lmin, lmax + 1]]) - # v = np.array([(qb.pop(k)[good]).mean()]) - # num_res += 1 qb_res[k] = v self.nbins_res = num_res @@ -6543,7 +6418,7 @@ def get_likelihood( self.log("Computing model spectrum", "debug") self.warn("Beam variation not implemented for case of no r fit") cbl = self.bin_cl_template(map_tag=map_tag) - cls_model = self.get_model_spectra(qb, cbl, delta=True, cls_noise=cls_noise) + cls_model = self.get_model_spectra(qb, cbl, cls_noise=cls_noise) else: qb = copy.deepcopy(qb) for spec in self.specs: @@ -6552,6 +6427,9 @@ def get_likelihood( if stag not in qb: continue qb[stag] = np.ones_like(qb[stag]) + qb_cmb = OrderedDict( + [(k, v) for k, v in qb.items() if k.startswith("cmb_")] + ) self.log("Computing r model spectrum", "debug") cls_shape_scalar = self.get_signal_shape( @@ -6563,28 +6441,26 @@ def get_likelihood( ) # load tensor and scalar terms separately + # include all model components here cbl_scalar = self.bin_cl_template(cls_shape_scalar, map_tag) cls_model_scalar = self.get_model_spectra( - qb, cbl_scalar, delta=True, cls_noise=cls_noise + qb, cbl_scalar, cls_noise=cls_noise ) - cbl_tensor = self.bin_cl_template(cls_shape_tensor, map_tag) - cls_model_tensor = self.get_model_spectra( - qb, cbl_tensor, delta=False, res=False + # include only CMB components for all additive terms below + cbl_tensor = self.bin_cl_template( + cls_shape_tensor, map_tag, component="cmb" ) + cls_model_tensor = self.get_model_spectra(qb_cmb, cbl_tensor) if beam_prior is not None: # load beam error term for tensor and scalar cbl_scalar_beam = self.bin_cl_template( - cls_shape_scalar, map_tag, beam_error=True - ) - cls_mod_scal_beam = self.get_model_spectra( - qb, cbl_scalar_beam, delta=True, res=False + cls_shape_scalar, map_tag, component="cmb", beam_error=True ) + cls_mod_scal_beam = self.get_model_spectra(qb_cmb, cbl_scalar_beam) cbl_tensor_beam = self.bin_cl_template( - cls_shape_tensor, map_tag, beam_error=True - ) - cls_mod_tens_beam = self.get_model_spectra( - qb, cbl_tensor_beam, delta=False, res=False + cls_shape_tensor, map_tag, component="cmb", beam_error=True ) + cls_mod_tens_beam = self.get_model_spectra(qb_cmb, cbl_tensor_beam) cbl = copy.deepcopy(cbl_scalar) cls_model = copy.deepcopy(cls_model_scalar) @@ -6719,7 +6595,6 @@ def log_like( beam_term += bc * cls_mod_scal_beam[ctag][xname][bn] dd[:] = cls_model_scalar[stag][xname] + beam_term - # compute noise model terms if res is None: clsm = cls_model0 diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index 5a9fab2b..facc37d6 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -38,12 +38,13 @@ def xfaster_run( signal_type="synfast", signal_subset="*", signal_transfer_type=None, - signal_transfer_type_dust=None, - signal_subset_dust=None, signal_type_sim=None, noise_type="stationary", noise_subset="*", noise_type_sim=None, + foreground_type=None, + foreground_subset="*", + foreground_transfer_type=None, foreground_type_sim=None, template_type=None, template_noise_type=None, @@ -60,6 +61,7 @@ def xfaster_run( bin_width_res=25, res_specs=None, foreground_fit=False, + beta_fit=False, bin_width_fg=30, # data options template_alpha_tags=None, @@ -95,15 +97,18 @@ def xfaster_run( save_iters=False, return_cls=False, fix_bb_transfer=False, + fix_fg_transfer=False, null_first_cmb=False, qb_file_sim=None, signal_spec=None, signal_transfer_spec=None, - signal_transfer_spec_dust=None, - separate_dust_Fl=False, + foreground_spec=None, + foreground_transfer_spec=None, model_r=None, - ref_freq=359.7, + freq_ref=359.7, beta_ref=1.54, + freq_ref_transfer=None, + beta_ref_transfer=None, delta_beta_prior=0.5, like_profiles=False, like_profile_sigma=3.0, @@ -175,14 +180,6 @@ def xfaster_run( signal_transfer_type : str The variant of signal sims to use for computing the transfer function. If not set, defaults to ``signal_type``. - signal_subset_dust: str - If not provided, use signal_subset values - Add more. - signal_transfer_type_dust : string - The variant of signal sims to use for transfer function for dust, - or, to add to signal_transfer type for raw (uncleaned) spectra, - multiplied by alpha90**2 or alpha150**2 - If provided, will compute separate transfer functions for dust and cmb. signal_type_sim : str The variant of signal sims to use for sim_index data maps. This enables having a different noise sim ensemble to use for @@ -197,6 +194,16 @@ def xfaster_run( The variant of noise sims to use for sim_index fake data map. This enables having a different noise sim ensemble to use for sim_index run than the ensemble from which the noise is computed. + foreground_type : str + The variant of foreground sims to use for the foreground component of + the covariance model. If not supplied, and ``foreground_fit`` is False, + no foreground component is included. + foreground_subset : str + The subset of the foreground sims to include. Must be a glob-parseable + string. + foreground_transfer_type : string + The variant of foreground sims to use for the transfer function for a + foreground component. If not set, defaults to ``foreground_type``. foreground_type_sim : str Tag for directory (foreground_) where foreground sims are stored that should be added to the signal and noise sims @@ -246,11 +253,19 @@ def xfaster_run( residuals simultaneously. If not supplied, this defaults to EEBB for polarized maps, or TT for unpolarized maps. foreground_fit : bool - Include foreground residuals in the estimator. + Include foreground bandpowers in the estimator. If ``foreground_type`` + is set, uses the foreground model in the appropriate sims directory, and + includes a separate transfer function calculation for the foreground + components of the covariance. Otherwise, uses a simple power-law model + with the same transfer function as the main signal component. + beta_fit : bool + If True, include a fit for ``delta_beta`` in the estimator. Otherwise, + only fit for foreground amplitudes. bin_width_fg : int or array_like of 6 ints - Width of each ell bin for each of the six output foreground spectra - (TT, EE, BB, TE, EB, TB). EE/BB bins should be the same - in order to handle mixing correctly. + Width of each ell bin for each of the six output foreground spectra (TT, + EE, BB, TE, EB, TB). EE/BB bins should be the same in order to handle + mixing correctly. Ignored if neither ``foreground_fit`` nor + ``foreground_type`` is set. template_alpha_tags : list of strings List of map tags from which foreground template maps should be subtracted. These should be the original map tags, not @@ -352,6 +367,10 @@ def xfaster_run( fix_bb_transfer : bool If True, after transfer functions have been calculated, impose that the BB transfer function is exactly equal to the EE transfer function. + fix_fg_transfer : bool + If True, the transfer function for the foreground component is fixed to + an interpolation of the CMB transfer function onto the foreground + binning. null_first_cmb : bool If True, keep first CMB bandpowers fixed to input shape (qb=1). qb_file_sim : str @@ -359,28 +378,40 @@ def xfaster_run( to correct the noise ensemble by an appropriate set of residual ``qb`` values. signal_spec : str - The spectrum data file to use for estimating bandpowers. If not - supplied, will search for ``spec_signal_.dat`` in the signal - sim directory. + The spectrum data file to use for estimating signal component + bandpowers. If not supplied, will search for + ``spec_signal_.dat`` in the signal sim directory. signal_transfer_spec : str - The spectrum data file used to generate signal sims. If not - supplied, will search for ``spec_signal_.dat`` in the - transfer signal sim directory. Used for computing transfer functions. - signal_transfer_spec_dust : string - The spectrum data file used to generate signal sims for dust transfer - function. If not supplied, will search for - `spec_signal_.dat` in the - dust transfer signal sim directory. Used for computing transfer functions. + The spectrum data file used to generate signal sims for computing the + signal transfer function. If not supplied, will search for + ``spec_signal_.dat`` in the transfer signal sim + directory. + foreground_spec : string + The spectrum data file to use for estimating foreground component + bandpowers. If not supplied, will search for + ``spec_foreground_.dat`` in the foreground sim + directory. + foreground_transfer_spec : string + The spectrum data file used to generate foreground sims for computing + the foreground transfer function. If not supplied, will search for + ``spec_foreground_.dat`` in the foreground + transfer sim directory. model_r : float The ``r`` value to use to compute a spectrum for estimating bandpowers. Overrides ``signal_spec``. - ref_freq : float + freq_ref : float In GHz, reference frequency for dust model. Dust bandpowers output will be at this reference frequency. beta_ref : float The spectral index of the dust model. This is a fixed value, with an additive deviation from this value fit for in foreground fitting mode. + freq_ref_transfer : float + Reference frequency for the dust model used for computing the foreground + transfer function. Defaults to ``freq_ref`` if not given. + beta_ref_transfer : float + Reference spectral index for the dust model used for computing the + foreground transfer function. Defaults to ``beta_ref`` if not given. delta_beta_prior : float The width of the prior on the additive change from beta_ref. If you don't want the code to fit for a spectral index different @@ -509,6 +540,7 @@ def xfaster_run( data_subset=data_subset, signal_subset=signal_subset, noise_subset=noise_subset, + foreground_subset=foreground_subset, data_type=data_type, noise_type=noise_type, noise_type_sim=noise_type_sim, @@ -516,11 +548,11 @@ def xfaster_run( signal_type=signal_type, signal_type_sim=signal_type_sim if sim_data_r is None else "r", signal_transfer_type=signal_transfer_type, - signal_transfer_type_dust=signal_transfer_type_dust, - signal_subset_dust=signal_subset_dust, + foreground_type=foreground_type, + foreground_transfer_type=foreground_transfer_type, + foreground_type_sim=foreground_type_sim, data_root2=data_root2, data_subset2=data_subset2, - foreground_type_sim=foreground_type_sim, template_type=template_type, template_noise_type=template_noise_type, template_type_sim=template_type_sim, @@ -555,6 +587,7 @@ def xfaster_run( res_specs=res_specs, bin_width_res=bin_width_res, foreground_fit=foreground_fit, + beta_fit=beta_fit, bin_width_fg=bin_width_fg, ) config_vars.update(bin_opts, "Binning Options") @@ -590,18 +623,27 @@ def xfaster_run( ) config_vars.update(kernel_opts, "Beam and Kernel Options") + if freq_ref_transfer is None: + freq_ref_transfer = freq_ref + if beta_ref_transfer is None: + beta_ref_transfer = beta_ref + spec_opts = dict( multi_map=multi_map, converge_criteria=converge_criteria, iter_max=iter_max, save_iters=save_iters, fix_bb_transfer=fix_bb_transfer, + fix_fg_transfer=fix_fg_transfer, signal_spec=signal_spec, signal_transfer_spec=signal_transfer_spec, - signal_transfer_spec_dust=signal_transfer_spec_dust, + foreground_spec=foreground_spec, + foreground_transfer_spec=foreground_transfer_spec, model_r=model_r, - ref_freq=ref_freq, + freq_ref=freq_ref, beta_ref=beta_ref, + freq_ref_transfer=freq_ref_transfer, + beta_ref_transfer=beta_ref_transfer, delta_beta_prior=delta_beta_prior, cond_noise=cond_noise, cond_criteria=cond_criteria, @@ -619,18 +661,22 @@ def xfaster_run( spec_opts.pop("multi_map") spec_opts.pop("signal_spec") spec_opts.pop("signal_transfer_spec") - spec_opts.pop("signal_transfer_spec_dust") + spec_opts.pop("foreground_spec") + spec_opts.pop("foreground_transfer_spec") + spec_opts.pop("freq_ref") + spec_opts.pop("beta_ref") + spec_opts.pop("freq_ref_transfer") + spec_opts.pop("beta_ref_transfer") spec_opts.pop("model_r") spec_opts.pop("qb_file") bandpwr_opts = spec_opts.copy() bandpwr_opts.pop("fix_bb_transfer") + bandpwr_opts.pop("fix_fg_transfer") spec_opts.pop("file_tag") transfer_opts = spec_opts.copy() transfer_opts.pop("cond_noise") transfer_opts.pop("cond_criteria") - transfer_opts.pop("ref_freq") - transfer_opts.pop("beta_ref") transfer_opts.pop("delta_beta_prior") transfer_opts.pop("null_first_cmb") transfer_opts.pop("return_cls") @@ -683,73 +729,26 @@ def xfaster_run( X.log("Computing kernels...", "notice") X.get_kernels(window_lmax=window_lmax) - X.log("Computing sim ensemble averages for transfer function...", - "notice") + X.log("Computing sim ensemble averages for transfer function...", "notice") # Do all the sims at once to also get the S+N sim ensemble average do_noise = signal_transfer_type in [signal_type, None] - - X.get_masked_sims( - transfer=True, - do_noise=do_noise, - qb_file=qb_file_sim) + X.get_masked_sims(transfer=True, do_noise=do_noise, qb_file=qb_file_sim) X.log("Computing beam window functions...", "notice") X.get_beams(pixwin=pixwin) - # TODO: load template_alpha in get_signal_shape - template_alpha90 = template_alpha["90"] - template_alpha150 = template_alpha["150a"] - X.log("Loading spectrum shape for transfer function...", "notice") - # load signal shapes for dust, cmb, and cmbdust - X.get_signal_shape( - #filename=signal_transfer_spec, transfer=True, - filename=[signal_transfer_spec, signal_transfer_spec_dust], - #tbeb=False, + # load signal shapes for dust and cmb + X.get_signal_shape( + filename=signal_transfer_spec, + filename_fg=foreground_transfer_spec, + freq_ref=freq_ref_transfer, + beta_ref=beta_ref_transfer, transfer=True, - transfer_dust=signal_transfer_type_dust is not None, - template_alpha90=template_alpha90, - template_alpha150=template_alpha150, - ) + ) - #X.log("Computing transfer functions for CMB...", "task") - #X.get_transfer( - # fix_bb_transfer=fix_bb_transfer, - # **transfer_opts) - dust = False - cmbdust = False - if (signal_transfer_spec_dust is not None): - if separate_dust_Fl: - dust = True - X.log("Computing transfer functions for CMB...", "task") - X.get_transfer( - fix_bb_transfer=fix_bb_transfer, - **transfer_opts, - ) - X.log("Computing transfer functions for Dust...", "task") - X.get_transfer( - fix_bb_transfer=fix_bb_transfer, - dust=dust, - template_alpha90=template_alpha90, - template_alpha150=template_alpha150, - **transfer_opts, - ) - else: - cmbdust=True - X.log("Computing transfer functions for Dust+CMB...", "task") - X.get_transfer( - cmbdust=cmbdust, - fix_bb_transfer=fix_bb_transfer, - template_alpha90=template_alpha90, - template_alpha150=template_alpha150, - **transfer_opts, - ) - else: - X.log("Computing transfer function...", "task") - X.get_transfer( - fix_bb_transfer=fix_bb_transfer, - **transfer_opts - ) + X.log("Computing transfer functions...", "notice") + X.get_transfer(**transfer_opts) X.log("Computing sim ensemble averages...", "notice") X.get_masked_sims(qb_file=qb_file_sim) @@ -769,13 +768,17 @@ def xfaster_run( X.get_signal_shape(flat=True) else: X.log("Loading spectrum shape for bandpowers...", "notice") - X.get_signal_shape(filename=signal_spec, r=model_r) + X.get_signal_shape( + filename=signal_spec, + r=model_r, + filename_fg=foreground_spec, + freq_ref=freq_ref, + beta_ref=beta_ref, + ) if multi_map: X.log("Computing multi-map bandpowers...", "notice") - qb, inv_fish = X.get_bandpowers(return_qb=True, - dust=dust, cmbdust=cmbdust, - **bandpwr_opts) + qb, inv_fish = X.get_bandpowers(return_qb=True, **bandpwr_opts) if likelihood: X.log("Computing multi-map likelihood...", "notice") @@ -1051,17 +1054,14 @@ def add_arg(P, name, argtype=None, default=None, short=None, help=None, **kwargs add_arg(G, "mask_type") add_arg(G, "signal_type") add_arg(G, "signal_subset") - #add_arg(G, "signal_transfer_type") - add_arg(G, "signal_transfer_type", - help="Signal sim variant for CMB transfer functions", - ) - add_arg(G, "signal_transfer_type_dust", - help="Signal sim variant for dust transfer functions", - ) + add_arg(G, "signal_transfer_type") add_arg(G, "signal_type_sim") add_arg(G, "noise_type") add_arg(G, "noise_subset") add_arg(G, "noise_type_sim") + add_arg(G, "foreground_type") + add_arg(G, "foreground_subset") + add_arg(G, "foreground_transfer_type") add_arg(G, "foreground_type_sim") add_arg(G, "template_type") add_arg(G, "template_noise_type") @@ -1086,6 +1086,7 @@ def add_arg(P, name, argtype=None, default=None, short=None, help=None, **kwargs metavar="SPEC", ) add_arg(G, "foreground_fit") + add_arg(G, "beta_fit") add_arg(G, "bin_width_fg", nargs="+") # data options @@ -1134,19 +1135,19 @@ def add_arg(P, name, argtype=None, default=None, short=None, help=None, **kwargs add_arg(G, "save_iters") add_arg(G, "return_cls") add_arg(G, "fix_bb_transfer") + add_arg(G, "fix_fg_transfer") add_arg(G, "null_first_cmb") add_arg(G, "qb_file_sim") E = G.add_mutually_exclusive_group() add_arg(E, "signal_spec") add_arg(E, "model_r") add_arg(G, "signal_transfer_spec") - add_arg(G, "signal_transfer_spec_dust", - help="Power spectrum used to create transfer signal simulations for dust. " - "Defaults to the spec_signal_{signal_transfer_type_dust}.dat file found " - "in the signal directory.", - ) - add_arg(G, "ref_freq") + add_arg(G, "foreground_spec") + add_arg(G, "foreground_transfer_spec") + add_arg(G, "freq_ref", argtype=float) add_arg(G, "beta_ref", argtype=float) + add_arg(G, "freq_ref_transfer", argtype=float) + add_arg(G, "beta_ref_transfer", argtype=float) add_arg(G, "delta_beta_prior", argtype=float) add_arg(G, "like_profiles") add_arg(G, "like_profile_sigma", argtype=float) @@ -1586,7 +1587,7 @@ def submit(self, group_by=None, verbose=True, **kwargs): group_by=group_by, serial=True, verbose=verbose, - **self.batch_args + **self.batch_args, ) self.reset() return job_ids From 3089e37788a29410fc4edc0d18181c1387088d23 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Tue, 12 Oct 2021 22:57:33 +1300 Subject: [PATCH 08/29] update tutorial notebook with foreground fitting changes --- docs/notebooks/XFaster_Tutorial.ipynb | 32 +++++++++++++++++---------- 1 file changed, 20 insertions(+), 12 deletions(-) diff --git a/docs/notebooks/XFaster_Tutorial.ipynb b/docs/notebooks/XFaster_Tutorial.ipynb index bc468d3e..33125504 100644 --- a/docs/notebooks/XFaster_Tutorial.ipynb +++ b/docs/notebooks/XFaster_Tutorial.ipynb @@ -13,7 +13,7 @@ "source": [ "This tutorial is also a Jupyter notebook, which can be found in [the example notebooks directory](https://github.com/annegambrel/xfaster/tree/main/example/notebooks). If you're running the notebook, rather than looking at the documentation produced from it, the links will not work. [Look at the docs for working links](https://annegambrel.github.io/xfaster/notebooks/XFaster_Tutorial.html).\n", "\n", - "The notebook reads intermediate output npz files from disk. To instead run the pieces starting from maps, first generate the example maps by running the script `xfaster/example/make_example_maps.py`, and set the checkpoint to something other than None." + "The notebook reads intermediate output npz files from disk. To instead run the pieces starting from maps, first generate the example maps by running the script `xfaster/example/make_example_maps.py`, and set the checkpoint to something other than `None`. This script generates sample signal (CMB), noise and foreground ensembles for the two SPIDER frequency bands, as well as two sets of simulated data maps (one with signal, noise and foreground components; and one with just signal and noise components)." ] }, { @@ -363,6 +363,7 @@ "Other options you might use here are \n", "\n", "* `template_type`, which points to files stored in `template_` for foreground template subtraction\n", + "* `foreground_type` and `foreground_transfer_type`, which points to files stored in `foreground_` and `foreground_`, to be used to construct the foreground signal component and its corresponding transfer function.\n", "* `data_root2`, which is used for null tests. This points to a second full map directory (data, signal and noise sims) for a set of data to be subtracted from the maps in data_root1.\n", "* `signal_type_sim`/`noise_type_sim`/`foreground_type_sim`/`template_type_sim` -- these are tags corresponding to directories `signal_`, `noise_`, `foreground_`, `template_` that are used when the `sim_data` argument is used in the `data` checkpoint. This mode chooses one of the sims for each each component in place of data maps. It adds the signal, noise, and/or foreground or template sims in alm-space. These options are not required to run in `sim_data` mode-- if they are not set, they default to `signal_type_sim=signal_type`, `noise_type_sim=noise_type`, `foreground_type_sim=None` and `template_type_sim=template_type`." ] @@ -417,7 +418,9 @@ "\n", "Residuals are fit per map, and by default are fit only for EE and BB, which are constrained to have the same fit parameter. To change this, use the option `res_specs`, which takes a list of the spectra you want to fit residuals for, ie. `[\"TT\", \"EE\", \"BB\"]` if you'd like to fit all of the spectra separately.\n", "\n", - "The last option you might wish to use is `weighted_bins`, which changes the default $\\chi_b(\\ell)$ binning operator from a tophat to one that weights by $\\ell(\\ell+1)$." + "Another option you might wish to use is `weighted_bins`, which changes the default $\\chi_b(\\ell)$ binning operator from a tophat to one that weights by $\\ell(\\ell+1)$.\n", + "\n", + "To enable fitting for a foreground component in the harmonic domain, set `foreground_fit=True`, and use the corresponding `bin_width_fg` to set the bin width for the foreground amplitude bins. Optionally, also set `beta_fit=True` to enable fitting for a frequency spectral index component for the foreground amplitude; otherwise, the frequency dependence is assumed to follow a single component dust model as defined in [scale_dust()](../api.rst#xfaster.spec_tools.scale_dust). If `foreground_fit` is True, use the `foreground_type` and `foreground_transfer_type` file options to enable the use of a particular foreground signal model and computation of a separate transfer function for the foreground component. If the file options are not set, the foreground model is assumed to be a simple power law model, as defined in [dust_model()](../api.rst#xfaster.spec_tools.dust_model), with the same transfer function as the CMB component." ] }, { @@ -570,7 +573,7 @@ "\n", "It only does `map2alm` once per map and caches the result for use in other cross spectra since this is the slowest step in the function.\n", "\n", - "Below, we call the function with `transfer=True`-- this tells the code to use the sims specified in `signal_type_transfer`, which defaults to being `signal_type` if the former is not specified." + "Below, we call the function with `transfer=True`-- this tells the code to use the sims specified in `signal_type_transfer`, which defaults to being `signal_type` if the former is not specified. This call also includes computing sims for the foreground component in the same way, if `foreground_type` and/or `foreground_transfer_type` are set." ] }, { @@ -589,7 +592,8 @@ "source": [ "The resulting outputs are:\n", "\n", - "* `cls_signal`: the average of the signal-only cross spectra, used to compute the transfer function\n", + "* `cls_signal`: the average of the signal-only cross spectra, used to compute the signal transfer function\n", + "* `cls_fg`: the average of the foreground-only cross spectra, used to compute the foreground transfer function\n", "* `cls_noise`: the average of the noise-only cross spectra, used as the noise model, $N_\\ell$\n", "* `cls_sim`: the average of the signal+noise spectra, where signal and noise maps are added in $a_{\\ell m}$s and thus the spectra include SxN terms.\n", "* `cls_med`: the median of the signal+noise spectra-- this is mainly used for debugging potential biases seen in the pipeline\n", @@ -734,7 +738,9 @@ "\n", "We have $K_{\\ell, \\ell'}$ and $B_{\\ell}$. For the transfer function calculation, we're going to set $F_\\ell$ to 1 so that we measure $q_b$s as the deviation from a uniform transfer function for our simulations. Binning, $\\chi_b$ has been chosen. All that's left is the full sky signal shape, $\\mathcal{C}_{\\ell'}^{XY (S)}$, loaded with [get_signal_shape()](../api.rst#xfaster.xfaster_class.XFaster.get_signal_shape). \n", "\n", - "For calculating the transfer function, this is just the shape spectrum that went into making our simulations. This spectrum can be specified by setting the argument `signal_transfer_spec` to a file containing the spectrum. If not provided, the code will look in the maps directory for the signal sims for a file labeled `spec_signal_.dat`. The file is expected to look like a CAMB output file, as demonstrated in `make_example_maps.py`, which writes such a file to the proper location in the signal sims directory." + "For calculating the transfer function, this is just the shape spectrum that went into making our simulations. This spectrum can be specified by setting the [xfaster_run()](../api.rst#xfaster.xfaster_exec.xfaster_run) argument `signal_transfer_spec` or `foreground_transfer_spec` to a file containing the spectrum for the signal or foreground component, respectively. If not provided, the code will look in the maps directory for the signal sims for a file labeled `spec_signal_.dat`, and similarly for foreground sims. The file is expected to look like a CAMB output file, as demonstrated in `make_example_maps.py`, which writes such a file to the proper location in the signal and foreground sims directories.\n", + "\n", + "For foreground fitting, this function also applies a frequency-dependent scaling to the input signal shape, using the [scale_dust()](../api.rst#xfaster.spec_tools.scale_dust) function. The arguments `freq_ref` and `beta_ref` are used to set the reference frequency and spectral index for the scaling function. Different reference points may be used for the transfer function computation, by setting the arguments `freq_ref_transfer` and `beta_ref_transfer`." ] }, { @@ -794,17 +800,17 @@ "The expression for the Fisher matrix does not change, other than the fact that its constituents are the same as detailed above for the transfer function.\n", "\n", "\n", - "Within the code, `get_transfer` basically has two steps within the function itself, which it performs per map. \n", + "Within the code, `get_transfer` basically has two steps within the function itself, which it performs per map and per signal component (either `\"cmb\"` or `\"fg\"`). \n", "\n", "1. Load up the $\\tilde{\\mathcal{C}}_{b\\ell}$: \n", "```python \n", - "cbl = self.bin_cl_template(map_tag=m0, transfer_run=True)\n", + "cbl = self.bin_cl_template(map_tag=m0, transfer_comp=\"cmb\")\n", "```\n", - "This uses the signal_shape internally that we calculated earlier, `m0` is the map, which is used to select the beam and kernel, and `transfer_run=True` sets the $F_\\ell$ term to 1.\n", + "This uses the signal_shape internally that we calculated earlier, `m0` is the map, which is used to select the beam and kernel, and `transfer_comp=\"cmb\"` sets the $F_\\ell$ term for the CMB component to 1.\n", "\n", "2. Run [fisher_iterate()](../api.rst#xfaster.xfaster_class.XFaster.fisher_iterate).\n", "```python\n", - "ret = self.fisher_iterate(cbl, m0, transfer_run=True,\n", + "ret = self.fisher_iterate(cbl, m0, transfer_comp=\"cmb\",\n", " iter_max=iter_max, converge_criteria=converge_criteria,\n", " save_iters=save_iters, ...)\n", "```\n", @@ -817,7 +823,9 @@ "```python\n", "qb['cmb_eb'] = np.sqrt(np.abs(qb['cmb_ee'] * qb['cmb_bb']))\n", "qb['cmb_tb'] = np.sqrt(np.abs(qb['cmb_tt'] * qb['cmb_bb']))\n", - "```" + "```\n", + "\n", + "When computing the transfer function for the foreground component, the TE transfer function is also handled similarly to the EB and TB terms." ] }, { @@ -987,7 +995,7 @@ "metadata": {}, "outputs": [], "source": [ - "cbl = X.bin_cl_template(map_tag=None, transfer_run=False, use_precalc=False)" + "cbl = X.bin_cl_template(map_tag=None, transfer_comp=None, use_precalc=False)" ] }, { @@ -1100,7 +1108,7 @@ "from collections import OrderedDict\n", "# construct a dummy qb array to compute input model\n", "qb = OrderedDict([(k, np.ones(len(v))) for k, v in X.bin_def.items()])\n", - "cls_model = X.get_model_spectra(qb, cbl, delta=True, cls_noise=X.cls_noise, cond_noise=1e-5)\n", + "cls_model = X.get_model_spectra(qb, cbl, cls_noise=X.cls_noise, cond_noise=1e-5)\n", "\n", "fig, axs = plt.subplots(2, 3, figsize=(11,7))\n", "axs = axs.ravel()\n", From 2177fdc9686da8e307f41d4e9473648e891eaee5 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Tue, 12 Oct 2021 23:02:36 +1300 Subject: [PATCH 09/29] example script for foreground fitting --- example/xfaster_example_fg.py | 69 +++++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) create mode 100644 example/xfaster_example_fg.py diff --git a/example/xfaster_example_fg.py b/example/xfaster_example_fg.py new file mode 100644 index 00000000..c5f20541 --- /dev/null +++ b/example/xfaster_example_fg.py @@ -0,0 +1,69 @@ +import xfaster as xf + +# use this script to run the xfaster algorithm with your own options +# like the example below + +submit = False # submit to queue +submit_opts = { + "wallt": 5, + "mem": 3, + "ppn": 1, + "omp_threads": 1, +} + +# Paths for the following keys are relative to where this script is run: +# config, data_root, data_root2, output_root +# If running this script from a different directory, make sure these +# keys point to the correct locations. + +xfaster_opts = { + # run options + "pol": True, + "pol_mask": True, + "bin_width": 25, + "lmax": 500, + "multi_map": True, + "likelihood": False, + "output_root": "outputs_example_fg", + "output_tag": "95x150", + # input files + "config": "config_example.ini", + "data_root": "maps_example", + "data_subset": "full/*95,full/*150", + "data_type": "cmbfg", + "noise_type": "gaussian", + "mask_type": "rectangle", + "signal_type": "synfast", + "foreground_type": "gaussian", + "data_root2": None, + "data_subset2": None, + # residual fitting + "residual_fit": True, + "bin_width_res": 100, + # foreground fitting + "foreground_fit": True, + "beta_fit": False, + "bin_width_fg": 40, + "beta_ref": 1.54, + "freq_ref": 359.7, + # spectrum + "ensemble_mean": False, + "tbeb": True, + "fix_fg_transfer": True, + "converge_criteria": 0.005, + "iter_max": 200, + "save_iters": False, + # likelihood + "like_lmin": 26, + "like_lmax": 250, + # beams + "pixwin": True, + "verbose": "info", + "checkpoint": "transfer", +} + +if submit: + xfaster_opts.update(**submit_opts) + xf.xfaster_submit(**xfaster_opts) +else: + xf.xfaster_run(**xfaster_opts) From 428f403b0c33fb8566c433ecc547f6fb2327911b Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Tue, 2 Nov 2021 22:04:36 +1300 Subject: [PATCH 10/29] keep track of model components with an attribute --- xfaster/xfaster_class.py | 140 ++++++++++++++++++--------------------- 1 file changed, 66 insertions(+), 74 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 97274d05..7206e9ec 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -3514,10 +3514,11 @@ def get_signal_shape( r : float If supplied and ``flat`` is False, a spectrum is computed using CAMB for the given ``r`` value. Overrides ``filename``. - component : 'scalar', 'tensor', 'fg' + component : 'scalar', 'tensor', 'cmb', 'fg' If 'scalar', and ``r`` is not None, return just the r=0 scalar terms in the signal model. If 'tensor', return just the tensor component - scaled by the input ``r`` value. If 'fg', return just fg terms + scaled by the input ``r`` value. If 'cmb' or 'fg', return just those + signal terms. flat : float If given, a spectrum that is flat in ell^2 Cl is returned, with amplitude given by the supplied value. Overrides all other options. @@ -3558,12 +3559,10 @@ def get_signal_shape( specs = list(self.specs) - comps = [] - if "cmb_tt" in self.bin_def or "cmb_ee" in self.bin_def: - comps += ["cmb"] - if not self.null_run and "fg_tt" in self.bin_def: - comps += ["fg"] - if component in ["cmb", "fg"]: + comps = list(self.signal_components) + if self.null_run: + comps.remove("fg") + if component in self.signal_components: comps = [component] if transfer: @@ -3977,6 +3976,9 @@ def get_bin_def( if self.pol and bin_width[1] != bin_width[2]: raise ValueError(bwerr) + comps = [] + signal_comps = [] + # Define bins nbins_cmb = 0 bin_def = OrderedDict() @@ -3985,6 +3987,8 @@ def get_bin_def( bins = np.append(bins, self.lmax + 1) bin_def["cmb_{}".format(spec)] = np.column_stack((bins[:-1], bins[1:])) nbins_cmb += len(bins) - 1 + comps += ["cmb"] + signal_comps += ["cmb"] self.log("Added {} CMB bins to bin_def".format(nbins_cmb), "debug") # Do the same for foreground bins @@ -4004,6 +4008,8 @@ def get_bin_def( nbins_fg += len(bins) - 1 if beta_fit: bin_def["delta_beta"] = np.array([[0, 0]]) + comps += ["fg"] + signal_comps += ["fg"] self.log( "Added {} foreground bins to bin_def".format(nbins_fg + int(beta_fit)), "debug", @@ -4035,6 +4041,7 @@ def get_bin_def( bin_def[btag] = np.column_stack((bins[:-1], bins[1:])) nbins_res += len(bins) - 1 + comps += ["res"] self.log("Added {} residual bins to bin_def".format(nbins_res), "debug") ell = np.arange(self.lmax + 1) @@ -4058,6 +4065,8 @@ def get_bin_def( self.specs = specs self.weighted_bins = weighted_bins self.bin_weights = bin_weights + self.components = comps + self.signal_components = signal_comps return self.bin_def @@ -4072,7 +4081,7 @@ def kernel_precalc(self, map_tag=None, component=None, transfer=False): If supplied, the kernels are computed only for the given map tag (or cross if map_tag is map_tag1:map_tag2). Otherwise, it is computed for all maps and crosses. - component : str + component : str or list If supplied, the matrix is computed only for the given component(s). Required if ``transfer=True``. transfer : bool @@ -4103,17 +4112,20 @@ def kernel_precalc(self, map_tag=None, component=None, transfer=False): lk = slice(0, lmax + 1) mll = OrderedDict() - all_comps = ["cmb", "fg"] if component is None: - component = all_comps - elif component not in all_comps: - raise ValueError( - "Invalid component, must be one of {}".format(", ".join(all_comps)) - ) + component = list(self.signal_components) + if not isinstance(component, list): + if component not in self.signal_components: + raise ValueError( + "Invalid component {}, must be one of {}".format( + component, ", ".join(self.signal_components) + ) + ) + component = [component] for stag in self.bin_def: comp, spec = stag.split("_", 1) - if comp not in ["cmb", "fg"]: + if comp not in self.signal_components: continue if transfer and comp not in component: @@ -4178,7 +4190,7 @@ def bin_cl_template( If supplied, the Cbl is computed only for the given map tag (or cross if map_tag is map_tag1:map_tag2). Otherwise, it is computed for all maps and crosses. - component : str + component : str or list If supplied, the Cbl is computed only for the given component(s). Required if ``transfer=True``. transfer : bool @@ -4221,10 +4233,9 @@ def bin_cl_template( specs = list(self.specs) if transfer: - assert component in [ - "cmb", - "fg", - ], "Argument `component` required if `transfer=True`" + assert ( + component in self.signal_components + ), "Argument `component` required if `transfer=True`" if "eb" in specs: specs.remove("eb") if "tb" in specs: @@ -4256,11 +4267,10 @@ def bin_cl_template( if component: comps = component if isinstance(component, list) else [component] else: - if "cmb_tt" in cls_shape or "cmb_ee" in cls_shape: - comps += ["cmb"] - if "fg_tt" in cls_shape or "fg_ee" in cls_shape: - comps += ["fg"] - if self.nbins_res > 0: + for comp in self.signal_components: + if any([k.startswith(comp + "_") for k in cls_shape]): + comps += [comp] + if "res" in self.components: comps += ["res"] cls_noise = self.cls_noise_null if self.null_run else self.cls_noise cls_res = self.cls_res_null if self.null_run else self.cls_res @@ -4430,10 +4440,9 @@ def get_model_spectra(self, qb, cbl, cls_noise=None, cond_noise=None): """ comps = [] - if any([k.startswith("cmb_") for k in qb]): - comps += ["cmb"] - if any([k.startswith("fg_") for k in qb]): - comps += ["fg"] + for comp in self.components: + if any([k.startswith(comp + "_") for k in qb]): + comps += [comp] delta_beta = 0.0 if "delta_beta" in qb: @@ -4441,9 +4450,6 @@ def get_model_spectra(self, qb, cbl, cls_noise=None, cond_noise=None): # in Fisher matrix, unless delta is True delta_beta = qb["delta_beta"][0] - if any([k.startswith("res_") for k in qb]): - comps += ["res"] - if cls_noise is not None: comps += ["noise"] @@ -4459,11 +4465,10 @@ def get_model_spectra(self, qb, cbl, cls_noise=None, cond_noise=None): specs = [] for spec in self.specs: - if "cmb_{}".format(spec) in cbl: - # Don't add entries that won't be filled in later - cls["total_{}".format(spec)] = OrderedDict() - elif "fg_{}".format(spec) in cbl: - cls["total_{}".format(spec)] = OrderedDict() + for comp in self.signal_components: + if "{}_{}".format(comp, spec) in cbl: + # Don't add entries that won't be filled in later + cls["total_{}".format(spec)] = OrderedDict() specs.append(spec) for comp in comps: @@ -4553,7 +4558,7 @@ def get_model_spectra(self, qb, cbl, cls_noise=None, cond_noise=None): cl1 /= 2.0 # compute model spectra - if comp in ["cmb", "fg"]: + if comp in self.signal_components: if xname not in cbl[stag]: continue cbl1 = cbl[stag][xname] @@ -4592,10 +4597,9 @@ def get_model_spectra(self, qb, cbl, cls_noise=None, cond_noise=None): cls.setdefault(stag, OrderedDict())[xname] = cl1 # add to total model - if not isinstance(cl1, dict): - ttag = "total_{}".format(spec) - cls[ttag].setdefault(xname, np.zeros_like(cl1)) - cls[ttag][xname] += cl1 + ttag = "total_{}".format(spec) + cls[ttag].setdefault(xname, np.zeros_like(cl1)) + cls[ttag][xname] += cl1 return cls @@ -4636,7 +4640,7 @@ def get_data_spectra(self, map_tag=None, transfer_comp=False, do_noise=True): map_pairs = pt.tag_pairs(map_tags) # select spectra - tbeb = "cmb_tb" in self.bin_def or "fg_tb" in self.bin_def + tbeb = "tb" in self.specs specs = list(self.specs) if transfer_comp or not tbeb: if "eb" in specs: @@ -4858,12 +4862,12 @@ def fisher_precalc( dim1 = pol_dim * num_maps comps = [] - if "cmb_tt" in cbl or "cmb_ee" in cbl: - comps += ["cmb"] - if "fg_tt" in cbl: - comps += ["fg"] - if "res_tt" in cbl or "res_ee" in cbl: - comps += ["res"] + sig_comps = [] + for comp in self.components: + if any([k.startswith(comp + "_") for k in cbl]): + comps += [comp] + if comp in self.signal_components: + sig_comps += [comp] specs = list(cls_input) @@ -4876,9 +4880,7 @@ def fisher_precalc( if windows: Dmat_obs = None Mmat = OrderedDict() - for comp in comps: - if comp not in ["cmb", "fg"]: - continue + for comp in sig_comps: for spec in specs: Mmat["{}_{}".format(comp, spec)] = OrderedDict() mll = getattr(self, "mll", None) @@ -4898,9 +4900,7 @@ def fisher_precalc( if likelihood: Dmat_obs_b[xname] = OrderedDict() elif windows: - for comp in comps: - if comp not in ["cmb", "fg"]: - continue + for comp in sig_comps: for spec in specs: Mmat["{}_{}".format(comp, spec)][xname] = OrderedDict() else: @@ -4911,9 +4911,7 @@ def fisher_precalc( # without bias subtraction for likelihood Dmat_obs_b[xname][spec] = cls_input[spec][xname] elif windows: - for comp in comps: - if comp not in ["cmb", "fg"]: - continue + for comp in sig_comps: stag = "{}_{}".format(comp, spec) Mmat[stag][xname][spec] = mll[stag][xname] if spec in ["ee", "bb"]: @@ -5207,9 +5205,7 @@ def fisher_calc( if not windows: Dmat_obs = pt.dict_to_dmat(Dmat_obs) # select only derivative terms for bins that are iterated - bin_def = OrderedDict( - [(k, v) for k, v in self.bin_def.items() if k in qb] - ) + bin_def = OrderedDict([(k, v) for k, v in self.bin_def.items() if k in qb]) dSdqb_mat1 = pt.dict_to_dsdqb_mat(dSdqb_mat1, bin_def) # apply ell limits @@ -5606,7 +5602,7 @@ def fisher_iterate( extra_tag=file_tag, ) - if "fg_tt" in self.bin_def: + if "fg" in self.components: out.update( map_freqs=self.map_freqs, freq_ref=self.freq_ref, @@ -5728,7 +5724,7 @@ def fisher_iterate( weighted_bins=self.weighted_bins, ) - if "fg_tt" in self.bin_def: + if "fg" in self.components: out.update( map_freqs=self.map_freqs, freq_ref=self.freq_ref, @@ -5768,7 +5764,7 @@ def fisher_iterate( self.log("Calculating final Fisher matrix without sample variance.", "info") qb_zeroed = copy.deepcopy(qb) qb_new_ns = copy.deepcopy(qb) - for comp in ["cmb", "fg"]: + for comp in self.signal_components: for spec in self.specs: stag = "{}_{}".format(comp, spec) if stag not in qb_zeroed: @@ -5964,11 +5960,7 @@ def get_transfer( Additionally the final output of ``fisher_iterate`` is stored in a dictionary called ``transfer_map`` for each map. """ - comps = [] - if "cmb_tt" in self.bin_def or "cmb_ee" in self.bin_def: - comps += ["cmb"] - if "fg_tt" in self.bin_def or "fg_ee" in self.bin_def: - comps += ["fg"] + comps = list(self.signal_components) opts = dict( converge_criteria=converge_criteria, @@ -5978,7 +5970,7 @@ def get_transfer( ) if "fg" in comps: - if self.cls_fg is None: + if getattr(self, "cls_fg") is None: fix_fg_transfer = True opts.update(fix_fg_transfer=fix_fg_transfer) @@ -6021,7 +6013,7 @@ def expand_transfer(qb_transfer, bin_def, check_only=False): if lq != lb: raise ValueError( "Found {} transfer bins for component {} map {}, expected {}".format( - lq, stag, m0, lb, + lq, stag, m0, lb ) ) if check_only: @@ -6253,7 +6245,7 @@ def get_bandpowers( ) # foreground fitting - if "fg_tt" in self.bin_def: + if "fg" in self.components: opts.update( freq_ref=self.freq_ref, beta_ref=self.beta_ref, @@ -6557,7 +6549,7 @@ def get_likelihood( else: qb = copy.deepcopy(qb) for spec in self.specs: - stags = ["cmb_{}".format(spec), "fg_{}".format(spec)] + stags = ["{}_{}".format(c, spec) for c in self.signal_components] for stag in stags: if stag not in qb: continue From bd73aec2d53819073e709dae07fe47284ae5d424 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Sun, 5 Dec 2021 02:45:58 +0200 Subject: [PATCH 11/29] fix bug in get likelihood --- xfaster/xfaster_class.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 7206e9ec..fdbbc5c6 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -776,6 +776,7 @@ def _get_sim_files( root1 = os.path.join(data_root, root) all_files = [] + print("DEBUG: ", root1, name, subset) for f in map_files: files = sorted( glob.glob( @@ -6569,7 +6570,8 @@ def get_likelihood( # load tensor and scalar terms separately # include all model components here - cbl_scalar = self.bin_cl_template(cls_shape_scalar, map_tag) + cbl_scalar = self.bin_cl_template( + cls_shape_scalar, map_tag, component="cmb") cls_model_scalar = self.get_model_spectra( qb, cbl_scalar, cls_noise=cls_noise ) From 92d8d8a9268d36d52d31edb6b8cf416b95d38bc1 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Fri, 7 Jan 2022 17:21:14 -0600 Subject: [PATCH 12/29] remove spurious print statement --- xfaster/xfaster_class.py | 1 - 1 file changed, 1 deletion(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index c0439617..320b6b9c 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -776,7 +776,6 @@ def _get_sim_files( root1 = os.path.join(data_root, root) all_files = [] - print("DEBUG: ", root1, name, subset) for f in map_files: files = sorted( glob.glob( From 4d78175e8db97daa6d3b25d763dccd5447e8910e Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Sun, 16 Jan 2022 21:23:35 -0600 Subject: [PATCH 13/29] make beta fit handling consistent with original code cleanup to simplify later merging --- docs/notebooks/XFaster_Tutorial.ipynb | 6 +-- example/xfaster_example_fg.py | 2 +- xfaster/parse_tools.py | 26 ++++++------ xfaster/xfaster_class.py | 58 +++++++++++++++------------ xfaster/xfaster_exec.py | 4 +- 5 files changed, 51 insertions(+), 45 deletions(-) diff --git a/docs/notebooks/XFaster_Tutorial.ipynb b/docs/notebooks/XFaster_Tutorial.ipynb index d4f36c58..d52003d1 100644 --- a/docs/notebooks/XFaster_Tutorial.ipynb +++ b/docs/notebooks/XFaster_Tutorial.ipynb @@ -806,9 +806,9 @@ "\n", "1. Load up the $\\tilde{\\mathcal{C}}_{b\\ell}$: \n", "```python \n", - "cbl = self.bin_cl_template(map_tag=m0, transfer_comp=\"cmb\")\n", + "cbl = self.bin_cl_template(map_tag=m0, transfer=True)\n", "```\n", - "This uses the signal_shape internally that we calculated earlier, `m0` is the map, which is used to select the beam and kernel, and `transfer_comp=\"cmb\"` sets the $F_\\ell$ term for the CMB component to 1.\n", + "This uses the signal_shape internally that we calculated earlier, `m0` is the map, which is used to select the beam and kernel, and `transfer=True` sets the $F_\\ell$ term for the CMB component to 1.\n", "\n", "2. Run [fisher_iterate()](../api.rst#xfaster.xfaster_class.XFaster.fisher_iterate).\n", "```python\n", @@ -997,7 +997,7 @@ "metadata": {}, "outputs": [], "source": [ - "cbl = X.bin_cl_template(map_tag=None, transfer_comp=None, use_precalc=False)" + "cbl = X.bin_cl_template(map_tag=None, transfer=False, use_precalc=False)" ] }, { diff --git a/example/xfaster_example_fg.py b/example/xfaster_example_fg.py index c5f20541..fdaeb0f3 100644 --- a/example/xfaster_example_fg.py +++ b/example/xfaster_example_fg.py @@ -59,7 +59,7 @@ # beams "pixwin": True, "verbose": "info", - "checkpoint": "transfer", + "checkpoint": None, } if submit: diff --git a/xfaster/parse_tools.py b/xfaster/parse_tools.py index cdcf8547..e2ed780b 100644 --- a/xfaster/parse_tools.py +++ b/xfaster/parse_tools.py @@ -424,7 +424,18 @@ def load_and_parse(filename, check_version=True): if "fix_bb_xfer" in data: data["fix_bb_transfer"] = data.pop("fix_bb_xfer") - if version >= 1: + # update data version in memory + data["data_version"] = dv + + if version == 2: + + if "reference_root" in data and "reference_type" not in data: + data["reference_type"] = None if data["reference_root"] is None else "sub" + + # update data version in memory + data["data_version"] = dv + + if version in [1, 2]: if "ref_freq" in data: data["freq_ref"] = data.pop("ref_freq") @@ -440,19 +451,6 @@ def load_and_parse(filename, check_version=True): if "cls_sim" in data and "cls_fg" not in data: data["cls_fg"] = None - # update data version in memory - data["data_version"] = dv - - if version == 2: - - if "reference_root" in data and "reference_type" not in data: - data["reference_type"] = None if data["reference_root"] is None else "sub" - - # update data version in memory - data["data_version"] = dv - - if version in [1, 2]: - if "map_files" in data: data["map_names"] = np.asarray( [os.path.relpath(f, data["map_root"]) for f in data["map_files"]] diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 99ce2b96..88e784ab 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -2987,7 +2987,6 @@ def get_masked_sims(self, transfer=False, do_noise=True, qb_file=None): shape=data_shape, shape_ref="cls_signal", ) - if ret is not None: if do_noise: if self.cls_noise is not None: @@ -3529,7 +3528,8 @@ def get_signal_shape( search for a spectrum stored in ``foreground_/spec_foreground_.dat``. Otherwise, if the filename is a relative path and not found, the - config directory is searched. + config directory is searched. If not supplied, a power law dust + model spectrum is assumed. freq_ref : float In GHz, reference frequency for dust model. Dust bandpowers output will be at this reference frequency. @@ -3557,15 +3557,18 @@ def get_signal_shape( Dictionary keyed by spectrum (cmb_tt, cmb_ee, ... , fg_tt, ...), each entry containing a vector of length 2 * lmax + 1. """ + lmax_kern = 2 * self.lmax specs = list(self.specs) comps = list(self.signal_components) - if self.null_run: + if self.null_run and "fg" in comps: comps.remove("fg") if component in self.signal_components: comps = [component] + elif component in ["scalar", "tensor"]: + comps = ["cmb"] if transfer: if "eb" in specs: @@ -4272,7 +4275,7 @@ def bin_cl_template( for comp in self.signal_components: if any([k.startswith(comp + "_") for k in cls_shape]): comps += [comp] - if "res" in self.components: + if "res" in self.components and not transfer: comps += ["res"] cls_noise = self.cls_noise_null if self.null_run else self.cls_noise cls_res = self.cls_res_null if self.null_run else self.cls_res @@ -4501,7 +4504,8 @@ def get_model_spectra(self, qb, cbl, cls_noise=None, cond_noise=None): elif comp == "fg": # beta scaling for foreground model - beta_scale = self.fg_scales[xname][1] ** delta_beta + # beta_scale = self.fg_scales[xname][1] ** delta_beta + beta_scale = 1 + self.fg_scales[xname][2] * delta_beta qbs = beta_scale * qb[stag] if spec in ["ee", "bb"]: qbm = beta_scale * qb[mstag] @@ -4702,7 +4706,7 @@ def get_data_spectra(self, map_tag=None, transfer_comp=False, do_noise=True): debias[spec][xname] = np.copy(self.cls_sim_null[spec][xname]) # Non-nulls are debiased by average of N sims - elif self.cls_noise is not None: + elif not self.null_run and self.cls_noise is not None: for spec in specs: nell[spec] = OrderedDict() debias[spec] = OrderedDict() @@ -4715,6 +4719,9 @@ def get_data_spectra(self, map_tag=None, transfer_comp=False, do_noise=True): nell[spec][xname] = np.copy(self.cls_noise[spec][xname]) debias[spec][xname] = np.copy(nell[spec][xname]) + else: + nell = debias = None + return obs, nell, debias def do_qb2cb(self, qb, inv_fish, wbl): @@ -5137,17 +5144,18 @@ def fisher_calc( # get foreground at pivot point spectral index _, beta_scale, log_beta_scale = self.fg_scales[xname] # separable beta correction - beta_scale = beta_scale ** delta_beta # non-linear - # beta_scale = (1 + log_beta_scale * delta_beta) # linear + # beta_scale = beta_scale ** delta_beta # non-linear + beta_scale = 1 + log_beta_scale * delta_beta # linear + beta_scale0 = log_beta_scale / beta_scale # scale foreground model by frequency scaling adjusted for beta for s1, sdat in dSdqb_mat1["fg"][xname].items(): for s2, sdat2 in sdat.items(): sdat2 *= beta_scale - # with linearized derivative term + # with linearized derivative term, evaluated at input beta dSdqb_mat1["delta_beta"][xname][s1] = OrderedDict( - {s1: cls_model["fg_{}".format(s1)][xname] * log_beta_scale} + {s1: cls_model["fg_{}".format(s1)][xname] * beta_scale0} ) # Set up Dmat -- if it's not well conditioned, add noise to the @@ -6022,7 +6030,7 @@ def expand_transfer(qb_transfer, bin_def, check_only=False): continue if spec == "bb" and fix_bb_transfer: v = transfer["{}_ee".format(comp)][m0] - elif spec in ["te", "tb", "eb"] and spec not in qb_transfer: + elif spec in ["te", "tb", "eb"] and stag not in qb_transfer: staga, stagb = ["{}_{}{}".format(comp, s, s) for s in spec] v = np.sqrt( np.abs(transfer[staga][m0] * transfer[stagb][m0]) @@ -6252,10 +6260,10 @@ def get_bandpowers( freq_ref=self.freq_ref, beta_ref=self.beta_ref, ) - if "delta_beta" in self.bin_def: - # priors on frequency spectral index - self.delta_beta_fix = 1.0e-8 - opts.update(delta_beta_prior=delta_beta_prior) + if "delta_beta" in self.bin_def: + # priors on frequency spectral index + self.delta_beta_fix = 1.0e-8 + opts.update(delta_beta_prior=delta_beta_prior) if self.template_cleaned: opts.update(template_alpha=self.template_alpha) @@ -6564,31 +6572,30 @@ def get_likelihood( cls_shape_scalar = self.get_signal_shape( r=1.0, save=False, component="scalar" ) + if "fg" in self.signal_components: + cls_shape_scalar.update( + {k: cls for k, cls in self.cls_shape.items() if k.startswith("fg_")} + ) cls_shape_tensor = self.get_signal_shape( r=1.0, save=False, component="tensor" ) # load tensor and scalar terms separately - # include all model components here - cbl_scalar = self.bin_cl_template( - cls_shape_scalar, map_tag, component="cmb") + cbl_scalar = self.bin_cl_template(cls_shape_scalar, map_tag) cls_model_scalar = self.get_model_spectra( qb, cbl_scalar, cls_noise=cls_noise ) - # include only CMB components for all additive terms below - cbl_tensor = self.bin_cl_template( - cls_shape_tensor, map_tag, component="cmb" - ) + cbl_tensor = self.bin_cl_template(cls_shape_tensor, map_tag) cls_model_tensor = self.get_model_spectra(qb_cmb, cbl_tensor) if beam_prior is not None: # load beam error term for tensor and scalar cbl_scalar_beam = self.bin_cl_template( - cls_shape_scalar, map_tag, component="cmb", beam_error=True + cls_shape_scalar, map_tag, beam_error=True ) - cls_mod_scal_beam = self.get_model_spectra(qb_cmb, cbl_scalar_beam) + cls_mod_scal_beam = self.get_model_spectra(qb, cbl_scalar_beam) cbl_tensor_beam = self.bin_cl_template( - cls_shape_tensor, map_tag, component="cmb", beam_error=True + cls_shape_tensor, map_tag, beam_error=True ) cls_mod_tens_beam = self.get_model_spectra(qb_cmb, cbl_tensor_beam) @@ -6690,6 +6697,7 @@ def log_like( beam_coeffs[xname] = d # compute new signal shape by scaling tensor component by r + # XXX handle beam fg terms here too if r is not None: for stag, d in cls_model0.items(): comp, spec = stag.split("_", 1) diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index 8fd2155c..eb2da599 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -657,11 +657,11 @@ def xfaster_run( spec_opts.pop("signal_transfer_spec") spec_opts.pop("foreground_spec") spec_opts.pop("foreground_transfer_spec") + spec_opts.pop("model_r") spec_opts.pop("freq_ref") spec_opts.pop("beta_ref") spec_opts.pop("freq_ref_transfer") spec_opts.pop("beta_ref_transfer") - spec_opts.pop("model_r") spec_opts.pop("qb_file") bandpwr_opts = spec_opts.copy() bandpwr_opts.pop("fix_bb_transfer") @@ -1579,7 +1579,7 @@ def submit(self, group_by=None, verbose=True, **kwargs): group_by=group_by, serial=True, verbose=verbose, - **self.batch_args, + **self.batch_args ) self.reset() return job_ids From 7e354abaf8987d5d897a5ae0c2b5b8ad403ac2d4 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 19 Jan 2022 00:28:44 -0600 Subject: [PATCH 14/29] bug fixes --- example/make_example_maps.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/example/make_example_maps.py b/example/make_example_maps.py index 982ee3c3..f11c4420 100644 --- a/example/make_example_maps.py +++ b/example/make_example_maps.py @@ -125,7 +125,7 @@ def sim_signal(isim, itag): def sim_noise(isim, itag): np.random.seed(nsim * (1 + itag) + isim) - noise = np.tile((3, 1), np.zeros(mask.size, dtype=float)) + noise = np.tile(np.zeros(mask.size, dtype=float), (3, 1)) noise[0][mask] = np.random.normal(scale=gnoise[itag], size=np.sum(mask)) noise[1][mask] = np.random.normal( scale=gnoise[itag] * np.sqrt(2), size=np.sum(mask) # pol @@ -155,7 +155,7 @@ def sim_fg(isim, itag): data_fg_file = os.path.join(data_fg_path, "map_{}.fits".format(tag)) data = None if not os.path.exists(data_file): - data = sim_sig(i, ti) + data_noise_scale * sim_noise(i, ti) + data = sim_signal(i, ti) + data_noise_scale * sim_noise(i, ti) write_map(data_file, data, mask) print("Wrote data signal+noise map {}".format(tag)) if do_fg and not os.path.exists(data_fg_file): @@ -172,7 +172,7 @@ def sim_fg(isim, itag): comps = [] if not os.path.exists(sig_file): - write_map(sig_file, sim_sig(i, ti), mask) + write_map(sig_file, sim_signal(i, ti), mask) comps += ["sig"] if not os.path.exists(noise_file): write_map(noise_file, sim_noise(i, ti), mask) From e133be946023d72a0a976018f93846245c7f2083 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 19 Jan 2022 00:30:29 -0600 Subject: [PATCH 15/29] frequency scaling in foreground window functions --- xfaster/xfaster_class.py | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 88e784ab..fdb387aa 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -4502,7 +4502,7 @@ def get_model_spectra(self, qb, cbl, cls_noise=None, cond_noise=None): if spec in ["ee", "bb"]: qbm = qb[mstag] - elif comp == "fg": + elif comp == "fg" and delta_beta != 0.0: # beta scaling for foreground model # beta_scale = self.fg_scales[xname][1] ** delta_beta beta_scale = 1 + self.fg_scales[xname][2] * delta_beta @@ -5280,7 +5280,7 @@ def fisher_calc( # include priors in likelihood if "delta_beta" in qb and delta_beta_prior is not None: - chi = (qb["delta_beta"] - self.delta_beta_fix) / delta_beta_prior + chi = (delta_beta - self.delta_beta_fix) / delta_beta_prior like -= chi ** 2 / 2.0 if null_first_cmb: @@ -5317,6 +5317,18 @@ def fisher_calc( wbl = OrderedDict() bin_index = pt.dict_to_index(qb) + # compute foreground scaling per frequency + if "fg" in self.signal_components: + fmat = OrderedDict() + for xname, (a, b, c) in self.fg_scales.items(): + fmat[xname] = OrderedDict() + amp = a + if "delta_beta" in qb: + amp *= 1.0 + delta_beta * c + # amp *= b ** delta_beta + fmat[xname] = OrderedDict([(s, amp) for s in self.specs]) + fmat = pt.dict_to_dmat(fmat, pol=self.pol) + # compute window functions for each spectrum for k, (left, right) in bin_index.items(): if k not in Mmat: @@ -5329,6 +5341,10 @@ def fisher_calc( # construct Mll' matrix marg = pt.dict_to_dmat(Mmat[k], pol=self.pol)[..., ell, :] + # apply frequency scaling + if comp == "fg": + marg = np.einsum("ij,ijlm->ijlm", fmat, marg) + # qb window function wbl1 = np.einsum("iil,ijkl,jilm->km", gmat, sarg, marg) / 2.0 / norm del marg From 303b1524b56ce2c21a438f4af28e763da91a468e Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Fri, 8 Jul 2022 23:33:13 -0500 Subject: [PATCH 16/29] fix filename --- xfaster/xfaster_class.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 0bf9edc8..268930dc 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -5500,7 +5500,7 @@ def fisher_iterate( if transfer_comp: null_first_cmb = False - save_name = "transfer_{}wbins" if self.weighted_bins else "transfer_{}" + save_name = "transfer_{}_wbins" if self.weighted_bins else "transfer_{}" save_name = save_name.format(transfer_comp) else: save_name = "bandpowers" From e42e672a67918b0c7ec7104dd389452616621cf2 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 20 Jul 2022 17:53:26 -0500 Subject: [PATCH 17/29] no need for backward compatibility in file parsing with new file handling --- xfaster/parse_tools.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/xfaster/parse_tools.py b/xfaster/parse_tools.py index 0f54ce8d..2cb1421b 100644 --- a/xfaster/parse_tools.py +++ b/xfaster/parse_tools.py @@ -443,17 +443,6 @@ def load_and_parse(filename, check_version=True): if "ref_freq" in data: data["freq_ref"] = data.pop("ref_freq") - if "foreground_type_sim" in data and "foreground_type" not in data: - for k in ["type", "root", "files", "root2", "files2"]: - data["foreground_" + k] = None - data["foreground_transfer_" + k] = None - data["num_foreground"] = 0 - data["num_foreground_transfer"] = 0 - data["foreground_subset"] = "*" - - if "cls_sim" in data and "cls_fg" not in data: - data["cls_fg"] = None - if "map_files" in data: data["map_names"] = np.asarray( [os.path.relpath(f, data["map_root"]) for f in data["map_files"]] From b8665d48d8b10e2bedce84ec969eb4d7fb117231 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 20 Jul 2022 19:08:21 -0500 Subject: [PATCH 18/29] no foreground component in nulls --- xfaster/xfaster_class.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index b56ce168..d6fe3ada 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -2901,9 +2901,6 @@ def get_masked_sims( "num_foreground", ] save_attrs += ["cls_fg"] - if null_run: - file_attrs += ["foreground_root2", "foreground_files2"] - save_attrs += ["cls_fg_null"] save_attrs += file_attrs From 5716f79afa076ffdb2f98078f376603fffef1bab Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 20 Jul 2022 19:20:31 -0500 Subject: [PATCH 19/29] missing variables --- xfaster/xfaster_class.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index d6fe3ada..f1561aca 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -2953,6 +2953,10 @@ def get_masked_sims( noise_files2 = self.noise_files2 if null_run else None num_noise = self.num_noise + if do_fg: + foreground_files = self.foreground_files + num_foreground = self.num_foreground + # process signal, foreground, noise, and S+N cls_sig = OrderedDict() cls_fg = OrderedDict() if do_fg else None From 0fec528b347b6552130fe8d41304880b8355b826 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 20 Jul 2022 19:49:27 -0500 Subject: [PATCH 20/29] simplify transfer component handling --- docs/notebooks/XFaster_Tutorial.ipynb | 2 +- xfaster/xfaster_class.py | 151 +++++++++++++------------- xfaster/xfaster_exec.py | 5 +- 3 files changed, 80 insertions(+), 78 deletions(-) diff --git a/docs/notebooks/XFaster_Tutorial.ipynb b/docs/notebooks/XFaster_Tutorial.ipynb index 8e6d11d9..4ce89c73 100644 --- a/docs/notebooks/XFaster_Tutorial.ipynb +++ b/docs/notebooks/XFaster_Tutorial.ipynb @@ -792,7 +792,7 @@ "\n", "2. Run [fisher_iterate()](../api.rst#xfaster.xfaster_class.XFaster.fisher_iterate).\n", "```python\n", - "ret = self.fisher_iterate(cbl, m0, transfer_comp=\"cmb\",\n", + "ret = self.fisher_iterate(cbl, m0, transfer=True,\n", " iter_max=iter_max, converge_criteria=converge_criteria,\n", " save_iters=save_iters, ...)\n", "```\n", diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index f1561aca..aec31dcf 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -2892,7 +2892,9 @@ def get_masked_sims( do_fg = foreground_type is not None if do_fg: - opts.update(foreground_type=foreground_type, foreground_subset=foreground_subset) + opts.update( + foreground_type=foreground_type, foreground_subset=foreground_subset + ) file_attrs += [ "foreground_type", "foreground_subset", @@ -2916,7 +2918,10 @@ def get_masked_sims( # don't force rerun the signal sims if previously loaded sims are of the # same type if hasattr(self, "signal_type") and self.signal_type == signal_type: - if not hasattr(self, "foreground_type") or self.foreground_type == foreground_type: + if ( + not hasattr(self, "foreground_type") + or self.foreground_type == foreground_type + ): self.force_rerun[cp] = False ret = self.load_data( @@ -4056,7 +4061,7 @@ def get_bin_def( return self.bin_def - def kernel_precalc(self, map_tag=None, component=None, transfer=False): + def kernel_precalc(self, map_tag=None, transfer=False): """ Compute the mixing kernels M_ll' = K_ll' * F_l' * B_l'^2. Called by ``bin_cl_template`` to pre-compute kernel terms. @@ -4067,12 +4072,11 @@ def kernel_precalc(self, map_tag=None, component=None, transfer=False): If supplied, the kernels are computed only for the given map tag (or cross if map_tag is map_tag1:map_tag2). Otherwise, it is computed for all maps and crosses. - component : str or list - If supplied, the matrix is computed only for the given component(s). - Required if ``transfer=True``. - transfer : bool - If True, assumes a unity transfer function for all bins for the - given component(s). + transfer : bool or str + If supplied as a signal component name, assumes a unity transfer + function for all bins for the given component. If True, the + ``"cmb"`` component is assumed. If False, the pre-computed + transfer function is applied. Returns ------- @@ -4098,23 +4102,23 @@ def kernel_precalc(self, map_tag=None, component=None, transfer=False): lk = slice(0, lmax + 1) mll = OrderedDict() - if component is None: - component = list(self.signal_components) - if not isinstance(component, list): - if component not in self.signal_components: + comps = [] + if transfer: + if transfer is True: + transfer = "cmb" + if transfer not in self.signal_components: raise ValueError( - "Invalid component {}, must be one of {}".format( + "Invalid transfer component {}, must be one of {}".format( component, ", ".join(self.signal_components) ) ) - component = [component] + comps = [transfer] + else: + comps = list(self.signal_components) for stag in self.bin_def: comp, spec = stag.split("_", 1) - if comp not in self.signal_components: - continue - - if transfer and comp not in component: + if comp not in comps: continue mll[stag] = OrderedDict() @@ -4156,7 +4160,6 @@ def bin_cl_template( self, cls_shape=None, map_tag=None, - component=None, transfer=False, beam_error=False, use_precalc=True, @@ -4176,14 +4179,12 @@ def bin_cl_template( If supplied, the Cbl is computed only for the given map tag (or cross if map_tag is map_tag1:map_tag2). Otherwise, it is computed for all maps and crosses. - component : str or list - If supplied, the Cbl is computed only for the given component(s). - Required if ``transfer=True``. - transfer : bool - If True, assumes a unity transfer function for all bins for the - given component, and the output Cbl is used to compute the transfer - functions for that component, which are then loaded when this method - is called with ``transfer=False``. + transfer : bool or str + If supplied, assumes a unity transfer function for all bins for the + given signal component, and the output Cbl is used to compute the + transfer functions for that component, which are then loaded when + this method is called with ``transfer=False``. If True, the + ``"cmb"`` component is assumed. beam_error : bool If True, use beam error envelope instead of beam to get cbls that are 1 sigma beam error envelope offset of signal terms. @@ -4219,23 +4220,26 @@ def bin_cl_template( specs = list(self.specs) if transfer: - assert ( - component in self.signal_components - ), "Argument `component` required if `transfer=True`" + if transfer is True: + transfer = "cmb" + if transfer not in self.signal_components: + raise ValueError( + "Invalid transfer component {}, must be one of {}".format( + component, ", ".join(self.signal_components) + ) + ) if "eb" in specs: specs.remove("eb") if "tb" in specs: specs.remove("tb") - if component == "fg" and "te" in specs: + if transfer == "fg" and "te" in specs: specs.remove("te") lmax = self.lmax lmax_kern = lmax # 2 * self.lmax if getattr(self, "mll", None) is None or not use_precalc: - mll = self.kernel_precalc( - map_tag=map_tag, component=component, transfer=transfer - ) + mll = self.kernel_precalc(map_tag=map_tag, transfer=transfer) if use_precalc: self.mll = mll else: @@ -4250,13 +4254,13 @@ def bin_cl_template( cbl = OrderedDict() comps = [] - if component: - comps = component if isinstance(component, list) else [component] + if transfer: + comps = [transfer] else: for comp in self.signal_components: if any([k.startswith(comp + "_") for k in cls_shape]): comps += [comp] - if "res" in self.components and not transfer: + if "res" in self.components: comps += ["res"] cls_noise = self.cls_noise_null if self.null_run else self.cls_noise cls_res = self.cls_res_null if self.null_run else self.cls_res @@ -4591,7 +4595,7 @@ def get_model_spectra(self, qb, cbl, cls_noise=None, cond_noise=None): return cls - def get_data_spectra(self, map_tag=None, transfer_comp=False, do_noise=True): + def get_data_spectra(self, map_tag=None, transfer=False, do_noise=True): """ Return data and noise spectra for the given map tag(s). Data spectra and signal/noise sim spectra must have been precomputed or loaded from @@ -4602,7 +4606,7 @@ def get_data_spectra(self, map_tag=None, transfer_comp=False, do_noise=True): map_tag : str If None, all map-map cross-spectra are included in the outputs. Otherwise, only the autospectra of the given map are included. - transfer_comp : str + transfer : bool or str If supplied, the data cls are the average of the simulations for the given component, and noise cls are ignored. If False, the data cls are either ``cls_data_null`` (for null tests) or ``cls_data``. See @@ -4616,10 +4620,10 @@ def get_data_spectra(self, map_tag=None, transfer_comp=False, do_noise=True): obs : OrderedDict Dictionary of data cross spectra nell : OrderedDict - Dictionary of noise cross spectra, or None if ``transfer_comp`` is set. + Dictionary of noise cross spectra, or None if ``transfer`` is set. debias : OrderedDict Dictionary of debiased data cross spectra, or None if - ``transfer_comp`` is set. + ``transfer`` is set. """ # select map pairs if map_tag is not None: @@ -4631,23 +4635,23 @@ def get_data_spectra(self, map_tag=None, transfer_comp=False, do_noise=True): # select spectra tbeb = "tb" in self.specs specs = list(self.specs) - if transfer_comp or not tbeb: + if transfer or not tbeb: if "eb" in specs: specs.remove("eb") if "tb" in specs: specs.remove("tb") - if transfer_comp == "fg" and "te" in specs: + if transfer == "fg" and "te" in specs: specs.remove("te") # obs depends on what you're computing - if transfer_comp: - if transfer_comp == "cmb": + if transfer: + if transfer is True or transfer == "cmb": obs_quant = self.cls_signal - elif transfer_comp == "fg": + elif transfer == "fg": obs_quant = self.cls_fg else: raise ValueError( - "Unknown transfer function component {}".format(transfer_comp) + "Unknown transfer function component {}".format(transfer) ) elif self.null_run: if self.reference_type is not None: @@ -4669,7 +4673,7 @@ def get_data_spectra(self, map_tag=None, transfer_comp=False, do_noise=True): if not do_noise: return obs - elif transfer_comp: + elif transfer: return obs, None, None nell = OrderedDict() @@ -5396,7 +5400,7 @@ def fisher_iterate( iter_max=200, converge_criteria=0.005, qb_start=None, - transfer_comp=None, + transfer=False, save_iters=False, null_first_cmb=False, delta_beta_prior=None, @@ -5430,7 +5434,7 @@ def fisher_iterate( qb_start : OrderedDict Initial guess at ``qb`` bandpower amplitudes. If None, unity is assumed for all bins. - transfer_comp : str + transfer : bool or str If set, the input Cls passed to ``fisher_calc`` are the average of the component simulations, and noise cls are ignored. If False, the input Cls are either ``cls_data_null`` @@ -5485,7 +5489,7 @@ def fisher_iterate( cls_shape : shape spectrum iters : number of iterations completed - If ``transfer_comp`` is not set, this dictionary also contains:: + If ``transfer`` is not set, this dictionary also contains:: qb_transfer : transfer function amplitudes transfer : ell-by-ell transfer function @@ -5495,16 +5499,19 @@ def fisher_iterate( Notes ----- This method stores outputs to files with name 'transfer' for transfer - function runs (if ``transfer_comp`` is set), otherwise with name + function runs (if ``transfer`` is set), otherwise with name 'bandpowers'. Outputs from each individual iteration, containing only the quantities that change with each step, are stored in separate files with the same name (and different index). """ - if transfer_comp: + if transfer: null_first_cmb = False save_name = "transfer_{}_wbins" if self.weighted_bins else "transfer_{}" - save_name = save_name.format(transfer_comp) + if transfer is True: + transfer = "cmb" + assert transfer in self.signal_components + save_name = save_name.format(transfer) else: save_name = "bandpowers" @@ -5516,12 +5523,12 @@ def fisher_iterate( qb = OrderedDict() for k, v in self.bin_def.items(): # transfer function for signal component - if transfer_comp: + if transfer: if "eb" in k or "tb" in k: continue - if transfer_comp == "fg" and "te" in k: + if transfer == "fg" and "te" in k: continue - if not k.startswith(transfer_comp): + if not k.startswith(transfer): continue qb[k] = np.ones(len(v)) continue @@ -5539,9 +5546,7 @@ def fisher_iterate( else: qb = qb_start - obs, nell, debias = self.get_data_spectra( - map_tag=map_tag, transfer_comp=transfer_comp - ) + obs, nell, debias = self.get_data_spectra(map_tag=map_tag, transfer=transfer) bin_index = pt.dict_to_index(self.bin_def) @@ -5635,7 +5640,7 @@ def fisher_iterate( beta_err=beta_err, ) - self.save_data(save_name, bp_opts=not transfer_comp, **out) + self.save_data(save_name, bp_opts=not transfer, **out) (nans,) = np.where(np.isnan(qb_new_arr)) if len(nans): @@ -5658,7 +5663,7 @@ def fisher_iterate( ) if np.nanmax(np.abs(fqb)) < converge_criteria: - if not transfer_comp: + if not transfer: # Calculate final fisher matrix without conditioning self.log("Calculating final Fisher matrix.", "info") _, inv_fish = self.fisher_calc( @@ -5687,15 +5692,13 @@ def fisher_iterate( else: msg = "{} {} did not converge in {} iterations".format( "Multi-map" if map_tag is None else "Map {}".format(map_tag), - "{} transfer function".format(transfer_comp) - if transfer_comp - else "spectrum", + "{} transfer function".format(transfer) if transfer else "spectrum", iter_max, ) # Check the slope of the last ten fqb_maxpoints. # If there's not a downward trend, adjust conditioning # criteria to help convergence. - if len(prev_fqb) <= 10 or transfer_comp: + if len(prev_fqb) <= 10 or transfer: continue m, b = np.polyfit(np.arange(10), prev_fqb[-10:], 1) if m > 0: # Not converging @@ -5770,7 +5773,7 @@ def fisher_iterate( Dmat_obs=self.Dmat_obs, ) - if not transfer_comp: + if not transfer: out.update( qb_transfer=self.qb_transfer, transfer=self.transfer, @@ -5778,7 +5781,7 @@ def fisher_iterate( if self.template_type is not None: out.update(template_alpha=self.template_alpha) - if success and not transfer_comp: + if success and not transfer: # do one more fisher calc that doesn't include sample variance # set qb=very close to 0. 0 causes singular matrix problems. # don't do this for noise residual bins @@ -5918,11 +5921,7 @@ def fisher_iterate( self.warn(msg) return self.save_data( - save_name, - map_tag=map_tag, - bp_opts=not transfer_comp, - extra_tag=file_tag, - **out + save_name, map_tag=map_tag, bp_opts=not transfer, extra_tag=file_tag, **out ) def get_transfer( @@ -6103,11 +6102,11 @@ def expand_transfer(qb_transfer, bin_def, check_only=False): "info", ) self.clear_precalc() - cbl = self.bin_cl_template(map_tag=m0, component=comp, transfer=True) + cbl = self.bin_cl_template(map_tag=m0, transfer=comp) ret = self.fisher_iterate( cbl, m0, - transfer_comp=comp, + transfer=comp, iter_max=iter_max, converge_criteria=converge_criteria, save_iters=save_iters, diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index 4ea31153..1faf0266 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -509,7 +509,10 @@ def xfaster_run( all_opts["foreground_transfer_type"] = foreground_transfer_type # make sure sims get rerun correctly - if signal_transfer_type == signal_type and foreground_transfer_type == foreground_type: + if ( + signal_transfer_type == signal_type + and foreground_transfer_type == foreground_type + ): # if signal types match, then sims are run before computing the # transfer function, so need to set the correct checkpoint to rerun if checkpoint == "sims": From 619c99e1a010cf4d1b23935d92f0bac08527dcbb Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 20 Jul 2022 21:21:57 -0500 Subject: [PATCH 21/29] cleanup --- xfaster/xfaster_class.py | 31 ++++++++++++++----------------- 1 file changed, 14 insertions(+), 17 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index aec31dcf..e8fc44a3 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -2890,7 +2890,7 @@ def get_masked_sims( file_attrs += ["noise_root2", "noise_files2"] save_attrs += ["cls_noise_null", "cls_res_null"] - do_fg = foreground_type is not None + do_fg = not null_run and foreground_type is not None if do_fg: opts.update( foreground_type=foreground_type, foreground_subset=foreground_subset @@ -2959,8 +2959,8 @@ def get_masked_sims( num_noise = self.num_noise if do_fg: - foreground_files = self.foreground_files - num_foreground = self.num_foreground + fg_files = self.foreground_files + num_fg = self.num_foreground # process signal, foreground, noise, and S+N cls_sig = OrderedDict() @@ -2986,7 +2986,7 @@ def get_masked_sims( nsim_min = min([num_signal, num_noise]) else: nsim_min = num_signal - nsim_max = max([num_signal, num_foreground, num_noise]) + nsim_max = max([num_signal, num_fg, num_noise]) cls_all = np.zeros( [nsim_max, len(map_pairs.items()), len(self.specs), self.lmax + 1] ) @@ -3085,13 +3085,9 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): if null_run: cls_null1t = np.copy(cls_null1_sig) - if do_fg and isim < num_foreground: - fimap_alms, _ = process_index( - foreground_files, None, idx, isim, fg_cache - ) - fjmap_alms, _ = process_index( - foreground_files, None, jdx, isim, fg_cache - ) + if do_fg and isim < num_fg: + fimap_alms, _ = process_index(fg_files, None, idx, isim, fg_cache) + fjmap_alms, _ = process_index(fg_files, None, jdx, isim, fg_cache) cls1_fg = self.alm2cl(fimap_alms, fjmap_alms) if do_noise and isim < num_noise: @@ -3149,7 +3145,7 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): quants += [[cls_sig, cls1_sig]] if null_run: quants += [[cls_null_sig, cls_null1_sig]] - if do_fg and isim < num_foreground: + if do_fg and isim < num_fg: quants += [[cls_fg, cls1_fg]] if do_noise and isim < num_noise: @@ -3195,7 +3191,8 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): cls_null_med[spec][xname] = cls_null_med_arr[xind][s] self.cls_signal = cls_sig - self.cls_fg = cls_fg + if do_fg: + self.cls_fg = cls_fg self.cls_signal_null = cls_null_sig self.cls_noise = cls_noise self.cls_noise_null = cls_null_noise @@ -4109,7 +4106,7 @@ def kernel_precalc(self, map_tag=None, transfer=False): if transfer not in self.signal_components: raise ValueError( "Invalid transfer component {}, must be one of {}".format( - component, ", ".join(self.signal_components) + transfer, ", ".join(self.signal_components) ) ) comps = [transfer] @@ -4225,7 +4222,7 @@ def bin_cl_template( if transfer not in self.signal_components: raise ValueError( "Invalid transfer component {}, must be one of {}".format( - component, ", ".join(self.signal_components) + transfer, ", ".join(self.signal_components) ) ) if "eb" in specs: @@ -4622,8 +4619,8 @@ def get_data_spectra(self, map_tag=None, transfer=False, do_noise=True): nell : OrderedDict Dictionary of noise cross spectra, or None if ``transfer`` is set. debias : OrderedDict - Dictionary of debiased data cross spectra, or None if - ``transfer`` is set. + Dictionary of debiased data cross spectra, or None if ``transfer`` + is set. """ # select map pairs if map_tag is not None: From 4e0e97bc84afc2dca3de55d6f21fd57cd86bf279 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Thu, 21 Jul 2022 00:40:51 -0500 Subject: [PATCH 22/29] bug fixes --- xfaster/parse_tools.py | 7 +++++-- xfaster/xfaster_class.py | 10 +++++++--- xfaster/xfaster_exec.py | 4 ++-- 3 files changed, 14 insertions(+), 7 deletions(-) diff --git a/xfaster/parse_tools.py b/xfaster/parse_tools.py index 2cb1421b..87a08577 100644 --- a/xfaster/parse_tools.py +++ b/xfaster/parse_tools.py @@ -560,18 +560,21 @@ def replace_root(k, v): v = data[k] if isinstance(v, str): data[k] = replace_root(k, v) + elif isinstance(v, list) and isinstance(v[0], str): + data[k] = [replace_root(k, vv) for vv in v] elif isinstance(v, np.ndarray) and isinstance(v.ravel()[0], str): varr = [replace_root(k, vv) for vv in v.ravel()] data[k] = np.array(varr).reshape(v.shape) elif isinstance(v, dict): v1 = list(v.values())[0] - if not isinstance(v1, np.ndarray): + if not (isinstance(v1, np.ndarray) and isinstance(v1.ravel()[0], str)): continue - if not isinstance(v1.ravel()[0], str): + if not (isinstance(v1, list) and isinstance(v1[0], str)): continue if not inplace: v = v.copy() for kk, vv in v.items(): + vv = np.asarray(vv) varr = [replace_root(k, vvv) for vvv in vv.ravel()] v[kk] = np.array(varr).reshape(vv.shape) if not inplace: diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index e8fc44a3..9e388a23 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -1638,14 +1638,14 @@ def save_data(self, name, from_attrs=[], file_attrs=[], **data): if len(file_attrs): data = data.copy() - file_data = {k: data[k] for k in file_attrs} + file_data = {k: copy.deepcopy(data[k]) for k in file_attrs} # strip data roots for storing to disk pt.fix_data_roots( file_data, mode="save", root=self.data_root, root2=getattr(self, "data_root2", None), - inplace=False, + inplace=True, ) data.update(file_data) @@ -2049,7 +2049,7 @@ def process_index(idx): self.log("Mask moments 4: {}".format(self.w4), "debug") self.log("G matrix: {}".format(self.gmat), "debug") - return self.save_data(save_name, from_attrs=save_attrs) + return self.save_data(save_name, from_attrs=save_attrs, file_attrs=file_attrs) def get_noise_residuals(self, filename): """ @@ -2946,6 +2946,10 @@ def get_masked_sims( ret.update(self._get_sim_files("signal", signal_type, signal_subset)) if do_noise: ret.update(self._get_sim_files("noise", noise_type, noise_subset)) + if do_fg: + ret.update( + self._get_sim_files("foreground", foreground_type, foreground_subset) + ) for k, v in ret.items(): setattr(self, k, v) diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index 1faf0266..cb118a8c 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -715,8 +715,8 @@ def xfaster_run( shape_opts.pop("freq_ref_transfer") shape_opts.pop("beta_ref_transfer") shape_transfer_opts = shape_opts.copy() - shape_transfer_opts["filename"] = filename_transfer - shape_transfer_opts["filename_fg"] = filename_transfer_fg + shape_transfer_opts["filename"] = signal_transfer_spec + shape_transfer_opts["filename_fg"] = foreground_transfer_spec shape_transfer_opts["freq_ref"] = freq_ref_transfer shape_transfer_opts["beta_ref"] = beta_ref_transfer From 74fc9841e379109b6aa08948f64de7443285fb06 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Thu, 21 Jul 2022 00:52:05 -0500 Subject: [PATCH 23/29] bug fix --- xfaster/xfaster_class.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 9e388a23..aac25204 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -3674,7 +3674,7 @@ def get_signal_shape( tbeb_flat = np.ones_like(cls_shape["fg_ee"]) * tbeb_flat cls_shape["fg_eb"] = np.copy(tbeb_flat) cls_shape["fg_tb"] = np.copy(tbeb_flat) - else: + elif not tbeb and not transfer: cls_shape["fg_eb"] = np.zeros_like(ell, dtype=float) cls_shape["fg_tb"] = np.zeros_like(ell, dtype=float) From 6f163622fed19975e65693795e24be222f345f40 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Thu, 6 Oct 2022 15:17:31 -0500 Subject: [PATCH 24/29] fix nsim_min handling --- xfaster/xfaster_class.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 0163afc8..19f5d19a 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -2965,6 +2965,8 @@ def get_masked_sims( if do_fg: fg_files = self.foreground_files num_fg = self.num_foreground + else: + num_fg = 0 # process signal, foreground, noise, and S+N cls_sig = OrderedDict() @@ -2986,10 +2988,7 @@ def get_masked_sims( if null_run: cls_null_res[k] = OrderedDict() - if num_noise != 0: - nsim_min = min([num_signal, num_noise]) - else: - nsim_min = num_signal + nsim_min = min([x for x in [num_signal, num_noise, num_fg] if x != 0]) nsim_max = max([num_signal, num_fg, num_noise]) cls_all = np.zeros( [nsim_max, len(map_pairs.items()), len(self.specs), self.lmax + 1] From 5c605ed5126aa7b8705612dea017697d717c6102 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 26 Jul 2023 18:26:19 -0500 Subject: [PATCH 25/29] handle foreground_subset correctly --- xfaster/xfaster_exec.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index 385a2a5e..b61dc345 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -670,6 +670,7 @@ def xfaster_run( noise_subset=noise_subset, foreground_type=foreground_type, foreground_transfer_type=foreground_transfer_type, + foreground_subset=foreground_subset, qb_file=qb_file_sim, ) config_vars.update(sim_opts, "Simulation Ensemble Options") @@ -1551,6 +1552,7 @@ def add_job(self, **kwargs): elif a in [ "noise_subset", "signal_subset", + "foreground_subset", "data_subset", "data_subset2", ]: From cb5b95cb60e08499142575bc5e262bac64ed4761 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Mon, 31 Jul 2023 11:44:30 -0700 Subject: [PATCH 26/29] correct filename variable --- xfaster/xfaster_class.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 083dce30..42f5adc3 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -3666,8 +3666,8 @@ def get_signal_shape( "debug", ) else: - if not os.path.exists(filename): - raise OSError("Missing foreground model file {}".format(filename)) + if not os.path.exists(filename_fg): + raise OSError("Missing foreground model file {}".format(filename_fg)) cls_fg = st.load_camb_cl(filename_fg, lmax=lmax_kern, pol=None) cls_shape.update({"fg_" + s: cls for s, cls in zip(specs, cls_fg)}) From adc4379e9aeede9e84143e7edd45252bcde6dcf3 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Tue, 19 Sep 2023 13:56:44 +0200 Subject: [PATCH 27/29] black --- xfaster/xfaster_class.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 67c91669..1a3cc29d 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -3663,7 +3663,9 @@ def get_signal_shape( ) else: if not os.path.exists(filename_fg): - raise OSError("Missing foreground model file {}".format(filename_fg)) + raise OSError( + "Missing foreground model file {}".format(filename_fg) + ) cls_fg = st.load_camb_cl(filename_fg, lmax=lmax_kern, pol=None) cls_shape.update({"fg_" + s: cls for s, cls in zip(specs, cls_fg)}) From e5ef0e5c39774517ac36fc1f74e3c977a9fe22c5 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Tue, 19 Sep 2023 14:23:13 +0200 Subject: [PATCH 28/29] black --- xfaster/xfaster_class.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 1a3cc29d..62f09377 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -6138,11 +6138,8 @@ def expand_transfer(qb_transfer, bin_def, check_only=False): "insufficient number of signal sims" ) except Exception as e: - msg = ( - "Unable to adjust negative bins for map {}: {}".format( - m0, str(e) - ) - ) + msg = "Unable to adjust negative bins for map {}: {}" + msg = msg.format(m0, str(e)) success = False # Set TE/EB/TB qb transfers to geometric means of components, if necessary From 54ed485d483f97190917215ce523249be6f4e0fa Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 4 Oct 2023 12:13:22 -0500 Subject: [PATCH 29/29] remove vestigial MPI and PBS support from batch_tools --- xfaster/batch_tools.py | 168 +++++++++++++--------------------------- xfaster/xfaster_exec.py | 15 +--- 2 files changed, 54 insertions(+), 129 deletions(-) diff --git a/xfaster/batch_tools.py b/xfaster/batch_tools.py index 6255b558..8defd96a 100644 --- a/xfaster/batch_tools.py +++ b/xfaster/batch_tools.py @@ -21,15 +21,7 @@ def get_job_logfile(): logfile : str Path to log file. """ - if os.getenv("PBS_O_WORKDIR"): - if os.getenv("PBS_ENVIRONMENT") != "PBS_INTERACTIVE": - workdir = os.getenv("PBS_O_WORKDIR") - jobname = os.getenv("PBS_JOBNAME") - jobid = os.getenv("PBS_JOBID").split(".", 1)[0] - logfile = os.path.join(workdir, "{}.u{}".format(jobname, jobid)) - else: - logfile = None - elif os.getenv("SLURM_SUBMIT_DIR"): + if os.getenv("SLURM_SUBMIT_DIR"): workdir = os.getenv("SLURM_SUBMIT_DIR") jobname = os.getenv("SLURM_JOB_NAME") jobid = os.getenv("SLURM_JOB_ID").split(".", 1)[0] @@ -37,7 +29,6 @@ def get_job_logfile(): logfile = None else: logfile = os.path.join(workdir, "{}-{}.log".format(jobname, jobid)) - # TODO generate different logs for multiple processes in same job? else: logfile = None return logfile @@ -89,11 +80,9 @@ def batch_sub( workdir=None, batch_args=[], omp_threads=None, - mpi_procs=None, - mpi_args="", env_script=None, env=None, - nice=0, + nice=None, echo=True, delete=True, submit=True, @@ -101,14 +90,16 @@ def batch_sub( debug=False, exclude=None, verbose=False, + srun=None, + srun_args="", ): """ - Create and submit a SLURM or PBS job. + Create and submit a SLURM job. Arguments --------- cmd : string or list of strings - A command sequence to run via SLURM or PBS. + A command sequence to run via SLURM. The command will be inserted into a qsub submission script with all of the options specified in the remaining arguments. name : string, optional @@ -118,12 +109,11 @@ def batch_sub( Or pass a string (eg '4gb') to use directly. nodes : int or string, optional Number of nodes to use in job - If a string, will be passed as-is to PBS -l node= resource If using SLURM and a string, will overwrite node_list if None node_list : string or list of strings List of nodes that can be used for job. SLURM-only. ppn : int, optional - Numper of processes per node + Number of processes per node cput : string or float or datetime.timedelta, optional Amount of CPU time requested. String values should be in the format HH:MM:SS, e.g. '10:00:00'. @@ -133,9 +123,9 @@ def batch_sub( String values should be in the format HH:MM:SS, e.g. '10:00:00'. Numerical values are interpreted as a number of hours. output : string, optional - PBS standard output filename. + Scheduler standard output filename. error : string, optional - PBS error output filename. + Scheduler error output filename. queue : string, optional The name of the queue to which to submit jobs dep_afterok : string or list of strings @@ -146,26 +136,17 @@ def batch_sub( This is where the output and error files will be created by default. Default: current directory. batch_args : string or list of strings, optional - Any additional arguments to pass to slurm/pbs. + Any additional arguments to pass to slurm. omp_threads : int, optional Number of OpenMP threads to use per process - mpi_procs : int - Number of MPI processes to use. - ``mpirun`` calls will be added to all lines of cmd as needed. - If cmd contains ``mpirun`` or ``mpiexec``, this does nothing. - mpi_args : string - Additional command line arguments for inserted ``mpirun`` commands. - If cmd contains ``mpirun`` or ``mpiexec``, this does nothing. env_script : string, optional Path to script to source during job script preamble For loading modules, setting environment variables, etc env : dict, optional Dictionary of environment variables to set in job script nice : int, optional - Adjust scheduling priority (SLURM only). Range from -5000 (highest - priority) to 5000 (lowest priority). - Note: actual submitted --nice value is 5000 higher, since negative - values require special privilege. + Adjust scheduling priority (SLURM only). + Unprivileged users can only increase niceness (lower job priority). echo : bool, optional Whether to use bash "set -x" in job script to echo commands to stdout. delete : bool, optional @@ -174,13 +155,19 @@ def batch_sub( If True (default) submit the job script once create. Will override the default option when False, to keep the script scheduler : string, optional - Which scheduler system to write a script for. One of "pbs" or "slurm" + Which scheduler system to write a script for. Only "slurm" is currently + supported. debug : bool, optional If True, print the contents of the job script to stdout for debugging. exclude : string or list of strings List of nodes that will be excluded for job. SLURM-only. verbose : bool, optional Print the working directory, and the job ID if submitted successfully. + srun : bool, optional + Invoke jobs with srun. Defaults to True when using SLURM, set to False + to disable. Use `srun_args` to pass extra arguments to this command. + srun_args : string, optional + If `srun=True`, additional arguments to pass to the `srun` command. Returns ------- @@ -203,8 +190,6 @@ def batch_sub( if mem is not None and not isinstance(mem, str): if mem < 0: mem = None - elif scheduler == "pbs": - mem = "{:d}mb".format(int(np.ceil(mem * 1024.0))) elif scheduler == "slurm": mem = "{:d}".format(int(np.ceil(mem * 1024.0))) if isinstance(dep_afterok, str): @@ -220,32 +205,13 @@ def batch_sub( if scheduler == "slurm" and node_list is None: node_list = nodes nodes = 1 + if srun is None: + srun = scheduler == "slurm" job_script = ["#!/usr/bin/env bash"] # TODO can maybe replace manual option with some automatic detection - if scheduler == "pbs": - # create PBS header - if name: - job_script += ["#PBS -N {:s}".format(name)] - if mem: - job_script += ["#PBS -l mem={:s}".format(mem)] - if nodes and ppn: - job_script += ["#PBS -l nodes={}:ppn={}".format(nodes, ppn)] - if cput: - job_script += ["#PBS -l cput={:s}".format(format_time(cput))] - if wallt: - job_script += ["#PBS -l walltime={:s}".format(format_time(wallt))] - if output: - job_script += ["#PBS -o {:s}".format(output)] - if error: - job_script += ["#PBS -e {:s}".format(error)] - if queue: - job_script += ["#PBS -q {:s}".format(queue)] - if dep_afterok: - job_script += ["#PBS -W depend=afterok:{}".format(":".join(dep_afterok))] - - elif scheduler == "slurm": + if scheduler == "slurm": # create slurm header if name: job_script += ["#SBATCH --job-name={:s}".format(name)] @@ -276,7 +242,6 @@ def batch_sub( if wallt: job_script += ["#SBATCH --time={:s}".format(format_time(wallt))] if nice is not None: - nice += 5000 job_script += ["#SBATCH --nice={}".format(nice)] if output: job_script += ["#SBATCH --output={:s}".format(output)] @@ -299,20 +264,27 @@ def batch_sub( if env: for k, v in env.items(): job_script += ["export {}={}".format(k, v)] - if scheduler == "pbs": - job_script += ["cd $PBS_O_WORKDIR"] - elif scheduler == "slurm": + if scheduler == "slurm": job_script += ["cd $SLURM_SUBMIT_DIR"] if omp_threads: job_script += ["export OMP_NUM_THREADS={}".format(omp_threads)] + if srun: + job_script += ["export SRUN_CPUS_PER_TASK={}".format(omp_threads)] # finally, add the command string to script - if mpi_procs is not None: - if "mpirun" not in cmd and "mpiexec" not in cmd: - mpi = "mpiexec -n {:d} {:s} ".format(mpi_procs, mpi_args) - cmd = "\n".join( - [(mpi + line) if line != "wait" else line for line in cmd.split("\n")] - ) + cmd_prefix = None + if srun: + if scheduler != "slurm": + raise ValueError('srun=True requires scheduler="slurm"') + if not cmd.startswith("srun"): + cmd_prefix = "srun {:s} ".format(srun_args) + if cmd_prefix is not None: + cmd = "\n".join( + [ + (cmd_prefix + line) if line != "wait" else line + for line in cmd.split("\n") + ] + ) job_script += [cmd] job_script = "\n".join(job_script) @@ -334,9 +306,7 @@ def batch_sub( # create and submit script prefix = "{}_".format(name if name else "job") - if scheduler == "pbs": - suffix = ".qsub" - else: + if scheduler == "slurm": suffix = ".slurm" with tempfile.NamedTemporaryFile( prefix=prefix, suffix=suffix, mode="w", dir=workdir, delete=delete @@ -345,12 +315,7 @@ def batch_sub( f.flush() os.chmod(f.name, stat.S_IRUSR | stat.S_IWUSR | stat.S_IXUSR | stat.S_IRGRP) if submit: - if scheduler == "pbs": - ret = sp.check_output(["qsub"] + batch_args + [f.name]).decode("UTF-8") - jobid = ret.split("\n")[0] # parse jobid - if not re.match("[0-9]+\.[\w]+", jobid): - raise RuntimeError("qsub error:\n{}".format(ret)) - elif scheduler == "slurm": + if scheduler == "slurm": ret = sp.check_output(["sbatch"] + batch_args + [f.name]).decode( "UTF-8" ) @@ -382,7 +347,7 @@ def batch_sub( def batch_group(cmds, group_by=1, serial=False, *args, **kwargs): """ - Create and submit SLURM or PBS job scripts for a group of similar commands. The + Create and submit SLURM job scripts for a group of similar commands. The commands can be grouped together into larger single jobs that run them in parallel on multiple processors on a node. @@ -403,7 +368,7 @@ def batch_group(cmds, group_by=1, serial=False, *args, **kwargs): Eg. on scinet use ``group_by=8`` to efficiently use whole nodes. serial : bool, optional Set to ``True`` to run cmds sequentially, rather than starting them all - in parallel. This will also work with MPI/OpenMP parallel jobs. + in parallel. This will also work with OpenMP parallel jobs. args, kwargs : arb Additional arguments passed to batch_sub @@ -452,10 +417,8 @@ def __init__( Standardized way to add job submission arguments to a script. Keyword arguments are used to fix values for parameters not needed by a - script. The corresponding command line arguments will not be added. For - example, if not using MPI, pass ``mpi_procs=None`` and then there will - be no ``--mpi-procs`` argument on the command line. See the - ``.opt_list`` attribute for a complete list. + script. The corresponding command line arguments will not be added. See + the ``.opt_list`` attribute for a complete list. Arguments --------- @@ -476,7 +439,7 @@ def __init__( Example ------- >>> jp = sa.batch.JobArgumentParser("some_serial_job", mem=4, time=1.5, - mpi_procs=None, omp_threads=1, workdir="/path/to/logs") + omp_threads=1, workdir="/path/to/logs") >>> AP = argparse.ArgumentParser(description="Do some job") >>> AP.add_argument(...) # other non-job arguments >>> jp.add_arguments(AP) @@ -519,7 +482,7 @@ def _opt_list(self): dict( type=int, default=None, - help="Processes per node. Default based on group, omp_threads, and mpi_procs", + help="Processes per node. Default based on group and omp_threads", ), ), ( @@ -531,14 +494,6 @@ def _opt_list(self): help="Nodes to exclude from jobs", ), ), - ( - "pbs", - dict( - action="store_true", - default=False, - help="Create PBS (rather than SLURM) job scripts", - ), - ), ( "use_cput", dict( @@ -559,8 +514,8 @@ def _opt_list(self): "nice", dict( type=int, - default=0, - help="Priority from -5000 (hi) to 5000 (lo). SLURM only", + default=None, + help="Adjust job priority. SLURM only", ), ), ( @@ -584,12 +539,6 @@ def _opt_list(self): type=int, default=1, help="Number of OpenMP threads per process" ), ), - ( - "mpi_procs", - dict( - type=int, default=1, help="Number of MPI processes (per node)" - ), - ), ( "group", dict( @@ -603,15 +552,7 @@ def _opt_list(self): dict( action="store_true", default=False, - help="Run grouped commands serially. Works for MPI", - ), - ), - ( - "procs_scale", - dict( - action="store_true", - default=False, - help="Scale time and memory by number of processes", + help="Run grouped commands serially.", ), ), ] @@ -698,8 +639,7 @@ def set_job_opts(self, args, load_defaults=True, **kwargs): if args["ppn"] is None: args["ppn"] = args["group"] - scale = 1.0 if not args["procs_scale"] else float(args["mpi_procs"]) - mem_scale = scale * (args["group"] if not args["serial"] else 1.0) + mem_scale = args["group"] if not args["serial"] else 1.0 mem = self.mem * mem_scale if self.mem is not None else None self.job_opts = dict( @@ -710,7 +650,6 @@ def set_job_opts(self, args, load_defaults=True, **kwargs): ppn=args["ppn"], queue=args["queue"], omp_threads=args["omp_threads"], - mpi_procs=args["mpi_procs"], env_script=args["env_script"], nice=args["nice"], delete=args["test"], @@ -719,12 +658,9 @@ def set_job_opts(self, args, load_defaults=True, **kwargs): group_by=args["group"], serial=args["serial"], ) - if args["pbs"]: - self.job_opts["scheduler"] = "pbs" - else: - self.job_opts["scheduler"] = "slurm" + self.job_opts["scheduler"] = "slurm" - time = self.time / args["cpu_speed"] / scale + time = self.time / args["cpu_speed"] if time < 0.25: time = 0.25 if args["use_cput"]: diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index e28976cc..bd1f20cc 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -1643,9 +1643,7 @@ def set_job_options( """ - # XFaster runs faster with OMP and does not use MPI - mpi = False - mpi_procs = None + # XFaster runs faster with OMP if omp_threads is None: omp_threads = ppn @@ -1670,15 +1668,7 @@ def set_job_options( if mem is not None: mem *= ppn - if not mpi: - mpi_procs = None - elif mpi_procs is None and ppn is not None: - mpi_procs = ppn - - if pbs: - scheduler = "pbs" - else: - scheduler = "slurm" + scheduler = "slurm" self.batch_args = dict( workdir=workdir, @@ -1690,7 +1680,6 @@ def set_job_options( queue=queue, env_script=env_script, omp_threads=omp_threads, - mpi_procs=mpi_procs, nice=nice, delete=False, submit=not test,