From 959582d54f6fb29881a57867e972078b2b9c888b Mon Sep 17 00:00:00 2001 From: Anne Gambrel Date: Wed, 19 May 2021 15:35:42 -0400 Subject: [PATCH 01/44] Create new branch fg_like to include updates of the branch xfaster/fg_like in labah repo. This commit is done by Vy Luu tluu@princeton.edu --- xfaster/parse_tools.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/xfaster/parse_tools.py b/xfaster/parse_tools.py index c8b98141..b3125552 100644 --- a/xfaster/parse_tools.py +++ b/xfaster/parse_tools.py @@ -206,6 +206,8 @@ def parse_data(data, field): if version == 1: if "foreground_type" in data: data["foreground_type_sim"] = data.pop("foreground_type") + if "clean_type" in data: + data["data_type"] = data.pop("clean_type") return data[field] From 0174c9f6f3a36baf65212b82d3defabc0426f595 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Sun, 30 May 2021 21:07:06 -0400 Subject: [PATCH 02/44] add function bin_cl_template_tmat to xfaster_class --- xfaster/xfaster_class.py | 362 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 362 insertions(+) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index cc7ebe56..891e286f 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -4584,6 +4584,364 @@ def bin_things(comp, d, md): bin_things(comp, d, md) return cbl + + def bin_cl_template_tmat( + self, + cls_shape, + map_tag=None, + beam_error=False, + fg_ell_ind=0, + ): + """ + *** Note: obtained from labah fg_like + Compute the Cbl matrix from the input shape spectrum. + Use transfer matrix instead of transfer function + + This method requires kernels to have been precomputed. + + Arguments + --------- + cls_shape : array_like + The shape spectrum to use. This can be computed using + `get_signal_shape` or otherwise. + map_tag : str + 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. + 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. + fg_ell_ind : float + If binning foreground shape, offset the ell index from the reference + by this amount. + + Returns + ------- + cbl : dict of arrays (num_bins, 2, lmax + 1) + The Cbl matrix, indexed by component and spectrum, then by map + cross, e.g. `cbl['cmb_tt']['map1:map2']` + E/B mixing terms are stored in elements `cbl[:, 1, :]`, + and unmixed terms are stored in elements `cbl[:, 0, :]`. + """ + from spider_analysis import nsi + + T = nsi.TransferMatrix() + + map_pairs = None + if map_tag is not None: + if map_tag in self.map_pairs: + map_pairs = {map_tag: self.map_pairs[map_tag]} + map_tags = list(set(self.map_pairs[map_tag])) + else: + map_tags = [map_tag] + else: + map_tags = self.map_tags + + if map_pairs is None: + map_pairs = pt.tag_pairs(map_tags) + + specs = list(self.specs) + + transfer_bins = np.arange( + 8, 8 + 25 * 13, 25 + ) # change this once transfer matrix adds lowest bin! + nbins_transfer = len(transfer_bins) - 1 + lmin = min(transfer_bins) + lmax = max(transfer_bins) + lmax_kern = lmax + + ls = slice(2, lmax + 1) + cbl = OrderedDict() + + comps = [] + if "cmb_tt" in cls_shape or "cmb_ee" in cls_shape: + comps += ["cmb"] + if "fg" 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_noise0 = self.cls_noise0_null if self.null_run else self.cls_noise0 + cls_noise1 = self.cls_noise1_null if self.null_run else self.cls_noise1 + cls_sxn0 = self.cls_sxn0_null if self.null_run else self.cls_sxn0 + cls_sxn1 = self.cls_sxn1_null if self.null_run else self.cls_sxn1 + cls_nxs0 = self.cls_nxs0_null if self.null_run else self.cls_nxs0 + cls_nxs1 = self.cls_nxs1_null if self.null_run else self.cls_nxs1 + + ell = np.arange(lmax_kern + 1) + lfac = ell * (ell + 1) / 2.0 / np.pi + if self.weighted_bins: + lfac_binned = stats.binned_statistic(ell, lfac, bins=transfer_bins)[0] + + def binup(d, left, right): + return d[..., left:right].sum(axis=-1) + + def bin_things(comp, d, md, d_b1, d_b2, d_b3, md_b1, md_b2, md_b3): + if "res" in comp: + return + for si, spec in enumerate(specs): + stag = "{}_{}".format(comp, spec) + cbl.setdefault(stag, OrderedDict()) + mstag = None + if spec in ["ee", "bb"]: + mstag = stag + "_mix" + cbl.setdefault(mstag, OrderedDict()) + bd = self.bin_def[stag] + for xi, (xname, (tag1, tag2)) in enumerate(map_pairs.items()): + if beam_error: + cbl[stag].setdefault(xname, OrderedDict()) + cbl[stag][xname]["b1"] = np.zeros((len(bd), lmax + 1)) + cbl[stag][xname]["b2"] = np.zeros((len(bd), lmax + 1)) + cbl[stag][xname]["b3"] = np.zeros((len(bd), lmax + 1)) + else: + cbl[stag][xname] = np.zeros((len(bd), lmax + 1)) + if spec in ["ee", "bb"]: + if beam_error: + cbl[mstag].setdefault(xname, OrderedDict()) + cbl[mstag][xname]["b1"] = np.zeros((len(bd), lmax + 1)) + cbl[mstag][xname]["b2"] = np.zeros((len(bd), lmax + 1)) + cbl[mstag][xname]["b3"] = np.zeros((len(bd), lmax + 1)) + else: + cbl[mstag][xname] = np.zeros((len(bd), lmax + 1)) + + # integrate per bin + for idx, (left, right) in enumerate(bd): + if beam_error: + cbl[stag][xname]["b1"][idx, ls] = binup( + d_b1[:, si, xi], left, right + ) + cbl[stag][xname]["b2"][idx, ls] = binup( + d_b2[:, si, xi], left, right + ) + cbl[stag][xname]["b3"][idx, ls] = binup( + d_b3[:, si, xi], left, right + ) + else: + cbl[stag][xname][idx, ls] = binup(d[:, si, xi], left, right) + if spec in ["ee", "bb"]: + if beam_error: + cbl[mstag][xname]["b1"][idx, ls] = binup( + md_b1[:, si - 1, xi], left, right + ) + cbl[mstag][xname]["b2"][idx, ls] = binup( + md_b2[:, si - 1, xi], left, right + ) + cbl[mstag][xname]["b3"][idx, ls] = binup( + md_b3[:, si - 1, xi], left, right + ) + else: + cbl[mstag][xname][idx, ls] = binup( + md[:, si - 1, xi], left, right + ) + + for comp in comps: + # convert to matrices to do multiplication to speed things up, + # except for res is weird so don't do it for that. + # need n_xname x n_spec x ell + nspec = len(specs) + nxmap = len(map_pairs.items()) + if comp == "fg" and fg_ell_ind != 0: + s_arr = (ell / 80.0) ** fg_ell_ind + s_arr[0] = 0 + if not beam_error: + # don't create a new object in memory each time + # use last one's space to save runtime + self.d = np.multiply(self.d_fg, s_arr, out=getattr(self, "d", None)) + self.md = np.multiply( + self.md_fg, s_arr, out=getattr(self, "md", None) + ) + bin_things( + comp, self.d, self.md, None, None, None, None, None, None + ) + else: + self.d_b1 = np.multiply( + self.d_fg_b1, s_arr, out=getattr(self, "d_b1", None) + ) + self.d_b2 = np.multiply( + self.d_fg_b2, s_arr, out=getattr(self, "d_b2", None) + ) + self.d_b3 = np.multiply( + self.d_fg_b3, s_arr, out=getattr(self, "d_b3", None) + ) + self.md_b1 = np.multiply( + self.md_fg_b1, s_arr, out=getattr(self, "md_b1", None) + ) + self.md_b2 = np.multiply( + self.md_fg_b2, s_arr, out=getattr(self, "md_b2", None) + ) + self.md_b3 = np.multiply( + self.md_fg_b3, s_arr, out=getattr(self, "md_b3", None) + ) + bin_things( + comp, + None, + None, + self.d_b1, + self.d_b2, + self.d_b3, + self.md_b1, + self.md_b2, + self.md_b3, + ) + else: + k_arr = np.zeros([nspec, nxmap, lmax - 1, lmax_kern + 1]) + mk_arr = np.zeros([2, nxmap, lmax - 1, lmax_kern + 1]) + # because transfer matrix includes flbl2, combine those, just store + # xfermat * binned shape + fbs_arr = np.zeros([nspec, nxmap, lmax_kern + 1]) + + for xi, (xname, (tag1, tag2)) in enumerate(map_pairs.items()): + for si, spec in enumerate(specs): + stag = "{}_{}".format(comp, spec) + mstag = None + if comp != "res" and spec in ["ee", "bb"]: + mstag = stag + "_mix" + + if "res" in comp: + s0, s1 = spec + res_tags = { + "s0m0": "res_{}_{}".format(s0 * 2, tag1), + "s0m1": "res_{}_{}".format(s0 * 2, tag2), + "s1m0": "res_{}_{}".format(s1 * 2, tag1), + "s1m1": "res_{}_{}".format(s1 * 2, tag2), + } + bd = [[0, lmax + 1]] + # if any component of XY spec is in residual bin + # def, use that bin def + for k, v in res_tags.items(): + spec0 = v.split("_")[1] + if v not in self.bin_def: + if spec0 in ["ee", "bb"]: + v = v.replace(spec0, "eebb") + if v in self.bin_def: + bd = self.bin_def[v] + else: + bd = self.bin_def[v] + for comp in [ + "res0_nxn", + "res1_nxn", + "res0_sxn", + "res1_sxn", + "res0_nxs", + "res1_nxs", + "res", + ]: + stag = "{}_{}".format(comp, spec) + cbl.setdefault(stag, OrderedDict()) + cbl[stag][xname] = np.zeros((len(bd), lmax + 1)) + if comp == "res0_nxn": + cl1 = cls_noise0[spec][xname] + elif comp == "res1_nxn": + cl1 = cls_noise1[spec][xname] + elif comp == "res0_sxn": + cl1 = cls_sxn0[spec][xname] + elif comp == "res1_sxn": + cl1 = cls_sxn1[spec][xname] + elif comp == "res0_nxs": + cl1 = cls_nxs0[spec][xname] + elif comp == "res1_nxs": + cl1 = cls_sxn1[spec][xname] + elif comp == "res": + cl1 = cls_noise[spec][xname] + for idx, (left, right) in enumerate(bd): + lls = slice(left, right) + cbl[stag][xname][idx, lls] = np.copy(cl1[lls]) + + continue + + xfermat = T.get( + self.map_reobs_freqs[tag1], self.map_reobs_freqs[tag1] + ) + + # use correct shape spectrum + if comp == "fg": + # single foreground spectrum + s_arr = ( + cls_shape["fg"][: lmax_kern + 1] + * (ell / 80.0) ** fg_ell_ind + ) + s_arr[0] = 0 + if self.weighted_bins: + s_arr = stats.binned_statistic( + ell, lfac * s_arr, bins=transfer_bins + )[0] + s_arr /= lfac_binned + else: + s_arr = stats.binned_statistic( + ell, s_arr, bins=transfer_bins + )[0] + + # repeat for all specs + s_arr = np.tile(s_arr, len(specs)) + else: + s_arr = np.zeros(nbins_transfer * len(specs)) + for si, spec in enumerate(specs): + s_spec = cls_shape["cmb_{}".format(spec)][: lmax_kern + 1] + if self.weighted_bins: + s_spec = stats.binned_statistic( + ell, lfac * s_spec, bins=transfer_bins + )[0] + s_spec /= lfac_binned + else: + s_spec = stats.binned_statistic( + ell, s_spec, bins=transfer_bins + )[0] + s_arr[ + si * nbins_transfer : (si + 1) * nbins_transfer + ] = s_spec + binned_fbs = np.matmul(np.linalg.inv(xfermat), s_arr) + # fbs_arr = np.zeros([nspec, nxmap, lmax_kern + 1]) + for si, spec in enumerate(specs): + for ib, b0 in enumerate(transfer_bins[:-1]): + fbs_arr[si][xi][b0 : transfer_bins[ib + 1]] = binned_fbs[ib] + + # get cross spectrum kernel terms + if spec == "tt": + k_arr[si, xi] = self.kern[xname][ls, : lmax_kern + 1] + elif spec in ["ee", "bb"]: + k_arr[si, xi] = self.pkern[xname][ls, : lmax_kern + 1] + mk_arr[si - 1, xi] = self.mkern[xname][ls, : lmax_kern + 1] + elif spec in ["te", "tb"]: + k_arr[si, xi] = self.xkern[xname][ls, : lmax_kern + 1] + elif spec == "eb": + k_arr[si, xi] = ( + self.pkern[xname][ls] - self.mkern[xname][ls] + )[:, : lmax_kern + 1] + # need last 3 dims of kernel to match other arrays + k_arr = np.transpose(k_arr, axes=[2, 0, 1, 3]) + mk_arr = np.transpose(mk_arr, axes=[2, 0, 1, 3]) + if not beam_error: + d = k_arr * fbs_arr + md = mk_arr * fbs_arr[1:3] + if comp == "fg": + self.d_fg = np.copy(d) + self.md_fg = np.copy(md) + d_b1 = None + d_b2 = None + d_b3 = None + md_b1 = None + md_b2 = None + md_b3 = None + else: + d = None + md = None + """ + d_b1 = k_arr * b1_arr * f_arr * s_arr + d_b2 = k_arr * b2_arr * f_arr * s_arr + d_b3 = k_arr * b3_arr * f_arr * s_arr + md_b1 = mk_arr * b1_arr[1:3] * f_arr[1:3] * s_arr_md + md_b2 = mk_arr * b2_arr[1:3] * f_arr[1:3] * s_arr_md + md_b3 = mk_arr * b3_arr[1:3] * f_arr[1:3] * s_arr_md + if comp == "fg": + self.d_fg_b1 = d_b1 + self.d_fg_b2 = d_b2 + self.d_fg_b3 = d_b3 + self.md_fg_b1 = md_b1 + self.md_fg_b2 = md_b2 + self.md_fg_b3 = md_b3 + """ + bin_things(comp, d, md, d_b1, d_b2, d_b3, md_b1, md_b2, md_b3) + return cbl def get_model_spectra( self, qb, cbl, delta=True, res=True, cls_noise=None, cond_noise=None @@ -6415,6 +6773,7 @@ def get_likelihood( converge_criteria=0.01, reset_backend=None, file_tag=None, + use_xfer_mat=True # labah ): """ Explore the likelihood, optionally with an MCMC sampler. @@ -6467,6 +6826,9 @@ def get_likelihood( set to True if the checkpoint has been forced to be rerun. file_tag : string If supplied, appended to the likelihood filename. + use_xfer_mat : bool + If True, use transfer matrix to construct model spectrum. If false, use + transfer function XFaster computes """ for x in [ From ed57c34ce4b0da326bdbef1e17bb35df3ce07db5 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Sun, 30 May 2021 22:00:08 -0400 Subject: [PATCH 03/44] finish adding stuffs from xfaster_class.py from old fg_like branch --- xfaster/xfaster_class.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 891e286f..f79a8274 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -6773,7 +6773,7 @@ def get_likelihood( converge_criteria=0.01, reset_backend=None, file_tag=None, - use_xfer_mat=True # labah + use_xfer_mat=True ): """ Explore the likelihood, optionally with an MCMC sampler. @@ -6927,8 +6927,16 @@ def get_likelihood( weighted_bins=self.weighted_bins, lmin=lmin, lmax=lmax, + use_xfer_mat=use_xfer_mat, ) - + + if use_xfer_mat: + bin_cl_template = self.bin_cl_template_tmat + self.lmax = 308 + self.lmin = 8 + else: + bin_cl_template = self.bin_cl_template + if mcmc and reset_backend is None: ret = self.load_data( save_name, From 259b603af5c1935cff108cf1306b56dc19f696d3 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Mon, 31 May 2021 14:40:21 -0400 Subject: [PATCH 04/44] add use_xfer_mat to xfaster_exec.py --- xfaster/xfaster_exec.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index e90ce51b..03e6ccd1 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -110,6 +110,7 @@ def xfaster_run( do_fake_signal=True, do_fake_noise=True, save_fake_data=False, + use_xfer_mat=False, ): """ Main function for running the XFaster algorithm. @@ -332,6 +333,9 @@ def xfaster_run( If true, use sim_index to set noise seed. If false, always use seed 0. save_fake_data : bool If true, save data_xcorr file to disk for fake data. + use_xfer_mat : bool + If True, use NSI transfer matrix instead of transfer function computed + by XFaster. """ from . import __version__ as version @@ -508,6 +512,7 @@ def xfaster_run( num_walkers=mcmc_walkers, converge_criteria=like_converge_criteria, file_tag=like_tag, + use_xfer_mat=use_xfer_mat, ) config_vars.update(like_opts, "Likelihood Estimation Options") config_vars.remove_option("XFaster General", "like_lmin") @@ -876,6 +881,12 @@ def add_arg(P, name, argtype=None, default=None, short=None, help=None, **kwargs help="Prior on dust ell index different from ref " "(should be [0, sig] for [mean 0, width sig] gaussian)", ) + add_arg( + G, + "use_xfer_mat", + help="Use NSI transfer matrix instead of transfer function" + " XFaster computes", + ) add_arg( G, "output_root", From 6239944c936ec1ee80fbb6cbf1c64b51dea89fa9 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Mon, 31 May 2021 20:39:47 -0400 Subject: [PATCH 05/44] fix typo line 2677 qrb to rqb --- xfaster/xfaster_class.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index f79a8274..502d5125 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -938,7 +938,7 @@ def get_files( signal_type : string The variant of signal simulation to use, typically identified by the input spectrum model used to generate it, e.g 'r0p03'. - signal_type_sim : string + signal_type_sim :string The variant of signal 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 signal is computed. @@ -2674,7 +2674,7 @@ def process_index(files, files2, idx, idx2=None, cache=None, qbf=None): ) ) rqb[bad] = 1 - mod = xft.expand_qb(qrb, rb0) + mod = xft.expand_qb(rqb, rb0) for rf in rfields[srb[1]]: if self.map_tags[idx] == srb[2]: @@ -6826,7 +6826,7 @@ def get_likelihood( set to True if the checkpoint has been forced to be rerun. file_tag : string If supplied, appended to the likelihood filename. - use_xfer_mat : bool + use_xfer_mat : bool If True, use transfer matrix to construct model spectrum. If false, use transfer function XFaster computes """ From 5e170b946a12591d2cef8de8f845c995ab0988e7 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Tue, 1 Jun 2021 22:33:08 -0400 Subject: [PATCH 06/44] modifying bin_cl_template_tmat in xfaster_class to work with the current data structure - ongoing[C --- xfaster/xfaster_class.py | 290 +++++++++++++++++++-------------------- 1 file changed, 140 insertions(+), 150 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 502d5125..e4b1e0ec 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -4587,7 +4587,7 @@ def bin_things(comp, d, md): def bin_cl_template_tmat( self, - cls_shape, + cls_shape=None, map_tag=None, beam_error=False, fg_ell_ind=0, @@ -4623,6 +4623,9 @@ def bin_cl_template_tmat( E/B mixing terms are stored in elements `cbl[:, 1, :]`, and unmixed terms are stored in elements `cbl[:, 0, :]`. """ + if cls_shape is None: + cls_shape = self.cls_shape + from spider_analysis import nsi T = nsi.TransferMatrix() @@ -4650,7 +4653,12 @@ def bin_cl_template_tmat( lmax = max(transfer_bins) lmax_kern = lmax + if beam_error: + beam_error = self.get_beam_errors() + beam_keys = ["b1", "b2", "b3"] + ls = slice(2, lmax + 1) + lk = slice(0, lmax_kern + 1) cbl = OrderedDict() comps = [] @@ -4673,10 +4681,10 @@ def bin_cl_template_tmat( if self.weighted_bins: lfac_binned = stats.binned_statistic(ell, lfac, bins=transfer_bins)[0] - def binup(d, left, right): - return d[..., left:right].sum(axis=-1) + def binup(d, left, right, weights): + return (d[..., left:right] * weights).sum(axis=-1) - def bin_things(comp, d, md, d_b1, d_b2, d_b3, md_b1, md_b2, md_b3): + def bin_things(comp, d, md): if "res" in comp: return for si, spec in enumerate(specs): @@ -4687,51 +4695,50 @@ def bin_things(comp, d, md, d_b1, d_b2, d_b3, md_b1, md_b2, md_b3): mstag = stag + "_mix" cbl.setdefault(mstag, OrderedDict()) bd = self.bin_def[stag] + bw = self.bin_weights[stag] for xi, (xname, (tag1, tag2)) in enumerate(map_pairs.items()): if beam_error: - cbl[stag].setdefault(xname, OrderedDict()) - cbl[stag][xname]["b1"] = np.zeros((len(bd), lmax + 1)) - cbl[stag][xname]["b2"] = np.zeros((len(bd), lmax + 1)) - cbl[stag][xname]["b3"] = np.zeros((len(bd), lmax + 1)) + cbl[mstag][xname] = OrderedDict( + [(k, np.zeros((len(bd), lmax + 1))) for k in beam_keys] + ) else: cbl[stag][xname] = np.zeros((len(bd), lmax + 1)) + if spec in ["ee", "bb"]: if beam_error: - cbl[mstag].setdefault(xname, OrderedDict()) - cbl[mstag][xname]["b1"] = np.zeros((len(bd), lmax + 1)) - cbl[mstag][xname]["b2"] = np.zeros((len(bd), lmax + 1)) - cbl[mstag][xname]["b3"] = np.zeros((len(bd), lmax + 1)) + cbl[mstag][xname] = OrderedDict( + [(k, np.zeros((len(bd), lmax + 1))) for k in beam_keys] + ) else: cbl[mstag][xname] = np.zeros((len(bd), lmax + 1)) # integrate per bin - for idx, (left, right) in enumerate(bd): + for idx, ((left, right), weights) in enumerate(zip(bd, bw)): if beam_error: - cbl[stag][xname]["b1"][idx, ls] = binup( - d_b1[:, si, xi], left, right - ) - cbl[stag][xname]["b2"][idx, ls] = binup( - d_b2[:, si, xi], left, right - ) - cbl[stag][xname]["b3"][idx, ls] = binup( - d_b3[:, si, xi], left, right - ) + for k in beam_keys: + cbl[stag][xname][k][idx, ls] = binup( + d[k][si, xi], left, right, weights + ) + #cbl[stag][xname]["b1"][idx, ls] = binup( + # d_b1[:, si, xi], left, right + #) else: - cbl[stag][xname][idx, ls] = binup(d[:, si, xi], left, right) + cbl[stag][xname][idx, ls] = binup( + d[si, xi], left, right, weights + ) if spec in ["ee", "bb"]: if beam_error: - cbl[mstag][xname]["b1"][idx, ls] = binup( - md_b1[:, si - 1, xi], left, right - ) - cbl[mstag][xname]["b2"][idx, ls] = binup( - md_b2[:, si - 1, xi], left, right - ) - cbl[mstag][xname]["b3"][idx, ls] = binup( - md_b3[:, si - 1, xi], left, right - ) + for k in beam_keys: + cbl[mstag][xname][k][idx, ls] = binup( + md[k][si - 1, xi], left, right, weights + ) + + #cbl[mstag][xname]["b1"][idx, ls] = binup( + # md_b1[:, si - 1, xi], left, right + #) else: cbl[mstag][xname][idx, ls] = binup( - md[:, si - 1, xi], left, right + md[si - 1, xi], left, right, weights ) for comp in comps: @@ -4750,53 +4757,45 @@ def bin_things(comp, d, md, d_b1, d_b2, d_b3, md_b1, md_b2, md_b3): self.md = np.multiply( self.md_fg, s_arr, out=getattr(self, "md", None) ) - bin_things( - comp, self.d, self.md, None, None, None, None, None, None - ) + # Question + #bin_things( + # comp, self.d, self.md, None, None, None, None, None, None + #) else: - self.d_b1 = np.multiply( - self.d_fg_b1, s_arr, out=getattr(self, "d_b1", None) - ) - self.d_b2 = np.multiply( - self.d_fg_b2, s_arr, out=getattr(self, "d_b2", None) - ) - self.d_b3 = np.multiply( - self.d_fg_b3, s_arr, out=getattr(self, "d_b3", None) - ) - self.md_b1 = np.multiply( - self.md_fg_b1, s_arr, out=getattr(self, "md_b1", None) - ) - self.md_b2 = np.multiply( - self.md_fg_b2, s_arr, out=getattr(self, "md_b2", None) - ) - self.md_b3 = np.multiply( - self.md_fg_b3, s_arr, out=getattr(self, "md_b3", None) - ) - bin_things( - comp, - None, - None, - self.d_b1, - self.d_b2, - self.d_b3, - self.md_b1, - self.md_b2, - self.md_b3, - ) + for k in beam_keys: + if not hasattr(self, "d"): + self.d = OrderedDict([(k, None) for k in beam_keys]) + self.md = OrderedDict([(k, None) for k in beam_keys]) + self.d[k] = np.multiply(self.d_fg[k], s_arr, out=self.d[k]) + self.md[k] = np.multiply(self.md_fg[k], s_arr, out=self.md[k]) + #self.d_b1 = np.multiply( + # self.d_fg_b1, s_arr, out=getattr(self, "d_b1", None) + #) + bin_things(comp, self.d, self.md) else: - k_arr = np.zeros([nspec, nxmap, lmax - 1, lmax_kern + 1]) - mk_arr = np.zeros([2, nxmap, lmax - 1, lmax_kern + 1]) + #k_arr = np.zeros([nspec, nxmap, self.lmax - 1, lmax_kern + 1]) + #mk_arr = np.zeros([2, nxmap, lmax - 1, lmax_kern + 1]) + kshape = [nspec, nxmap, self.lmax - 1, lmax_kern + 1] + mkshape = [2] + kshape[1:] + k_arr = np.zeros(kshape) + mk_arr = np.zeros(mkshape) + # because transfer matrix includes flbl2, combine those, just store # xfermat * binned shape - fbs_arr = np.zeros([nspec, nxmap, lmax_kern + 1]) - - for xi, (xname, (tag1, tag2)) in enumerate(map_pairs.items()): - for si, spec in enumerate(specs): - stag = "{}_{}".format(comp, spec) - mstag = None - if comp != "res" and spec in ["ee", "bb"]: - mstag = stag + "_mix" + + #fbs_arr = np.zeros([nspec, nxmap, lmax_kern + 1]) + shape = [nspec, nxmap, 1, lmax_kern + 1] + s_arr = np.zeros(shape) + if beam_error: + b_arr = {k: np.zeros(shape) for k in beam_keys} + for si, spec in enumerate(specs): + stag = "{}_{}".format(comp, spec) + mstag = None + if comp != "res" and spec in ["ee", "bb"]: + mstag = stag + "_mix" + + for xi, (xname, (tag1, tag2)) in enumerate(map_pairs.items()): if "res" in comp: s0, s1 = spec res_tags = { @@ -4817,96 +4816,87 @@ def bin_things(comp, d, md, d_b1, d_b2, d_b3, md_b1, md_b2, md_b3): bd = self.bin_def[v] else: bd = self.bin_def[v] - for comp in [ - "res0_nxn", - "res1_nxn", - "res0_sxn", - "res1_sxn", - "res0_nxs", - "res1_nxs", - "res", + for comp, cls in [ + ("res0_nxn", cls_noise0), + ("res1_nxn", cls_noise1), + ("res0_sxn", cls_sxn0), + ("res1_sxn", cls_sxn1), + ("res0_nxs", cls_nxs0), + ("res1_nxs", cls_nxs1), + ("res", cls_noise), ]: stag = "{}_{}".format(comp, spec) cbl.setdefault(stag, OrderedDict()) cbl[stag][xname] = np.zeros((len(bd), lmax + 1)) - if comp == "res0_nxn": - cl1 = cls_noise0[spec][xname] - elif comp == "res1_nxn": - cl1 = cls_noise1[spec][xname] - elif comp == "res0_sxn": - cl1 = cls_sxn0[spec][xname] - elif comp == "res1_sxn": - cl1 = cls_sxn1[spec][xname] - elif comp == "res0_nxs": - cl1 = cls_nxs0[spec][xname] - elif comp == "res1_nxs": - cl1 = cls_sxn1[spec][xname] - elif comp == "res": - cl1 = cls_noise[spec][xname] + cl1 = cls[spec][xname] for idx, (left, right) in enumerate(bd): lls = slice(left, right) cbl[stag][xname][idx, lls] = np.copy(cl1[lls]) continue - - xfermat = T.get( - self.map_reobs_freqs[tag1], self.map_reobs_freqs[tag1] - ) - - # use correct shape spectrum - if comp == "fg": - # single foreground spectrum - s_arr = ( - cls_shape["fg"][: lmax_kern + 1] - * (ell / 80.0) ** fg_ell_ind + + if beam_error: + b_arr["b1"][si, xi] = beam_error[spec][tag1] + b_arr["b2"][si, xi] = beam_error[spec][tag2] + b_arr["b3"][si, xi] = ( + b_arr["b1"][si, xi] * b_arr["b2"][si, xi] + ) + + xfermat = T.get( + self.map_reobs_freqs[tag1], self.map_reobs_freqs[tag1] ) - s_arr[0] = 0 - if self.weighted_bins: - s_arr = stats.binned_statistic( - ell, lfac * s_arr, bins=transfer_bins - )[0] - s_arr /= lfac_binned - else: - s_arr = stats.binned_statistic( - ell, s_arr, bins=transfer_bins - )[0] - # repeat for all specs - s_arr = np.tile(s_arr, len(specs)) - else: - s_arr = np.zeros(nbins_transfer * len(specs)) - for si, spec in enumerate(specs): - s_spec = cls_shape["cmb_{}".format(spec)][: lmax_kern + 1] + # use correct shape spectrum + if comp == "fg": + # single foreground spectrum + s_arr = cls_shape["fg"][lk] * (ell / 80.0) ** fg_ell_ind + s_arr[0] = 0 if self.weighted_bins: - s_spec = stats.binned_statistic( - ell, lfac * s_spec, bins=transfer_bins + s_arr = stats.binned_statistic( + ell, lfac * s_arr, bins=transfer_bins )[0] - s_spec /= lfac_binned + s_arr /= lfac_binned else: - s_spec = stats.binned_statistic( - ell, s_spec, bins=transfer_bins + s_arr = stats.binned_statistic( + ell, s_arr, bins=transfer_bins )[0] - s_arr[ - si * nbins_transfer : (si + 1) * nbins_transfer - ] = s_spec - binned_fbs = np.matmul(np.linalg.inv(xfermat), s_arr) - # fbs_arr = np.zeros([nspec, nxmap, lmax_kern + 1]) - for si, spec in enumerate(specs): - for ib, b0 in enumerate(transfer_bins[:-1]): - fbs_arr[si][xi][b0 : transfer_bins[ib + 1]] = binned_fbs[ib] - - # get cross spectrum kernel terms - if spec == "tt": - k_arr[si, xi] = self.kern[xname][ls, : lmax_kern + 1] - elif spec in ["ee", "bb"]: - k_arr[si, xi] = self.pkern[xname][ls, : lmax_kern + 1] - mk_arr[si - 1, xi] = self.mkern[xname][ls, : lmax_kern + 1] - elif spec in ["te", "tb"]: - k_arr[si, xi] = self.xkern[xname][ls, : lmax_kern + 1] - elif spec == "eb": - k_arr[si, xi] = ( - self.pkern[xname][ls] - self.mkern[xname][ls] - )[:, : lmax_kern + 1] + # FROM HERE ON ... + # repeat for all specs + s_arr = np.tile(s_arr, len(specs)) + else: + s_arr = np.zeros(nbins_transfer * len(specs)) + for si, spec in enumerate(specs): + s_spec = cls_shape["cmb_{}".format(spec)][: lmax_kern + 1] + if self.weighted_bins: + s_spec = stats.binned_statistic( + ell, lfac * s_spec, bins=transfer_bins + )[0] + s_spec /= lfac_binned + else: + s_spec = stats.binned_statistic( + ell, s_spec, bins=transfer_bins + )[0] + s_arr[ + si * nbins_transfer : (si + 1) * nbins_transfer + ] = s_spec + binned_fbs = np.matmul(np.linalg.inv(xfermat), s_arr) + # fbs_arr = np.zeros([nspec, nxmap, lmax_kern + 1]) + for si, spec in enumerate(specs): + for ib, b0 in enumerate(transfer_bins[:-1]): + fbs_arr[si][xi][b0 : transfer_bins[ib + 1]] = binned_fbs[ib] + + # get cross spectrum kernel terms + if spec == "tt": + k_arr[si, xi] = self.kern[xname][ls, : lmax_kern + 1] + elif spec in ["ee", "bb"]: + k_arr[si, xi] = self.pkern[xname][ls, : lmax_kern + 1] + mk_arr[si - 1, xi] = self.mkern[xname][ls, : lmax_kern + 1] + elif spec in ["te", "tb"]: + k_arr[si, xi] = self.xkern[xname][ls, : lmax_kern + 1] + elif spec == "eb": + k_arr[si, xi] = ( + self.pkern[xname][ls] - self.mkern[xname][ls] + )[:, : lmax_kern + 1] # need last 3 dims of kernel to match other arrays k_arr = np.transpose(k_arr, axes=[2, 0, 1, 3]) mk_arr = np.transpose(mk_arr, axes=[2, 0, 1, 3]) From 0acf9b78d8be4a482746db345e33f8f06a02340f Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Wed, 2 Jun 2021 23:12:02 -0400 Subject: [PATCH 07/44] making data structure in bin_cl_template_tmat compliable to the new package setup --- xfaster/xfaster_class.py | 164 +++++++++++++++++++-------------------- 1 file changed, 82 insertions(+), 82 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index e4b1e0ec..0938c1d2 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -4626,7 +4626,7 @@ def bin_cl_template_tmat( if cls_shape is None: cls_shape = self.cls_shape - from spider_analysis import nsi + from spider_analysis import si T = nsi.TransferMatrix() @@ -4681,8 +4681,8 @@ def bin_cl_template_tmat( if self.weighted_bins: lfac_binned = stats.binned_statistic(ell, lfac, bins=transfer_bins)[0] - def binup(d, left, right, weights): - return (d[..., left:right] * weights).sum(axis=-1) + def binup(d, left, right): + return (d[..., left:right]).sum(axis=-1) def bin_things(comp, d, md): if "res" in comp: @@ -4695,7 +4695,7 @@ def bin_things(comp, d, md): mstag = stag + "_mix" cbl.setdefault(mstag, OrderedDict()) bd = self.bin_def[stag] - bw = self.bin_weights[stag] + #bw = self.bin_weights[stag] for xi, (xname, (tag1, tag2)) in enumerate(map_pairs.items()): if beam_error: cbl[mstag][xname] = OrderedDict( @@ -4713,24 +4713,24 @@ def bin_things(comp, d, md): cbl[mstag][xname] = np.zeros((len(bd), lmax + 1)) # integrate per bin - for idx, ((left, right), weights) in enumerate(zip(bd, bw)): + for idx, (left, right) in enumerate(bd): if beam_error: for k in beam_keys: cbl[stag][xname][k][idx, ls] = binup( - d[k][si, xi], left, right, weights + d[k][si, xi], left, right #, weight ) #cbl[stag][xname]["b1"][idx, ls] = binup( # d_b1[:, si, xi], left, right #) else: cbl[stag][xname][idx, ls] = binup( - d[si, xi], left, right, weights + d[si, xi], left, right #, weights ) if spec in ["ee", "bb"]: if beam_error: for k in beam_keys: cbl[mstag][xname][k][idx, ls] = binup( - md[k][si - 1, xi], left, right, weights + md[k][si - 1, xi], left, right#, weights ) #cbl[mstag][xname]["b1"][idx, ls] = binup( @@ -4738,7 +4738,7 @@ def bin_things(comp, d, md): #) else: cbl[mstag][xname][idx, ls] = binup( - md[si - 1, xi], left, right, weights + md[si - 1, xi], left, right#, weights ) for comp in comps: @@ -4757,7 +4757,6 @@ def bin_things(comp, d, md): self.md = np.multiply( self.md_fg, s_arr, out=getattr(self, "md", None) ) - # Question #bin_things( # comp, self.d, self.md, None, None, None, None, None, None #) @@ -4773,8 +4772,6 @@ def bin_things(comp, d, md): #) bin_things(comp, self.d, self.md) else: - #k_arr = np.zeros([nspec, nxmap, self.lmax - 1, lmax_kern + 1]) - #mk_arr = np.zeros([2, nxmap, lmax - 1, lmax_kern + 1]) kshape = [nspec, nxmap, self.lmax - 1, lmax_kern + 1] mkshape = [2] + kshape[1:] k_arr = np.zeros(kshape) @@ -4783,19 +4780,19 @@ def bin_things(comp, d, md): # because transfer matrix includes flbl2, combine those, just store # xfermat * binned shape - #fbs_arr = np.zeros([nspec, nxmap, lmax_kern + 1]) - shape = [nspec, nxmap, 1, lmax_kern + 1] - s_arr = np.zeros(shape) + fbs_arr = np.zeros([nspec, nxmap, lmax_kern + 1]) + + # for compute d of non-fg with beam error if beam_error: b_arr = {k: np.zeros(shape) for k in beam_keys} - for si, spec in enumerate(specs): - stag = "{}_{}".format(comp, spec) - mstag = None - if comp != "res" and spec in ["ee", "bb"]: - mstag = stag + "_mix" - - for xi, (xname, (tag1, tag2)) in enumerate(map_pairs.items()): + for xi, (xname, (tag1, tag2)) in enumerate(map_pairs.items()): + for si, spec in enumerate(specs): + stag = "{}_{}".format(comp, spec) + mstag = None + if comp != "res" and spec in ["ee", "bb"]: + mstag = stag + "_mix" + if "res" in comp: s0, s1 = spec res_tags = { @@ -4841,62 +4838,63 @@ def bin_things(comp, d, md): b_arr["b3"][si, xi] = ( b_arr["b1"][si, xi] * b_arr["b2"][si, xi] ) - - xfermat = T.get( - self.map_reobs_freqs[tag1], self.map_reobs_freqs[tag1] - ) - # use correct shape spectrum - if comp == "fg": - # single foreground spectrum - s_arr = cls_shape["fg"][lk] * (ell / 80.0) ** fg_ell_ind - s_arr[0] = 0 + xfermat = T.get( + self.map_reobs_freqs[tag1], self.map_reobs_freqs[tag1] + ) + + # use correct shape spectrum + if comp == "fg": + # single foreground spectrum + s_arr = cls_shape["fg"][lk] * (ell / 80.0) ** fg_ell_ind + s_arr[0] = 0 + if self.weighted_bins: + s_arr = stats.binned_statistic( + ell, lfac * s_arr, bins=transfer_bins + )[0] + s_arr /= lfac_binned + else: + s_arr = stats.binned_statistic( + ell, s_arr, bins=transfer_bins + )[0] + # repeat for all specs + s_arr = np.tile(s_arr, len(specs)) + else: + s_arr = np.zeros(nbins_transfer * len(specs)) + for si, spec in enumerate(specs): + s_spec = cls_shape["cmb_{}".format(spec)][: lk] if self.weighted_bins: - s_arr = stats.binned_statistic( - ell, lfac * s_arr, bins=transfer_bins + s_spec = stats.binned_statistic( + ell, lfac * s_spec, bins=transfer_bins )[0] - s_arr /= lfac_binned + s_spec /= lfac_binned else: - s_arr = stats.binned_statistic( - ell, s_arr, bins=transfer_bins + s_spec = stats.binned_statistic( + ell, s_spec, bins=transfer_bins )[0] - # FROM HERE ON ... - # repeat for all specs - s_arr = np.tile(s_arr, len(specs)) - else: - s_arr = np.zeros(nbins_transfer * len(specs)) - for si, spec in enumerate(specs): - s_spec = cls_shape["cmb_{}".format(spec)][: lmax_kern + 1] - if self.weighted_bins: - s_spec = stats.binned_statistic( - ell, lfac * s_spec, bins=transfer_bins - )[0] - s_spec /= lfac_binned - else: - s_spec = stats.binned_statistic( - ell, s_spec, bins=transfer_bins - )[0] - s_arr[ - si * nbins_transfer : (si + 1) * nbins_transfer - ] = s_spec - binned_fbs = np.matmul(np.linalg.inv(xfermat), s_arr) - # fbs_arr = np.zeros([nspec, nxmap, lmax_kern + 1]) - for si, spec in enumerate(specs): - for ib, b0 in enumerate(transfer_bins[:-1]): - fbs_arr[si][xi][b0 : transfer_bins[ib + 1]] = binned_fbs[ib] - - # get cross spectrum kernel terms - if spec == "tt": - k_arr[si, xi] = self.kern[xname][ls, : lmax_kern + 1] - elif spec in ["ee", "bb"]: - k_arr[si, xi] = self.pkern[xname][ls, : lmax_kern + 1] - mk_arr[si - 1, xi] = self.mkern[xname][ls, : lmax_kern + 1] - elif spec in ["te", "tb"]: - k_arr[si, xi] = self.xkern[xname][ls, : lmax_kern + 1] - elif spec == "eb": - k_arr[si, xi] = ( - self.pkern[xname][ls] - self.mkern[xname][ls] - )[:, : lmax_kern + 1] + s_arr[ + si * nbins_transfer : (si + 1) * nbins_transfer + ] = s_spec + binned_fbs = np.matmul(np.linalg.inv(xfermat), s_arr) + # fbs_arr = np.zeros([nspec, nxmap, lmax_kern + 1]) + for si, spec in enumerate(specs): + for ib, b0 in enumerate(transfer_bins[:-1]): + fbs_arr[si][xi][b0 : transfer_bins[ib + 1]] = binned_fbs[ib] + + # get cross spectrum kernel terms + if spec == "tt": + k_arr[si, xi] = self.kern[xname][ls, : lmax_kern + 1] + elif spec in ["ee", "bb"]: + k_arr[si, xi] = self.pkern[xname][ls, : lmax_kern + 1] + mk_arr[si - 1, xi] = self.mkern[xname][ls, : lmax_kern + 1] + elif spec in ["te", "tb"]: + k_arr[si, xi] = self.xkern[xname][ls, : lmax_kern + 1] + elif spec == "eb": + k_arr[si, xi] = ( + self.pkern[xname][ls] - self.mkern[xname][ls] + )[:, : lmax_kern + 1] + + # need last 3 dims of kernel to match other arrays k_arr = np.transpose(k_arr, axes=[2, 0, 1, 3]) mk_arr = np.transpose(mk_arr, axes=[2, 0, 1, 3]) @@ -4913,6 +4911,8 @@ def bin_things(comp, d, md): md_b2 = None md_b3 = None else: + # why not include error here + d = None md = None """ @@ -6988,7 +6988,7 @@ def get_likelihood( if r_prior is None: 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) + cbl = bin_cl_template(map_tag=map_tag) cls_model = self.get_model_spectra(qb, cbl, delta=True, cls_noise=cls_noise) else: qb = copy.deepcopy(qb) @@ -7009,23 +7009,23 @@ def get_likelihood( ) # load tensor and scalar terms separately - cbl_scalar = self.bin_cl_template(cls_shape_scalar, map_tag) + cbl_scalar = bin_cl_template(cls_shape_scalar, map_tag) cls_model_scalar = self.get_model_spectra( qb, cbl_scalar, delta=True, cls_noise=cls_noise ) - cbl_tensor = self.bin_cl_template(cls_shape_tensor, map_tag) + cbl_tensor = bin_cl_template(cls_shape_tensor, map_tag) cls_model_tensor = self.get_model_spectra( qb, cbl_tensor, delta=False, res=False ) if beam_prior is not None: # load beam error term for tensor and scalar - cbl_scalar_beam = self.bin_cl_template( + cbl_scalar_beam = 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 ) - cbl_tensor_beam = self.bin_cl_template( + cbl_tensor_beam = bin_cl_template( cls_shape_tensor, map_tag, beam_error=True ) cls_mod_tens_beam = self.get_model_spectra( @@ -7037,9 +7037,9 @@ def get_likelihood( cls_shape_dust = self.get_signal_shape(save=False, component="fg") # if dust_ellind_prior is None: # # can preload shape since not varying ell index - cbl_fg = self.bin_cl_template(cls_shape_dust, map_tag=map_tag) + cbl_fg = bin_cl_template(cls_shape_dust, map_tag=map_tag) if beam_prior is not None: - cbl_fg_beam = self.bin_cl_template( + cbl_fg_beam = bin_cl_template( cls_shape_dust, map_tag, beam_error=True ) @@ -7205,11 +7205,11 @@ def log_like( else: qb["delta_beta"][:] = betad if dust_ellind is not None: - cbl_fg0 = self.bin_cl_template( + cbl_fg0 = bin_cl_template( cls_shape_dust, map_tag=map_tag, fg_ell_ind=dust_ellind ) if beam is not None: - cbl_fg_beam0 = self.bin_cl_template( + cbl_fg_beam0 = bin_cl_template( cls_shape_dust, map_tag, fg_ell_ind=dust_ellind, From 7d718bda5f350a6ef4cee11fe3a35eb135dd2e89 Mon Sep 17 00:00:00 2001 From: "Anne E. Gambrel" Date: Thu, 10 Jun 2021 21:44:38 +0200 Subject: [PATCH 08/44] Load up shape spectra properly if including foregrounds --- xfaster/xfaster_class.py | 171 ++++++++++++++++++++------------------- 1 file changed, 87 insertions(+), 84 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 0938c1d2..c4ce8456 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -3791,9 +3791,6 @@ def get_signal_shape( if not self.null_run and "fg_tt" in self.bin_def: specs.append("fg") - if component == "fg": - specs = ["fg"] - nspecs = len(specs) if save: @@ -3921,6 +3918,9 @@ def get_signal_shape( if not masked: cls_shape[csk] *= 1.0e-12 + if component == "fg": + cls_shape = {"fg": cls_shape["fg"]} + if save: self.cls_shape = cls_shape opts["cls_shape"] = cls_shape @@ -4584,7 +4584,7 @@ def bin_things(comp, d, md): bin_things(comp, d, md) return cbl - + def bin_cl_template_tmat( self, cls_shape=None, @@ -4695,7 +4695,7 @@ def bin_things(comp, d, md): mstag = stag + "_mix" cbl.setdefault(mstag, OrderedDict()) bd = self.bin_def[stag] - #bw = self.bin_weights[stag] + # bw = self.bin_weights[stag] for xi, (xname, (tag1, tag2)) in enumerate(map_pairs.items()): if beam_error: cbl[mstag][xname] = OrderedDict( @@ -4703,7 +4703,7 @@ def bin_things(comp, d, md): ) else: cbl[stag][xname] = np.zeros((len(bd), lmax + 1)) - + if spec in ["ee", "bb"]: if beam_error: cbl[mstag][xname] = OrderedDict( @@ -4717,28 +4717,28 @@ def bin_things(comp, d, md): if beam_error: for k in beam_keys: cbl[stag][xname][k][idx, ls] = binup( - d[k][si, xi], left, right #, weight + d[k][si, xi], left, right # , weight ) - #cbl[stag][xname]["b1"][idx, ls] = binup( + # cbl[stag][xname]["b1"][idx, ls] = binup( # d_b1[:, si, xi], left, right - #) + # ) else: cbl[stag][xname][idx, ls] = binup( - d[si, xi], left, right #, weights + d[si, xi], left, right # , weights ) if spec in ["ee", "bb"]: if beam_error: for k in beam_keys: cbl[mstag][xname][k][idx, ls] = binup( - md[k][si - 1, xi], left, right#, weights + md[k][si - 1, xi], left, right # , weights ) - - #cbl[mstag][xname]["b1"][idx, ls] = binup( + + # cbl[mstag][xname]["b1"][idx, ls] = binup( # md_b1[:, si - 1, xi], left, right - #) + # ) else: cbl[mstag][xname][idx, ls] = binup( - md[si - 1, xi], left, right#, weights + md[si - 1, xi], left, right # , weights ) for comp in comps: @@ -4757,9 +4757,9 @@ def bin_things(comp, d, md): self.md = np.multiply( self.md_fg, s_arr, out=getattr(self, "md", None) ) - #bin_things( + # bin_things( # comp, self.d, self.md, None, None, None, None, None, None - #) + # ) else: for k in beam_keys: if not hasattr(self, "d"): @@ -4767,9 +4767,9 @@ def bin_things(comp, d, md): self.md = OrderedDict([(k, None) for k in beam_keys]) self.d[k] = np.multiply(self.d_fg[k], s_arr, out=self.d[k]) self.md[k] = np.multiply(self.md_fg[k], s_arr, out=self.md[k]) - #self.d_b1 = np.multiply( + # self.d_b1 = np.multiply( # self.d_fg_b1, s_arr, out=getattr(self, "d_b1", None) - #) + # ) bin_things(comp, self.d, self.md) else: kshape = [nspec, nxmap, self.lmax - 1, lmax_kern + 1] @@ -4779,10 +4779,10 @@ def bin_things(comp, d, md): # because transfer matrix includes flbl2, combine those, just store # xfermat * binned shape - + fbs_arr = np.zeros([nspec, nxmap, lmax_kern + 1]) - - # for compute d of non-fg with beam error + + # for compute d of non-fg with beam error if beam_error: b_arr = {k: np.zeros(shape) for k in beam_keys} @@ -4792,7 +4792,7 @@ def bin_things(comp, d, md): mstag = None if comp != "res" and spec in ["ee", "bb"]: mstag = stag + "_mix" - + if "res" in comp: s0, s1 = spec res_tags = { @@ -4831,7 +4831,7 @@ def bin_things(comp, d, md): cbl[stag][xname][idx, lls] = np.copy(cl1[lls]) continue - + if beam_error: b_arr["b1"][si, xi] = beam_error[spec][tag1] b_arr["b2"][si, xi] = beam_error[spec][tag2] @@ -4862,7 +4862,7 @@ def bin_things(comp, d, md): else: s_arr = np.zeros(nbins_transfer * len(specs)) for si, spec in enumerate(specs): - s_spec = cls_shape["cmb_{}".format(spec)][: lk] + s_spec = cls_shape["cmb_{}".format(spec)][:lk] if self.weighted_bins: s_spec = stats.binned_statistic( ell, lfac * s_spec, bins=transfer_bins @@ -4893,7 +4893,6 @@ def bin_things(comp, d, md): k_arr[si, xi] = ( self.pkern[xname][ls] - self.mkern[xname][ls] )[:, : lmax_kern + 1] - # need last 3 dims of kernel to match other arrays k_arr = np.transpose(k_arr, axes=[2, 0, 1, 3]) @@ -4912,7 +4911,7 @@ def bin_things(comp, d, md): md_b3 = None else: # why not include error here - + d = None md = None """ @@ -6392,7 +6391,8 @@ def fisher_iterate( # use max likelihood qb's qb2cb to convert qb->cb if "res" not in stag: cb_arr[iq] = np.einsum( - "ij,j->i", qb2cb[stag], qb1[stag])[ibin] + "ij,j->i", qb2cb[stag], qb1[stag] + )[ibin] try: like = self.fisher_calc( qb1, @@ -6763,62 +6763,62 @@ def get_likelihood( converge_criteria=0.01, reset_backend=None, file_tag=None, - use_xfer_mat=True + use_xfer_mat=True, ): """ - Explore the likelihood, optionally with an MCMC sampler. - - Arguments - --------- - qb : OrderedDict - Bandpower parameters previously computed by Fisher iteration. - inv_fish : array_like - Inverse Fisher matrix computed with the input qb's. - map_tag : string - If not None, then the likelihood is sampled using the spectra - corresponding to the given map, rather over all possible - combinations of map-map cross-spectra. The input qb's and inv_fish - must have been computed with the same option. - mcmc : bool - If True, sample the likelihood using an MCMC sampler. Remaining options - determine parameter space and sampler configuration. - r_prior : 2-list or None - Prior upper and lower bound on tensor to scalar ratio. If None, the - fiducial shape spectrum is assumed, and the r parameter space is not - varied. - alpha_prior : 2-list or None - Prior upper and lower bound on template coefficients. If None, the - alpha parameter space is not varied. - res_prior : 2-list or none - Prior upper and lower bound on residual qbs. If None, the - res parameter space is not varied. - beam_prior : 2-list or none - Prior mean and width of gaussian width on beam error (when - multiplied by beam error envelope). If None, the - beam parameter space is not varied. - betad_prior : 2-list or none - Prior mean and width of gaussian width on dust spectral index. - If None, the dust index parameter space is not varied. - dust_amp_prior : 2-list or none - Prior upper and lower bound on dust amplitude. - If None, the dust amp parameter space is not varied. - dust_ellind_prior : 2-list or none - Prior mean and width of Gaussian prior on difference in dust ell - power law index. If None, don't vary from reference if fitting dust - power spectrum model. - num_walkers : int - Number of unique walkers with which to sample the parameter space. - num_steps : int - Number of steps each walker should take in sampling the parameter space. - reset_backend : bool - If True, clear the backend buffer before sampling. If False, - samples are appended to the existing buffer. If not supplied, - set to True if the checkpoint has been forced to be rerun. - file_tag : string - If supplied, appended to the likelihood filename. - use_xfer_mat : bool - If True, use transfer matrix to construct model spectrum. If false, use - transfer function XFaster computes + Explore the likelihood, optionally with an MCMC sampler. + + Arguments + --------- + qb : OrderedDict + Bandpower parameters previously computed by Fisher iteration. + inv_fish : array_like + Inverse Fisher matrix computed with the input qb's. + map_tag : string + If not None, then the likelihood is sampled using the spectra + corresponding to the given map, rather over all possible + combinations of map-map cross-spectra. The input qb's and inv_fish + must have been computed with the same option. + mcmc : bool + If True, sample the likelihood using an MCMC sampler. Remaining options + determine parameter space and sampler configuration. + r_prior : 2-list or None + Prior upper and lower bound on tensor to scalar ratio. If None, the + fiducial shape spectrum is assumed, and the r parameter space is not + varied. + alpha_prior : 2-list or None + Prior upper and lower bound on template coefficients. If None, the + alpha parameter space is not varied. + res_prior : 2-list or none + Prior upper and lower bound on residual qbs. If None, the + res parameter space is not varied. + beam_prior : 2-list or none + Prior mean and width of gaussian width on beam error (when + multiplied by beam error envelope). If None, the + beam parameter space is not varied. + betad_prior : 2-list or none + Prior mean and width of gaussian width on dust spectral index. + If None, the dust index parameter space is not varied. + dust_amp_prior : 2-list or none + Prior upper and lower bound on dust amplitude. + If None, the dust amp parameter space is not varied. + dust_ellind_prior : 2-list or none + Prior mean and width of Gaussian prior on difference in dust ell + power law index. If None, don't vary from reference if fitting dust + power spectrum model. + num_walkers : int + Number of unique walkers with which to sample the parameter space. + num_steps : int + Number of steps each walker should take in sampling the parameter space. + reset_backend : bool + If True, clear the backend buffer before sampling. If False, + samples are appended to the existing buffer. If not supplied, + set to True if the checkpoint has been forced to be rerun. + file_tag : string + If supplied, appended to the likelihood filename. + use_xfer_mat : bool + If True, use transfer matrix to construct model spectrum. If false, use + transfer function XFaster computes """ for x in [ @@ -6895,6 +6895,9 @@ def get_likelihood( "dust_amp_prior": dust_amp_prior, "dust_ellind_prior": dust_ellind_prior, } + + self.log("Priors: {}".format(priors), "debug") + # priors on quantities that affect Dmat_obs or gmat (precalculated) obs_priors = [alpha_prior] @@ -6919,14 +6922,14 @@ def get_likelihood( lmax=lmax, use_xfer_mat=use_xfer_mat, ) - + if use_xfer_mat: bin_cl_template = self.bin_cl_template_tmat self.lmax = 308 self.lmin = 8 else: bin_cl_template = self.bin_cl_template - + if mcmc and reset_backend is None: ret = self.load_data( save_name, From c08c2958e119f33c188994b8cc96a2b519e62346 Mon Sep 17 00:00:00 2001 From: Anne Gambrel Date: Mon, 12 Jul 2021 15:56:09 -0600 Subject: [PATCH 09/44] Clean up bin_cl_template_tmat function --- xfaster/xfaster_class.py | 67 ++++++++++------------------------------ 1 file changed, 16 insertions(+), 51 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index ec6e1c33..3f11f438 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -4756,14 +4756,16 @@ def bin_cl_template_tmat( cbl : dict of arrays (num_bins, 2, lmax + 1) The Cbl matrix, indexed by component and spectrum, then by map cross, e.g. `cbl['cmb_tt']['map1:map2']` - E/B mixing terms are stored in elements `cbl[:, 1, :]`, - and unmixed terms are stored in elements `cbl[:, 0, :]`. + E/B mixing terms are stored in elements ``cbl['cmb_ee_mix']`` and + ``cbl['cmb_bb_mix']``, and unmixed terms are stored in elements + ``cbl['cmb_ee']`` and ``cbl['cmb_bb']``. """ if cls_shape is None: cls_shape = self.cls_shape from spider_analysis import si + # Need to update this to use ell-by-ell anafast transfer matrix T = nsi.TransferMatrix() map_pairs = None @@ -4781,6 +4783,7 @@ def bin_cl_template_tmat( specs = list(self.specs) + # This will need to change for ell-by-ell transfer matrix transfer_bins = np.arange( 8, 8 + 25 * 13, 25 ) # change this once transfer matrix adds lowest bin! @@ -4813,12 +4816,9 @@ def bin_cl_template_tmat( cls_nxs1 = self.cls_nxs1_null if self.null_run else self.cls_nxs1 ell = np.arange(lmax_kern + 1) - lfac = ell * (ell + 1) / 2.0 / np.pi - if self.weighted_bins: - lfac_binned = stats.binned_statistic(ell, lfac, bins=transfer_bins)[0] - def binup(d, left, right): - return (d[..., left:right]).sum(axis=-1) + def binup(d, left, right, weights): + return (d[..., left:right] * weights).sum(axis=-1) def bin_things(comp, d, md): if "res" in comp: @@ -4831,7 +4831,7 @@ def bin_things(comp, d, md): mstag = stag + "_mix" cbl.setdefault(mstag, OrderedDict()) bd = self.bin_def[stag] - # bw = self.bin_weights[stag] + bw = self.bin_weights[stag] for xi, (xname, (tag1, tag2)) in enumerate(map_pairs.items()): if beam_error: cbl[mstag][xname] = OrderedDict( @@ -4849,32 +4849,25 @@ def bin_things(comp, d, md): cbl[mstag][xname] = np.zeros((len(bd), lmax + 1)) # integrate per bin - for idx, (left, right) in enumerate(bd): + for idx, ((left, right), weights) in enumerate(zip(bd, bw)): if beam_error: for k in beam_keys: cbl[stag][xname][k][idx, ls] = binup( - d[k][si, xi], left, right # , weight + d[k][si, xi], left, right, weights ) - # cbl[stag][xname]["b1"][idx, ls] = binup( - # d_b1[:, si, xi], left, right - # ) else: cbl[stag][xname][idx, ls] = binup( - d[si, xi], left, right # , weights + d[si, xi], left, right, weights ) if spec in ["ee", "bb"]: if beam_error: for k in beam_keys: cbl[mstag][xname][k][idx, ls] = binup( - md[k][si - 1, xi], left, right # , weights + md[k][si - 1, xi], left, right, weights ) - - # cbl[mstag][xname]["b1"][idx, ls] = binup( - # md_b1[:, si - 1, xi], left, right - # ) else: cbl[mstag][xname][idx, ls] = binup( - md[si - 1, xi], left, right # , weights + md[si - 1, xi], left, right, weights ) for comp in comps: @@ -4882,7 +4875,7 @@ def bin_things(comp, d, md): # except for res is weird so don't do it for that. # need n_xname x n_spec x ell nspec = len(specs) - nxmap = len(map_pairs.items()) + nxmap = len(map_pairs) if comp == "fg" and fg_ell_ind != 0: s_arr = (ell / 80.0) ** fg_ell_ind s_arr[0] = 0 @@ -4892,10 +4885,6 @@ def bin_things(comp, d, md): self.d = np.multiply(self.d_fg, s_arr, out=getattr(self, "d", None)) self.md = np.multiply( self.md_fg, s_arr, out=getattr(self, "md", None) - ) - # bin_things( - # comp, self.d, self.md, None, None, None, None, None, None - # ) else: for k in beam_keys: if not hasattr(self, "d"): @@ -4903,10 +4892,8 @@ def bin_things(comp, d, md): self.md = OrderedDict([(k, None) for k in beam_keys]) self.d[k] = np.multiply(self.d_fg[k], s_arr, out=self.d[k]) self.md[k] = np.multiply(self.md_fg[k], s_arr, out=self.md[k]) - # self.d_b1 = np.multiply( - # self.d_fg_b1, s_arr, out=getattr(self, "d_b1", None) - # ) bin_things(comp, self.d, self.md) + else: kshape = [nspec, nxmap, self.lmax - 1, lmax_kern + 1] mkshape = [2] + kshape[1:] @@ -4918,10 +4905,6 @@ def bin_things(comp, d, md): fbs_arr = np.zeros([nspec, nxmap, lmax_kern + 1]) - # for compute d of non-fg with beam error - if beam_error: - b_arr = {k: np.zeros(shape) for k in beam_keys} - for xi, (xname, (tag1, tag2)) in enumerate(map_pairs.items()): for si, spec in enumerate(specs): stag = "{}_{}".format(comp, spec) @@ -4984,30 +4967,12 @@ def bin_things(comp, d, md): # single foreground spectrum s_arr = cls_shape["fg"][lk] * (ell / 80.0) ** fg_ell_ind s_arr[0] = 0 - if self.weighted_bins: - s_arr = stats.binned_statistic( - ell, lfac * s_arr, bins=transfer_bins - )[0] - s_arr /= lfac_binned - else: - s_arr = stats.binned_statistic( - ell, s_arr, bins=transfer_bins - )[0] # repeat for all specs s_arr = np.tile(s_arr, len(specs)) else: s_arr = np.zeros(nbins_transfer * len(specs)) for si, spec in enumerate(specs): s_spec = cls_shape["cmb_{}".format(spec)][:lk] - if self.weighted_bins: - s_spec = stats.binned_statistic( - ell, lfac * s_spec, bins=transfer_bins - )[0] - s_spec /= lfac_binned - else: - s_spec = stats.binned_statistic( - ell, s_spec, bins=transfer_bins - )[0] s_arr[ si * nbins_transfer : (si + 1) * nbins_transfer ] = s_spec @@ -5065,7 +5030,7 @@ def bin_things(comp, d, md): self.md_fg_b2 = md_b2 self.md_fg_b3 = md_b3 """ - bin_things(comp, d, md, d_b1, d_b2, d_b3, md_b1, md_b2, md_b3) + bin_things(comp, d, md) return cbl def get_model_spectra( From 491b126f9815d183e91bdaaddf757f6f57caf817 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Sun, 22 Aug 2021 21:37:32 +0200 Subject: [PATCH 10/44] adding anafast l-by-l transfer matrix --- xfaster/xfaster_class.py | 35 ++++++++++++++++++++++++----------- 1 file changed, 24 insertions(+), 11 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 3f11f438..bbdb9d96 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -4763,10 +4763,10 @@ def bin_cl_template_tmat( if cls_shape is None: cls_shape = self.cls_shape - from spider_analysis import si - + #from spider_analysis import nsi # Need to update this to use ell-by-ell anafast transfer matrix - T = nsi.TransferMatrix() + #T = nsi.TransferMatrix() + # updated to anafast ll xfermat, no need nsi. map_pairs = None if map_tag is not None: @@ -4784,9 +4784,11 @@ def bin_cl_template_tmat( specs = list(self.specs) # This will need to change for ell-by-ell transfer matrix - transfer_bins = np.arange( - 8, 8 + 25 * 13, 25 - ) # change this once transfer matrix adds lowest bin! + #transfer_bins = np.arange( + # 8, 8 + 25 * 13, 25 + #) # change this once transfer matrix adds lowest bin! + transfer_bins = np.arange(5, 301, 25) + nbins_transfer = len(transfer_bins) - 1 lmin = min(transfer_bins) lmax = max(transfer_bins) @@ -4883,8 +4885,7 @@ def bin_things(comp, d, md): # don't create a new object in memory each time # use last one's space to save runtime self.d = np.multiply(self.d_fg, s_arr, out=getattr(self, "d", None)) - self.md = np.multiply( - self.md_fg, s_arr, out=getattr(self, "md", None) + self.md = np.multiply(self.md_fg, s_arr, out=getattr(self, "md", None)) else: for k in beam_keys: if not hasattr(self, "d"): @@ -4958,9 +4959,21 @@ def bin_things(comp, d, md): b_arr["b1"][si, xi] * b_arr["b2"][si, xi] ) - xfermat = T.get( - self.map_reobs_freqs[tag1], self.map_reobs_freqs[tag1] - ) + # updating to anafast ll transfer matrix, no need nsi + #xfermat = T.get( + # self.map_reobs_freqs[tag1], self.map_reobs_freqs[tag1] + #) + + # loading anafast ll transfer matrix + ll_xfermat_dir = "/data/vyluu/anafast_matrix_blocks" + ll_xfermat_fn = "{:s}x{:s}_{:s}_to_{:s}_block.dat".format( + self.map_reobs_freqs[tag1], + self.map_reobs_freqs[tag2], + specs[0], + specs[1] + ) + xfermat = np.loadtxt(os.path.join(ll_xfermat_dir, ll_xfermat_fn)) + xfermat = xfermat[5:301] # use correct shape spectrum if comp == "fg": From 49db9cc35fb20c2090e0ffa90016a8b11805c78e Mon Sep 17 00:00:00 2001 From: Sherry Song Date: Mon, 12 Jul 2021 15:48:44 -0400 Subject: [PATCH 11/44] Adding gcorr_config_null.ini file to script/gcorr folder. --- scripts/gcorr/gcorr_config_null.ini | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 scripts/gcorr/gcorr_config_null.ini diff --git a/scripts/gcorr/gcorr_config_null.ini b/scripts/gcorr/gcorr_config_null.ini new file mode 100644 index 00000000..6917d90d --- /dev/null +++ b/scripts/gcorr/gcorr_config_null.ini @@ -0,0 +1,18 @@ +# Options specific to gcorr calculation +[gcorr_opts] +null = true +map_tags = 90,150a +data_subset = full +output_root = /home/xsong/nulltest_testing/gcorr_run +nsim = 500 + +# Options we can directly pass to XFaster +[xfaster_opts] +config = /home/xsong/nulltest_testing/config.ini +data_root = /projects/WCJONES/spider/XF_Nulls_Sept18/port +data_root2 = /projects/WCJONES/spider/XF_Nulls_Sept18/starboard +data_type = raw +signal_type = flatBBDl +noise_type = stationary +mask_type = pointsource_latlon +verbose = debug From 1c91f219be4eb2b9eaef774f0c3aeede1eaa96d1 Mon Sep 17 00:00:00 2001 From: Sherry Song Date: Wed, 14 Jul 2021 16:06:13 -0400 Subject: [PATCH 12/44] Modifying gcorr_config_null to be consistent with the example code package. --- scripts/gcorr/gcorr_config_null.ini | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/scripts/gcorr/gcorr_config_null.ini b/scripts/gcorr/gcorr_config_null.ini index 6917d90d..74b68598 100644 --- a/scripts/gcorr/gcorr_config_null.ini +++ b/scripts/gcorr/gcorr_config_null.ini @@ -1,18 +1,17 @@ # Options specific to gcorr calculation [gcorr_opts] null = true -map_tags = 90,150a +map_tags = 90, 150a data_subset = full -output_root = /home/xsong/nulltest_testing/gcorr_run -nsim = 500 +output_root = ../../example/null/gcorr_run +nsim = 100 # Options we can directly pass to XFaster [xfaster_opts] -config = /home/xsong/nulltest_testing/config.ini -data_root = /projects/WCJONES/spider/XF_Nulls_Sept18/port -data_root2 = /projects/WCJONES/spider/XF_Nulls_Sept18/starboard -data_type = raw -signal_type = flatBBDl -noise_type = stationary -mask_type = pointsource_latlon +config = ../../example/config_example.ini +data_root = ../../example/maps_example_half1 +data_root2 = ../../example/maps_example_half2 +signal_type = synfast +noise_type = gaussian #necessary for null gcorr +mask_type = rectangle verbose = debug From 16cf745360ef335e15eb98a5a1e369d5641af2e0 Mon Sep 17 00:00:00 2001 From: Sherry Song Date: Wed, 14 Jul 2021 16:10:29 -0400 Subject: [PATCH 13/44] Modifying gcorr_config_null. --- scripts/gcorr/gcorr_config_null.ini | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/scripts/gcorr/gcorr_config_null.ini b/scripts/gcorr/gcorr_config_null.ini index 74b68598..8717799c 100644 --- a/scripts/gcorr/gcorr_config_null.ini +++ b/scripts/gcorr/gcorr_config_null.ini @@ -9,8 +9,9 @@ nsim = 100 # Options we can directly pass to XFaster [xfaster_opts] config = ../../example/config_example.ini +#Make null map set using the same set as for the signal run. data_root = ../../example/maps_example_half1 -data_root2 = ../../example/maps_example_half2 +data_root2 = ../../example/maps_example_half2 #necessary for a null run signal_type = synfast noise_type = gaussian #necessary for null gcorr mask_type = rectangle From 94a3a5d9a0fc2d3b0d23408243e66226038027ac Mon Sep 17 00:00:00 2001 From: Anne Gambrel Date: Tue, 20 Jul 2021 13:41:46 -0600 Subject: [PATCH 14/44] Allow beam/beam error products to have unused fields (#13) * Allow beam/beam error products to have unused fields * Same fix for other config fields * Keep check for transfer function keys --- xfaster/xfaster_class.py | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index bbdb9d96..d89ce3e2 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -318,14 +318,20 @@ def load_map_config(self, filename): self.fwhm = { k: np.radians(cfg.getfloat("fwhm", k) / 60.0) for k in cfg["fwhm"] } - assert tagset >= set(self.fwhm), "Unknown tags in [fwhm]" + # Only use the FWHM for fields in tagset + for ftag in list(self.fwhm.keys()): + if ftag not in tagset: + self.fwhm.pop(ftag) else: self.fwhm = {} # beam fwhm error for each tag, if not supplied in beam_error_product if "fwhm_err" in cfg: self.fwhm_err = {k: cfg.getfloat("fwhm_err", k) for k in cfg["fwhm_err"]} - assert tagset >= set(self.fwhm_err), "Unknown tags in [fwhm_err]" + # Only use the FWHM for fields in tagset + for ftag in list(self.fwhm_err.keys()): + if ftag not in tagset: + self.fwhm_err.pop(ftag) else: self.fwhm_err = {} @@ -337,8 +343,10 @@ def load_map_config(self, filename): v = os.path.join(self.config_root, v) assert os.path.exists(v), "Missing beam product file {}".format(v) self.beam_product = pt.load_compat(v) - beam_set = set(self.beam_product) - assert tagset >= beam_set, "Unknown tags in beam product" + # Only use the fields in the beam product we need + for btag in list(self.beam_product.keys()): + if btag not in tagset: + self.beam_product.pop(btag) else: self.beam_product = {} @@ -348,8 +356,10 @@ def load_map_config(self, filename): v = os.path.join(self.config_root, v) assert os.path.exists(v), "Missing beam error product file {}".format(v) self.beam_error_product = pt.load_compat(v) - beam_set = set(self.beam_error_product) - assert tagset >= beam_set, "Unknown tags in beam error product" + # Only use the fields in the beam product we need + for btag in list(self.beam_error_product.keys()): + if btag not in tagset: + self.beam_error_product.pop(btag) else: self.beam_error_product = {} else: @@ -371,6 +381,10 @@ def load_map_config(self, filename): self.fit_transfer = { k: cfg.getboolean("transfer", k) for k in cfg["transfer"] } + # Remove any fields not needed + for ftag in list(self.fit_transfer.keys()): + if ftag not in tagset: + self.fit_transfer.pop(ftag) assert tagset == set(self.fit_transfer), "Missing tags in [transfer]" else: # assume true for all tags otherwise From b8ed01292afd73cad683c50f0408ba275b30226f Mon Sep 17 00:00:00 2001 From: Anne Gambrel Date: Mon, 16 Aug 2021 07:43:22 -0600 Subject: [PATCH 15/44] Fixing bugs introduced by mistake in public code (#15) * Fixing bug where first gcorr file loaded was used for all subsequent maps * Fix bug that made cleaned bandpowers always rerun. Template alpha tags were being converted to an array to compare shape to what's on disk, and then being compared in value to a dict, causing a force_rerun. * Have like_alpha_tags default to template_alpha_tags if not set. Now don't need to set both arguments, but could still only fit for a subset of alphas in the likelihood. * Break if likelihood is run with OMP threads > 1 * Pipe through likelihood r and template specs, sub_hm_noise options * Don't set EB shape spectrum to flat for likelihood. Unintended consequence of using class properties to get whether to set EB and TB spectra to flat spectrum-- don't want that for likelihood (using r instead of a spectrum file for the shape), otherwise EB changes with r, which we don't want. * upgrade pip before doing rest of docs * install g++ * install compilers before building anything * Revert "install compilers before building anything" This reverts commit 4d70f4683b527d17d729569aa96af5209e6b1e50. Co-authored-by: Anne E. Gambrel Co-authored-by: Sasha Rahlin --- .github/workflows/documentation_pr.yaml | 5 +- xfaster/xfaster_class.py | 119 +++++++++++++++++------- xfaster/xfaster_exec.py | 39 +++++++- 3 files changed, 123 insertions(+), 40 deletions(-) diff --git a/.github/workflows/documentation_pr.yaml b/.github/workflows/documentation_pr.yaml index 1d7b0079..96b59311 100644 --- a/.github/workflows/documentation_pr.yaml +++ b/.github/workflows/documentation_pr.yaml @@ -13,6 +13,9 @@ jobs: steps: - name: Checkout github repo uses: actions/checkout@v2 + - name: Upgrade pip + run: | + pip install --upgrade pip - name: Cache pip uses: actions/cache@v2 with: @@ -28,7 +31,7 @@ jobs: # Standard drop-in approach that should work for most people. - uses: ammaraskar/sphinx-action@master with: - pre-build-command: "apt-get update -y && apt-get install -y gfortran && apt-get install -y pandoc" + pre-build-command: "apt-get update -y && apt-get install -y g++ && apt-get install -y gfortran && apt-get install -y pandoc" docs-folder: "docs/" # Create an artifact of the html output. - uses: actions/upload-artifact@v2 diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index d89ce3e2..f8efe4e1 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -898,6 +898,7 @@ def get_files( foreground_type_sim=None, template_type=None, sub_planck=False, + sub_hm_noise=False, ): """ Find all files for the given data root. The data structure is:: @@ -1012,6 +1013,9 @@ def get_files( missions so no Planck autos are used. Useful for removing expected signal residuals from null tests. Maps are expected to be in reobs_planck directory + sub_hm_noise : bool + If True, subtract average of Planck ffp10 noise crosses to debias + template-cleaned spectra Returns ------- @@ -1125,18 +1129,22 @@ def get_template_files(fs, template_type): if not os.path.exists(nfile): raise OSError("Missing hm-{} template for {}".format(hm, f)) tfiles.append(nfile) - nfiles = sorted( - glob.glob( - f.replace(fs["map_root"], tnroot).replace( - ".fits", "_*.fits" + if sub_hm_noise: + nfiles = sorted( + glob.glob( + f.replace(fs["map_root"], tnroot).replace( + ".fits", "_*.fits" + ) ) ) - ) - if not len(nfiles): - raise OSError( - "Missing hm-{} template noise for {}".format(hm, f) - ) - tnfiles.append(nfiles) + if not len(nfiles): + raise OSError( + "Missing hm-{} template noise for {}".format(hm, f) + ) + tnfiles.append(nfiles) + else: + nfiles = [] + if num_template_noise is not None: if len(nfiles) != num_template_noise: raise OSError( @@ -1163,12 +1171,13 @@ def get_template_files(fs, template_type): ), "info", ) - self.log( - "Found {} template noise files in {}".format( - fs["num_template_noise"], fs["template_noise_root"] - ), - "info", - ) + if sub_hm_noise: + self.log( + "Found {} template noise files in {}".format( + fs["num_template_noise"], fs["template_noise_root"] + ), + "info", + ) self.log("Template files: {}".format(fs["template_files"]), "debug") fields = [ @@ -1800,6 +1809,12 @@ def force_rerun_children(): if value_ref is not None: for k in [field, attr]: vref = value_ref.pop(k, "undef") + # turn it into an array if possible to compare + try: + vref = pt.dict_to_arr(vref) + except: + pass + if not vref == "undef" and np.any(v != vref): self.warn( "{}: Field {} has value {}, expected {}".format( @@ -2023,8 +2038,8 @@ def get_mask_weights(self, apply_gcorr=False, reload_gcorr=False, gcorr_file=Non iteratively solving for the correction terms. gcorr_file : str If not None, path to gcorr file. Otherwise, use file labeled - mask_map__gcorr.npy in mask directory for signal, or - mask_map__gcorr_null.npy for nulls. + mask_map__gcorr.npz in mask directory for signal, or + mask_map__gcorr_null.npz for nulls. Notes ----- @@ -2063,7 +2078,7 @@ def get_mask_weights(self, apply_gcorr=False, reload_gcorr=False, gcorr_file=Non alt_name="data_xcorr", ) - def process_gcorr(gcorr_file): + def process_gcorr(gcorr_file_in): if not hasattr(self, "gcorr"): self.gcorr = None if apply_gcorr and self.gcorr is None: @@ -2076,14 +2091,18 @@ def process_gcorr(gcorr_file): if not reload_gcorr and tag in self.gcorr: continue - if gcorr_file is None: + if gcorr_file_in is None: if self.null_run: gcorr_file = mfile.replace(".fits", "_gcorr_null.npz") else: gcorr_file = mfile.replace(".fits", "_gcorr.npz") + else: + gcorr_file = gcorr_file_in + if not os.path.exists(gcorr_file): self.warn("G correction file {} not found".format(gcorr_file)) continue + gdata = pt.load_and_parse(gcorr_file) gcorr = gdata["gcorr"] for k, g in gcorr.items(): @@ -2246,6 +2265,7 @@ def get_masked_data( template_alpha=None, sub_planck=False, sub_hm_noise=True, + temp_specs=None, ): """ Compute cross spectra of the data maps. @@ -2281,6 +2301,8 @@ def get_masked_data( sub_hm_noise : bool If True, subtract average of Planck ffp10 noise crosses to debias template-cleaned spectra. + temp_specs : list + Which spectra to use for alpha in the likelihood. Notes ----- @@ -2306,6 +2328,9 @@ def get_masked_data( null_run = self.null_run map_files2 = self.map_files2 if null_run else None + if temp_specs is None: + temp_specs = ["ee", "bb", "te", "eb", "tb"] + # ensure dictionary if template_alpha is None: template_alpha = OrderedDict() @@ -2352,7 +2377,7 @@ def apply_template(): """ cls_clean = getattr(self, "cls_data_clean", OrderedDict()) - for spec in self.specs: + for spec in set(self.specs) & set(temp_specs): cls_clean[spec] = copy.deepcopy(self.cls_data[spec]) if spec not in self.cls_template: continue @@ -3932,7 +3957,7 @@ def get_signal_shape( specs.remove("eb") if "tb" in specs: specs.remove("tb") - tbeb = "eb" in specs and "tb" in specs + 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: @@ -4035,12 +4060,16 @@ def get_signal_shape( cls_shape[spec] = np.append([0, 0], d[: ltmp - 2]) # EB and TB flat l^2 * C_l - if self.pol and 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) + 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: + 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 @@ -6924,7 +6953,9 @@ def get_likelihood( converge_criteria=0.01, reset_backend=None, file_tag=None, - use_xfer_mat=True, + sub_hm_noise=False, + r_specs=["ee", "bb"], + temp_specs=["ee", "bb", "eb"], ): """ Explore the likelihood, optionally with an MCMC sampler. @@ -6994,11 +7025,22 @@ def get_likelihood( set to True if the checkpoint has been forced to be rerun. file_tag : string If supplied, appended to the likelihood filename. - use_xfer_mat : bool - If True, use transfer matrix to construct model spectrum. If false, use - transfer function XFaster computes + sub_hm_noise : bool + If True, subtract average of Planck ffp10 noise crosses to debias + template-cleaned spectra + r_specs : list + Which spectra to use in the r likelihood. + temp_specs : list + Which spectra to use for alpha in the likelihood. """ + # Check if OMP threads > 1. It actually slows the code way down to have + # more threads. Break instead. + assert os.getenv("OMP_NUM_THREADS") in [ + None, + "1", + ], "OMP threads must be set to 1 for likelihood" + for x in [ r_prior, alpha_prior, @@ -7011,6 +7053,9 @@ def get_likelihood( if x is not None: x[:] = [float(x[0]), float(x[1])] + r_specs = [x.lower() for x in r_specs] + temp_specs = [x.lower() for x in temp_specs] + save_name = "like_mcmc" if not mcmc: alpha_prior = None @@ -7320,7 +7365,11 @@ def log_like( if alpha is None: clsi = cls_input else: - self.get_masked_data(template_alpha=OrderedDict(zip(alpha_tags, alpha))) + self.get_masked_data( + template_alpha=OrderedDict(zip(alpha_tags, alpha)), + sub_hm_noise=sub_hm_noise, + temp_specs=temp_specs, + ) clsi = self.get_data_spectra(do_noise=False) if beam is not None: @@ -7341,7 +7390,7 @@ def log_like( if r is not None: for stag, d in cls_model0.items(): comp, spec = stag.split("_", 1) - if spec not in ["ee", "bb"] or comp not in ["cmb", "total"]: + if spec not in r_specs or comp not in ["cmb", "total"]: continue ctag = "cmb_{}".format(spec) for xname, dd in d.items(): @@ -7363,7 +7412,7 @@ def log_like( elif beam is not None: for stag, d in cls_model0.items(): comp, spec = stag.split("_", 1) - if spec not in ["ee", "bb"] or comp not in ["cmb", "total"]: + if spec not in r_specs or comp not in ["cmb", "total"]: continue ctag = "cmb_{}".format(spec) for xname, dd in d.items(): diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index 96d25ffa..52908310 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -70,6 +70,8 @@ def xfaster_run( like_profiles=False, like_profile_sigma=3.0, like_profile_points=100, + like_r_specs=["EE", "BB"], + like_temp_specs=["EE", "BB", "EB"], pixwin=True, save_iters=False, verbose="notice", @@ -91,7 +93,7 @@ def xfaster_run( reload_gcorr=False, gcorr_file=None, qb_file=None, - like_alpha_tags=["95", "150"], + like_alpha_tags=None, alpha_prior=[-np.inf, np.inf], r_prior=[-np.inf, np.inf], res_prior=None, @@ -248,6 +250,10 @@ def xfaster_run( Range in units of 1sigma over which to compute profile likelihoods like_profile_points : int Number of points to sample along the likelihood profile + like_r_specs : list + Which spectra to use in the r likelihood. + like_temp_specs : list + Which spectra to use for alpha in the likelihood. pixwin : bool If True, apply pixel window functions to beam windows. save_iters : bool @@ -318,7 +324,8 @@ def xfaster_run( by the residual qb values stored in qb_file. like_alpha_tags : list of strings List of map tags from which foreground template maps should be - subtracted and fit in the likelihood. + subtracted and fit in the likelihood. If None, defaults to + template_alpha_tags. alpha_prior: list of floats Flat prior edges for allowed alpha values in the likelihood. Set to None to not fit for alpha values in the likelihood. @@ -393,6 +400,9 @@ def xfaster_run( xfc.XFasterWarning, ) + if like_alpha_tags is None: + like_alpha_tags = template_alpha_tags + if len(template_alpha_tags) != len(template_alpha): raise ValueError( "template_alpha_tags and template_alpha must be the same length" @@ -440,6 +450,7 @@ def xfaster_run( foreground_type_sim=foreground_type_sim, template_type=template_type, sub_planck=sub_planck, + sub_hm_noise=sub_hm_noise, ) config_vars.update(file_opts, "File Options") @@ -537,6 +548,8 @@ def xfaster_run( mcmc=mcmc, lmin=like_lmin, lmax=like_lmax, + r_specs=like_r_specs, + temp_specs=like_temp_specs, null_first_cmb=null_first_cmb, alpha_tags=like_alpha_tags, alpha_prior=alpha_prior, @@ -550,11 +563,13 @@ def xfaster_run( num_walkers=mcmc_walkers, converge_criteria=like_converge_criteria, file_tag=like_tag, - use_xfer_mat=use_xfer_mat, + sub_hm_noise=sub_hm_noise, ) config_vars.update(like_opts, "Likelihood Estimation Options") config_vars.remove_option("XFaster General", "like_lmin") config_vars.remove_option("XFaster General", "like_lmax") + config_vars.remove_option("XFaster General", "like_r_specs") + config_vars.remove_option("XFaster General", "like_temp_specs") config_vars.remove_option("XFaster General", "mcmc_walkers") config_vars.remove_option("XFaster General", "like_converge_criteria") config_vars.remove_option("XFaster General", "like_tag") @@ -613,7 +628,7 @@ def xfaster_run( weighted_bins=weighted_bins, ) - if template_type is not None: + if template_type is not None and sub_hm_noise: X.log("Computing template noise ensemble averages...", "notice") X.get_masked_template_noise(template_type) @@ -1097,6 +1112,20 @@ def add_arg(P, name, argtype=None, default=None, short=None, help=None, **kwargs argtype=int, help="Maximum ell to include in the likelihood calculation", ) + add_arg( + G, + "like_r_specs", + nargs="+", + help="Which spectra to use in the r likelihood calculation", + choices=["TT", "EE", "BB", "TE", "EB", "TB"], + ) + add_arg( + G, + "like_temp_specs", + nargs="+", + help="Which spectra to use in the alpha likelihood calculation", + choices=["TT", "EE", "BB", "TE", "EB", "TB"], + ) add_arg( G, "like_profile_sigma", @@ -1439,6 +1468,8 @@ def add_job(self, **kwargs): "like_beam_tags", "beam_prior", "res_prior", + "like_r_specs", + "like_temp_specs", "betad_prior", "dust_amp_prior", "dust_ellind_prior", From 76f1be77e7c59a8916c9ca2aff930e654ae8f300 Mon Sep 17 00:00:00 2001 From: Anne Gambrel Date: Mon, 16 Aug 2021 07:52:31 -0600 Subject: [PATCH 16/44] Update main github action file to match fix for PR one --- .github/workflows/documentation.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/documentation.yaml b/.github/workflows/documentation.yaml index 8d833378..28586111 100644 --- a/.github/workflows/documentation.yaml +++ b/.github/workflows/documentation.yaml @@ -28,7 +28,7 @@ jobs: # Standard drop-in approach that should work for most people. - uses: ammaraskar/sphinx-action@master with: - pre-build-command: "apt-get update -y && apt-get install -y gfortran && apt-get install -y pandoc" + pre-build-command: "apt-get update -y && apt-get install -y g++ && apt-get install -y gfortran && apt-get install -y pandoc" docs-folder: "docs/" # Create an artifact of the html output. - uses: actions/upload-artifact@v2 From fc744cc9924cb409ad7b78b795ce80bd10ddd6aa Mon Sep 17 00:00:00 2001 From: Anne Gambrel Date: Mon, 16 Aug 2021 08:16:51 -0600 Subject: [PATCH 17/44] add pip upgrade --- .github/workflows/documentation.yaml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.github/workflows/documentation.yaml b/.github/workflows/documentation.yaml index 28586111..9be14874 100644 --- a/.github/workflows/documentation.yaml +++ b/.github/workflows/documentation.yaml @@ -13,6 +13,9 @@ jobs: steps: - name: Checkout github repo uses: actions/checkout@v2 + - name: Upgrade pip + run: | + pip install --upgrade pip - name: Cache pip uses: actions/cache@v2 with: From 79372d99a6033fcdb3a9d17e583b22c9524d6cf2 Mon Sep 17 00:00:00 2001 From: Anne Gambrel Date: Mon, 16 Aug 2021 08:45:43 -0600 Subject: [PATCH 18/44] specifying versions of some requirements to hopefully reduce python package dependency errors --- docs/requirements.txt | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/requirements.txt b/docs/requirements.txt index d23826b2..55da13d0 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,6 +1,6 @@ sphinx_rtd_theme -numpy -nbsphinx -ipykernel +numpy==1.18.4 +nbsphinx==0.8.3 +ipykernel==5.5.3 pandoc -e . From 903ca3681a93e6d413817f50c1db038b45a4b559 Mon Sep 17 00:00:00 2001 From: Anne Gambrel Date: Mon, 16 Aug 2021 09:21:03 -0600 Subject: [PATCH 19/44] again trying to fix CI --- .github/workflows/documentation.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/documentation.yaml b/.github/workflows/documentation.yaml index 9be14874..57295e2b 100644 --- a/.github/workflows/documentation.yaml +++ b/.github/workflows/documentation.yaml @@ -31,7 +31,7 @@ jobs: # Standard drop-in approach that should work for most people. - uses: ammaraskar/sphinx-action@master with: - pre-build-command: "apt-get update -y && apt-get install -y g++ && apt-get install -y gfortran && apt-get install -y pandoc" + pre-build-command: "apt-get update -y && apt-get install --reinstall -y g++ && apt-get install -y gfortran && apt-get install -y pandoc" docs-folder: "docs/" # Create an artifact of the html output. - uses: actions/upload-artifact@v2 From 064eaddb60156a31c9ffd7bcbf6769e921842f12 Mon Sep 17 00:00:00 2001 From: Anne Gambrel Date: Mon, 16 Aug 2021 09:31:05 -0600 Subject: [PATCH 20/44] again trying to fix CI --- .github/workflows/documentation.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/documentation.yaml b/.github/workflows/documentation.yaml index 57295e2b..2c88bcf3 100644 --- a/.github/workflows/documentation.yaml +++ b/.github/workflows/documentation.yaml @@ -31,7 +31,7 @@ jobs: # Standard drop-in approach that should work for most people. - uses: ammaraskar/sphinx-action@master with: - pre-build-command: "apt-get update -y && apt-get install --reinstall -y g++ && apt-get install -y gfortran && apt-get install -y pandoc" + pre-build-command: "apt-get update -y && apt-get install --reinstall build-essential && apt-get install -y g++ && apt-get install -y gfortran && apt-get install -y pandoc" docs-folder: "docs/" # Create an artifact of the html output. - uses: actions/upload-artifact@v2 From e66a3d1771f75777e1badddd5022e32ab73c0469 Mon Sep 17 00:00:00 2001 From: Anne Gambrel Date: Mon, 16 Aug 2021 10:31:54 -0600 Subject: [PATCH 21/44] again trying to fix CI 3 --- .github/workflows/documentation.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/documentation.yaml b/.github/workflows/documentation.yaml index 2c88bcf3..992cc7d9 100644 --- a/.github/workflows/documentation.yaml +++ b/.github/workflows/documentation.yaml @@ -21,7 +21,7 @@ jobs: with: path: ~/.cache/pip # Look to see if there is a cache hit for the corresponding requirements file - key: ${{ runner.os }}-pip-${{ hashFiles('requirements.txt') }} + key: ${{ runner.os }}-pip-${{ hashFiles('**/requirements.txt') }} restore-keys: | ${{ runner.os }}-pip- ${{ runner.os }}- @@ -31,7 +31,7 @@ jobs: # Standard drop-in approach that should work for most people. - uses: ammaraskar/sphinx-action@master with: - pre-build-command: "apt-get update -y && apt-get install --reinstall build-essential && apt-get install -y g++ && apt-get install -y gfortran && apt-get install -y pandoc" + pre-build-command: "apt-get update -y && apt-get install --reinstall build-essential && apt-get install -y g++ && apt-get install -y gfortran && apt-get install -y pandoc && pip install --upgrade pip" docs-folder: "docs/" # Create an artifact of the html output. - uses: actions/upload-artifact@v2 From 1d18a0a6665c7f7f6ad9ed7d628c3b2dd9e90db0 Mon Sep 17 00:00:00 2001 From: Anne Gambrel Date: Mon, 16 Aug 2021 10:45:31 -0600 Subject: [PATCH 22/44] again trying to fix CI 4 --- .github/workflows/documentation.yaml | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/.github/workflows/documentation.yaml b/.github/workflows/documentation.yaml index 992cc7d9..121565da 100644 --- a/.github/workflows/documentation.yaml +++ b/.github/workflows/documentation.yaml @@ -29,9 +29,10 @@ jobs: run: | pip install -e . # Standard drop-in approach that should work for most people. - - uses: ammaraskar/sphinx-action@master + - name: Build docs + uses: ammaraskar/sphinx-action@master with: - pre-build-command: "apt-get update -y && apt-get install --reinstall build-essential && apt-get install -y g++ && apt-get install -y gfortran && apt-get install -y pandoc && pip install --upgrade pip" + pre-build-command: "apt-get update -y && apt-get install --reinstall -y build-essential && apt-get install -y g++ && apt-get install -y gfortran && apt-get install -y pandoc && pip install --upgrade pip" docs-folder: "docs/" # Create an artifact of the html output. - uses: actions/upload-artifact@v2 From c3224e728c97d9f67c0b041e9b1aac7ad0a6612b Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Wed, 25 Aug 2021 17:55:28 +0200 Subject: [PATCH 23/44] cherry-pick recent updates in main branch to fg_like --- xfaster/xfaster_class.py | 20 +++++++++++++++----- xfaster/xfaster_exec.py | 1 + 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index f8efe4e1..58e410b5 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -4879,7 +4879,7 @@ def bin_things(comp, d, md): bw = self.bin_weights[stag] for xi, (xname, (tag1, tag2)) in enumerate(map_pairs.items()): if beam_error: - cbl[mstag][xname] = OrderedDict( + cbl[stag][xname] = OrderedDict( # change from mstag [(k, np.zeros((len(bd), lmax + 1))) for k in beam_keys] ) else: @@ -5008,7 +5008,7 @@ def bin_things(comp, d, md): #) # loading anafast ll transfer matrix - ll_xfermat_dir = "/data/vyluu/anafast_matrix_blocks" + ll_xfermat_dir = "/data/vyluu/anafast_matrix_blocks" # make this into param ll_xfermat_fn = "{:s}x{:s}_{:s}_to_{:s}_block.dat".format( self.map_reobs_freqs[tag1], self.map_reobs_freqs[tag2], @@ -6795,6 +6795,7 @@ def get_bandpowers( like_profile_sigma=3.0, like_profile_points=100, file_tag=None, + use_xfer_mat=True, ): """ Compute the maximum likelihood bandpowers of the data, assuming @@ -6903,7 +6904,12 @@ def get_bandpowers( self.clear_precalc() - cbl = self.bin_cl_template(map_tag=map_tag, transfer_run=False) + if use_xfer_mat: + cbl = self.bin_cl_template_tmat(map_tag=map_tag, transfer_run=False) + self.lmax = 300 #308 + self.lmin = 5 #8 + else: + cbl = self.bin_cl_template(map_tag=map_tag, transfer_run=False) ret = self.fisher_iterate( cbl, @@ -6953,6 +6959,7 @@ def get_likelihood( converge_criteria=0.01, reset_backend=None, file_tag=None, + use_xfer_mat=True, sub_hm_noise=False, r_specs=["ee", "bb"], temp_specs=["ee", "bb", "eb"], @@ -7025,6 +7032,9 @@ def get_likelihood( set to True if the checkpoint has been forced to be rerun. file_tag : string If supplied, appended to the likelihood filename. + use_xfer_mat : bool + If True, use transfer matrix to construct model spectrum. If false, use + transfer function XFaster computes sub_hm_noise : bool If True, subtract average of Planck ffp10 noise crosses to debias template-cleaned spectra @@ -7148,8 +7158,8 @@ def get_likelihood( if use_xfer_mat: bin_cl_template = self.bin_cl_template_tmat - self.lmax = 308 - self.lmin = 8 + self.lmax = 300 #308 + self.lmin = 5 #8 else: bin_cl_template = self.bin_cl_template diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index 52908310..905a261c 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -563,6 +563,7 @@ def xfaster_run( num_walkers=mcmc_walkers, converge_criteria=like_converge_criteria, file_tag=like_tag, + use_xfer_mat=use_xfer_mat, sub_hm_noise=sub_hm_noise, ) config_vars.update(like_opts, "Likelihood Estimation Options") From 520a42b73ba37f1dad232acfb875de42936dcd4b Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 1 Sep 2021 01:58:53 +1200 Subject: [PATCH 24/44] merge detritus --- docs/requirements.txt | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/requirements.txt b/docs/requirements.txt index 55da13d0..d23826b2 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,6 +1,6 @@ sphinx_rtd_theme -numpy==1.18.4 -nbsphinx==0.8.3 -ipykernel==5.5.3 +numpy +nbsphinx +ipykernel pandoc -e . From b35cb1c32db75a00b738ad5198bf6ac6dd260f10 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Thu, 2 Sep 2021 14:05:24 +1200 Subject: [PATCH 25/44] Rewire transfer matrix handling * ``transfer_matrix_root`` argument where transfer matrix data are stored. If set, short circuit the transfer function calculation and use the matrix to build Cbl's for both bandpowers and likelihoods. * ``get_transfer_matrix`` function currently loads matrix blocks from separate text files into a dictionary, and writes to an output checkpoint npz file. * check the ``transfer_matrix`` attribute in all the necessary functions, rather than carrying around an extra keyword argument * transfer function and beams are excluded from cached Mll' matrix if the transfer function has been loaded * transfer matrix is applied to the input cls_shape before multiplying by Mll' to construct Cbl's --- xfaster/xfaster_class.py | 496 ++++++++++----------------------------- xfaster/xfaster_exec.py | 37 +-- 2 files changed, 144 insertions(+), 389 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index cc1875bd..3526556f 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -134,6 +134,7 @@ class XFaster(object): "sims_transfer": ["transfer"], "shape_transfer": ["transfer"], "transfer": ["bandpowers"], + "transfer_matrix": ["bandpowers"], "sims": ["bandpowers"], "beams": ["transfer"], "data": ["bandpowers"], @@ -1395,7 +1396,8 @@ def get_filename( template_cleaned, reference_subtracted. bp_opts : bool If True, also check the following attributes (in addition to those - checked if ``data_opts`` is True): weighted_bins, return_cls. + checked if ``data_opts`` is True): weighted_bins, return_cls, + transfer_matrix. Returns ------- @@ -1435,6 +1437,8 @@ def get_filename( name += ["wbins"] if getattr(self, "return_cls", False): name += ["cl"] + if hasattr(self, "transfer_matrix"): + name += ["xfermat"] if map_tag is not None: name += ["map", map_tag] @@ -3830,30 +3834,30 @@ def kernel_precalc(self, 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: - transfer[spec][tag] = self.transfer[spec][tag] - lk = slice(0, lmax + 1) mll = OrderedDict() + use_transfer_matrix = hasattr(self, "transfer_matrix") + for spec in specs: mll[spec] = OrderedDict() if spec in ["ee", "bb"]: mspec = "{}_mix".format(spec) mll[mspec] = OrderedDict() - for xname, (m0, m1) in map_pairs.items(): - # beams - fb2 = self.beam_windows[spec][m0][lk] * self.beam_windows[spec][m1][lk] + # assume that beam and filtering transfer effects are included in + # the transfer matrix + if not use_transfer_matrix: + bw = self.beam_windows[spec] + tf = self.transfer[spec] - # transfer function - if not transfer_run: - fb2 *= np.sqrt(transfer[spec][m0][lk] * transfer[spec][m1][lk]) + for xname, (m0, m1) in map_pairs.items(): + if not use_transfer_matrix: + # beams + fb2 = bw[m0][lk] * bw[m1][lk] + # transfer function + if not transfer_run: + fb2 *= np.sqrt(tf[m0][lk] * tf[m1][lk]) # kernels if spec == "tt": @@ -3867,9 +3871,9 @@ def kernel_precalc(self, map_tag=None, transfer_run=False): k = self.pkern[xname][:, lk] - self.mkern[xname][:, lk] # store final product - mll[spec][xname] = k * fb2 + mll[spec][xname] = k if use_transfer_matrix else k * fb2 if spec in ["ee", "bb"]: - mll[mspec][xname] = mk * fb2 + mll[mspec][xname] = mk if use_transfer_matrix else mk * fb2 return mll @@ -4038,12 +4042,26 @@ def binup2(d, bd, bw): continue - # use correct shape spectrum - if comp == "fg": - # single foreground spectrum - s_arr[xi] = cls_shape["fg"][lk] + if hasattr(self, "transfer_matrix"): + # construct modified shape spectrum C^X_l = sum_{l',Y} J^XY_ll' * C^Y_l' + for spec2 in specs: + # use correct shape spectrum + if comp == "fg": + cl1 = cls_shape["fg"][lk] + else: + cl1 = cls_shape["cmb_{}".format(spec2)][lk] + + xfermat = self.transfer_matrix[xname][spec][spec2][lk, lk] + + # matrix multiplication over l', loop sum over Y + s_arr[xi] += np.einsum("ij,j->i", xfermat, cl1) else: - s_arr[xi] = cls_shape["cmb_{}".format(spec)][lk] + # use correct shape spectrum + if comp == "fg": + # single foreground spectrum + s_arr[xi] = cls_shape["fg"][lk] + else: + s_arr[xi] = cls_shape["cmb_{}".format(spec)][lk] # use correct beam error shape if beam_error: @@ -4096,331 +4114,6 @@ def binup2(d, bd, bw): return cbl - def bin_cl_template_tmat( - self, - cls_shape=None, - map_tag=None, - beam_error=False, - fg_ell_ind=0, - ): - """ - *** Note: obtained from labah fg_like - Compute the Cbl matrix from the input shape spectrum. - Use transfer matrix instead of transfer function - - This method requires kernels to have been precomputed. - - Arguments - --------- - cls_shape : array_like - The shape spectrum to use. This can be computed using - `get_signal_shape` or otherwise. - map_tag : str - 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. - 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. - fg_ell_ind : float - If binning foreground shape, offset the ell index from the reference - by this amount. - - Returns - ------- - cbl : dict of arrays (num_bins, 2, lmax + 1) - The Cbl matrix, indexed by component and spectrum, then by map - cross, e.g. `cbl['cmb_tt']['map1:map2']` - E/B mixing terms are stored in elements ``cbl['cmb_ee_mix']`` and - ``cbl['cmb_bb_mix']``, and unmixed terms are stored in elements - ``cbl['cmb_ee']`` and ``cbl['cmb_bb']``. - """ - if cls_shape is None: - cls_shape = self.cls_shape - - #from spider_analysis import nsi - # Need to update this to use ell-by-ell anafast transfer matrix - #T = nsi.TransferMatrix() - # updated to anafast ll xfermat, no need nsi. - - map_pairs = None - if map_tag is not None: - if map_tag in self.map_pairs: - map_pairs = {map_tag: self.map_pairs[map_tag]} - map_tags = list(set(self.map_pairs[map_tag])) - else: - map_tags = [map_tag] - else: - map_tags = self.map_tags - - if map_pairs is None: - map_pairs = pt.tag_pairs(map_tags) - - specs = list(self.specs) - - # This will need to change for ell-by-ell transfer matrix - #transfer_bins = np.arange( - # 8, 8 + 25 * 13, 25 - #) # change this once transfer matrix adds lowest bin! - transfer_bins = np.arange(5, 301, 25) - - nbins_transfer = len(transfer_bins) - 1 - lmin = min(transfer_bins) - lmax = max(transfer_bins) - lmax_kern = lmax - - if beam_error: - beam_error = self.get_beam_errors() - beam_keys = ["b1", "b2", "b3"] - - ls = slice(2, lmax + 1) - lk = slice(0, lmax_kern + 1) - cbl = OrderedDict() - - comps = [] - if "cmb_tt" in cls_shape or "cmb_ee" in cls_shape: - comps += ["cmb"] - if "fg" 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_noise0 = self.cls_noise0_null if self.null_run else self.cls_noise0 - cls_noise1 = self.cls_noise1_null if self.null_run else self.cls_noise1 - cls_sxn0 = self.cls_sxn0_null if self.null_run else self.cls_sxn0 - cls_sxn1 = self.cls_sxn1_null if self.null_run else self.cls_sxn1 - cls_nxs0 = self.cls_nxs0_null if self.null_run else self.cls_nxs0 - cls_nxs1 = self.cls_nxs1_null if self.null_run else self.cls_nxs1 - - ell = np.arange(lmax_kern + 1) - - def binup(d, left, right, weights): - return (d[..., left:right] * weights).sum(axis=-1) - - def bin_things(comp, d, md): - if "res" in comp: - return - for si, spec in enumerate(specs): - stag = "{}_{}".format(comp, spec) - cbl.setdefault(stag, OrderedDict()) - mstag = None - if spec in ["ee", "bb"]: - mstag = stag + "_mix" - cbl.setdefault(mstag, OrderedDict()) - bd = self.bin_def[stag] - bw = self.bin_weights[stag] - for xi, (xname, (tag1, tag2)) in enumerate(map_pairs.items()): - if beam_error: - cbl[stag][xname] = OrderedDict( # change from mstag - [(k, np.zeros((len(bd), lmax + 1))) for k in beam_keys] - ) - else: - cbl[stag][xname] = np.zeros((len(bd), lmax + 1)) - - if spec in ["ee", "bb"]: - if beam_error: - cbl[mstag][xname] = OrderedDict( - [(k, np.zeros((len(bd), lmax + 1))) for k in beam_keys] - ) - else: - cbl[mstag][xname] = np.zeros((len(bd), lmax + 1)) - - # integrate per bin - for idx, ((left, right), weights) in enumerate(zip(bd, bw)): - if beam_error: - for k in beam_keys: - cbl[stag][xname][k][idx, ls] = binup( - d[k][si, xi], left, right, weights - ) - else: - cbl[stag][xname][idx, ls] = binup( - d[si, xi], left, right, weights - ) - if spec in ["ee", "bb"]: - if beam_error: - for k in beam_keys: - cbl[mstag][xname][k][idx, ls] = binup( - md[k][si - 1, xi], left, right, weights - ) - else: - cbl[mstag][xname][idx, ls] = binup( - md[si - 1, xi], left, right, weights - ) - - for comp in comps: - # convert to matrices to do multiplication to speed things up, - # except for res is weird so don't do it for that. - # need n_xname x n_spec x ell - nspec = len(specs) - nxmap = len(map_pairs) - if comp == "fg" and fg_ell_ind != 0: - s_arr = (ell / 80.0) ** fg_ell_ind - s_arr[0] = 0 - if not beam_error: - # don't create a new object in memory each time - # use last one's space to save runtime - self.d = np.multiply(self.d_fg, s_arr, out=getattr(self, "d", None)) - self.md = np.multiply(self.md_fg, s_arr, out=getattr(self, "md", None)) - else: - for k in beam_keys: - if not hasattr(self, "d"): - self.d = OrderedDict([(k, None) for k in beam_keys]) - self.md = OrderedDict([(k, None) for k in beam_keys]) - self.d[k] = np.multiply(self.d_fg[k], s_arr, out=self.d[k]) - self.md[k] = np.multiply(self.md_fg[k], s_arr, out=self.md[k]) - bin_things(comp, self.d, self.md) - - else: - kshape = [nspec, nxmap, self.lmax - 1, lmax_kern + 1] - mkshape = [2] + kshape[1:] - k_arr = np.zeros(kshape) - mk_arr = np.zeros(mkshape) - - # because transfer matrix includes flbl2, combine those, just store - # xfermat * binned shape - - fbs_arr = np.zeros([nspec, nxmap, lmax_kern + 1]) - - for xi, (xname, (tag1, tag2)) in enumerate(map_pairs.items()): - for si, spec in enumerate(specs): - stag = "{}_{}".format(comp, spec) - mstag = None - if comp != "res" and spec in ["ee", "bb"]: - mstag = stag + "_mix" - - if "res" in comp: - s0, s1 = spec - res_tags = { - "s0m0": "res_{}_{}".format(s0 * 2, tag1), - "s0m1": "res_{}_{}".format(s0 * 2, tag2), - "s1m0": "res_{}_{}".format(s1 * 2, tag1), - "s1m1": "res_{}_{}".format(s1 * 2, tag2), - } - bd = [[0, lmax + 1]] - # if any component of XY spec is in residual bin - # def, use that bin def - for k, v in res_tags.items(): - spec0 = v.split("_")[1] - if v not in self.bin_def: - if spec0 in ["ee", "bb"]: - v = v.replace(spec0, "eebb") - if v in self.bin_def: - bd = self.bin_def[v] - else: - bd = self.bin_def[v] - for comp, cls in [ - ("res0_nxn", cls_noise0), - ("res1_nxn", cls_noise1), - ("res0_sxn", cls_sxn0), - ("res1_sxn", cls_sxn1), - ("res0_nxs", cls_nxs0), - ("res1_nxs", cls_nxs1), - ("res", cls_noise), - ]: - stag = "{}_{}".format(comp, spec) - cbl.setdefault(stag, OrderedDict()) - cbl[stag][xname] = np.zeros((len(bd), lmax + 1)) - cl1 = cls[spec][xname] - for idx, (left, right) in enumerate(bd): - lls = slice(left, right) - cbl[stag][xname][idx, lls] = np.copy(cl1[lls]) - - continue - - if beam_error: - b_arr["b1"][si, xi] = beam_error[spec][tag1] - b_arr["b2"][si, xi] = beam_error[spec][tag2] - b_arr["b3"][si, xi] = ( - b_arr["b1"][si, xi] * b_arr["b2"][si, xi] - ) - - # updating to anafast ll transfer matrix, no need nsi - #xfermat = T.get( - # self.map_reobs_freqs[tag1], self.map_reobs_freqs[tag1] - #) - - # loading anafast ll transfer matrix - ll_xfermat_dir = "/data/vyluu/anafast_matrix_blocks" # make this into param - ll_xfermat_fn = "{:s}x{:s}_{:s}_to_{:s}_block.dat".format( - self.map_reobs_freqs[tag1], - self.map_reobs_freqs[tag2], - specs[0], - specs[1] - ) - xfermat = np.loadtxt(os.path.join(ll_xfermat_dir, ll_xfermat_fn)) - xfermat = xfermat[5:301] - - # use correct shape spectrum - if comp == "fg": - # single foreground spectrum - s_arr = cls_shape["fg"][lk] * (ell / 80.0) ** fg_ell_ind - s_arr[0] = 0 - # repeat for all specs - s_arr = np.tile(s_arr, len(specs)) - else: - s_arr = np.zeros(nbins_transfer * len(specs)) - for si, spec in enumerate(specs): - s_spec = cls_shape["cmb_{}".format(spec)][:lk] - s_arr[ - si * nbins_transfer : (si + 1) * nbins_transfer - ] = s_spec - binned_fbs = np.matmul(np.linalg.inv(xfermat), s_arr) - # fbs_arr = np.zeros([nspec, nxmap, lmax_kern + 1]) - for si, spec in enumerate(specs): - for ib, b0 in enumerate(transfer_bins[:-1]): - fbs_arr[si][xi][b0 : transfer_bins[ib + 1]] = binned_fbs[ib] - - # get cross spectrum kernel terms - if spec == "tt": - k_arr[si, xi] = self.kern[xname][ls, : lmax_kern + 1] - elif spec in ["ee", "bb"]: - k_arr[si, xi] = self.pkern[xname][ls, : lmax_kern + 1] - mk_arr[si - 1, xi] = self.mkern[xname][ls, : lmax_kern + 1] - elif spec in ["te", "tb"]: - k_arr[si, xi] = self.xkern[xname][ls, : lmax_kern + 1] - elif spec == "eb": - k_arr[si, xi] = ( - self.pkern[xname][ls] - self.mkern[xname][ls] - )[:, : lmax_kern + 1] - - # need last 3 dims of kernel to match other arrays - k_arr = np.transpose(k_arr, axes=[2, 0, 1, 3]) - mk_arr = np.transpose(mk_arr, axes=[2, 0, 1, 3]) - if not beam_error: - d = k_arr * fbs_arr - md = mk_arr * fbs_arr[1:3] - if comp == "fg": - self.d_fg = np.copy(d) - self.md_fg = np.copy(md) - d_b1 = None - d_b2 = None - d_b3 = None - md_b1 = None - md_b2 = None - md_b3 = None - else: - # why not include error here - - d = None - md = None - """ - d_b1 = k_arr * b1_arr * f_arr * s_arr - d_b2 = k_arr * b2_arr * f_arr * s_arr - d_b3 = k_arr * b3_arr * f_arr * s_arr - md_b1 = mk_arr * b1_arr[1:3] * f_arr[1:3] * s_arr_md - md_b2 = mk_arr * b2_arr[1:3] * f_arr[1:3] * s_arr_md - md_b3 = mk_arr * b3_arr[1:3] * f_arr[1:3] * s_arr_md - if comp == "fg": - self.d_fg_b1 = d_b1 - self.d_fg_b2 = d_b2 - self.d_fg_b3 = d_b3 - self.md_fg_b1 = md_b1 - self.md_fg_b2 = md_b2 - self.md_fg_b3 = md_b3 - """ - bin_things(comp, d, md) - return cbl - def get_model_spectra( self, qb, cbl, delta=True, res=True, cls_noise=None, cond_noise=None ): @@ -6113,6 +5806,81 @@ def expand_transfer(qb_transfer, bin_def): return self.transfer + def get_transfer_matrix(self, file_root=None): + """ + Read in an ell-by-ell transfer matrix, constructed from ratios of + `anafast` spectra using an ensemble of simulations. This matrix should + be used in place of the standard XFaster transfer function calculation, + and includes both filtering and beam effects. + + Arguments + --------- + file_root : str + Directory where the transfer matrix blocks are stored. Expects + column text files with naming pattern: + ``x__to__block.dat`` + Required if the transfer matrix has not already been computed. + + Returns + ------- + transfer_matrix : OrderedDict + Dictionary of matrix blocks of shape (lmax + 1, lmax + 1), + keyed by map cross (xname), output spectrum, and input spectrum. + This dictionary is also cached in the ``transfer_matrix`` attribute. + + Notes + ----- + This function is called at the ``transfer_matrix`` checkpoint, and is + meant to be an alternative to the standard ``transfer`` sequence. + + This option uses a significant amount of RAM, especially for large + numbers of map crosses, since the entire matrix is stored in memory + for efficient computation. + """ + + tm_shape = ( + self.num_corr * (self.num_spec ** 2) * (self.lmax + 1), + self.lmax + 1, + ) + + save_name = "transfer_matrix" + save_attrs = ["transfer_file_root", "transfer_matrix"] + ret = self.load_data( + save_name, + "transfer_matrix", + fields=save_attrs, + to_attrs=True, + shape=tm_shape, + shape_ref="transfer_matrix", + ) + + if ret is not None: + return ret["transfer_matrix"] + + lk = slice(0, self.lmax + 1) + + transfer_matrix = OrderedDict() + for xname, (tag1, tag2) in self.map_pairs.items(): + dx = transfer_matrix.setdefault(xname, OrderedDict()) + + for spec_out in self.specs: + dsi = dx.setdefault(spec_out, OrderedDict()) + + for spec_in in self.specs: + fname = os.path.join( + file_root, + "{}x{}_{}_to_{}_block.dat".format( + tag1, tag2, spec_in, spec_out + ), + ) + dsi[spec_out] = np.loadtxt(fname)[lk, lk] + + self.transfer_file_root = file_root + self.transfer_matrix = transfer_matrix + + self.save_data(save_name, from_attrs=save_attrs) + return transfer_matrix + def get_bandpowers( self, map_tag=None, @@ -6131,7 +5899,6 @@ def get_bandpowers( like_profile_sigma=3.0, like_profile_points=100, file_tag=None, - use_xfer_mat=True, ): """ Compute the maximum likelihood bandpowers of the data, assuming @@ -6255,12 +6022,7 @@ def get_bandpowers( self.clear_precalc() - if use_xfer_mat: - cbl = self.bin_cl_template_tmat(map_tag=map_tag, transfer_run=False) - self.lmax = 300 #308 - self.lmin = 5 #8 - else: - cbl = self.bin_cl_template(map_tag=map_tag, transfer_run=False) + cbl = self.bin_cl_template(map_tag=map_tag, transfer_run=False) ret = self.fisher_iterate( cbl, @@ -6310,7 +6072,6 @@ def get_likelihood( subtract_template_noise=False, r_specs=["ee", "bb"], template_specs=["ee", "bb", "eb"], - use_xfer_mat=False, ): """ Explore the likelihood, optionally with an MCMC sampler. @@ -6385,9 +6146,6 @@ def get_likelihood( Which spectra to use in the r likelihood. template_specs : list Which spectra to use for alpha in the likelihood. - use_xfer_mat : bool - If True, use transfer matrix to construct model spectrum. If false, use - transfer function XFaster computes """ # Check if OMP threads > 1. It actually slows the code way down to have @@ -6462,9 +6220,6 @@ def get_likelihood( "res_prior": res_prior, "beam_prior": beam_prior, } - - self.log("Priors: {}".format(priors), "debug") - # priors on quantities that affect Dmat_obs or gmat (precalculated) obs_priors = [alpha_prior] @@ -6484,16 +6239,8 @@ def get_likelihood( weighted_bins=self.weighted_bins, lmin=lmin, lmax=lmax, - use_xfer_mat=use_xfer_mat, ) - if use_xfer_mat: - bin_cl_template = self.bin_cl_template_tmat - self.lmax = 300 #308 - self.lmin = 5 #8 - else: - bin_cl_template = self.bin_cl_template - if mcmc and reset_backend is None: ret = self.load_data( save_name, @@ -6555,7 +6302,7 @@ def get_likelihood( if r_prior is None: self.log("Computing model spectrum", "debug") self.warn("Beam variation not implemented for case of no r fit") - cbl = bin_cl_template(map_tag=map_tag) + cbl = self.bin_cl_template(map_tag=map_tag) cls_model = self.get_model_spectra(qb, cbl, delta=True, cls_noise=cls_noise) else: qb = copy.deepcopy(qb) @@ -6576,23 +6323,23 @@ def get_likelihood( ) # load tensor and scalar terms separately - cbl_scalar = bin_cl_template(cls_shape_scalar, map_tag) + 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 ) - cbl_tensor = bin_cl_template(cls_shape_tensor, map_tag) + 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 ) if beam_prior is not None: # load beam error term for tensor and scalar - cbl_scalar_beam = bin_cl_template( + 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 ) - cbl_tensor_beam = bin_cl_template( + cbl_tensor_beam = self.bin_cl_template( cls_shape_tensor, map_tag, beam_error=True ) cls_mod_tens_beam = self.get_model_spectra( @@ -6732,6 +6479,7 @@ 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 1deb01ed..2d4a8010 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -98,6 +98,7 @@ def xfaster_run( signal_spec=None, signal_transfer_spec=None, model_r=None, + transfer_matrix_root=None, ref_freq=359.7, beta_ref=1.54, delta_beta_prior=0.5, @@ -120,7 +121,6 @@ def xfaster_run( res_prior=None, like_beam_tags="all", beam_prior=None, - use_xfer_mat=False, ): """ Main function for running the XFaster algorithm. @@ -358,6 +358,10 @@ def xfaster_run( model_r : float The ``r`` value to use to compute a spectrum for estimating bandpowers. Overrides ``signal_spec``. + transfer_matrix_root : str + If set, use the pre-computed ell-by-ell transfer matrix stored in this + directory, instead of the standard transfer function computed by the + XFaster estimator. ref_freq : float In GHz, reference frequency for dust model. Dust bandpowers output will be at this reference frequency. @@ -418,9 +422,6 @@ def xfaster_run( Gaussian prior mean and number of strandard deviations for beam error. This Gaussian is applied as a prior in fitting for beam error in the likelihood. Set to None to not fit for beam error. - use_xfer_mat : bool - If True, use NSI transfer matrix instead of transfer function computed - by XFaster. """ from . import __version__ as version @@ -584,6 +585,7 @@ def xfaster_run( signal_spec=signal_spec, signal_transfer_spec=signal_transfer_spec, model_r=model_r, + transfer_matrix_root=transfer_matrix_root, ref_freq=ref_freq, beta_ref=beta_ref, delta_beta_prior=delta_beta_prior, @@ -605,6 +607,7 @@ def xfaster_run( spec_opts.pop("signal_transfer_spec") spec_opts.pop("model_r") spec_opts.pop("qb_file") + spec_opts.pop("transfer_matrix_root") bandpwr_opts = spec_opts.copy() bandpwr_opts.pop("fix_bb_transfer") spec_opts.pop("file_tag") @@ -639,7 +642,6 @@ def xfaster_run( converge_criteria=like_converge_criteria, file_tag=like_tag, subtract_template_noise=subtract_template_noise, - use_xfer_mat=use_xfer_mat, ) config_vars.update(like_opts, "Likelihood Estimation Options") config_vars.remove_option("XFaster General", "like_lmin") @@ -667,19 +669,24 @@ 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") - # 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.log("Computing beam window functions...", "notice") X.get_beams(pixwin=pixwin) - X.log("Loading spectrum shape for transfer function...", "notice") - X.get_signal_shape(filename=signal_transfer_spec, transfer=True) + if transfer_matrix_root: + X.log("Loading pre-computed transfer matrix...", "notice") + X.get_transfer_matrix(file_root=transfer_matrix_root) + + else: + 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.log("Loading spectrum shape for transfer function...", "notice") + X.get_signal_shape(filename=signal_transfer_spec, transfer=True) - X.log("Computing transfer functions...", "notice") - X.get_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) @@ -1062,6 +1069,7 @@ 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, "transfer_matrix_root") add_arg(G, "ref_freq") add_arg(G, "beta_ref", argtype=float) add_arg(G, "delta_beta_prior", argtype=float) @@ -1097,7 +1105,6 @@ def add_arg(P, name, argtype=None, default=None, short=None, help=None, **kwargs add_arg(G, "res_prior", nargs=2) add_arg(G, "like_beam_tags", nargs="+", metavar="TAG") add_arg(G, "beam_prior", nargs=2) - add_arg(G, "use_xfer_mat") # submit args if mode == "submit": From 197329f14bc49f051bc5268822993bafda7e5001 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Thu, 23 Sep 2021 02:39:42 +1200 Subject: [PATCH 26/44] 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 3526556f..d0d9a88d 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -5873,7 +5873,7 @@ def get_transfer_matrix(self, file_root=None): tag1, tag2, spec_in, spec_out ), ) - dsi[spec_out] = np.loadtxt(fname)[lk, lk] + dsi[spec_in] = np.loadtxt(fname)[lk, lk] self.transfer_file_root = file_root self.transfer_matrix = transfer_matrix From 67b076b3d0bc2411e5c03b617f2611c56286dfc8 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Wed, 13 Oct 2021 01:50:48 +0200 Subject: [PATCH 27/44] Debugging transfer_matrix branch --- xfaster/xfaster_class.py | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 3526556f..014f5a89 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -116,6 +116,7 @@ class XFaster(object): "sims_transfer", "shape_transfer", "transfer", + "transfer_matrix", "sims", "beams", "data", @@ -4052,7 +4053,7 @@ def binup2(d, bd, bw): cl1 = cls_shape["cmb_{}".format(spec2)][lk] xfermat = self.transfer_matrix[xname][spec][spec2][lk, lk] - + print(xfermat.shape, cl1.shape) # matrix multiplication over l', loop sum over Y s_arr[xi] += np.einsum("ij,j->i", xfermat, cl1) else: @@ -5867,13 +5868,21 @@ def get_transfer_matrix(self, file_root=None): dsi = dx.setdefault(spec_out, OrderedDict()) for spec_in in self.specs: - fname = os.path.join( - file_root, - "{}x{}_{}_to_{}_block.dat".format( - tag1, tag2, spec_in, spec_out - ), - ) - dsi[spec_out] = np.loadtxt(fname)[lk, lk] + if spec_in[0]==spec_in[1] and spec_in==spec_out: + fname = os.path.join( + file_root, + "{}x{}_{}_to_{}_block.dat".format( + tag1, tag2, spec_in.upper(), spec_out.upper() + ), + ) + if not os.path.isfile(fname): + fname = os.path.join( + file_root, + "{}x{}_{}_to_{}_block.dat".format( + tag2, tag1, spec_in.upper(), spec_out.upper() + ), + ) + dsi[spec_out] = np.loadtxt(fname)[lk, lk] self.transfer_file_root = file_root self.transfer_matrix = transfer_matrix From 61e94cf0597bc745ed3afbf098210d341dee953a Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Thu, 14 Oct 2021 17:42:31 +0200 Subject: [PATCH 28/44] modify the get_transfer_matrix function to approximate cross1 -> cross1 transfer matrix and provide zero arrays for auto1/cross1 -> auto2/cross2 when files don't exist --- xfaster/xfaster_class.py | 72 ++++++++++++++++++++++++++++++---------- 1 file changed, 55 insertions(+), 17 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 014f5a89..092a1687 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -1926,7 +1926,7 @@ def process_gcorr(gcorr_file_in): bd = gdata["bin_def"]["cmb_{}".format(k)] if len(bd0) < len(bd): bd = bd[: len(bd0)] - gcorr[k] = g[: len(bd)] + gcorr[k] = g[: len(bd)] if not np.all(bd0 == bd): self.warn( "G correction for map {} has incompatible bin def".format( @@ -4045,6 +4045,7 @@ def binup2(d, bd, bw): if hasattr(self, "transfer_matrix"): # construct modified shape spectrum C^X_l = sum_{l',Y} J^XY_ll' * C^Y_l' + print("\nbin_cl_template: ", specs) for spec2 in specs: # use correct shape spectrum if comp == "fg": @@ -5860,29 +5861,66 @@ def get_transfer_matrix(self, file_root=None): lk = slice(0, self.lmax + 1) + def load_matrix(spec_in_in, spec_out_in): + # there exist tag 90x150a but not 150ax90 + # TODO: better way to make this isensitive to upper/lower case + fname = os.path.join( + file_root, + "{}x{}_{}_to_{}_block.dat".format( + tag1, tag2, spec_in_in.upper(), spec_out_in.upper() + ), + ) + if os.path.isfile(fname): + return np.loadtxt(fname) + else: + fname = os.path.join( + file_root, + "{}x{}_{}_to_{}_block.dat".format( + tag2, tag1, spec_in_in.upper(), spec_out_in.upper() + ), + ) + if os.path.isfile(fname): + return np.loadtxt(fname) + else: + return None + + def get_cross_matrix(cross_spec): + s1, s2 = cross_spec + s1 = s1 * 2 + s2 = s2 * 2 + return np.sqrt(np.abs(np.multiply(load_matrix(s1, s1), + load_matrix(s2, s2)))) + transfer_matrix = OrderedDict() for xname, (tag1, tag2) in self.map_pairs.items(): dx = transfer_matrix.setdefault(xname, OrderedDict()) - + + #print("\nget transfer matrix: ", self.specs) for spec_out in self.specs: dsi = dx.setdefault(spec_out, OrderedDict()) for spec_in in self.specs: - if spec_in[0]==spec_in[1] and spec_in==spec_out: - fname = os.path.join( - file_root, - "{}x{}_{}_to_{}_block.dat".format( - tag1, tag2, spec_in.upper(), spec_out.upper() - ), - ) - if not os.path.isfile(fname): - fname = os.path.join( - file_root, - "{}x{}_{}_to_{}_block.dat".format( - tag2, tag1, spec_in.upper(), spec_out.upper() - ), - ) - dsi[spec_out] = np.loadtxt(fname)[lk, lk] + #print("\nspecs: ", xname, tag1, tag2, spec_out, spec_in) + # case auto1- => auto2- and + # auto1- => auto1- + tmp = load_matrix(spec_in, spec_out) + if tmp is None: + # cross- => auto- and + # auto1- => auto2- whose files dont exist + # TODO: this is floppy. fix this. + matrix_shape = dx[self.specs[0]][self.specs[0]].shape + if spec_out[0]==spec_out[1]: + tmp = np.zeros(matrix_shape) + else: + # cross1- => cross1- + tmp = get_cross_matrix(spec_out) + # auto- => cross- and + # cross1- => cross2- + if spec_in != spec_out: + tmp = np.zeros(matrix_shape) + + dsi[spec_in] = tmp[lk, lk] + del tmp self.transfer_file_root = file_root self.transfer_matrix = transfer_matrix From e73b8ef5ed19f6019daaa4b6b78bb6a667a40340 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Thu, 14 Oct 2021 19:43:22 +0200 Subject: [PATCH 29/44] comment out qb_transfer update in fisher_iterate function --- 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 092a1687..8fab43ce 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -5462,7 +5462,7 @@ def fisher_iterate( if not transfer_run: out.update( - qb_transfer=self.qb_transfer, + #qb_transfer=self.qb_transfer, transfer=self.transfer, ) if self.template_cleaned: From a86640c9b42d936aa684d9e987fced379bb57073 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Fri, 15 Oct 2021 23:12:04 +0200 Subject: [PATCH 30/44] comment out transfer and qb_transfer --- xfaster/xfaster_class.py | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index ada9595e..1725acff 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -4046,7 +4046,6 @@ def binup2(d, bd, bw): if hasattr(self, "transfer_matrix"): # construct modified shape spectrum C^X_l = sum_{l',Y} J^XY_ll' * C^Y_l' - print("\nbin_cl_template: ", specs) for spec2 in specs: # use correct shape spectrum if comp == "fg": @@ -5462,10 +5461,10 @@ def fisher_iterate( ) if not transfer_run: - out.update( - #qb_transfer=self.qb_transfer, - transfer=self.transfer, - ) + # out.update( + # qb_transfer=self.qb_transfer, + # transfer=self.transfer, + # ) if self.template_cleaned: out.update(template_alpha=self.template_alpha) @@ -5901,7 +5900,6 @@ def get_cross_matrix(cross_spec): dsi = dx.setdefault(spec_out, OrderedDict()) for spec_in in self.specs: -<<<<<<< HEAD #print("\nspecs: ", xname, tag1, tag2, spec_out, spec_in) # case auto1- => auto2- and # auto1- => auto1- @@ -5923,15 +5921,6 @@ def get_cross_matrix(cross_spec): dsi[spec_in] = tmp[lk, lk] del tmp -======= - fname = os.path.join( - file_root, - "{}x{}_{}_to_{}_block.dat".format( - tag1, tag2, spec_in, spec_out - ), - ) - dsi[spec_in] = np.loadtxt(fname)[lk, lk] ->>>>>>> 2739e163a6c0571b11971cd319f36c88f2adeab9 self.transfer_file_root = file_root self.transfer_matrix = transfer_matrix From efabd96aa40dbb67ce7e451b56dd61ab2df845a6 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Mon, 18 Oct 2021 04:42:18 +1300 Subject: [PATCH 31/44] simplify transfer matrix file handling don't store matrices filled with zeros for components that are null --- xfaster/xfaster_class.py | 112 +++++++++++++++++---------------------- 1 file changed, 48 insertions(+), 64 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 1725acff..725ca691 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -1927,7 +1927,7 @@ def process_gcorr(gcorr_file_in): bd = gdata["bin_def"]["cmb_{}".format(k)] if len(bd0) < len(bd): bd = bd[: len(bd0)] - gcorr[k] = g[: len(bd)] + gcorr[k] = g[: len(bd)] if not np.all(bd0 == bd): self.warn( "G correction for map {} has incompatible bin def".format( @@ -4045,18 +4045,18 @@ def binup2(d, bd, bw): continue if hasattr(self, "transfer_matrix"): + xfermat = self.transfer_matrix[xname].get(spec, {}) + # construct modified shape spectrum C^X_l = sum_{l',Y} J^XY_ll' * C^Y_l' - for spec2 in specs: + for spec2 in set(specs) & set(xfermat): # use correct shape spectrum if comp == "fg": cl1 = cls_shape["fg"][lk] else: cl1 = cls_shape["cmb_{}".format(spec2)][lk] - xfermat = self.transfer_matrix[xname][spec][spec2][lk, lk] - print(xfermat.shape, cl1.shape) # matrix multiplication over l', loop sum over Y - s_arr[xi] += np.einsum("ij,j->i", xfermat, cl1) + s_arr[xi] += np.einsum("ij,j->i", xfermat[spec2][lk, lk], cl1) else: # use correct shape spectrum if comp == "fg": @@ -5461,10 +5461,11 @@ def fisher_iterate( ) if not transfer_run: - # out.update( - # qb_transfer=self.qb_transfer, - # transfer=self.transfer, - # ) + if not hasattr(self, "transfer_matrix"): + out.update( + qb_transfer=self.qb_transfer, + transfer=self.transfer, + ) if self.template_cleaned: out.update(template_alpha=self.template_alpha) @@ -5859,68 +5860,51 @@ def get_transfer_matrix(self, file_root=None): if ret is not None: return ret["transfer_matrix"] - lk = slice(0, self.lmax + 1) - - def load_matrix(spec_in_in, spec_out_in): - # there exist tag 90x150a but not 150ax90 - # TODO: better way to make this isensitive to upper/lower case - fname = os.path.join( - file_root, - "{}x{}_{}_to_{}_block.dat".format( - tag1, tag2, spec_in_in.upper(), spec_out_in.upper() - ), - ) - if os.path.isfile(fname): - return np.loadtxt(fname) - else: - fname = os.path.join( - file_root, - "{}x{}_{}_to_{}_block.dat".format( - tag2, tag1, spec_in_in.upper(), spec_out_in.upper() - ), - ) - if os.path.isfile(fname): - return np.loadtxt(fname) - else: - return None - - def get_cross_matrix(cross_spec): - s1, s2 = cross_spec - s1 = s1 * 2 - s2 = s2 * 2 - return np.sqrt(np.abs(np.multiply(load_matrix(s1, s1), - load_matrix(s2, s2)))) - transfer_matrix = OrderedDict() for xname, (tag1, tag2) in self.map_pairs.items(): dx = transfer_matrix.setdefault(xname, OrderedDict()) - - #print("\nget transfer matrix: ", self.specs) + + ftag = None + for t in ["{}x{}".format(tag1, tag2), "{}x{}".format(tag2, tag1)]: + fname = os.path.join(file_root, "{}_*_to_*_block.dat".format(t)) + if len(glob.glob(fname)): + ftag = t + break + else: + self.warn("Missing transfer matrix for {}".format(xname)) + continue + + # file name format string + fname = os.path.join(file_root, "{}_{{}}_to_{{}}_block.dat".format(ftag)) + for spec_out in self.specs: dsi = dx.setdefault(spec_out, OrderedDict()) for spec_in in self.specs: - #print("\nspecs: ", xname, tag1, tag2, spec_out, spec_in) - # case auto1- => auto2- and - # auto1- => auto1- - tmp = load_matrix(spec_in, spec_out) - if tmp is None: - # cross- => auto- and - # auto1- => auto2- whose files dont exist - # TODO: this is floppy. fix this. - matrix_shape = dx[self.specs[0]][self.specs[0]].shape - if spec_out[0]==spec_out[1]: - tmp = np.zeros(matrix_shape) - else: - # cross1- => cross1- - tmp = get_cross_matrix(spec_out) - # auto- => cross- and - # cross1- => cross2- - if spec_in != spec_out: - tmp = np.zeros(matrix_shape) - - dsi[spec_in] = tmp[lk, lk] - del tmp + if spec_in in ["tt", "ee", "bb"] and spec_out in ["tt", "ee", "bb"]: + # auto- => auto- + fname1 = fname.format(spec_in.upper(), spec_out.upper()) + if not os.path.exists(fname1): + # some are not computed + continue + mat = np.loadtxt(fname1) + + elif spec_in == spec_out: + # cross1- => cross1- + s1, s2 = spec_in.upper() + mat = np.loadtxt(fname.format(s1 * 2, s1 * 2)) + mat *= np.loadtxt(fname.format(s2 * 2, s2 * 2)) + mat = np.sqrt(np.abs(mat)) + + else: + # cross- => auto- + continue + + # populate matrix up to lmax + dsi[spec_in] = np.zeros((self.lmax + 1, self.lmax + 1)) + nx, ny = min([s, self.lmax + 1]) for s in mat.shape + dsi[spec_in][:nx, :ny] = mat[:nx, :ny] + del mat self.transfer_file_root = file_root self.transfer_matrix = transfer_matrix From 5d0b0a090c47aa0fc6ca801cf819ee975ede8231 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Sun, 17 Oct 2021 23:54:29 +0200 Subject: [PATCH 32/44] adding [] to a list syntax --- 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 13fb1161..0055df4b 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -1925,6 +1925,8 @@ def process_gcorr(gcorr_file_in): # check bins match g bd0 = self.bin_def["cmb_{}".format(k)] bd = gdata["bin_def"]["cmb_{}".format(k)] + print('bd0:', bd0) + print('bd: ', bd) if len(bd0) < len(bd): bd = bd[: len(bd0)] gcorr[k] = g[: len(bd)] @@ -5902,7 +5904,7 @@ def get_transfer_matrix(self, file_root=None): # populate matrix up to lmax dsi[spec_in] = np.zeros((self.lmax + 1, self.lmax + 1)) - nx, ny = min([s, self.lmax + 1]) for s in mat.shape + nx, ny = [min([s, self.lmax + 1]) for s in mat.shape] dsi[spec_in][:nx, :ny] = mat[:nx, :ny] del mat From ab3de721811d26e41d41289d90a7acb6539c6793 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Sun, 17 Oct 2021 23:55:20 +0200 Subject: [PATCH 33/44] remove debug printing --- xfaster/xfaster_class.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 0055df4b..6adcd5b8 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -1925,8 +1925,6 @@ def process_gcorr(gcorr_file_in): # check bins match g bd0 = self.bin_def["cmb_{}".format(k)] bd = gdata["bin_def"]["cmb_{}".format(k)] - print('bd0:', bd0) - print('bd: ', bd) if len(bd0) < len(bd): bd = bd[: len(bd0)] gcorr[k] = g[: len(bd)] From bf53e5a8d744cfe5f9c735a3cc232df1cd8e1e1c Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Mon, 18 Oct 2021 21:51:52 +1300 Subject: [PATCH 34/44] replace Mll altogether with transfer matrix, rather than introducing another sum over Kll further simplify file handling --- xfaster/xfaster_class.py | 120 ++++++++++++++++++--------------------- xfaster/xfaster_exec.py | 12 ++-- 2 files changed, 61 insertions(+), 71 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 6adcd5b8..72c6d0fb 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -3803,7 +3803,8 @@ def get_bin_def( def kernel_precalc(self, 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. + ``bin_cl_template`` to pre-compute kernel terms. Constructs M_ll' + directly from the ``transfer_matrix`` attribute, if present. Arguments --------- @@ -3840,6 +3841,10 @@ def kernel_precalc(self, map_tag=None, transfer_run=False): mll = OrderedDict() use_transfer_matrix = hasattr(self, "transfer_matrix") + if use_transfer_matrix and transfer_run: + raise ValueError( + "Transfer matrix cannot be used to compute transfer functions." + ) for spec in specs: mll[spec] = OrderedDict() @@ -3847,19 +3852,26 @@ def kernel_precalc(self, map_tag=None, transfer_run=False): mspec = "{}_mix".format(spec) mll[mspec] = OrderedDict() - # assume that beam and filtering transfer effects are included in - # the transfer matrix if not use_transfer_matrix: bw = self.beam_windows[spec] tf = self.transfer[spec] for xname, (m0, m1) in map_pairs.items(): - if not use_transfer_matrix: - # beams - fb2 = bw[m0][lk] * bw[m1][lk] - # transfer function - if not transfer_run: - fb2 *= np.sqrt(tf[m0][lk] * tf[m1][lk]) + # extract transfer matrix terms + if use_transfer_matrix: + tm = self.transfer_matrix[xname][spec] + mll[spec][xname] = tm[spec][:, lk] + if spec in ["ee", "bb"]: + spec2 = "bb" if spec == "ee" else "ee" + mll[mspec][xname] = tm[spec2][:, lk] + continue + + # beams + fb2 = bw[m0][lk] * bw[m1][lk] + + # transfer function + if not transfer_run: + fb2 *= np.sqrt(tf[m0][lk] * tf[m1][lk]) # kernels if spec == "tt": @@ -3873,9 +3885,9 @@ def kernel_precalc(self, map_tag=None, transfer_run=False): k = self.pkern[xname][:, lk] - self.mkern[xname][:, lk] # store final product - mll[spec][xname] = k if use_transfer_matrix else k * fb2 + mll[spec][xname] = k * fb2 if spec in ["ee", "bb"]: - mll[mspec][xname] = mk if use_transfer_matrix else mk * fb2 + mll[mspec][xname] = mk * fb2 return mll @@ -4044,26 +4056,12 @@ def binup2(d, bd, bw): continue - if hasattr(self, "transfer_matrix"): - xfermat = self.transfer_matrix[xname].get(spec, {}) - - # construct modified shape spectrum C^X_l = sum_{l',Y} J^XY_ll' * C^Y_l' - for spec2 in set(specs) & set(xfermat): - # use correct shape spectrum - if comp == "fg": - cl1 = cls_shape["fg"][lk] - else: - cl1 = cls_shape["cmb_{}".format(spec2)][lk] - - # matrix multiplication over l', loop sum over Y - s_arr[xi] += np.einsum("ij,j->i", xfermat[spec2][lk, lk], cl1) + # use correct shape spectrum + if comp == "fg": + # single foreground spectrum + s_arr[xi] = cls_shape["fg"][lk] else: - # use correct shape spectrum - if comp == "fg": - # single foreground spectrum - s_arr[xi] = cls_shape["fg"][lk] - else: - s_arr[xi] = cls_shape["cmb_{}".format(spec)][lk] + s_arr[xi] = cls_shape["cmb_{}".format(spec)][lk] # use correct beam error shape if beam_error: @@ -5814,7 +5812,7 @@ def get_transfer_matrix(self, file_root=None): Read in an ell-by-ell transfer matrix, constructed from ratios of `anafast` spectra using an ensemble of simulations. This matrix should be used in place of the standard XFaster transfer function calculation, - and includes both filtering and beam effects. + and includes masking, filtering and beam effects. Arguments --------- @@ -5835,19 +5833,15 @@ def get_transfer_matrix(self, file_root=None): ----- This function is called at the ``transfer_matrix`` checkpoint, and is meant to be an alternative to the standard ``transfer`` sequence. - - This option uses a significant amount of RAM, especially for large - numbers of map crosses, since the entire matrix is stored in memory - for efficient computation. """ tm_shape = ( - self.num_corr * (self.num_spec ** 2) * (self.lmax + 1), + self.num_corr * (self.num_spec + 2) * (self.lmax + 1), self.lmax + 1, ) save_name = "transfer_matrix" - save_attrs = ["transfer_file_root", "transfer_matrix"] + save_attrs = ["transfer_matrix_root", "transfer_matrix"] ret = self.load_data( save_name, "transfer_matrix", @@ -5864,15 +5858,17 @@ def get_transfer_matrix(self, file_root=None): for xname, (tag1, tag2) in self.map_pairs.items(): dx = transfer_matrix.setdefault(xname, OrderedDict()) - ftag = None - for t in ["{}x{}".format(tag1, tag2), "{}x{}".format(tag2, tag1)]: - fname = os.path.join(file_root, "{}_*_to_*_block.dat".format(t)) - if len(glob.glob(fname)): - ftag = t - break + if tag1 == tag2: + ftag = "{}x{}".format(tag1, tag2) else: - self.warn("Missing transfer matrix for {}".format(xname)) - continue + for t in ["{}x{}".format(tag1, tag2), "{}x{}".format(tag2, tag1)]: + fname = os.path.join(file_root, "{}_*_to_*_block.dat".format(t)) + if len(glob.glob(fname)): + ftag = t + break + else: + self.warn("Missing transfer matrix for {}".format(xname)) + continue # file name format string fname = os.path.join(file_root, "{}_{{}}_to_{{}}_block.dat".format(ftag)) @@ -5880,33 +5876,27 @@ def get_transfer_matrix(self, file_root=None): for spec_out in self.specs: dsi = dx.setdefault(spec_out, OrderedDict()) - for spec_in in self.specs: - if spec_in in ["tt", "ee", "bb"] and spec_out in ["tt", "ee", "bb"]: + specs_in = ["ee", "bb"] if spec_out in ["ee", "bb"] else [spec_out] + for spec_in in specs_in: + if spec_out in ["tt", "ee", "bb"]: # auto- => auto- fname1 = fname.format(spec_in.upper(), spec_out.upper()) - if not os.path.exists(fname1): - # some are not computed - continue mat = np.loadtxt(fname1) + # transpose for input axis on second dimension ? + # mat = mat.T - elif spec_in == spec_out: - # cross1- => cross1- - s1, s2 = spec_in.upper() - mat = np.loadtxt(fname.format(s1 * 2, s1 * 2)) - mat *= np.loadtxt(fname.format(s2 * 2, s2 * 2)) - mat = np.sqrt(np.abs(mat)) + # populate matrix up to lmax + dsi[spec_in] = np.zeros((self.lmax + 1, self.lmax + 1)) + nx, ny = [min([s, self.lmax + 1]) for s in mat.shape] + dsi[spec_in][:nx, :ny] = mat[:nx, :ny] + del mat else: - # cross- => auto- - continue - - # populate matrix up to lmax - dsi[spec_in] = np.zeros((self.lmax + 1, self.lmax + 1)) - nx, ny = [min([s, self.lmax + 1]) for s in mat.shape] - dsi[spec_in][:nx, :ny] = mat[:nx, :ny] - del mat + # cross1- => cross1- + s1, s2 = [s + s for s in spec_in] + dsi[spec_in] = np.sqrt(np.abs(dx[s1][s1] * dx[s2][s2])) - self.transfer_file_root = file_root + self.transfer_matrix_root = file_root self.transfer_matrix = transfer_matrix self.save_data(save_name, from_attrs=save_attrs) diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index 2d4a8010..47d6d43f 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -666,17 +666,17 @@ def xfaster_run( apply_gcorr=apply_gcorr, reload_gcorr=reload_gcorr, gcorr_file=gcorr_file ) - X.log("Computing kernels...", "notice") - X.get_kernels(window_lmax=window_lmax) - - X.log("Computing beam window functions...", "notice") - X.get_beams(pixwin=pixwin) - if transfer_matrix_root: X.log("Loading pre-computed transfer matrix...", "notice") X.get_transfer_matrix(file_root=transfer_matrix_root) else: + X.log("Computing kernels...", "notice") + X.get_kernels(window_lmax=window_lmax) + + X.log("Computing beam window functions...", "notice") + X.get_beams(pixwin=pixwin) + 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] From ade66be2d3d11e6acb05382e4ed83d29ba2dab07 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Tue, 19 Oct 2021 04:18:09 +0200 Subject: [PATCH 35/44] Transpose the transfer matrix because the second axis is input --- 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 72c6d0fb..342620fd 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -5882,8 +5882,8 @@ def get_transfer_matrix(self, file_root=None): # auto- => auto- fname1 = fname.format(spec_in.upper(), spec_out.upper()) mat = np.loadtxt(fname1) - # transpose for input axis on second dimension ? - # mat = mat.T + # transpose for input axis on second dimension. + mat = mat.T # populate matrix up to lmax dsi[spec_in] = np.zeros((self.lmax + 1, self.lmax + 1)) From 4ba398ffb21398eaa912866b4e39e657c3c3c012 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Thu, 21 Oct 2021 03:27:29 +1300 Subject: [PATCH 36/44] script for building a dummy transfer matrix for debugging --- example/make_example_transfer_matrix.py | 31 +++++++++++++++++++++++++ 1 file changed, 31 insertions(+) create mode 100644 example/make_example_transfer_matrix.py diff --git a/example/make_example_transfer_matrix.py b/example/make_example_transfer_matrix.py new file mode 100644 index 00000000..67766c9c --- /dev/null +++ b/example/make_example_transfer_matrix.py @@ -0,0 +1,31 @@ +import numpy as np +import xfaster as xf + +kerns = xf.load_and_parse("outputs_example/95x150/kernels_95x150.npz") +beams = xf.load_and_parse("outputs_example/95x150/beams_95x150.npz") +tf = np.loadtxt("maps_example/transfer_example.txt") +tf.shape +lmax = 500 +lk = slice(0, lmax + 1) +bw = beams["beam_windows"] + +for xname in kerns["kern"].keys(): + tag1, tag2 = xname.split(":") + fname = "transfer_matrix/{}x{}_{{}}_to_{{}}_block.dat".format(tag1, tag2) + for spec in ["tt", "ee", "bb"]: + fb2 = tf[lk] * bw[spec][tag1][lk] * bw[spec][tag2][lk] + + if spec == "tt": + k = kerns["kern"][xname][..., lk] + else: + k = kerns["pkern"][xname][..., lk] + mk = kerns["mkern"][xname][..., lk] + + mat = k * fb2 + su = spec.upper() + np.savetxt(fname.format(su, su), mat.T) + + if spec != "tt": + mat = mk * fb2 + su2 = "BB" if spec == "ee" else "EE" + np.savetxt(fname.format(su, su2), mat.T) From 14554ffe2df0a1dff1a91c95b4d0d427e47377b7 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Thu, 21 Oct 2021 03:32:17 +1300 Subject: [PATCH 37/44] makedirs --- example/make_example_transfer_matrix.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/example/make_example_transfer_matrix.py b/example/make_example_transfer_matrix.py index 67766c9c..895be41f 100644 --- a/example/make_example_transfer_matrix.py +++ b/example/make_example_transfer_matrix.py @@ -1,5 +1,6 @@ import numpy as np import xfaster as xf +import os kerns = xf.load_and_parse("outputs_example/95x150/kernels_95x150.npz") beams = xf.load_and_parse("outputs_example/95x150/beams_95x150.npz") @@ -9,6 +10,8 @@ lk = slice(0, lmax + 1) bw = beams["beam_windows"] +os.makedirs("transfer_matrix") + for xname in kerns["kern"].keys(): tag1, tag2 = xname.split(":") fname = "transfer_matrix/{}x{}_{{}}_to_{{}}_block.dat".format(tag1, tag2) From 5401d10f5e6ce9a30d2235e857ea66e831901253 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Thu, 21 Oct 2021 19:02:48 +0200 Subject: [PATCH 38/44] transfer_matrix branch has bug when transfer_matrix_root=None. debug --- 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 342620fd..c20c932c 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -3852,7 +3852,7 @@ def kernel_precalc(self, map_tag=None, transfer_run=False): mspec = "{}_mix".format(spec) mll[mspec] = OrderedDict() - if not use_transfer_matrix: + if not use_transfer_matrix and not transfer_run: bw = self.beam_windows[spec] tf = self.transfer[spec] From 431cb8c246aef30484514d4cf89a994985ae8dfb Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Thu, 21 Oct 2021 19:10:29 +0200 Subject: [PATCH 39/44] continue debug .. --- xfaster/xfaster_class.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index c20c932c..4af1ae8c 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -3852,9 +3852,10 @@ def kernel_precalc(self, map_tag=None, transfer_run=False): mspec = "{}_mix".format(spec) mll[mspec] = OrderedDict() - if not use_transfer_matrix and not transfer_run: + if not use_transfer_matrix: bw = self.beam_windows[spec] - tf = self.transfer[spec] + if not transfer_run: + tf = self.transfer[spec] for xname, (m0, m1) in map_pairs.items(): # extract transfer matrix terms From 41c2c3644b539453f2fa7bae64b2f8fc2318da56 Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Mon, 25 Oct 2021 19:59:24 +0200 Subject: [PATCH 40/44] remove transpose of the anafast transfer matrix after test on simple (cmb + noise) sims --- 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 4af1ae8c..28270769 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -5884,7 +5884,7 @@ def get_transfer_matrix(self, file_root=None): fname1 = fname.format(spec_in.upper(), spec_out.upper()) mat = np.loadtxt(fname1) # transpose for input axis on second dimension. - mat = mat.T + #mat = mat.T # populate matrix up to lmax dsi[spec_in] = np.zeros((self.lmax + 1, self.lmax + 1)) From 0bbc1d68eb5538dc4c9992aaf7bf4c99bac7e61f Mon Sep 17 00:00:00 2001 From: Vy Luu Date: Tue, 26 Oct 2021 18:44:45 +0200 Subject: [PATCH 41/44] move ref_freq and beta_ref from get_bandpower to get_transfer for code to work in foreground_fit=True case --- xfaster/xfaster_class.py | 19 +++++++++++++++---- xfaster/xfaster_exec.py | 4 ++-- 2 files changed, 17 insertions(+), 6 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index 28270769..4f2ee05d 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -5621,6 +5621,8 @@ def get_transfer( iter_max=200, save_iters=False, fix_bb_transfer=False, + ref_freq=359.7, + beta_ref=.154, ): """ Compute the transfer function from signal simulations created using @@ -5719,6 +5721,15 @@ def expand_transfer(qb_transfer, bin_def): else: self.transfer = expand_transfer(ret["qb_transfer"], ret["bin_def"]) return self.transfer + + # foreground fitting, required in fisher_iterate: + # TODO: JUST TO PASS CODE ERROR; REPEATED IN GET_GETBANDPOWER, CLEAN UP + 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 self.qb_transfer = OrderedDict() for spec in self.specs: @@ -6011,11 +6022,11 @@ 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 + # reference frequency and spectral index, move to get_transfer + #self.ref_freq = ref_freq + #self.beta_ref = beta_ref # priors on frequency spectral index - self.delta_beta_fix = 1.0e-8 + #self.delta_beta_fix = 1.0e-8 opts.update( delta_beta_prior=delta_beta_prior, ref_freq=self.ref_freq, diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index 47d6d43f..7e567429 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -615,8 +615,8 @@ def xfaster_run( 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("ref_freq") + #transfer_opts.pop("beta_ref") transfer_opts.pop("delta_beta_prior") transfer_opts.pop("null_first_cmb") transfer_opts.pop("return_cls") From 7a300936dd02761fac28a974371ce1d8e4bdfc59 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Tue, 2 Nov 2021 21:41:03 +1300 Subject: [PATCH 42/44] correct bug fix for get_transfer with foreground fitting enabled --- xfaster/xfaster_class.py | 27 +++++++-------------------- xfaster/xfaster_exec.py | 4 ++-- 2 files changed, 9 insertions(+), 22 deletions(-) diff --git a/xfaster/xfaster_class.py b/xfaster/xfaster_class.py index a74a0723..69d1467e 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -3993,7 +3993,7 @@ def kernel_precalc(self, map_tag=None, transfer_run=False): mspec = "{}_mix".format(spec) mll[mspec] = OrderedDict() - if not use_transfer_matrix: + if not use_transfer_matrix: bw = self.beam_windows[spec] if not transfer_run: tf = self.transfer[spec] @@ -5462,7 +5462,7 @@ def fisher_iterate( extra_tag=file_tag, ) - if "fg_tt" in self.bin_def: + if "fg_tt" in self.bin_def and not transfer_run: out.update( beta_fit=beta_fit, beta_err=beta_err, @@ -5579,7 +5579,7 @@ def fisher_iterate( weighted_bins=self.weighted_bins, ) - if "fg_tt" in self.bin_def: + if "fg_tt" in self.bin_def and not transfer_run: out.update( delta_beta_prior=delta_beta_prior, beta_fit=beta_fit, @@ -5762,8 +5762,6 @@ def get_transfer( iter_max=200, save_iters=False, fix_bb_transfer=False, - ref_freq=359.7, - beta_ref=.154, ): """ Compute the transfer function from signal simulations created using @@ -5862,15 +5860,6 @@ def expand_transfer(qb_transfer, bin_def): else: self.transfer = expand_transfer(ret["qb_transfer"], ret["bin_def"]) return self.transfer - - # foreground fitting, required in fisher_iterate: - # TODO: JUST TO PASS CODE ERROR; REPEATED IN GET_GETBANDPOWER, CLEAN UP - 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 self.qb_transfer = OrderedDict() for spec in self.specs: @@ -6035,8 +6024,6 @@ def get_transfer_matrix(self, file_root=None): # auto- => auto- fname1 = fname.format(spec_in.upper(), spec_out.upper()) mat = np.loadtxt(fname1) - # transpose for input axis on second dimension. - #mat = mat.T # populate matrix up to lmax dsi[spec_in] = np.zeros((self.lmax + 1, self.lmax + 1)) @@ -6163,11 +6150,11 @@ def get_bandpowers( # foreground fitting if "fg_tt" in self.bin_def: - # reference frequency and spectral index, move to get_transfer - #self.ref_freq = ref_freq - #self.beta_ref = beta_ref + # 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 + self.delta_beta_fix = 1.0e-8 opts.update( delta_beta_prior=delta_beta_prior, ref_freq=self.ref_freq, diff --git a/xfaster/xfaster_exec.py b/xfaster/xfaster_exec.py index ae25d9ff..35c3b698 100644 --- a/xfaster/xfaster_exec.py +++ b/xfaster/xfaster_exec.py @@ -609,8 +609,8 @@ def xfaster_run( 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("ref_freq") + transfer_opts.pop("beta_ref") transfer_opts.pop("delta_beta_prior") transfer_opts.pop("null_first_cmb") transfer_opts.pop("return_cls") From 8307efe7c264e1a5ba43e350f914c3ce1913d342 Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Fri, 18 Feb 2022 17:39:27 -0600 Subject: [PATCH 43/44] use correct tags in kernel_precalc --- 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 a4314138..69c2fda3 100644 --- a/xfaster/xfaster_class.py +++ b/xfaster/xfaster_class.py @@ -4071,10 +4071,10 @@ def kernel_precalc(self, map_tag=None, transfer=False): # extract transfer matrix terms if use_transfer_matrix: tm = self.transfer_matrix[xname][spec] - mll[spec][xname] = tm[spec][:, lk] + mll[stag][xname] = tm[spec][:, lk] if spec in ["ee", "bb"]: spec2 = "bb" if spec == "ee" else "ee" - mll[mspec][xname] = tm[spec2][:, lk] + mll[mstag][xname] = tm[spec2][:, lk] continue # beams From 5ad3ce538851b6e1b1fc53db01d4d75e83bd13eb Mon Sep 17 00:00:00 2001 From: Sasha Rahlin Date: Wed, 20 Jul 2022 23:48:20 -0500 Subject: [PATCH 44/44] fix working directory in github action --- .github/workflows/pull_request.yaml | 5 +++-- .github/workflows/release.yaml | 5 +++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/.github/workflows/pull_request.yaml b/.github/workflows/pull_request.yaml index c1d71eff..411e5ee8 100644 --- a/.github/workflows/pull_request.yaml +++ b/.github/workflows/pull_request.yaml @@ -22,6 +22,7 @@ jobs: run: | python -m pip install -e . - name: Run tutorial to test + working-directory: docs/notebooks run: | - jupyter nbconvert --to script docs/notebooks/XFaster_Tutorial.ipynb - python docs/notebooks/XFaster_Tutorial.py + jupyter nbconvert --to script XFaster_Tutorial.ipynb + python XFaster_Tutorial.py diff --git a/.github/workflows/release.yaml b/.github/workflows/release.yaml index dc66c866..f909f17c 100644 --- a/.github/workflows/release.yaml +++ b/.github/workflows/release.yaml @@ -22,9 +22,10 @@ jobs: python -m build python -m pip install dist/xfaster*.whl - name: Run tutorial to test + working-directory: docs/notebooks run: | - jupyter nbconvert --to script docs/notebooks/XFaster_Tutorial.ipynb - python docs/notebooks/XFaster_Tutorial.py + jupyter nbconvert --to script XFaster_Tutorial.ipynb + python XFaster_Tutorial.py - name: Upload to PyPI env: TWINE_USERNAME: __token__