diff --git a/docs/notebooks/XFaster_Tutorial.ipynb b/docs/notebooks/XFaster_Tutorial.ipynb index 614972e7..58ea0fb7 100644 --- a/docs/notebooks/XFaster_Tutorial.ipynb +++ b/docs/notebooks/XFaster_Tutorial.ipynb @@ -404,7 +404,7 @@ "\n", "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, and set `lmin_fg` and `lmax_fg` to optionally limit the bandwidth over which foregrounds are to be fit. 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, 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." + "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, and set `lmin_fg` and `lmax_fg` to optionally limit the bandwidth over which foregrounds are to be fit. 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." ] }, { @@ -541,7 +541,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now we want to get the ensemble average of all of our signal and noise simulations, which we're going to use to calculate the filter transfer function and the noise shape, respectively. This is done with [get_masked_sims()](../api.rst#xfaster.xfaster_class.XFaster.get_masked_sims). The method will also compute the signal cross noise terms, which are used for null tests, where they can contribute significantly to the expected residuals that are subtracted from the data." + "Now we want to get the ensemble average of all of our signal (and/or optionally foreground) and noise simulations, which we're going to use to calculate the filter transfer function and the noise shape, respectively. This is done with [get_masked_sims()](../api.rst#xfaster.xfaster_class.XFaster.get_masked_sims). The method will also compute the signal cross noise terms, which are used for null tests, where they can contribute significantly to the expected residuals that are subtracted from the data." ] }, { @@ -574,7 +574,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", @@ -719,9 +720,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 [xfaster_run()](../api.rst#xfaster.xfaster_exec.xfaster_run) argument `signal_transfer_spec` to a file containing the spectrum for the signal component. 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.\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 [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." + "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`." ] }, { @@ -781,13 +782,13 @@ "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=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=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=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", @@ -804,7 +805,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." ] }, { diff --git a/example/xfaster_example_fg.py b/example/xfaster_example_fg.py index 6fbb2b7e..79ed1a36 100644 --- a/example/xfaster_example_fg.py +++ b/example/xfaster_example_fg.py @@ -34,6 +34,7 @@ "noise_type": "gaussian", "mask_type": "rectangle", "signal_type": "synfast", + "foreground_type": "gaussian", "data_root2": None, "data_subset2": None, # residual fitting @@ -50,6 +51,7 @@ # spectrum "ensemble_mean": False, "tbeb": True, + "fix_fg_transfer": True, "converge_criteria": 0.005, "iter_max": 200, "save_iters": False, diff --git a/xfaster/parse_tools.py b/xfaster/parse_tools.py index 2a218461..b58c0c39 100644 --- a/xfaster/parse_tools.py +++ b/xfaster/parse_tools.py @@ -315,6 +315,9 @@ def load_and_parse(filename, check_version=True): ]: data.pop(k, None) + # 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: diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 6ce451b9..a9196df1 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -2825,6 +2825,8 @@ def get_masked_sims( signal_subset="*", noise_type=None, noise_subset="*", + foreground_type=None, + foreground_subset="*", qb_file=None, ): """ @@ -2864,6 +2866,15 @@ def get_masked_sims( 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_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_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]*'``. qb_file : string Pointer to a bandpowers.npz file in the output directory, used to correct the ensemble mean noise spectrum by the appropriate @@ -2876,6 +2887,8 @@ def get_masked_sims( cls_signal, cls_signal_null: Mean signal spectra + cls_fg: + Mean foreground spectra cls_noise, cls_noise_null: Mean noise spectra cls_sim, cls_sim_null: @@ -2930,6 +2943,20 @@ def get_masked_sims( file_attrs += ["noise_root2", "noise_files2"] save_attrs += ["cls_noise_null", "cls_res_null"] + do_fg = not null_run and foreground_type is not None + if do_fg: + opts.update( + foreground_type=foreground_type, foreground_subset=foreground_subset + ) + file_attrs += [ + "foreground_type", + "foreground_subset", + "foreground_root", + "foreground_files", + "num_foreground", + ] + save_attrs += ["cls_fg"] + save_attrs += file_attrs alt_name = None @@ -2937,12 +2964,18 @@ def get_masked_sims( if do_noise: alt_name = save_name save_name = "{}_{}".format(save_name, noise_type) + if do_fg: + save_name = "{}_{}".format(save_name, foreground_type) cp = "sims_transfer" if transfer else "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: - self.force_rerun[cp] = False + if ( + not hasattr(self, "foreground_type") + or self.foreground_type == foreground_type + ): + self.force_rerun[cp] = False ret = self.load_data( save_name, @@ -2966,6 +2999,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) @@ -2980,8 +3017,15 @@ def get_masked_sims( else: num_noise = 0 - # process signal, noise, and S+N + 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() + cls_fg = OrderedDict() if do_fg else None cls_noise = OrderedDict() if do_noise else None cls_tot = OrderedDict() cls_med = OrderedDict() @@ -3002,11 +3046,8 @@ 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_max = max([num_signal, num_noise]) + 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] ) @@ -3061,7 +3102,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: @@ -3074,10 +3115,12 @@ def process_index(files, files2, idx, idx2=None, cache=None, rls=None): return cache[idx] sig_cache = dict() + fg_cache = dict() noise_cache = dict() for isim in range(nsim_max): sig_cache.clear() + fg_cache.clear() noise_cache.clear() for xind, (xname, (idx, jdx)) in enumerate(map_pairs.items()): self.log( @@ -3103,6 +3146,11 @@ 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_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: cls1_res = OrderedDict() if null_run: @@ -3158,6 +3206,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 do_fg and isim < num_fg: + quants += [[cls_fg, cls1_fg]] if do_noise and isim < num_noise: quants += [[cls_noise, cls1_noise]] @@ -3202,6 +3252,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 + if do_fg: + self.cls_fg = cls_fg self.cls_noise = cls_noise self.cls_sim = cls_tot self.cls_med = cls_med @@ -3513,9 +3565,12 @@ def get_signal_shape( 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 the filename is a - relative path and not found, the config directory is searched. If - not supplied, a power law dust model spectrum is assumed. + 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. 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. @@ -3531,6 +3586,9 @@ def get_signal_shape( and ``r`` is None and ``flat`` is False, will search for a signal 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. @@ -3546,7 +3604,7 @@ def get_signal_shape( specs = list(self.specs) comps = list(self.signal_components) - if (self.null_run or transfer) and "fg" in comps: + if self.null_run and "fg" in comps: comps.remove("fg") if component in self.signal_components: comps = [component] @@ -3650,6 +3708,14 @@ def get_signal_shape( 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_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) @@ -3660,8 +3726,6 @@ def get_signal_shape( "debug", ) else: - if not os.path.exists(filename_fg): - filename_fg = os.path.join(self.config_root, filename_fg) if not os.path.exists(filename_fg): raise OSError( "Missing foreground model file {}".format(filename_fg) @@ -3676,7 +3740,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) @@ -4075,9 +4139,11 @@ def kernel_precalc(self, map_tag=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. - transfer : bool - If True, set transfer function to 1 to solve for transfer function - qbs. + 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 ------- @@ -4103,8 +4169,17 @@ def kernel_precalc(self, map_tag=None, transfer=False): lk = slice(0, lmax + 1) mll = OrderedDict() + comps = [] if transfer: - comps = ["cmb"] + if transfer is True: + transfer = "cmb" + if transfer not in self.signal_components: + raise ValueError( + "Invalid transfer component {}, must be one of {}".format( + transfer, ", ".join(self.signal_components) + ) + ) + comps = [transfer] else: comps = list(self.signal_components) @@ -4171,10 +4246,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. - transfer : 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 = 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. @@ -4210,10 +4287,20 @@ def bin_cl_template( specs = list(self.specs) if transfer: + if transfer is True: + transfer = "cmb" + if transfer not in self.signal_components: + raise ValueError( + "Invalid transfer component {}, must be one of {}".format( + transfer, ", ".join(self.signal_components) + ) + ) if "eb" in specs: specs.remove("eb") if "tb" in specs: specs.remove("tb") + if transfer == "fg" and "te" in specs: + specs.remove("te") lmax = self.lmax lmax_kern = lmax # 2 * self.lmax @@ -4234,15 +4321,16 @@ def bin_cl_template( cbl = OrderedDict() comps = [] - for comp in self.signal_components: - if any([k.startswith(comp + "_") for k in cls_shape]): - comps += [comp] - if "fg" in comps and transfer: - comps.remove("fg") - 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 + 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: + 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) @@ -4585,10 +4673,10 @@ def get_data_spectra(self, map_tag=None, transfer=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 : 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 : 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 ``get_masked_data`` for how these are computed. The input noise is similarly either ``cls_noise_null`` or ``cls_noise``. do_noise : bool @@ -4599,10 +4687,10 @@ def get_data_spectra(self, map_tag=None, transfer=False, do_noise=True): obs : OrderedDict Dictionary of data cross spectra nell : OrderedDict - Dictionary of noise cross spectra, or None if ``transfer`` is True. + Dictionary of noise cross spectra, or None if ``transfer`` is set. debias : OrderedDict Dictionary of debiased data cross spectra, or None if ``transfer`` - is True. + is set. """ # select map pairs if map_tag is not None: @@ -4613,14 +4701,25 @@ def get_data_spectra(self, map_tag=None, transfer=False, do_noise=True): # select spectra tbeb = "tb" in self.specs + specs = list(self.specs) if transfer or not tbeb: - specs = self.specs[:4] - else: - specs = self.specs + if "eb" in specs: + specs.remove("eb") + if "tb" in specs: + specs.remove("tb") + if transfer == "fg" and "te" in specs: + specs.remove("te") # obs depends on what you're computing if transfer: - obs_quant = self.cls_signal + if transfer is True or transfer == "cmb": + obs_quant = self.cls_signal + elif transfer == "fg": + obs_quant = self.cls_fg + else: + raise ValueError( + "Unknown transfer function component {}".format(transfer) + ) elif self.null_run: if self.reference_type is not None: obs_quant = self.cls_data_sub_null @@ -5402,9 +5501,9 @@ def fisher_iterate( qb_start : OrderedDict Initial guess at ``qb`` bandpower amplitudes. If None, unity is assumed for all bins. - transfer : bool - If True, the input Cls passed to ``fisher_calc`` are the average - of the signal simulations, and noise cls are ignored. + 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`` (for null tests) or ``cls_data``. See ``get_masked_data`` for how these are computed. The input noise is similarly either @@ -5460,7 +5559,7 @@ def fisher_iterate( cls_shape : shape spectrum iters : number of iterations completed - If ``transfer`` is False, this dictionary also contains:: + If ``transfer`` is not set, this dictionary also contains:: qb_transfer : transfer function amplitudes transfer : ell-by-ell transfer function @@ -5470,7 +5569,7 @@ def fisher_iterate( Notes ----- This method stores outputs to files with name 'transfer' for transfer - function runs (if ``transfer = True``), 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). @@ -5478,7 +5577,11 @@ def fisher_iterate( if transfer: null_first_cmb = False - save_name = "transfer_wbins" if self.weighted_bins else "transfer" + save_name = "transfer_{}_wbins" if self.weighted_bins else "transfer_{}" + if transfer is True: + transfer = "cmb" + assert transfer in self.signal_components + save_name = save_name.format(transfer) else: save_name = "bandpowers" @@ -5489,9 +5592,16 @@ def fisher_iterate( if qb_start is None: qb = OrderedDict() for k, v in self.bin_def.items(): + # transfer function for signal component if transfer: - if "cmb" not in k or "eb" in k or "tb" in k: + if "eb" in k or "tb" in k: + continue + if transfer == "fg" and "te" in k: continue + if not k.startswith(transfer): + 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 @@ -5588,7 +5698,7 @@ def fisher_iterate( extra_tag=file_tag, ) - if "fg" in self.components and not transfer: + if "fg" in self.components: out.update( map_freqs=self.map_freqs, freq_ref=self.freq_ref, @@ -5652,7 +5762,7 @@ 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 else "spectrum", + "{} transfer function".format(transfer) if transfer else "spectrum", iter_max, ) # Check the slope of the last ten fqb_maxpoints. @@ -5708,7 +5818,7 @@ def fisher_iterate( weighted_bins=self.weighted_bins, ) - if "fg" in self.components and not transfer: + if "fg" in self.components: out.update( map_freqs=self.map_freqs, freq_ref=self.freq_ref, @@ -5891,6 +6001,7 @@ def get_transfer( iter_max=200, save_iters=False, fix_bb_transfer=False, + fix_fg_transfer=False, ): """ Compute the transfer function from signal simulations created using the @@ -5913,6 +6024,10 @@ 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. + 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 ------- @@ -5945,6 +6060,11 @@ def get_transfer( weighted_bins=self.weighted_bins, ) + if "fg" in comps: + if getattr(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) @@ -5960,7 +6080,7 @@ def get_transfer( def expand_transfer(qb_transfer, bin_def, check_only=False): if not check_only: transfer = OrderedDict() - if not check_only and "fg" in comps: + 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_") qb_transfer[ftag] = OrderedDict() @@ -5998,7 +6118,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 ["tb", "eb"] and stag 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]) @@ -6012,8 +6132,6 @@ def expand_transfer(qb_transfer, bin_def, check_only=False): if ret is not None: try: check_only = ret.get("transfer", None) is not None - if "fg" in comps: - check_only &= any(k.startswith("fg_") for k in ret["qb_transfer"]) xfer = expand_transfer(ret["qb_transfer"], ret["bin_def"], check_only) except: ret = None @@ -6023,74 +6141,88 @@ def expand_transfer(qb_transfer, bin_def, check_only=False): return self.transfer self.qb_transfer = OrderedDict() - for spec in self.specs: - self.qb_transfer["cmb_" + spec] = 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) - ) - self.log("Setting map {} transfer to unity".format(m0), "info") - continue + for comp in comps: + for spec in self.specs: + self.qb_transfer["{}_{}".format(comp, spec)] = OrderedDict() - 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=True) - ret = self.fisher_iterate( - cbl, - m0, - transfer=True, - 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( - "Transfer function amplitude {} < 0" - "for {} bin {} of map {}".format(v, k, negbin, m0) + 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" ) - # 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 + continue + + self.log( + "Computing {} transfer function for map {}/{}".format( + comp, im0 + 1, self.num_maps + ), + "info", + ) + self.clear_precalc() + cbl = self.bin_cl_template(map_tag=m0, transfer=comp) + ret = self.fisher_iterate( + cbl, + m0, + transfer=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" - ) - except Exception as e: - msg = "Unable to adjust negative bins for map {}: {}".format( - m0, str(e) + "Transfer function amplitude {} < 0" + "for {} bin {} of map {}".format(v, k, negbin, m0) ) - 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"])) + # 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 {}: {}" + msg = msg.format(m0, str(e)) + success = False + + # 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 + for stag, qbdat in qb.items(): + self.qb_transfer[stag][m0] = qbdat if success: self.transfer = expand_transfer(self.qb_transfer, self.bin_def) diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index e639e27a..4e9a7d29 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -67,6 +67,9 @@ def xfaster_run( signal_type="synfast", signal_subset="*", signal_transfer_type=None, + foreground_type=None, + foreground_subset="*", + foreground_transfer_type=None, noise_type="stationary", noise_subset="*", qb_file_sim=None, @@ -74,9 +77,12 @@ def xfaster_run( signal_spec=None, signal_transfer_spec=None, foreground_spec=None, + foreground_transfer_spec=None, model_r=None, freq_ref=359.7, beta_ref=1.54, + freq_ref_transfer=None, + beta_ref_transfer=None, # data options data_type="raw", template_type=None, @@ -113,6 +119,7 @@ def xfaster_run( return_cls=False, qb_only=False, fix_bb_transfer=False, + fix_fg_transfer=False, null_first_cmb=False, delta_beta_prior=None, like_profiles=False, @@ -217,14 +224,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. lmin_fg : int Minimum ell to use for defining foreground bins. If not set, defaults to ``lmin``. @@ -263,6 +275,16 @@ def xfaster_run( noise_subset : str The subset of the noise sims to include. Must be a glob-parseable string. + 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``. qb_file_sim : str If not None, pointer to a bandpowers.npz file in the output directory, to correct the noise ensemble by an appropriate set of residual ``qb`` @@ -278,8 +300,15 @@ def xfaster_run( directory. foreground_spec : string The spectrum data file to use for estimating foreground component - bandpowers. If not supplied, use a simple power law model for dust + bandpowers. If not supplied, will search for + ``spec_foreground_.dat`` in the foreground sim + directory. Otherwise, use a simple power law model for dust foregrounds. + 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``. @@ -290,6 +319,12 @@ def xfaster_run( 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. data_type : str The data type to use template_type : str @@ -410,6 +445,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). delta_beta_prior : float @@ -479,8 +518,15 @@ def xfaster_run( signal_transfer_type = signal_type all_opts["signal_transfer_type"] = signal_transfer_type + if foreground_transfer_type is None: + foreground_transfer_type = foreground_type + all_opts["foreground_transfer_type"] = foreground_transfer_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 checkpoint == "sims": @@ -640,32 +686,58 @@ def xfaster_run( signal_subset=signal_subset, noise_type=noise_type, 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") config_vars.remove_option("XFaster General", "qb_file_sim") sim_opts.pop("signal_transfer_type") + sim_opts.pop("foreground_transfer_type") sim_transfer_opts = sim_opts.copy() if signal_transfer_type != signal_type: sim_transfer_opts["signal_type"] = signal_transfer_type sim_transfer_opts.pop("noise_type") sim_transfer_opts.pop("noise_subset") sim_transfer_opts.pop("qb_file") + if foreground_transfer_type != foreground_type: + sim_transfer_opts["foreground_type"] = foreground_transfer_type + sim_transfer_opts.pop("noise_type", None) + sim_transfer_opts.pop("noise_subset", None) + sim_transfer_opts.pop("qb_file", None) + + if freq_ref_transfer is None: + freq_ref_transfer = freq_ref + if beta_ref_transfer is None: + beta_ref_transfer = beta_ref shape_opts = dict( filename=signal_spec, filename_transfer=signal_transfer_spec, filename_fg=foreground_spec, + filename_transfer_fg=foreground_transfer_spec, r=model_r, freq_ref=freq_ref, beta_ref=beta_ref, + freq_ref_transfer=freq_ref_transfer, + beta_ref_transfer=beta_ref_transfer, ) config_vars.update(shape_opts, "Shape Spectrum Options") config_vars.remove_option("XFaster General", "signal_spec") config_vars.remove_option("XFaster General", "signal_transfer_spec") config_vars.remove_option("XFaster General", "foreground_spec") + config_vars.remove_option("XFaster General", "foreground_transfer_spec") config_vars.remove_option("XFaster General", "model_r") shape_opts.pop("filename_transfer") + shape_opts.pop("filename_transfer_fg") + shape_opts.pop("freq_ref_transfer") + shape_opts.pop("beta_ref_transfer") + shape_transfer_opts = shape_opts.copy() + 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 spec_opts = dict( multi_map=multi_map, @@ -673,6 +745,7 @@ def xfaster_run( iter_max=iter_max, save_iters=save_iters, fix_bb_transfer=fix_bb_transfer, + fix_fg_transfer=fix_fg_transfer, delta_beta_prior=delta_beta_prior, cond_noise=cond_noise, cond_criteria=cond_criteria, @@ -689,6 +762,7 @@ def xfaster_run( spec_opts.pop("multi_map") 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() @@ -756,7 +830,7 @@ def xfaster_run( X.get_masked_sims(transfer=True, **sim_transfer_opts) X.log("Loading spectrum shape for transfer function...", "notice") - X.get_signal_shape(filename=signal_transfer_spec, transfer=True) + X.get_signal_shape(transfer=True, **shape_transfer_opts) X.log("Computing transfer functions...", "notice") X.get_transfer(**transfer_opts) @@ -1198,6 +1272,9 @@ def add_arg( add_arg(G, "signal_transfer_type") add_arg(G, "noise_type") add_arg(G, "noise_subset") + add_arg(G, "foreground_type") + add_arg(G, "foreground_subset") + add_arg(G, "foreground_transfer_type") add_arg(G, "qb_file_sim") # shape spectrum options @@ -1207,8 +1284,11 @@ def add_arg( add_arg(E, "model_r") add_arg(G, "signal_transfer_spec") 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) # data options G = PP.add_argument_group("data options") @@ -1256,6 +1336,7 @@ def add_arg( add_arg(G, "return_cls") add_arg(G, "qb_only") add_arg(G, "fix_bb_transfer") + add_arg(G, "fix_fg_transfer") add_arg(G, "null_first_cmb") add_arg(G, "delta_beta_prior", argtype=float) add_arg(G, "like_profiles") @@ -1520,6 +1601,7 @@ def add_job(self, **kwargs): elif a in [ "noise_subset", "signal_subset", + "foreground_subset", "data_subset", "data_subset2", ]: