From 0c7874e1c5d95c93943c592bb9be0cbb6735a8e0 Mon Sep 17 00:00:00 2001 From: Maximilian H Abitbol Date: Wed, 11 Nov 2020 16:41:29 +0000 Subject: [PATCH 01/20] hybrid stuff i think --- bbpower/__init__.py | 5 ----- bbpower/compsep.py | 36 +++++++++++++----------------------- 2 files changed, 13 insertions(+), 28 deletions(-) diff --git a/bbpower/__init__.py b/bbpower/__init__.py index 2635b16..3d118fb 100644 --- a/bbpower/__init__.py +++ b/bbpower/__init__.py @@ -1,8 +1,3 @@ from bbpipe import PipelineStage -from .mask_preproc import BBMaskPreproc -from .maps_preproc import BBMapsPreproc -from .power_specter import BBPowerSpecter -from .power_summarizer import BBPowerSummarizer -from .covfefe import BBCovFeFe from .compsep import BBCompSep from .plotter import BBPlotter diff --git a/bbpower/compsep.py b/bbpower/compsep.py index b4c54ce..7cb7cc1 100644 --- a/bbpower/compsep.py +++ b/bbpower/compsep.py @@ -30,6 +30,7 @@ def setup_compsep(self): self.load_cmb() self.fg_model = FGModel(self.config) self.params = ParameterManager(self.config) + self.hybridparams = np.load(self.config['resid_seds']) if self.use_handl: self.prepare_h_and_l() return @@ -70,16 +71,13 @@ def parse_sacc_file(self): """ Reads the data in the sacc file included the power spectra, bandpasses, and window functions. """ - #Decide if you're using H&L self.use_handl = self.config['likelihood_type'] == 'h&l' - #Read data self.s = sacc.Sacc.load_fits(self.get_input('cells_coadded')) if self.use_handl: s_fid = sacc.Sacc.load_fits(self.get_input('cells_fiducial')) s_noi = sacc.Sacc.load_fits(self.get_input('cells_noise')) - # Keep only desired correlations self.pols = self.config['pol_channels'] corr_all = ['cl_ee', 'cl_eb', 'cl_be', 'cl_bb'] corr_keep = [] @@ -94,7 +92,6 @@ def parse_sacc_file(self): s_fid.remove_selection(c) s_noi.remove_selection(c) - # Scale cuts self.s.remove_selection(ell__gt=self.config['l_max']) self.s.remove_selection(ell__lt=self.config['l_min']) if self.use_handl: @@ -111,7 +108,6 @@ def parse_sacc_file(self): self.ncross = (self.nmaps * (self.nmaps + 1)) // 2 self.pol_order=dict(zip(self.pols,range(self.npol))) - #Collect bandpasses self.bpss = [] for i_t, tn in enumerate(tr_names): t = self.s.tracers[tn] @@ -123,24 +119,18 @@ def parse_sacc_file(self): bnu = t.bandpass self.bpss.append(Bandpass(nu, dnu, bnu, i_t+1, self.config)) - # Get ell sampling - # Example power spectrum self.ell_b, _ = self.s.get_ell_cl('cl_' + 2 * self.pols[0].lower(), tr_names[0], tr_names[0]) - #Avoid l<2 win0 = self.s.data[0]['window'] mask_w = win0.values > 1 self.bpw_l = win0.values[mask_w] self.n_ell = len(self.bpw_l) self.n_bpws = len(self.ell_b) - # D_ell factor self.dl2cl = 2 * np.pi / (self.bpw_l * (self.bpw_l + 1)) self.windows = np.zeros([self.ncross, self.n_bpws, self.n_ell]) - #Get power spectra and covariances if not (self.s.covariance.covmat.shape[-1] == len(self.s.mean) == self.n_bpws * self.ncross): raise ValueError("C_ell vector's size is wrong") - ###cv = self.s.precision.getCovarianceMatrix() v2d = np.zeros([self.n_bpws, self.ncross]) if self.use_handl: @@ -151,7 +141,6 @@ def parse_sacc_file(self): self.vector_indices = self.vector_to_matrix(np.arange(self.ncross, dtype=int)).astype(int) self.indx = [] - #Parse into the right ordering for b1, b2, p1, p2, m1, m2, ind_vec in self._freq_pol_iterator(): t1 = tr_names[b1] t2 = tr_names[b2] @@ -196,7 +185,6 @@ def load_cmb(self): mask = (self.cmb_ells <= self.bpw_l.max()) & (self.cmb_ells > 1) self.cmb_ells = self.cmb_ells[mask] - # TODO: this is a patch nell = len(self.cmb_ells) self.cmb_tens = np.zeros([self.npol, self.npol, nell]) self.cmb_lens = np.zeros([self.npol, self.npol, nell]) @@ -221,15 +209,23 @@ def integrate_seds(self, params): sed_params = [params[comp['names_sed_dict'][k]] for k in comp['sed'].params] rot_matrices.append([]) - def sed(nu): - return comp['sed'].eval(nu, *sed_params) + + if self.config.get("diff"): + def sed(nu): + return comp['sed'].diff(nu, *sed_params) + else: + def sed(nu): + return comp['sed'].eval(nu, *sed_params) for tn in range(self.nfreqs): sed_b, rot = self.bpss[tn].convolve_sed(sed, params) fg_scaling[i_c, tn] = sed_b * units rot_matrices[i_c].append(rot) + + if self.config.get("diff"): + fg_scaling[i_c] = np.dot(self.hybridparams['Q'], fg_scaling[i_c]) - return fg_scaling.T,rot_matrices + return fg_scaling.T, rot_matrices def evaluate_power_spectra(self, params): fg_pspectra = np.zeros([self.fg_model.n_components, @@ -268,28 +264,23 @@ def model(self, params): fg_scaling, rot_m = self.integrate_seds(params) # [nfreq, ncomp], [ncomp,nfreq,[matrix]] fg_cell = self.evaluate_power_spectra(params) # [ncomp,ncomp,npol,npol,nell] - # Add all components scaled in frequency (and HWP-rotated if needed) cls_array_fg = np.zeros([self.nfreqs,self.nfreqs,self.n_ell,self.npol,self.npol]) fg_cell = np.transpose(fg_cell, axes = [0,1,4,2,3]) # [ncomp,ncomp,nell,npol,npol] cmb_cell = np.transpose(cmb_cell, axes = [2,0,1]) # [nell,npol,npol] for f1 in range(self.nfreqs): - for f2 in range(f1,self.nfreqs): # Note that we only need to fill in half of the frequencies + for f2 in range(f1,self.nfreqs): cls=cmb_cell.copy() - # Loop over component pairs for c1 in range(self.fg_model.n_components): mat1=rot_m[c1][f1] a1=fg_scaling[f1,c1] for c2 in range(self.fg_model.n_components): mat2=rot_m[c2][f2] a2=fg_scaling[f2,c2] - # Rotate if needed clrot=rotate_cells_mat(mat2,mat1,fg_cell[c1,c2]) - # Scale in frequency and add cls += clrot*a1*a2 cls_array_fg[f1,f2]=cls - # Window convolution cls_array_list = np.zeros([self.n_bpws, self.nfreqs, self.npol, self.nfreqs, self.npol]) for f1 in range(self.nfreqs): for p1 in range(self.npol): @@ -304,7 +295,6 @@ def model(self, params): if m1!=m2: cls_array_list[:, f2, p2, f1, p1] = clband - # Polarization angle rotation for f1 in range(self.nfreqs): for f2 in range(self.nfreqs): cls_array_list[:,f1,:,f2,:] = rotate_cells(self.bpss[f2], self.bpss[f1], From 14a52bb18571ccb0dbe6e4eb40efda7546552e9f Mon Sep 17 00:00:00 2001 From: Maximilian H Abitbol Date: Mon, 8 Feb 2021 23:01:42 +0000 Subject: [PATCH 02/20] moments option --- bbpower/compsep.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/bbpower/compsep.py b/bbpower/compsep.py index 7cb7cc1..ad0841d 100644 --- a/bbpower/compsep.py +++ b/bbpower/compsep.py @@ -30,7 +30,8 @@ def setup_compsep(self): self.load_cmb() self.fg_model = FGModel(self.config) self.params = ParameterManager(self.config) - self.hybridparams = np.load(self.config['resid_seds']) + if self.config.get("diff"): + self.hybridparams = np.load(self.config['resid_seds']) if self.use_handl: self.prepare_h_and_l() return From 32ae7b22a6b54ee22b1cb1cafbb4b42b920b7734 Mon Sep 17 00:00:00 2001 From: Maximilian H Abitbol Date: Tue, 9 Feb 2021 15:06:44 +0000 Subject: [PATCH 03/20] testing token --- tester | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 tester diff --git a/tester b/tester new file mode 100644 index 0000000..e69de29 From 67663754e906e2cd6ebcadb22a22bb19d8e061bc Mon Sep 17 00:00:00 2001 From: Maximilian H Abitbol Date: Tue, 9 Feb 2021 18:26:50 +0000 Subject: [PATCH 04/20] auto queue over realizations --- autoqbaseline.sh | 29 +++++++++++++++++++++++++++++ autoqresiduals.sh | 25 +++++++++++++++++++++++++ 2 files changed, 54 insertions(+) create mode 100755 autoqbaseline.sh create mode 100755 autoqresiduals.sh diff --git a/autoqbaseline.sh b/autoqbaseline.sh new file mode 100755 index 0000000..0c7c926 --- /dev/null +++ b/autoqbaseline.sh @@ -0,0 +1,29 @@ +#!/bin/bash + +for k in {0..3} +do + for j in {0..19} + do + output="sim0$k/output$j/" + mkdir -p $output + + datadir="/mnt/zfsusers/mabitbol/BBHybrid/data/sim0$k/" + coadd=$datadir"cls_coadd_base$j.fits" + noise=$datadir"cls_noise_base$j.fits" + fid=$datadir"cls_fid_base$j.fits" + chain=$output"param_chains.npz" + configcopy=$output"config_copy.yml" + config="config.yml" + + addqueue -q cmb -c "2 days" -m 1 -s -n 1x12 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config + done +done + +#--cells_coadded=/mnt/zfsusers/mabitbol/BBHybrid/data/sim00/cls_coadd_base0.fits +#--cells_noise=/mnt/zfsusers/mabitbol/BBHybrid/data/sim00/cls_noise_base0.fits +#--cells_fiducial=/mnt/zfsusers/mabitbol/BBHybrid/data/sim00/cls_fid_base0.fits +#--param_chains=./moments/baseline0chi2/output/param_chains.npz +#--config_copy=./moments/baseline0chi2/output/config_copy.yml +#--config=./moments/baseline0chi2/config.yml + +#python3 -m bbpower BBCompSep --cells_coadded=/mnt/zfsusers/mabitbol/BBHybrid/data/sim00/cls_coadd_base0.fits --cells_noise=/mnt/zfsusers/mabitbol/BBHybrid/data/sim00/cls_noise_base0.fits --cells_fiducial=/mnt/zfsusers/mabitbol/BBHybrid/data/sim00/cls_fid_base0.fits --param_chains=./moments/baseline0chi2/output/param_chains.npz --config_copy=./moments/baseline0chi2/output/config_copy.yml --config=./moments/baseline0chi2/config.yml diff --git a/autoqresiduals.sh b/autoqresiduals.sh new file mode 100755 index 0000000..e164c96 --- /dev/null +++ b/autoqresiduals.sh @@ -0,0 +1,25 @@ +#!/bin/bash + +for k in {0..3} +do + for j in {0..19} + do + output="sim0$k/output$j/" + mkdir -p $output + + cp config.yml $output + config=$output"config.yml" + sed -i "s/sim0X/sim0$k/" $config + sed -i "s/paramsY/params$j/" $config + + datadir="/mnt/zfsusers/mabitbol/BBHybrid/data/sim0$k/" + coadd=$datadir"cls_coadd$j.fits" + noise=$datadir"cls_noise$j.fits" + fid=$datadir"cls_fid$j.fits" + chain=$output"param_chains.npz" + configcopy=$output"config_copy.yml" + + addqueue -q cmb -c "2 days" -m 1 -s -n 1x12 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config + done +done + From 1e8ae8c301d817b05801be67e8554198eb960c6f Mon Sep 17 00:00:00 2001 From: Maximilian H Abitbol Date: Tue, 9 Feb 2021 18:41:41 +0000 Subject: [PATCH 05/20] auto queue --- autoqbaseline.sh | 18 ++++++------------ autoqresiduals.sh | 8 +++++--- 2 files changed, 11 insertions(+), 15 deletions(-) diff --git a/autoqbaseline.sh b/autoqbaseline.sh index 0c7c926..22eebde 100755 --- a/autoqbaseline.sh +++ b/autoqbaseline.sh @@ -1,10 +1,12 @@ #!/bin/bash +mdir="baseline_noiseless/" + for k in {0..3} do - for j in {0..19} + for j in {0..20} do - output="sim0$k/output$j/" + output=$mdir"sim0$k/output$j/" mkdir -p $output datadir="/mnt/zfsusers/mabitbol/BBHybrid/data/sim0$k/" @@ -13,17 +15,9 @@ do fid=$datadir"cls_fid_base$j.fits" chain=$output"param_chains.npz" configcopy=$output"config_copy.yml" - config="config.yml" + config=$mdir"config.yml" - addqueue -q cmb -c "2 days" -m 1 -s -n 1x12 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config + addqueue -q cmb -c "2 days" -m 1 -s -n 1x8 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config done done -#--cells_coadded=/mnt/zfsusers/mabitbol/BBHybrid/data/sim00/cls_coadd_base0.fits -#--cells_noise=/mnt/zfsusers/mabitbol/BBHybrid/data/sim00/cls_noise_base0.fits -#--cells_fiducial=/mnt/zfsusers/mabitbol/BBHybrid/data/sim00/cls_fid_base0.fits -#--param_chains=./moments/baseline0chi2/output/param_chains.npz -#--config_copy=./moments/baseline0chi2/output/config_copy.yml -#--config=./moments/baseline0chi2/config.yml - -#python3 -m bbpower BBCompSep --cells_coadded=/mnt/zfsusers/mabitbol/BBHybrid/data/sim00/cls_coadd_base0.fits --cells_noise=/mnt/zfsusers/mabitbol/BBHybrid/data/sim00/cls_noise_base0.fits --cells_fiducial=/mnt/zfsusers/mabitbol/BBHybrid/data/sim00/cls_fid_base0.fits --param_chains=./moments/baseline0chi2/output/param_chains.npz --config_copy=./moments/baseline0chi2/output/config_copy.yml --config=./moments/baseline0chi2/config.yml diff --git a/autoqresiduals.sh b/autoqresiduals.sh index e164c96..7378e09 100755 --- a/autoqresiduals.sh +++ b/autoqresiduals.sh @@ -1,10 +1,12 @@ #!/bin/bash +mdir="residual_noiseless/" + for k in {0..3} do - for j in {0..19} + for j in {0..20} do - output="sim0$k/output$j/" + output=$mdir"sim0$k/output$j/" mkdir -p $output cp config.yml $output @@ -19,7 +21,7 @@ do chain=$output"param_chains.npz" configcopy=$output"config_copy.yml" - addqueue -q cmb -c "2 days" -m 1 -s -n 1x12 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config + addqueue -q cmb -c "2 days" -m 1 -s -n 1x8 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config done done From 2462723ee5c672709ea784a801177ef3aca8a533 Mon Sep 17 00:00:00 2001 From: Maximilian H Abitbol Date: Tue, 9 Feb 2021 19:57:21 +0000 Subject: [PATCH 06/20] plotting --- autoqresiduals.sh | 5 +++-- plotting/labels.py | 38 +++++++++++++++++++++++++++++++++ plotting/saveclean.py | 49 +++++++++++++++++++++++++++++++++++++++++++ tester | 0 4 files changed, 90 insertions(+), 2 deletions(-) create mode 100644 plotting/labels.py create mode 100755 plotting/saveclean.py delete mode 100644 tester diff --git a/autoqresiduals.sh b/autoqresiduals.sh index 7378e09..250f1ee 100755 --- a/autoqresiduals.sh +++ b/autoqresiduals.sh @@ -1,6 +1,7 @@ #!/bin/bash mdir="residual_noiseless/" +cf="residual_noiseless/config.yml" for k in {0..3} do @@ -9,7 +10,7 @@ do output=$mdir"sim0$k/output$j/" mkdir -p $output - cp config.yml $output + cp $cf $output config=$output"config.yml" sed -i "s/sim0X/sim0$k/" $config sed -i "s/paramsY/params$j/" $config @@ -21,7 +22,7 @@ do chain=$output"param_chains.npz" configcopy=$output"config_copy.yml" - addqueue -q cmb -c "2 days" -m 1 -s -n 1x8 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config + addqueue -q berg -c "2 days" -m 1 -s -n 1x8 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config done done diff --git a/plotting/labels.py b/plotting/labels.py new file mode 100644 index 0000000..18e8473 --- /dev/null +++ b/plotting/labels.py @@ -0,0 +1,38 @@ +alllabels = {'A_lens':'A_{lens}', 'r_tensor':'r', 'beta_d':'\\beta_d', 'epsilon_ds':'\epsilon', 'birefringence':'\\theta', + 'alpha_d_ee':'\\alpha^{EE}_d', 'amp_d_ee':'A^{EE}_d', 'alpha_d_bb':'\\alpha^{BB}_d', + 'amp_d_bb':'A^{BB}_d', 'beta_s':'\\beta_s', 'alpha_s_ee':'\\alpha^{EE}_s', + 'amp_s_ee':'A^{EE}_s', 'alpha_s_bb':'\\alpha^{BB}_s', 'amp_s_bb':'A^{BB}_s', + 'amp_d_eb':'A^{EB}_d', 'alpha_d_eb':'\\alpha^{EB}_d', 'amp_s_eb':'A^{EB}_s', + 'alpha_s_eb':'\\alpha^{EB}_s', 'decorr_amp_d':'\Delta_d', 'decorr_amp_s':'\Delta_s', + 'gain_1':'G_1', 'gain_2':'G_2', 'gain_3':'G_3', 'gain_4':'G_4', 'gain_5':'G_5', 'gain_6':'G_6', + 'gain_7':'G_7', 'gain_8':'G_8', 'gain_9':'G_9', 'gain_10':'{G_10}', 'gain_11':'{G_11}', 'gain_12':'{G_12}', + 'shift_1':'\Delta\\nu_1', 'shift_2':'\Delta\\nu_2', 'shift_3':'\Delta\\nu_3', + 'shift_4':'\Delta\\nu_4', 'shift_5':'\Delta\\nu_5', 'shift_6':'\Delta\\nu_6', + 'shift_7':'\Delta\\nu_7', 'shift_8':'\Delta\\nu_8', 'shift_9':'\Delta\\nu_9', + 'shift_10':'\Delta\\nu_{10}', 'shift_11':'\Delta\\nu_{11}', 'shift_12':'\Delta\\nu_{12}', + 'angle_1':'\phi_1', 'angle_2':'\phi_2', 'angle_3':'\phi_3', 'angle_4':'\phi_4', + 'angle_5':'\phi_5', 'angle_6':'\phi_6', 'dphi1_1':'\Delta\phi_1', 'dphi1_2':'\Delta\phi_2', + 'dphi1_3':'\Delta\phi_3', 'dphi1_4':'\Delta\phi_4', 'dphi1_5':'\Delta\phi_5', 'dphi1_6':'\Delta\phi_6'} + + +allranges={'$r$': [-1.0, 1.0], + '$A_{lens}$': [0.0, 2.0], + '$\\beta_d$': [0.1, 10.0], + '$A^{EE}_d$': [0.0, None], + '$\\alpha^{EE}_d$': [-3.0, 1.0], + '$A^{BB}_d$': [0.0, None], + '$\\alpha^{BB}_d$': [-3.0, 1.0], + '$A^{EB}_d$': [-100.0, 100.0], + '$\\alpha^{EB}_d$': [-3.0, 1.0], + '$\\epsilon$': [-1.0, 1.0], + '$\\Delta_d$': [0.9, 1.1], + '$\\beta_s$': [-10.0, 0.0], + '$A^{EE}_s$': [0.0, None], + '$\\alpha^{EE}_s$': [-3.0, 1.0], + '$A^{BB}_s$': [0.0, None], + '$\\alpha^{BB}_s$': [-3.0, 1.0], + '$A^{EB}_s$': [-100.0, 100.0], + '$\\alpha^{EB}_s$': [-3.0, 1.0], + '$\\theta$': [-30., 30.], + '$\\Delta_s$': [0.9, 1.1]} + diff --git a/plotting/saveclean.py b/plotting/saveclean.py new file mode 100755 index 0000000..8680378 --- /dev/null +++ b/plotting/saveclean.py @@ -0,0 +1,49 @@ +import numpy as np +import emcee +from labels import alllabels +import glob + +allfs = ['/mnt/zfsusers/mabitbol/BBPower/moments/hybrid0chi2p3/output/'] + + +def save_cleaned_chains(fdir): + outf = fdir+'chain_info.txt' + reader = emcee.backends.HDFBackend(fdir+'param_chains.npz.h5') + x = np.load(fdir+'param_chains.npz') + labels = ['$%s$' %alllabels[k] for k in x['names']] + + tau = reader.get_autocorr_time(tol=0) + burnin = int(10. * np.mean(tau)) + thin = int(0.5 * np.mean(tau)) + + samples = reader.get_chain(discard=burnin, flat=True, thin=thin) + chains = reader.get_chain(discard=burnin, flat=False) + + N = chains.shape[0] + M = chains.shape[1] + chain_mean = np.mean(chains, axis=0) + chain_var = np.var(chains, axis=0) + samp_mean = np.mean(chains, axis=(0,1)) + B = N / (M-1) * np.sum( (chain_mean-samp_mean)**2, axis=0 ) + W = (1./M) * np.sum(chain_var, axis=0) + Vbar = (N-1)/N * W + (M+1)/(M*N) * B + GR = Vbar/W + + with open(outf, 'w') as of: + inds = int((N/np.mean(tau))) + print("chains: ", chains.shape, file=of) + print("tau: ", np.mean(tau), np.min(tau), np.max(tau), np.std(tau), file=of) + print("burnin: %d, thin: %d" %(burnin, thin), file=of) + print("independent samps per chain: %d" %inds, file=of) + print("GR", file=of) + print(GR, file=of) + if np.any(GR > 1.1): + print("FAILED GR", file=of) + if inds < 50: + print("POTENTIALLY BAD TAU", file=of) + np.savez(fdir+'cleaned_chains', samples=samples, names=x['names'], labels=labels) + return + +for af in allfs: + save_cleaned_chains(af) + diff --git a/tester b/tester deleted file mode 100644 index e69de29..0000000 From 67c7f3520e3edde5bd288a1b98a92b94716b8b0a Mon Sep 17 00:00:00 2001 From: Maximilian H Abitbol Date: Tue, 9 Feb 2021 22:43:57 +0000 Subject: [PATCH 07/20] cleaning and plotting --- plotting/plotter.py | 37 ++++++++++++++++ plotting/saveclean.py | 99 ++++++++++++++++++++++++++++++++----------- 2 files changed, 111 insertions(+), 25 deletions(-) create mode 100644 plotting/plotter.py diff --git a/plotting/plotter.py b/plotting/plotter.py new file mode 100644 index 0000000..c7d48d4 --- /dev/null +++ b/plotting/plotter.py @@ -0,0 +1,37 @@ +import numpy as np +from getdist import plots, MCSamples +from labels import alllabels, allranges +import matplotlib as mpl + + +allfs = ['/mnt/zfsusers/mabitbol/BBPower/baseline_noiseless/sim00/output0/'] + + +def getdist_clean(fdir): + sampler = np.load(fdir+'cleaned_chains.npz') + vals = [allranges[k] for k in sampler['labels']] + ranges = dict(zip(sampler['labels'], vals)) + gd_samps = MCSamples(samples=sampler['samples'], + names=sampler['labels'], + labels=[p.strip('$') for p in sampler['labels']], + ranges=ranges) + return gd_samps + +for af in allfs: + sname = af.split('/')[-3] + af.split('/')[-2] + outf = af+'parameter_info.txt' + + samps = getdist_clean(af) + samps.getTable(limit=1).write(outf) + + z = samps.paramNames.list() + g = plots.get_subplot_plotter(width_inch=16) + g.settings.axes_fontsize=14 + g.settings.axes_labelsize=14 + g.settings.line_styles = 'tab10' + g.triangle_plot([samps], z, shaded=True, title_limit=1) + mpl.pyplot.savefig('/mnt/zfsusers/mabitbol/notebooks/data_and_figures/'+sname+'_triangle') + + + + diff --git a/plotting/saveclean.py b/plotting/saveclean.py index 8680378..70bd623 100755 --- a/plotting/saveclean.py +++ b/plotting/saveclean.py @@ -3,7 +3,50 @@ from labels import alllabels import glob -allfs = ['/mnt/zfsusers/mabitbol/BBPower/moments/hybrid0chi2p3/output/'] +allfs = glob.glob('/mnt/zfsusers/mabitbol/BBPower/baseline_noiseless/sim0*/output*/') +allfs.sort() + +allfs = ['/mnt/zfsusers/mabitbol/BBPower/baseline_noiseless/sim00/output0/'] + +def next_pow_two(n): + i = 1 + while i < n: + i = i << 1 + return i + +def autocorr_func_1d(x, norm=True): + x = np.atleast_1d(x) + if len(x.shape) != 1: + raise ValueError("invalid dimensions for 1D autocorrelation function") + n = next_pow_two(len(x)) + + # Compute the FFT and then (from that) the auto-correlation function + f = np.fft.fft(x - np.mean(x), n=2 * n) + acf = np.fft.ifft(f * np.conjugate(f))[: len(x)].real + acf /= 4 * n + + # Optionally normalize + if norm: + acf /= acf[0] + + return acf + +# Automated windowing procedure following Sokal (1989) +def auto_window(taus, c): + m = np.arange(len(taus)) < c * taus + if np.any(m): + return np.argmin(m) + return len(taus) - 1 + +def autocorr_new(y, c=5.0): + f = np.zeros(y.shape[1]) + for yy in y: + f += autocorr_func_1d(yy) + f /= len(y) + taus = 2.0 * np.cumsum(f) - 1.0 + window = auto_window(taus, c) + return taus[window] + def save_cleaned_chains(fdir): @@ -12,34 +55,39 @@ def save_cleaned_chains(fdir): x = np.load(fdir+'param_chains.npz') labels = ['$%s$' %alllabels[k] for k in x['names']] - tau = reader.get_autocorr_time(tol=0) - burnin = int(10. * np.mean(tau)) - thin = int(0.5 * np.mean(tau)) + tau0 = reader.get_autocorr_time(tol=0) + burnin0 = int(5. * np.max(tau0)) + samps = reader.get_chain(discard=burnin0) + + nm = samps.shape[0] + nc = samps.shape[1] + npar = samps.shape[2] + + tau = np.zeros(npar) + for k in range(npar): + tau[k] = autocorr_new(samps[:, :, k].T) + + burnin = int(5. * np.max(tau)) + thin = int(0.5 * np.max(tau)) samples = reader.get_chain(discard=burnin, flat=True, thin=thin) - chains = reader.get_chain(discard=burnin, flat=False) - - N = chains.shape[0] - M = chains.shape[1] - chain_mean = np.mean(chains, axis=0) - chain_var = np.var(chains, axis=0) - samp_mean = np.mean(chains, axis=(0,1)) - B = N / (M-1) * np.sum( (chain_mean-samp_mean)**2, axis=0 ) - W = (1./M) * np.sum(chain_var, axis=0) - Vbar = (N-1)/N * W + (M+1)/(M*N) * B - GR = Vbar/W - + N = samples.shape[0] + with open(outf, 'w') as of: - inds = int((N/np.mean(tau))) - print("chains: ", chains.shape, file=of) - print("tau: ", np.mean(tau), np.min(tau), np.max(tau), np.std(tau), file=of) + inds = nm/tau + print("N: ", N, file=of) print("burnin: %d, thin: %d" %(burnin, thin), file=of) - print("independent samps per chain: %d" %inds, file=of) - print("GR", file=of) - print(GR, file=of) - if np.any(GR > 1.1): - print("FAILED GR", file=of) - if inds < 50: + print("mean tau: ", np.mean(tau), file=of) + print("mean number independent samps per chain: %d" %np.mean(inds), file=of) + print("percent Mean chain uncertainty: ", 100./np.sqrt(N), file=of) + print("\n", file=of) + print("number independent samps per chain:", file=of) + print(inds, file=of) + print("tau:", file=of) + print(tau, file=of) + print("percent parameter uncertainty: ", file=of) + print(100./np.sqrt(inds), file=of) + if np.any(inds < 50): print("POTENTIALLY BAD TAU", file=of) np.savez(fdir+'cleaned_chains', samples=samples, names=x['names'], labels=labels) return @@ -47,3 +95,4 @@ def save_cleaned_chains(fdir): for af in allfs: save_cleaned_chains(af) + From 238e0f0610db396e843a4cf4501159ae384bcb61 Mon Sep 17 00:00:00 2001 From: Maximilian H Abitbol Date: Wed, 10 Feb 2021 00:22:23 +0000 Subject: [PATCH 08/20] plotting and cleaned up some trash in compsep --- bbpower/compsep.py | 27 ++++++++++++--------------- plotting/plotter.py | 10 +++++----- plotting/saveclean.py | 21 +++++++++++++-------- 3 files changed, 30 insertions(+), 28 deletions(-) diff --git a/bbpower/compsep.py b/bbpower/compsep.py index ad0841d..d4875db 100644 --- a/bbpower/compsep.py +++ b/bbpower/compsep.py @@ -368,10 +368,8 @@ def emcee_sampler(self): Sample the model with MCMC. """ import emcee - from multiprocessing import Pool fname_temp = self.get_output('param_chains')+'.h5' - backend = emcee.backends.HDFBackend(fname_temp) nwalkers = self.config['nwalkers'] @@ -381,19 +379,21 @@ def emcee_sampler(self): if not found_file: backend.reset(nwalkers,ndim) - pos = [self.params.p0 + 1.e-3*np.random.randn(ndim) for i in range(nwalkers)] + p0= self.minimizer() + pos = [p0 + 1.e-3*np.random.randn(ndim) for i in range(nwalkers)] nsteps_use = n_iters else: print("Restarting from previous run") pos = None nsteps_use = max(n_iters-len(backend.get_chain()), 0) - - with Pool() as pool: - sampler = emcee.EnsembleSampler(nwalkers, ndim, self.lnprob,backend=backend) - if nsteps_use > 0: - sampler.run_mcmc(pos, nsteps_use, store=True, progress=True); - return sampler + import time + start = time.time() + sampler = emcee.EnsembleSampler(nwalkers, ndim, self.lnprob, backend=backend) + if nsteps_use > 0: + sampler.run_mcmc(pos, nsteps_use, store=True, progress=True); + end = time.time() + return sampler, end-start def minimizer(self): """ @@ -435,19 +435,16 @@ def timing(self, n_eval=300): for i in range(n_eval): lik = self.lnprob(self.params.p0) end = time.time() - return end-start, (end-start)/n_eval def run(self): from shutil import copyfile - copyfile(self.get_input('config'), self.get_output('config_copy')) self.setup_compsep() if self.config.get('sampler')=='emcee': - sampler = self.emcee_sampler() + sampler, timing = self.emcee_sampler() np.savez(self.get_output('param_chains'), - chain=sampler.chain, - names=self.params.p_free_names) - print("Finished sampling") + names=self.params.p_free_names, time=timing) + print("Finished sampling: ", timing) elif self.config.get('sampler')=='fisher': fisher = self.fisher() cov = np.linalg.inv(fisher) diff --git a/plotting/plotter.py b/plotting/plotter.py index c7d48d4..29d637f 100644 --- a/plotting/plotter.py +++ b/plotting/plotter.py @@ -3,9 +3,11 @@ from labels import alllabels, allranges import matplotlib as mpl +#allfs = glob.glob('/mnt/zfsusers/mabitbol/BBPower/residual_noiseless/sim0*/output*/') +#allfs.sort() -allfs = ['/mnt/zfsusers/mabitbol/BBPower/baseline_noiseless/sim00/output0/'] - +allfs = glob.glob('/mnt/zfsusers/mabitbol/BBPower/baseline_noiseless/sim0*/output*/') +allfs.sort() def getdist_clean(fdir): sampler = np.load(fdir+'cleaned_chains.npz') @@ -19,7 +21,7 @@ def getdist_clean(fdir): for af in allfs: sname = af.split('/')[-3] + af.split('/')[-2] - outf = af+'parameter_info.txt' + outf = af+'results.txt' samps = getdist_clean(af) samps.getTable(limit=1).write(outf) @@ -33,5 +35,3 @@ def getdist_clean(fdir): mpl.pyplot.savefig('/mnt/zfsusers/mabitbol/notebooks/data_and_figures/'+sname+'_triangle') - - diff --git a/plotting/saveclean.py b/plotting/saveclean.py index 70bd623..8f856dd 100755 --- a/plotting/saveclean.py +++ b/plotting/saveclean.py @@ -1,12 +1,23 @@ +import os.path import numpy as np +import glob import emcee from labels import alllabels -import glob + +#allfs = glob.glob('/mnt/zfsusers/mabitbol/BBPower/residual_noiseless/sim0*/output*/') +#allfs.sort() allfs = glob.glob('/mnt/zfsusers/mabitbol/BBPower/baseline_noiseless/sim0*/output*/') allfs.sort() -allfs = ['/mnt/zfsusers/mabitbol/BBPower/baseline_noiseless/sim00/output0/'] +dones = [] +for af in allfs: + if os.path.isfile(af+'cleaned_chains.npz'): + dones.append(af) +for df in dones: + allfs.remove(df) +print(len(allfs)) + def next_pow_two(n): i = 1 @@ -20,18 +31,14 @@ def autocorr_func_1d(x, norm=True): raise ValueError("invalid dimensions for 1D autocorrelation function") n = next_pow_two(len(x)) - # Compute the FFT and then (from that) the auto-correlation function f = np.fft.fft(x - np.mean(x), n=2 * n) acf = np.fft.ifft(f * np.conjugate(f))[: len(x)].real acf /= 4 * n - # Optionally normalize if norm: acf /= acf[0] - return acf -# Automated windowing procedure following Sokal (1989) def auto_window(taus, c): m = np.arange(len(taus)) < c * taus if np.any(m): @@ -47,8 +54,6 @@ def autocorr_new(y, c=5.0): window = auto_window(taus, c) return taus[window] - - def save_cleaned_chains(fdir): outf = fdir+'chain_info.txt' reader = emcee.backends.HDFBackend(fdir+'param_chains.npz.h5') From 7509cd9cac86c8cef9525f6d50c5791e6de44264 Mon Sep 17 00:00:00 2001 From: Maximilian H Abitbol Date: Wed, 10 Mar 2021 15:31:08 +0000 Subject: [PATCH 09/20] autoq update --- autoqbaseline.sh | 31 +++++++++++++++---- autoqresiduals.sh | 17 +++++----- bbpower/compsep.py | 2 +- plotting/combined.py | 72 +++++++++++++++++++++++++++++++++++++++++++ plotting/plotter.py | 14 ++++----- plotting/saveclean.py | 71 +++++------------------------------------- 6 files changed, 121 insertions(+), 86 deletions(-) create mode 100755 plotting/combined.py diff --git a/autoqbaseline.sh b/autoqbaseline.sh index 22eebde..c176f6b 100755 --- a/autoqbaseline.sh +++ b/autoqbaseline.sh @@ -1,23 +1,42 @@ #!/bin/bash -mdir="baseline_noiseless/" +mdir="baseline_full/" +for k in {0..3} +do + for j in {0000..0020} + do + output=$mdir"sim0$k/output$j/" + mkdir -p $output + + datadir="/mnt/zfsusers/mabitbol/BBHybrid/data/sim0$k/" + coadd=$datadir"cls_coadd_baseline_$j.fits" + noise=$datadir"cls_noise_baseline_$j.fits" + fid=$datadir"cls_fid_baseline_$j.fits" + chain=$output"param_chains.npz" + configcopy=$output"config_copy.yml" + config=$mdir"config.yml" + + addqueue -q cmb -c "2 hours" -m 1 -s -n 1x8 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config + done +done +mdir="baseline_masked/" for k in {0..3} do - for j in {0..20} + for j in {0000..0020} do output=$mdir"sim0$k/output$j/" mkdir -p $output datadir="/mnt/zfsusers/mabitbol/BBHybrid/data/sim0$k/" - coadd=$datadir"cls_coadd_base$j.fits" - noise=$datadir"cls_noise_base$j.fits" - fid=$datadir"cls_fid_base$j.fits" + coadd=$datadir"cls_coadd_baseline_masked_$j.fits" + noise=$datadir"cls_noise_baseline_masked_$j.fits" + fid=$datadir"cls_fid_baseline_masked_$j.fits" chain=$output"param_chains.npz" configcopy=$output"config_copy.yml" config=$mdir"config.yml" - addqueue -q cmb -c "2 days" -m 1 -s -n 1x8 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config + addqueue -q cmb -c "2 hours" -m 1 -s -n 1x8 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config done done diff --git a/autoqresiduals.sh b/autoqresiduals.sh index 250f1ee..6d2260d 100755 --- a/autoqresiduals.sh +++ b/autoqresiduals.sh @@ -1,11 +1,10 @@ #!/bin/bash -mdir="residual_noiseless/" -cf="residual_noiseless/config.yml" - +mdir="moments/hybrid_masked_10p/" +cf=$mdir"config.yml" for k in {0..3} do - for j in {0..20} + for j in {0000..0020} do output=$mdir"sim0$k/output$j/" mkdir -p $output @@ -13,16 +12,16 @@ do cp $cf $output config=$output"config.yml" sed -i "s/sim0X/sim0$k/" $config - sed -i "s/paramsY/params$j/" $config + sed -i "s/paramsY/params_$j/" $config datadir="/mnt/zfsusers/mabitbol/BBHybrid/data/sim0$k/" - coadd=$datadir"cls_coadd$j.fits" - noise=$datadir"cls_noise$j.fits" - fid=$datadir"cls_fid$j.fits" + coadd=$datadir"cls_coadd_residual_masked_$j.fits" + noise=$datadir"cls_noise_residual_masked_$j.fits" + fid=$datadir"cls_fid_residual_masked_$j.fits" chain=$output"param_chains.npz" configcopy=$output"config_copy.yml" - addqueue -q berg -c "2 days" -m 1 -s -n 1x8 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config + addqueue -q cmb -c "4 hours" -m 1 -s -n 1x4 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config done done diff --git a/bbpower/compsep.py b/bbpower/compsep.py index d4875db..c45a997 100644 --- a/bbpower/compsep.py +++ b/bbpower/compsep.py @@ -379,7 +379,7 @@ def emcee_sampler(self): if not found_file: backend.reset(nwalkers,ndim) - p0= self.minimizer() + p0 = self.minimizer() pos = [p0 + 1.e-3*np.random.randn(ndim) for i in range(nwalkers)] nsteps_use = n_iters else: diff --git a/plotting/combined.py b/plotting/combined.py new file mode 100755 index 0000000..9eb790d --- /dev/null +++ b/plotting/combined.py @@ -0,0 +1,72 @@ +import os.path +import numpy as np +import glob +import emcee +from labels import alllabels + +import matplotlib as mpl +from getdist import plots, MCSamples +from labels import alllabels, allranges + +npf = 'hybrid_masked_10p' +allfs = glob.glob(f'/mnt/zfsusers/mabitbol/BBPower/{npf}/sim0*/output*/') +allfs.sort() + + +def save_cleaned_chains(fdir): + reader = emcee.backends.HDFBackend(f'{fdir}param_chains.npz.h5') + x = np.load(f'{fdir}param_chains.npz') + labels = ['$%s$' %alllabels[k] for k in x['names']] + + try: + tau = reader.get_autocorr_time(tol=50) + taust = 'good tau' + except: + taust = 'POTENTIALLY BAD TAU' + tau = reader.get_autocorr_time(tol=0) + burnin = int(5. * np.max(tau)) + thin = int(0.5 * np.max(tau)) + + samples = reader.get_chain(discard=burnin, flat=True, thin=thin) + N = samples.shape[0] + + with open(f'{fdir}chain_info.txt', 'w') as of: + print("N: ", N, file=of) + print("burnin: %d, thin: %d" %(burnin, thin), file=of) + print("tau:", file=of) + print(tau, file=of) + print(taust, file=of) + np.savez(f'{fdir}cleaned_chains', samples=samples, names=x['names'], labels=labels) + return + +for af in allfs: + save_cleaned_chains(af) + + +def getdist_clean(fdir): + sampler = np.load(f'{fdir}cleaned_chains.npz') + vals = [allranges[k] for k in sampler['labels']] + ranges = dict(zip(sampler['labels'], vals)) + gd_samps = MCSamples(samples=sampler['samples'], + names=sampler['labels'], + labels=[p.strip('$') for p in sampler['labels']], + ranges=ranges) + return gd_samps + +for af in allfs: + sname = af.split('/')[-4] + af.split('/')[-3] + af.split('/')[-2] + #sname = af.split('/')[-3] + outf = f'{af}results.txt' + + samps = getdist_clean(af) + samps.getTable(limit=1).write(outf) + + z = samps.paramNames.list() + g = plots.get_subplot_plotter(width_inch=16) + g.settings.axes_fontsize=14 + g.settings.axes_labelsize=14 + g.settings.line_styles = 'tab10' + g.triangle_plot([samps], z, shaded=True, title_limit=1) + mpl.pyplot.savefig(f'/mnt/zfsusers/mabitbol/notebooks/data_and_figures/{sname}_triangle') + mpl.pyplot.close() + diff --git a/plotting/plotter.py b/plotting/plotter.py index 29d637f..12dfa5a 100644 --- a/plotting/plotter.py +++ b/plotting/plotter.py @@ -1,12 +1,10 @@ import numpy as np +import matplotlib as mpl +import glob from getdist import plots, MCSamples from labels import alllabels, allranges -import matplotlib as mpl - -#allfs = glob.glob('/mnt/zfsusers/mabitbol/BBPower/residual_noiseless/sim0*/output*/') -#allfs.sort() -allfs = glob.glob('/mnt/zfsusers/mabitbol/BBPower/baseline_noiseless/sim0*/output*/') +allfs = glob.glob('/mnt/zfsusers/mabitbol/BBPower/baseline_full/sim0*/output*/') allfs.sort() def getdist_clean(fdir): @@ -20,7 +18,8 @@ def getdist_clean(fdir): return gd_samps for af in allfs: - sname = af.split('/')[-3] + af.split('/')[-2] + sname = af.split('/')[-4] + af.split('/')[-3] + af.split('/')[-2] + #sname = af.split('/')[-3] outf = af+'results.txt' samps = getdist_clean(af) @@ -32,6 +31,7 @@ def getdist_clean(fdir): g.settings.axes_labelsize=14 g.settings.line_styles = 'tab10' g.triangle_plot([samps], z, shaded=True, title_limit=1) - mpl.pyplot.savefig('/mnt/zfsusers/mabitbol/notebooks/data_and_figures/'+sname+'_triangle') + mpl.pyplot.savefig(f'/mnt/zfsusers/mabitbol/notebooks/data_and_figures/{sname}_triangle') + mpl.pyplot.close() diff --git a/plotting/saveclean.py b/plotting/saveclean.py index 8f856dd..e1fae51 100755 --- a/plotting/saveclean.py +++ b/plotting/saveclean.py @@ -4,10 +4,7 @@ import emcee from labels import alllabels -#allfs = glob.glob('/mnt/zfsusers/mabitbol/BBPower/residual_noiseless/sim0*/output*/') -#allfs.sort() - -allfs = glob.glob('/mnt/zfsusers/mabitbol/BBPower/baseline_noiseless/sim0*/output*/') +allfs = glob.glob('/mnt/zfsusers/mabitbol/BBPower/baseline_full/sim0*/output*/') allfs.sort() dones = [] @@ -16,43 +13,7 @@ dones.append(af) for df in dones: allfs.remove(df) -print(len(allfs)) - - -def next_pow_two(n): - i = 1 - while i < n: - i = i << 1 - return i - -def autocorr_func_1d(x, norm=True): - x = np.atleast_1d(x) - if len(x.shape) != 1: - raise ValueError("invalid dimensions for 1D autocorrelation function") - n = next_pow_two(len(x)) - - f = np.fft.fft(x - np.mean(x), n=2 * n) - acf = np.fft.ifft(f * np.conjugate(f))[: len(x)].real - acf /= 4 * n - if norm: - acf /= acf[0] - return acf - -def auto_window(taus, c): - m = np.arange(len(taus)) < c * taus - if np.any(m): - return np.argmin(m) - return len(taus) - 1 - -def autocorr_new(y, c=5.0): - f = np.zeros(y.shape[1]) - for yy in y: - f += autocorr_func_1d(yy) - f /= len(y) - taus = 2.0 * np.cumsum(f) - 1.0 - window = auto_window(taus, c) - return taus[window] def save_cleaned_chains(fdir): outf = fdir+'chain_info.txt' @@ -60,18 +21,12 @@ def save_cleaned_chains(fdir): x = np.load(fdir+'param_chains.npz') labels = ['$%s$' %alllabels[k] for k in x['names']] - tau0 = reader.get_autocorr_time(tol=0) - burnin0 = int(5. * np.max(tau0)) - samps = reader.get_chain(discard=burnin0) - - nm = samps.shape[0] - nc = samps.shape[1] - npar = samps.shape[2] - - tau = np.zeros(npar) - for k in range(npar): - tau[k] = autocorr_new(samps[:, :, k].T) - + try: + tau = reader.get_autocorr_time(tol=50) + taust = 'good tau' + except: + taust = 'POTENTIALLY BAD TAU' + tau = reader.get_autocorr_time(tol=0) burnin = int(5. * np.max(tau)) thin = int(0.5 * np.max(tau)) @@ -79,21 +34,11 @@ def save_cleaned_chains(fdir): N = samples.shape[0] with open(outf, 'w') as of: - inds = nm/tau print("N: ", N, file=of) print("burnin: %d, thin: %d" %(burnin, thin), file=of) - print("mean tau: ", np.mean(tau), file=of) - print("mean number independent samps per chain: %d" %np.mean(inds), file=of) - print("percent Mean chain uncertainty: ", 100./np.sqrt(N), file=of) - print("\n", file=of) - print("number independent samps per chain:", file=of) - print(inds, file=of) print("tau:", file=of) print(tau, file=of) - print("percent parameter uncertainty: ", file=of) - print(100./np.sqrt(inds), file=of) - if np.any(inds < 50): - print("POTENTIALLY BAD TAU", file=of) + print(taust, file=of) np.savez(fdir+'cleaned_chains', samples=samples, names=x['names'], labels=labels) return From 95d96f22cbc53abf083c8ab6d0e26e704daa3f51 Mon Sep 17 00:00:00 2001 From: Maximilian H Abitbol Date: Tue, 22 Jun 2021 16:53:31 +0100 Subject: [PATCH 10/20] update --- autoqresiduals.sh | 4 ++-- plotting/combined.py | 11 +++++++---- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/autoqresiduals.sh b/autoqresiduals.sh index 6d2260d..0bf101a 100755 --- a/autoqresiduals.sh +++ b/autoqresiduals.sh @@ -1,6 +1,6 @@ #!/bin/bash -mdir="moments/hybrid_masked_10p/" +mdir="moments/hybrid_masked_alphap/" cf=$mdir"config.yml" for k in {0..3} do @@ -21,7 +21,7 @@ do chain=$output"param_chains.npz" configcopy=$output"config_copy.yml" - addqueue -q cmb -c "4 hours" -m 1 -s -n 1x4 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config + addqueue -q redwood -c "4 hours" -m 1 -s -n 1x1 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config done done diff --git a/plotting/combined.py b/plotting/combined.py index 9eb790d..99a98b0 100755 --- a/plotting/combined.py +++ b/plotting/combined.py @@ -8,8 +8,8 @@ from getdist import plots, MCSamples from labels import alllabels, allranges -npf = 'hybrid_masked_10p' -allfs = glob.glob(f'/mnt/zfsusers/mabitbol/BBPower/{npf}/sim0*/output*/') +npf = 'hybrid_masked_alphap' +allfs = glob.glob(f'/mnt/zfsusers/mabitbol/BBPower/moments/{npf}/sim0*/output*/') allfs.sort() @@ -47,10 +47,13 @@ def getdist_clean(fdir): sampler = np.load(f'{fdir}cleaned_chains.npz') vals = [allranges[k] for k in sampler['labels']] ranges = dict(zip(sampler['labels'], vals)) + #gd_samps = MCSamples(samples=sampler['samples'], + # names=sampler['labels'], + # labels=[p.strip('$') for p in sampler['labels']], + # ranges=ranges) gd_samps = MCSamples(samples=sampler['samples'], names=sampler['labels'], - labels=[p.strip('$') for p in sampler['labels']], - ranges=ranges) + labels=[p.strip('$') for p in sampler['labels']]) return gd_samps for af in allfs: From 30c6ca0a1f03cea804cfbc4a5f6e95928a38995b Mon Sep 17 00:00:00 2001 From: Maximilian H Abitbol Date: Tue, 22 Jun 2021 16:59:35 +0100 Subject: [PATCH 11/20] config files --- moments/config.yml | 52 ++++++++++++++++++++++++++++++++++++++++++++ moments/rdir.sh | 11 ++++++++++ moments/settings.yml | 17 +++++++++++++++ 3 files changed, 80 insertions(+) create mode 100644 moments/config.yml create mode 100755 moments/rdir.sh create mode 100644 moments/settings.yml diff --git a/moments/config.yml b/moments/config.yml new file mode 100644 index 0000000..832ddd5 --- /dev/null +++ b/moments/config.yml @@ -0,0 +1,52 @@ +global: + nside: 256 + compute_dell: True + +BBCompSep: + sampler: 'emcee' + nwalkers: 128 + n_iters: 10000 + likelihood_type: 'chi2' + pol_channels: ['B'] + l_min: 30 + l_max: 300 + diff: False + + cmb_model: + cmb_templates: ["./examples/data/camb_lens_nobb.dat", + "./examples/data/camb_lens_r1.dat"] + params: + r_tensor: ['r_tensor', 'tophat', [-0.1, 0.00, 0.1]] + A_lens: ['A_lens', 'tophat', [0.00, 1.0, 2.00]] + + resid_seds: "/mnt/zfsusers/mabitbol/BBHybrid/data/sim0X/hybrid_params0.npz" + fg_model: + component_1: + name: Dust + sed: Dust + cl: + BB: ClPowerLaw + sed_parameters: + beta_d: ['beta_d', 'Gaussian', [1.60, 0.1]] + temp_d: ['temp', 'fixed', [19.6]] + nu0_d: ['nu0', 'fixed', [353.]] + cl_parameters: + BB: + amp_d_bb: ['amp', 'tophat', [0., 1., 10.]] + alpha_d_bb: ['alpha', 'tophat', [-3., -0.5, 1.]] + l0_d_bb: ['ell0', 'fixed', [80.]] + + component_2: + name: Synchrotron + sed: Synchrotron + cl: + BB: ClPowerLaw + sed_parameters: + beta_s: ['beta_pl', 'Gaussian', [-3.0, 0.3]] + nu0_s: ['nu0', 'fixed', [23.]] + cl_parameters: + BB: + amp_s_bb: ['amp', 'tophat', [0., 1., 10.]] + alpha_s_bb: ['alpha', 'tophat', [-3., -0.5, 1.]] + l0_s_bb: ['ell0', 'fixed', [80.]] + diff --git a/moments/rdir.sh b/moments/rdir.sh new file mode 100755 index 0000000..33b5f61 --- /dev/null +++ b/moments/rdir.sh @@ -0,0 +1,11 @@ +#!/bin/bash + +sname=$1"/" +fdir="/mnt/zfsusers/mabitbol/BBPower/moments/" +ndir=$fdir$sname +mkdir $ndir +cp config.yml $ndir +mkdir $ndir"output/" +sed "s/xxx/$sname" settings.yml > $ndir"settings.yml" +#sed "s/xxx/$sname" $fdir"settings.yml" > $ndir"settings.yml" + diff --git a/moments/settings.yml b/moments/settings.yml new file mode 100644 index 0000000..98d6f98 --- /dev/null +++ b/moments/settings.yml @@ -0,0 +1,17 @@ +modules: bbpower +launcher: local +stages: + - name: BBCompSep + nprocess: 1 + +inputs: + cells_coadded: /mnt/zfsusers/mabitbol/BBHybrid/data/sim02/cls_coadd0.fits + cells_fiducial: /mnt/zfsusers/mabitbol/BBHybrid/data/sim02/cls_fid0.fits + cells_noise: /mnt/zfsusers/mabitbol/BBHybrid/data/sim02/cls_noise0.fits + +resume: False +config: ./moments/xxx/config.yml +output_dir: ./moments/xxx/output +log_dir: ./moments/xxx/output +pipeline_log: ./moments/xxx/output/log.txt + From fbfb8b998c1c2084d0972976793ac81b0f5cf0d2 Mon Sep 17 00:00:00 2001 From: Maximilian H Abitbol Date: Tue, 22 Jun 2021 17:05:18 +0100 Subject: [PATCH 12/20] config --- moments/config.yml | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/moments/config.yml b/moments/config.yml index 832ddd5..0b394c9 100644 --- a/moments/config.yml +++ b/moments/config.yml @@ -10,7 +10,7 @@ BBCompSep: pol_channels: ['B'] l_min: 30 l_max: 300 - diff: False + diff: True cmb_model: cmb_templates: ["./examples/data/camb_lens_nobb.dat", @@ -19,7 +19,7 @@ BBCompSep: r_tensor: ['r_tensor', 'tophat', [-0.1, 0.00, 0.1]] A_lens: ['A_lens', 'tophat', [0.00, 1.0, 2.00]] - resid_seds: "/mnt/zfsusers/mabitbol/BBHybrid/data/sim0X/hybrid_params0.npz" + resid_seds: "/mnt/zfsusers/mabitbol/BBHybrid/data/sim0X/masked_hybrid_paramsY.npz" fg_model: component_1: name: Dust @@ -27,12 +27,12 @@ BBCompSep: cl: BB: ClPowerLaw sed_parameters: - beta_d: ['beta_d', 'Gaussian', [1.60, 0.1]] + beta_d: ['beta_d', 'Gaussian', [1.60, 0.016]] temp_d: ['temp', 'fixed', [19.6]] nu0_d: ['nu0', 'fixed', [353.]] cl_parameters: BB: - amp_d_bb: ['amp', 'tophat', [0., 1., 10.]] + amp_d_bb: ['amp', 'tophat', [-10., 1., 10.]] alpha_d_bb: ['alpha', 'tophat', [-3., -0.5, 1.]] l0_d_bb: ['ell0', 'fixed', [80.]] @@ -42,11 +42,11 @@ BBCompSep: cl: BB: ClPowerLaw sed_parameters: - beta_s: ['beta_pl', 'Gaussian', [-3.0, 0.3]] + beta_s: ['beta_pl', 'Gaussian', [-3.0, 0.03]] nu0_s: ['nu0', 'fixed', [23.]] cl_parameters: BB: - amp_s_bb: ['amp', 'tophat', [0., 1., 10.]] + amp_s_bb: ['amp', 'tophat', [-10., 1., 10.]] alpha_s_bb: ['alpha', 'tophat', [-3., -0.5, 1.]] l0_s_bb: ['ell0', 'fixed', [80.]] From a5c413e6a65ff5ccbde8f8207b86314afe3a6371 Mon Sep 17 00:00:00 2001 From: susannaaz Date: Thu, 2 Dec 2021 12:14:52 +0000 Subject: [PATCH 13/20] off-diag coadd correction --- autoqbaseline.sh | 54 ++++++++++++++++++++++--------------- autoqresiduals.sh | 16 ++++++----- bbpower/__init__.py | 2 ++ bbpower/compsep.py | 42 +++++++++++++++++++++++++++-- bbpower/power_summarizer.py | 4 ++- plotting/combined.py | 16 ++++++++--- plotting/labels.py | 6 ++--- 7 files changed, 103 insertions(+), 37 deletions(-) diff --git a/autoqbaseline.sh b/autoqbaseline.sh index c176f6b..820e00c 100755 --- a/autoqbaseline.sh +++ b/autoqbaseline.sh @@ -1,34 +1,44 @@ #!/bin/bash -mdir="baseline_full/" -for k in {0..3} -do - for j in {0000..0020} - do - output=$mdir"sim0$k/output$j/" - mkdir -p $output - - datadir="/mnt/zfsusers/mabitbol/BBHybrid/data/sim0$k/" - coadd=$datadir"cls_coadd_baseline_$j.fits" - noise=$datadir"cls_noise_baseline_$j.fits" - fid=$datadir"cls_fid_baseline_$j.fits" - chain=$output"param_chains.npz" - configcopy=$output"config_copy.yml" - config=$mdir"config.yml" - - addqueue -q cmb -c "2 hours" -m 1 -s -n 1x8 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config - done -done +##mdir="baseline_full/" +##for k in {0..3} +##do +## for j in 0000 #{0000..0020} +## do +## output=$mdir"sim0$k/output$j/" +## mkdir -p $output +## cp config.yml $output +## cp config.yml $mdir +## +## datadir="/mnt/zfsusers/susanna/BBHybrid/data/sim0$k/" +## coadd=$datadir"cls_coadd_baseline_$j.fits" +## noise=$datadir"cls_noise_baseline_$j.fits" +## fid=$datadir"cls_fid_baseline_$j.fits" +## chain=$output"param_chains.npz" +## configcopy=$output"config_copy.yml" +## config=$mdir"config.yml" +## +## addqueue -q berg -c "2 hours" -m 1 -s -n 1x8 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config +## done +##done mdir="baseline_masked/" for k in {0..3} do - for j in {0000..0020} + for j in 0000 #{0000..0020} do output=$mdir"sim0$k/output$j/" mkdir -p $output - datadir="/mnt/zfsusers/mabitbol/BBHybrid/data/sim0$k/" + #cp test_hybrid/config.yml $output + #cp test_hybrid/config.yml $mdir + #cp baseline_masked/config_widerprior.yml $output + #cp baseline_masked/config.yml $output + cp test_BBSims_Hybrid/config.yaml ${mdir}config.yml + cp test_BBSims_Hybrid/config.yaml ${output}config.yml + + datadir="/mnt/zfsusers/susanna/BBHybrid/data/sim0$k/" + #datadir="/mnt/zfsusers/susanna/BBHybrid/data/sims_DA/sim0$k/" coadd=$datadir"cls_coadd_baseline_masked_$j.fits" noise=$datadir"cls_noise_baseline_masked_$j.fits" fid=$datadir"cls_fid_baseline_masked_$j.fits" @@ -36,7 +46,7 @@ do configcopy=$output"config_copy.yml" config=$mdir"config.yml" - addqueue -q cmb -c "2 hours" -m 1 -s -n 1x8 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config + addqueue -s -q cmb -c '2 hours' -m 4 -s -n 1x12 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config done done diff --git a/autoqresiduals.sh b/autoqresiduals.sh index 0bf101a..cb62bbf 100755 --- a/autoqresiduals.sh +++ b/autoqresiduals.sh @@ -1,27 +1,31 @@ #!/bin/bash -mdir="moments/hybrid_masked_alphap/" +mdir="test_hybrid/hybrid_masked_alphap/" cf=$mdir"config.yml" for k in {0..3} do - for j in {0000..0020} + for j in 0000 #{0000..0020} do - output=$mdir"sim0$k/output$j/" + #output=$mdir"sim0$k/output$j/" + output=$mdir"sim0$k/output$j/" mkdir -p $output + cp test_BBSims_Hybrid/config_hybrid.yaml ${mdir}config.yml + #cp test_BBSims_Hybrid/config_hybrid.yaml ${output}config.yml + cp $cf $output config=$output"config.yml" sed -i "s/sim0X/sim0$k/" $config sed -i "s/paramsY/params_$j/" $config - datadir="/mnt/zfsusers/mabitbol/BBHybrid/data/sim0$k/" + datadir="/mnt/zfsusers/susanna/BBHybrid/data/sim0$k/" + #datadir="/mnt/zfsusers/susanna/BBHybrid/data/sim00_noiselessMA/" coadd=$datadir"cls_coadd_residual_masked_$j.fits" noise=$datadir"cls_noise_residual_masked_$j.fits" fid=$datadir"cls_fid_residual_masked_$j.fits" chain=$output"param_chains.npz" configcopy=$output"config_copy.yml" - addqueue -q redwood -c "4 hours" -m 1 -s -n 1x1 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config + addqueue -s -q berg -c '2 hours' -m 4 -s -n 1x12 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config done done - diff --git a/bbpower/__init__.py b/bbpower/__init__.py index 3d118fb..4ae7ebc 100644 --- a/bbpower/__init__.py +++ b/bbpower/__init__.py @@ -1,3 +1,5 @@ from bbpipe import PipelineStage +from .power_specter import BBPowerSpecter +from .power_summarizer import BBPowerSummarizer from .compsep import BBCompSep from .plotter import BBPlotter diff --git a/bbpower/compsep.py b/bbpower/compsep.py index c45a997..a5cd8d6 100644 --- a/bbpower/compsep.py +++ b/bbpower/compsep.py @@ -351,17 +351,41 @@ def lnprob(self, par): """ Likelihood with priors. """ + ##prior = self.params.lnprior(par) + ##if not np.isfinite(prior): + ## return -np.inf + ## + ##params = self.params.build_params(par) + ##if self.use_handl: + ## dx = self.h_and_l_dx(params) + ##else: + ## dx = self.chi_sq_dx(params) + ##like = -0.5 * np.einsum('i, ij, j',dx,self.invcov,dx) + ##return prior + like prior = self.params.lnprior(par) if not np.isfinite(prior): return -np.inf + return prior + self.lnlike(par) + + def lnlike(self, par): + """ + Likelihood without priors. + """ params = self.params.build_params(par) if self.use_handl: dx = self.h_and_l_dx(params) + if np.any(np.isinf(dx)): + return -np.inf else: dx = self.chi_sq_dx(params) - like = -0.5 * np.einsum('i, ij, j',dx,self.invcov,dx) - return prior + like + # invcov = 567, 567 + # dx = 567 + #print(self.invcov.shape) + #print(dx.shape) + like = -0.5 * np.dot(dx, np.dot(self.invcov, dx)) + return like + def emcee_sampler(self): """ @@ -404,6 +428,20 @@ def chi2(par): c2=-2*self.lnprob(par) return c2 res=minimize(chi2, self.params.p0, method="Powell") + ## Plot Cls data vs model + ##bbdata_cls = self.bbdata # n_bpws, nfreqs, nfreqs + ##names=self.params.p_free_names + ##parameters = dict(zip(list(names),res.x)) + ##model_cls = self.model(parameters) # n_bpws, nfreqs, nfreqs + ##for i in range(self.nfreqs): + ## for j in range(i, self.nfreqs): + ## import matplotlib.pyplot as plt + ## plt.figure() + ## plt.plot(self.ell_b, bbdata_cls[:,i,j], label='data') + ## plt.plot(self.ell_b, model_cls[:,i,j], label='theory') + ## plt.legend() + ## plt.savefig(f'model_vs_data_{i}_{j}.png',bbox_inches='tight') + return res.x def fisher(self): diff --git a/bbpower/power_summarizer.py b/bbpower/power_summarizer.py index 18d9a23..0fa249a 100644 --- a/bbpower/power_summarizer.py +++ b/bbpower/power_summarizer.py @@ -266,7 +266,9 @@ def parse_splits_sacc_file(self, s, get_saccs=False, with_windows=False): weights_total) # Off-diagonal coadding - spectra_xcorr = np.mean(spectra[np.triu_indices(self.nsplits,1)],axis=0) + triu_mean = np.mean(spectra[np.triu_indices(self.nsplits, 1)], axis=0) + tril_mean = np.mean(spectra[np.tril_indices(self.nsplits, -1)], axis=0) + spectra_xcorr = 0.5*(tril_mean+triu_mean) # Noise power spectra spectra_noise = spectra_total - spectra_xcorr diff --git a/plotting/combined.py b/plotting/combined.py index 99a98b0..6ee94bf 100755 --- a/plotting/combined.py +++ b/plotting/combined.py @@ -8,8 +8,15 @@ from getdist import plots, MCSamples from labels import alllabels, allranges -npf = 'hybrid_masked_alphap' -allfs = glob.glob(f'/mnt/zfsusers/mabitbol/BBPower/moments/{npf}/sim0*/output*/') +hybrid = 0 +baseline = 1 + +if hybrid: + npf = 'hybrid_masked_alphap' + allfs = glob.glob(f'/mnt/zfsusers/susanna/BBPower/test_hybrid/{npf}/sim0*/output*/') + #/mnt/zfsusers/susanna/BBHybrid/data/sim00/ +elif baseline: + allfs = glob.glob('/mnt/zfsusers/susanna/BBPower/baseline_masked/sim0*/output*/') allfs.sort() @@ -70,6 +77,9 @@ def getdist_clean(fdir): g.settings.axes_labelsize=14 g.settings.line_styles = 'tab10' g.triangle_plot([samps], z, shaded=True, title_limit=1) - mpl.pyplot.savefig(f'/mnt/zfsusers/mabitbol/notebooks/data_and_figures/{sname}_triangle') + if hybrid: + mpl.pyplot.savefig(f'/mnt/zfsusers/susanna/BBPower/test_hybrid/plots/{sname}_triangle') + elif baseline: + mpl.pyplot.savefig(f'/mnt/zfsusers/susanna/BBPower/baseline_masked/plots/{sname}_triangle') mpl.pyplot.close() diff --git a/plotting/labels.py b/plotting/labels.py index 18e8473..59e9915 100644 --- a/plotting/labels.py +++ b/plotting/labels.py @@ -20,8 +20,8 @@ '$\\beta_d$': [0.1, 10.0], '$A^{EE}_d$': [0.0, None], '$\\alpha^{EE}_d$': [-3.0, 1.0], - '$A^{BB}_d$': [0.0, None], - '$\\alpha^{BB}_d$': [-3.0, 1.0], + '$A^{BB}_d$': [-10.0, None], + '$\\alpha^{BB}_d$': [-3.0, 3.0], '$A^{EB}_d$': [-100.0, 100.0], '$\\alpha^{EB}_d$': [-3.0, 1.0], '$\\epsilon$': [-1.0, 1.0], @@ -29,7 +29,7 @@ '$\\beta_s$': [-10.0, 0.0], '$A^{EE}_s$': [0.0, None], '$\\alpha^{EE}_s$': [-3.0, 1.0], - '$A^{BB}_s$': [0.0, None], + '$A^{BB}_s$': [-10.0, None], '$\\alpha^{BB}_s$': [-3.0, 1.0], '$A^{EB}_s$': [-100.0, 100.0], '$\\alpha^{EB}_s$': [-3.0, 1.0], From 97f5895f4b2663a2f647ca9527c7cdcacd688702 Mon Sep 17 00:00:00 2001 From: susannaaz Date: Wed, 9 Mar 2022 11:55:25 +0000 Subject: [PATCH 14/20] Q to R modifications --- autoqbaseline.sh | 2 +- bbpower/compsep.py | 128 ++++++++++++++++++++++++------------------- bbpower/fg_model.py | 10 ++++ plotting/combined.py | 16 ++++-- 4 files changed, 94 insertions(+), 62 deletions(-) diff --git a/autoqbaseline.sh b/autoqbaseline.sh index 820e00c..1f531f7 100755 --- a/autoqbaseline.sh +++ b/autoqbaseline.sh @@ -23,7 +23,7 @@ ##done mdir="baseline_masked/" -for k in {0..3} +for k in 0 #{0..3} do for j in 0000 #{0000..0020} do diff --git a/bbpower/compsep.py b/bbpower/compsep.py index a5cd8d6..88d801d 100644 --- a/bbpower/compsep.py +++ b/bbpower/compsep.py @@ -32,12 +32,18 @@ def setup_compsep(self): self.params = ParameterManager(self.config) if self.config.get("diff"): self.hybridparams = np.load(self.config['resid_seds']) + hyb_params = self.hybridparams['params_cent'] + hyb_sigma = self.hybridparams['params_sigma'] + # Beta_s: [centre, sigma] + self.params.p_free_priors[6][2] = [hyb_params[0], hyb_sigma[0]] + # Beta_d: [centre, sigma] + self.params.p_free_priors[2][2] = [hyb_params[1], hyb_sigma[1]] if self.use_handl: self.prepare_h_and_l() return def matrix_to_vector(self, mat): - return mat[..., self.index_ut[0], self.index_ut[1]] + return mat[..., self.index_ut_modes[0], self.index_ut_modes[1]] def vector_to_matrix(self, vec): if vec.ndim == 1: @@ -167,12 +173,47 @@ def parse_sacc_file(self): cv2d[:, ind_vec, :, ind_vecb] = self.s.covariance.covmat[ind_a][:, ind_b] #Store data - self.bbdata = self.vector_to_matrix(v2d) + self.bbdata = self.vector_to_matrix(v2d) #(nbpw, nf, nf) + if self.config.get("diff"): + from scipy.linalg import eig, inv, pinv + self.hybridparams = np.load(self.config['resid_seds']) + Qmat = self.hybridparams['Q'] + eval, evec = eig(Qmat) + mask = eval > 0.1 #take away rows equivalent to eigenvalue=0 + self.R = np.linalg.inv(evec)[mask, :] + self.bbdata = np.einsum("aj,bk,ijk->iab", self.R, self.R, self.bbdata) + cov_all = np.zeros([self.n_bpws, self.nmaps, self.nmaps, self.n_bpws, self.nmaps, self.nmaps]) + for ix in range(self.ncross): + i, j = np.array(self.index_ut)[:,ix] + for jx in range(self.ncross): + m, n = np.array(self.index_ut)[:,jx] + cov_all[:, i, j, :, m, n] = cv2d[:, ix, :, jx] + cov_all[:, j, i, :, m, n] = cv2d[:, ix, :, jx] + cov_all[:, i, j, :, n, m] = cv2d[:, ix, :, jx] + cov_all[:, j, i, :, n, m] = cv2d[:, ix, :, jx] + ccov = np.einsum("aj,bk,gm,dn,ijklmn->iablgd", self.R, self.R, self.R, self.R, cov_all) + self.nmodes = 4 + self.nmaps_modes = self.nmodes * self.npol + self.index_ut_modes = np.triu_indices(self.nmodes) + self.ncross_modes = (self.nmaps_modes * (self.nmaps_modes + 1)) // 2 + cv2d = np.zeros([self.n_bpws, self.ncross_modes, self.n_bpws, self.ncross_modes]) + for ix in range(self.ncross_modes): + i, j = np.array(self.index_ut_modes)[:,ix] + for jx in range(self.ncross_modes): + m, n = np.array(self.index_ut_modes)[:,jx] + cv2d[:, ix, :, jx] = ccov[:, i, j, :, m, n] + else: + self.nmodes = self.nfreqs + self.nmaps_modes = self.nmaps + self.index_ut_modes = self.index_ut + self.ncross_modes = self.ncross + self.bbcovar = cv2d.reshape([self.n_bpws * self.ncross_modes, self.n_bpws * self.ncross_modes]) + if self.use_handl: self.bbnoise = self.vector_to_matrix(v2d_noi) self.bbfiducial = self.vector_to_matrix(v2d_fid) - self.bbcovar = cv2d.reshape([self.n_bpws * self.ncross, self.n_bpws * self.ncross]) - self.invcov = np.linalg.solve(self.bbcovar, np.identity(len(self.bbcovar))) + + self.invcov = np.linalg.inv(self.bbcovar) return def load_cmb(self): @@ -210,21 +251,18 @@ def integrate_seds(self, params): sed_params = [params[comp['names_sed_dict'][k]] for k in comp['sed'].params] rot_matrices.append([]) - + if self.config.get("diff"): def sed(nu): - return comp['sed'].diff(nu, *sed_params) + return comp['sed'].diff(nu, *sed_params) else: def sed(nu): return comp['sed'].eval(nu, *sed_params) - + for tn in range(self.nfreqs): sed_b, rot = self.bpss[tn].convolve_sed(sed, params) fg_scaling[i_c, tn] = sed_b * units rot_matrices[i_c].append(rot) - - if self.config.get("diff"): - fg_scaling[i_c] = np.dot(self.hybridparams['Q'], fg_scaling[i_c]) return fg_scaling.T, rot_matrices @@ -264,14 +302,14 @@ def model(self, params): self.cmb_scal) * self.dl2cl # [npol,npol,nell] fg_scaling, rot_m = self.integrate_seds(params) # [nfreq, ncomp], [ncomp,nfreq,[matrix]] fg_cell = self.evaluate_power_spectra(params) # [ncomp,ncomp,npol,npol,nell] - + cls_array_fg = np.zeros([self.nfreqs,self.nfreqs,self.n_ell,self.npol,self.npol]) fg_cell = np.transpose(fg_cell, axes = [0,1,4,2,3]) # [ncomp,ncomp,nell,npol,npol] cmb_cell = np.transpose(cmb_cell, axes = [2,0,1]) # [nell,npol,npol] for f1 in range(self.nfreqs): for f2 in range(f1,self.nfreqs): cls=cmb_cell.copy() - + for c1 in range(self.fg_model.n_components): mat1=rot_m[c1][f1] a1=fg_scaling[f1,c1] @@ -295,14 +333,16 @@ def model(self, params): cls_array_list[:, f1, p1, f2, p2] = clband if m1!=m2: cls_array_list[:, f2, p2, f1, p1] = clband - + for f1 in range(self.nfreqs): for f2 in range(self.nfreqs): cls_array_list[:,f1,:,f2,:] = rotate_cells(self.bpss[f2], self.bpss[f1], cls_array_list[:,f1,:,f2,:], params) - - return cls_array_list.reshape([self.n_bpws, self.nmaps, self.nmaps]) + cls_array_list = cls_array_list.reshape([self.n_bpws, self.nmaps, self.nmaps]) + if self.config.get("diff"): + cls_array_list = np.einsum("aj,bk,ijk->iab", self.R, self.R, cls_array_list) + return cls_array_list def chi_sq_dx(self, params): """ @@ -318,54 +358,51 @@ def prepare_h_and_l(self): return def h_and_l_dx(self, params): - """ - Hamimeche and Lewis likelihood. - Taken from Cobaya written by H, L and Torrado - See: https://github.com/CobayaSampler/cobaya/blob/master/cobaya/likelihoods/_cmblikes_prototype/_cmblikes_prototype.py + """ + Hamimeche and Lewis likelihood. + Taken from Cobaya written by H, L and Torrado + See: https://github.com/CobayaSampler/cobaya/blob/master/cobaya/likelihoods/_cmblikes_prototype/_cmblikes_prototype.py """ model_cls = self.model(params) dx_vec = [] for k in range(model_cls.shape[0]): C = model_cls[k] + self.bbnoise[k] X = self.h_and_l(C, self.observed_cls[k], self.Cfl_sqrt[k]) + if np.any(np.isinf(X)): + return [np.inf] dx = self.matrix_to_vector(X).flatten() dx_vec = np.concatenate([dx_vec, dx]) return dx_vec - + def h_and_l(self, C, Chat, Cfl_sqrt): - diag, U = np.linalg.eigh(C) + try: + diag, U = np.linalg.eigh(C) + except: + return [np.inf] rot = U.T.dot(Chat).dot(U) roots = np.sqrt(diag) for i, root in enumerate(roots): rot[i, :] /= root rot[:, i] /= root U.dot(rot.dot(U.T), rot) - diag, rot = np.linalg.eigh(rot) - diag = np.sign(diag - 1) * np.sqrt(2 * np.maximum(0, diag - np.log(diag) - 1)) + try: + diag, rot = np.linalg.eigh(rot) + except: + return [np.inf] + diag = (np.sign(diag - 1) * + np.sqrt(2 * np.maximum(0, diag - np.log(diag) - 1))) Cfl_sqrt.dot(rot, U) for i, d in enumerate(diag): rot[:, i] = U[:, i] * d return rot.dot(U.T) - + def lnprob(self, par): """ Likelihood with priors. """ - ##prior = self.params.lnprior(par) - ##if not np.isfinite(prior): - ## return -np.inf - ## - ##params = self.params.build_params(par) - ##if self.use_handl: - ## dx = self.h_and_l_dx(params) - ##else: - ## dx = self.chi_sq_dx(params) - ##like = -0.5 * np.einsum('i, ij, j',dx,self.invcov,dx) - ##return prior + like prior = self.params.lnprior(par) if not np.isfinite(prior): return -np.inf - return prior + self.lnlike(par) def lnlike(self, par): @@ -379,14 +416,9 @@ def lnlike(self, par): return -np.inf else: dx = self.chi_sq_dx(params) - # invcov = 567, 567 - # dx = 567 - #print(self.invcov.shape) - #print(dx.shape) like = -0.5 * np.dot(dx, np.dot(self.invcov, dx)) return like - def emcee_sampler(self): """ Sample the model with MCMC. @@ -428,20 +460,6 @@ def chi2(par): c2=-2*self.lnprob(par) return c2 res=minimize(chi2, self.params.p0, method="Powell") - ## Plot Cls data vs model - ##bbdata_cls = self.bbdata # n_bpws, nfreqs, nfreqs - ##names=self.params.p_free_names - ##parameters = dict(zip(list(names),res.x)) - ##model_cls = self.model(parameters) # n_bpws, nfreqs, nfreqs - ##for i in range(self.nfreqs): - ## for j in range(i, self.nfreqs): - ## import matplotlib.pyplot as plt - ## plt.figure() - ## plt.plot(self.ell_b, bbdata_cls[:,i,j], label='data') - ## plt.plot(self.ell_b, model_cls[:,i,j], label='theory') - ## plt.legend() - ## plt.savefig(f'model_vs_data_{i}_{j}.png',bbox_inches='tight') - return res.x def fisher(self): diff --git a/bbpower/fg_model.py b/bbpower/fg_model.py index 99e8fb6..1fae6db 100644 --- a/bbpower/fg_model.py +++ b/bbpower/fg_model.py @@ -85,6 +85,16 @@ def load_foregrounds(self, config): val = None params_fgl[k][l[0]]=val + ###Added------------------------------------------------------------------------ + ### Moment parameters + ##comp['names_moments_dict'] = {} + ##d = component.get('moments') + ##if d and config['fg_model'].get('use_moments'): + ## comp['moments_pameters'] = component['moments'] + ## for k, l in component['moments'].items(): + ## comp['names_moments_dict'][l[0]] = k + ###------------------------------------------------------------------------ + # Set Cl functions comp['cl'] = {} for k,c in component['cl'].items(): diff --git a/plotting/combined.py b/plotting/combined.py index 6ee94bf..978b384 100755 --- a/plotting/combined.py +++ b/plotting/combined.py @@ -8,13 +8,15 @@ from getdist import plots, MCSamples from labels import alllabels, allranges -hybrid = 0 -baseline = 1 +hybrid = 1 +baseline = 0 if hybrid: - npf = 'hybrid_masked_alphap' - allfs = glob.glob(f'/mnt/zfsusers/susanna/BBPower/test_hybrid/{npf}/sim0*/output*/') - #/mnt/zfsusers/susanna/BBHybrid/data/sim00/ + #npf = 'hybrid_masked_alphap' + #allfs = glob.glob(f'/mnt/zfsusers/susanna/BBPower/test_hybrid/{npf}/sim0*/output*/') + npf = 'hybrid_masked_outs' + #allfs = glob.glob(f'/mnt/extraspace/susanna/BBHybrid/{npf}/gaussian_priorsfromBfore/sim0*/output*/') + allfs = glob.glob(f'/mnt/extraspace/susanna/BBHybrid/{npf}/realistic/d1s1/output000_pinv/') elif baseline: allfs = glob.glob('/mnt/zfsusers/susanna/BBPower/baseline_masked/sim0*/output*/') allfs.sort() @@ -78,7 +80,9 @@ def getdist_clean(fdir): g.settings.line_styles = 'tab10' g.triangle_plot([samps], z, shaded=True, title_limit=1) if hybrid: - mpl.pyplot.savefig(f'/mnt/zfsusers/susanna/BBPower/test_hybrid/plots/{sname}_triangle') + #mpl.pyplot.savefig(f'/mnt/zfsusers/susanna/BBPower/test_hybrid/plots/{sname}_triangle') + #mpl.pyplot.savefig(f'/mnt/zfsusers/susanna/BBPower/plotting/{sname}_triangle_betasB4') + mpl.pyplot.savefig(f'/mnt/zfsusers/susanna/BBPower/plotting/{sname}_triangle_d1s1pinv') elif baseline: mpl.pyplot.savefig(f'/mnt/zfsusers/susanna/BBPower/baseline_masked/plots/{sname}_triangle') mpl.pyplot.close() From 874ead2d3879946f2a30ec03881d7103575a4e46 Mon Sep 17 00:00:00 2001 From: susannaaz Date: Tue, 23 Aug 2022 14:08:39 +0100 Subject: [PATCH 15/20] Latest hybrid compsep --- autoqbaseline.sh | 60 ++++---- autoqresiduals.sh | 68 ++++++--- bbpower/compsep.py | 293 +++++++++++++++++++++++++++++++----- bbpower/fg_model.py | 105 +++++++------ bbpower/param_manager.py | 42 ++++-- bbpower/power_specter.py | 282 +++++++++++++++++++--------------- bbpower/power_summarizer.py | 1 + plotting/labels.py | 10 +- 8 files changed, 588 insertions(+), 273 deletions(-) diff --git a/autoqbaseline.sh b/autoqbaseline.sh index 1f531f7..a65ed02 100755 --- a/autoqbaseline.sh +++ b/autoqbaseline.sh @@ -21,32 +21,40 @@ ## addqueue -q berg -c "2 hours" -m 1 -s -n 1x8 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config ## done ##done - -mdir="baseline_masked/" -for k in 0 #{0..3} +ellmin=30 #30 #2 #10 +#simt=d1s1_maskpysm_Bonly +#simt=d1s1_maskpysm_Bonly_r0.01 +#simt=d1s1_maskpysm_Bonly_r0.01_whitenoiONLY +simt=d1s1_maskpysm_Bonly_whitenoiONLY +fsky=0.3 #0.6 #0.8 +#mdir="/mnt/extraspace/susanna/BBHybrid/baseline_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}/" +#mdir="/mnt/extraspace/susanna/BBHybrid/baseline_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_r0.01/" +#mdir="/mnt/extraspace/susanna/BBHybrid/baseline_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_r0.01_whitenoiONLY/" +mdir="/mnt/extraspace/susanna/BBHybrid/baseline_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_whitenoiONLY/" +cf=$mdir"config.yml" +for j in {0000..0020} #0000 do - for j in 0000 #{0000..0020} - do - output=$mdir"sim0$k/output$j/" - mkdir -p $output - - #cp test_hybrid/config.yml $output - #cp test_hybrid/config.yml $mdir - #cp baseline_masked/config_widerprior.yml $output - #cp baseline_masked/config.yml $output - cp test_BBSims_Hybrid/config.yaml ${mdir}config.yml - cp test_BBSims_Hybrid/config.yaml ${output}config.yml - - datadir="/mnt/zfsusers/susanna/BBHybrid/data/sim0$k/" - #datadir="/mnt/zfsusers/susanna/BBHybrid/data/sims_DA/sim0$k/" - coadd=$datadir"cls_coadd_baseline_masked_$j.fits" - noise=$datadir"cls_noise_baseline_masked_$j.fits" - fid=$datadir"cls_fid_baseline_masked_$j.fits" - chain=$output"param_chains.npz" - configcopy=$output"config_copy.yml" - config=$mdir"config.yml" - - addqueue -s -q cmb -c '2 hours' -m 4 -s -n 1x12 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config - done + output=$mdir"output$j/" + mkdir -p $output + + #cp test_BBSims_Hybrid/config_baseline_d1s1.yaml ${mdir}config.yml + cp test_BBSims_Hybrid/config_baseline_fsky.yaml ${mdir}config.yml + + cp $cf $output + config=$output"config.yml" + sed -i "s/KK/${ellmin}/" $config + + echo $config + + datadir="/mnt/zfsusers/susanna/BBHybrid/data/sim${simt}/fsky${fsky}/" + coadd=$datadir"cls_coadd_baseline_masked_$j.fits" + noise=$datadir"cls_noise_baseline_masked_$j.fits" + fid=$datadir"cls_fid_baseline_masked_$j.fits" + chain=$output"param_chains.npz" + configcopy=$output"config_copy.yml" + + addqueue -s -q cmb -c 'baseline fsky0.3 1 hour' -m 4 -s -n 1x12 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config + done + diff --git a/autoqresiduals.sh b/autoqresiduals.sh index cb62bbf..be1d22b 100755 --- a/autoqresiduals.sh +++ b/autoqresiduals.sh @@ -1,31 +1,49 @@ #!/bin/bash -mdir="test_hybrid/hybrid_masked_alphap/" -cf=$mdir"config.yml" -for k in {0..3} -do - for j in 0000 #{0000..0020} - do - #output=$mdir"sim0$k/output$j/" - output=$mdir"sim0$k/output$j/" - mkdir -p $output +ellmin=30 #2 #10 #30 +#simt=d1s1_maskpysm_Bonly +#simt=d1s1_maskpysm_Bonly_r0.01 +#simt=d1s1_maskpysm_Bonly_r0.01_whitenoiONLY +simt=d1s1_maskpysm_Bonly_whitenoiONLY + +fsky=0.3 #0.3 0.6 #0.8 +zz=masked_ #full_ - cp test_BBSims_Hybrid/config_hybrid.yaml ${mdir}config.yml - #cp test_BBSims_Hybrid/config_hybrid.yaml ${output}config.yml - - cp $cf $output - config=$output"config.yml" - sed -i "s/sim0X/sim0$k/" $config - sed -i "s/paramsY/params_$j/" $config +#mdir="/mnt/extraspace/susanna/BBHybrid/hybrid_masked_outs/mask_planck/realistic/d1s1/fsky${fsky}/" +#mdir="/mnt/extraspace/susanna/BBHybrid/hybrid_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_epsilonfixed/" +#mdir="/mnt/extraspace/susanna/BBHybrid/hybrid_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}/" +#mdir="/mnt/extraspace/susanna/BBHybrid/hybrid_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_r0.01/" +#mdir="/mnt/extraspace/susanna/BBHybrid/hybrid_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_r0.01_whitenoiONLY/" +mdir="/mnt/extraspace/susanna/BBHybrid/hybrid_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_whitenoiONLY/" +cf=$mdir"config.yml" - datadir="/mnt/zfsusers/susanna/BBHybrid/data/sim0$k/" - #datadir="/mnt/zfsusers/susanna/BBHybrid/data/sim00_noiselessMA/" - coadd=$datadir"cls_coadd_residual_masked_$j.fits" - noise=$datadir"cls_noise_residual_masked_$j.fits" - fid=$datadir"cls_fid_residual_masked_$j.fits" - chain=$output"param_chains.npz" - configcopy=$output"config_copy.yml" +for j in {0000..0020} +do + output=$mdir"output$j/" + echo $output + mkdir -p $output + + #cp test_BBSims_Hybrid/config_hybrid.yaml ${mdir}config.yml + #cp test_BBSims_Hybrid/config_hybrid_fixedepsilon.yaml ${mdir}config.yml + cp test_BBSims_Hybrid/config_hybrid_ellmin_fsky.yaml ${mdir}config.yml + + cp $cf $output + config=$output"config.yml" + sed -i "s/simJ/sim${simt}/" $config + sed -i "s/fskyX/fsky${fsky}/" $config + sed -i "s/paramsY/params_$j/" $config + sed -i "s/Zhybrid/${zz}hybrid/" $config - addqueue -s -q berg -c '2 hours' -m 4 -s -n 1x12 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config - done + sed -i "s/KK/${ellmin}/" $config + + datadir="/mnt/zfsusers/susanna/BBHybrid/data/sim${simt}/fsky${fsky}/" + #datadir="/mnt/zfsusers/susanna/BBHybrid/data/simd1s1_maskpysm_Bonly/fsky${fsky}/" + coadd=$datadir"cls_coadd_residual_masked_$j.fits" + noise=$datadir"cls_noise_residual_masked_$j.fits" + fid=$datadir"cls_fid_residual_masked_$j.fits" + chain=$output"param_chains.npz" + configcopy=$output"config_copy.yml" + + addqueue -s -q cmb -c 'res fsky0.3 1 hour' -m 4 -s -n 1x12 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config done + diff --git a/bbpower/compsep.py b/bbpower/compsep.py index 88d801d..8c78e00 100644 --- a/bbpower/compsep.py +++ b/bbpower/compsep.py @@ -17,20 +17,28 @@ class BBCompSep(PipelineStage): The foreground model parameters are defined in the config.yml file. """ name = "BBCompSep" - inputs = [('cells_coadded', FitsFile),('cells_noise', FitsFile),('cells_fiducial', FitsFile)] - outputs = [('param_chains', NpzFile), ('config_copy', YamlFile)] + inputs = [('cells_coadded', FitsFile), + ('cells_noise', FitsFile), + ('cells_fiducial', FitsFile), + ('cells_coadded_cov', FitsFile)] # Read covariance from separate file + outputs = [('param_chains', NpzFile), + ('config_copy', YamlFile)] config_options={'likelihood_type':'h&l', 'n_iters':32, 'nwalkers':16, 'r_init':1.e-3, - 'sampler':'emcee'} + 'sampler':'emcee', 'bands': 'all'} def setup_compsep(self): """ Pre-load the data, CMB BB power spectrum, and foreground models. """ self.parse_sacc_file() + if self.config['fg_model'].get('use_moments'): + print('Running moments') + self.precompute_w3j() self.load_cmb() self.fg_model = FGModel(self.config) self.params = ParameterManager(self.config) if self.config.get("diff"): + print('Running hybrid') self.hybridparams = np.load(self.config['resid_seds']) hyb_params = self.hybridparams['params_cent'] hyb_sigma = self.hybridparams['params_sigma'] @@ -38,10 +46,44 @@ def setup_compsep(self): self.params.p_free_priors[6][2] = [hyb_params[0], hyb_sigma[0]] # Beta_d: [centre, sigma] self.params.p_free_priors[2][2] = [hyb_params[1], hyb_sigma[1]] + if hyb_sigma[0]==np.nan: + hyb_sigma[0]=0.1 + if hyb_sigma[1]==np.nan: + hyb_sigma[1]=0.1 if self.use_handl: self.prepare_h_and_l() return + def get_moments_lmax(self): + return self.config['fg_model'].get('moments_lmax', 384) + + def precompute_w3j(self): + from pyshtools.utils import Wigner3j + + lmax = self.get_moments_lmax() + ells_w3j = np.arange(0, lmax) + w3j = np.zeros_like(ells_w3j, dtype=float) + self.big_w3j = np.zeros((lmax, lmax, lmax)) + for ell1 in ells_w3j[1:]: + for ell2 in ells_w3j[1:]: + w3j_array, ellmin, ellmax = Wigner3j(ell1, ell2, 0, 0, 0) + w3j_array = w3j_array[:ellmax - ellmin + 1] + # make the w3j_array the same shape as the w3j + if len(w3j_array) < len(ells_w3j): + reference = np.zeros(len(w3j)) + reference[:w3j_array.shape[0]] = w3j_array + w3j_array = reference + + w3j_array = np.concatenate([w3j_array[-ellmin:], + w3j_array[:-ellmin]]) + w3j_array = w3j_array[:len(ells_w3j)] + w3j_array[:ellmin] = 0 + + self.big_w3j[:, ell1, ell2] = w3j_array + + self.big_w3j = self.big_w3j**2 + + def matrix_to_vector(self, mat): return mat[..., self.index_ut_modes[0], self.index_ut_modes[1]] @@ -76,15 +118,20 @@ def _freq_pol_iterator(self): def parse_sacc_file(self): """ - Reads the data in the sacc file included the power spectra, bandpasses, and window functions. + Reads the data in the sacc file included the power spectra, + bandpasses, and window functions. """ + # Decide if you're using H&L self.use_handl = self.config['likelihood_type'] == 'h&l' + # Read data self.s = sacc.Sacc.load_fits(self.get_input('cells_coadded')) + self.s_cov = sacc.Sacc.load_fits(self.get_input('cells_coadded_cov')) if self.use_handl: s_fid = sacc.Sacc.load_fits(self.get_input('cells_fiducial')) s_noi = sacc.Sacc.load_fits(self.get_input('cells_noise')) + # Keep only desired correlations self.pols = self.config['pol_channels'] corr_all = ['cl_ee', 'cl_eb', 'cl_be', 'cl_bb'] corr_keep = [] @@ -95,19 +142,26 @@ def parse_sacc_file(self): for c in corr_all: if c not in corr_keep: self.s.remove_selection(c) + self.s_cov.remove_selection(c) if self.use_handl: s_fid.remove_selection(c) s_noi.remove_selection(c) + # Scale cuts self.s.remove_selection(ell__gt=self.config['l_max']) self.s.remove_selection(ell__lt=self.config['l_min']) + self.s_cov.remove_selection(ell__gt=self.config['l_max']) + self.s_cov.remove_selection(ell__lt=self.config['l_min']) if self.use_handl: s_fid.remove_selection(ell__gt=self.config['l_max']) s_fid.remove_selection(ell__lt=self.config['l_min']) s_noi.remove_selection(ell__gt=self.config['l_max']) s_noi.remove_selection(ell__lt=self.config['l_min']) - tr_names = sorted(list(self.s.tracers.keys())) + if self.config['bands'] == 'all': + tr_names = sorted(list(self.s.tracers.keys())) + else: + tr_names = self.config['bands'] self.nfreqs = len(tr_names) self.npol = len(self.pols) self.nmaps = self.nfreqs * self.npol @@ -115,6 +169,7 @@ def parse_sacc_file(self): self.ncross = (self.nmaps * (self.nmaps + 1)) // 2 self.pol_order=dict(zip(self.pols,range(self.npol))) + # Collect bandpasses self.bpss = [] for i_t, tn in enumerate(tr_names): t = self.s.tracers[tn] @@ -126,19 +181,24 @@ def parse_sacc_file(self): bnu = t.bandpass self.bpss.append(Bandpass(nu, dnu, bnu, i_t+1, self.config)) + # Get ell sampling + # Example power spectrum self.ell_b, _ = self.s.get_ell_cl('cl_' + 2 * self.pols[0].lower(), tr_names[0], tr_names[0]) + # Avoid l<2 win0 = self.s.data[0]['window'] mask_w = win0.values > 1 self.bpw_l = win0.values[mask_w] self.n_ell = len(self.bpw_l) self.n_bpws = len(self.ell_b) + # D_ell factor self.dl2cl = 2 * np.pi / (self.bpw_l * (self.bpw_l + 1)) self.windows = np.zeros([self.ncross, self.n_bpws, self.n_ell]) - if not (self.s.covariance.covmat.shape[-1] == len(self.s.mean) == self.n_bpws * self.ncross): - raise ValueError("C_ell vector's size is wrong") - + # Get power spectra and covariances + if self.config['bands'] == 'all': + if not (self.s_cov.covariance.covmat.shape[-1] == len(self.s.mean) == self.n_bpws * self.ncross): + raise ValueError("C_ell vector's size is wrong") v2d = np.zeros([self.n_bpws, self.ncross]) if self.use_handl: v2d_noi = np.zeros([self.n_bpws, self.ncross]) @@ -148,6 +208,7 @@ def parse_sacc_file(self): self.vector_indices = self.vector_to_matrix(np.arange(self.ncross, dtype=int)).astype(int) self.indx = [] + # Parse into the right ordering for b1, b2, p1, p2, m1, m2, ind_vec in self._freq_pol_iterator(): t1 = tr_names[b1] t2 = tr_names[b2] @@ -170,7 +231,7 @@ def parse_sacc_file(self): pol2b = self.pols[p2b].lower() cl_typb = f'cl_{pol1b}{pol2b}' ind_b = self.s.indices(cl_typb, (t1b, t2b)) - cv2d[:, ind_vec, :, ind_vecb] = self.s.covariance.covmat[ind_a][:, ind_b] + cv2d[:, ind_vec, :, ind_vecb] = self.s_cov.covariance.covmat[ind_a][:, ind_b] #Store data self.bbdata = self.vector_to_matrix(v2d) #(nbpw, nf, nf) @@ -200,22 +261,24 @@ def parse_sacc_file(self): for ix in range(self.ncross_modes): i, j = np.array(self.index_ut_modes)[:,ix] for jx in range(self.ncross_modes): - m, n = np.array(self.index_ut_modes)[:,jx] + m, n = np.array(self.index_ut_modes)[:,jx] + #import traceback + #traceback.print_stack() cv2d[:, ix, :, jx] = ccov[:, i, j, :, m, n] else: self.nmodes = self.nfreqs self.nmaps_modes = self.nmaps self.index_ut_modes = self.index_ut self.ncross_modes = self.ncross - self.bbcovar = cv2d.reshape([self.n_bpws * self.ncross_modes, self.n_bpws * self.ncross_modes]) - if self.use_handl: self.bbnoise = self.vector_to_matrix(v2d_noi) self.bbfiducial = self.vector_to_matrix(v2d_fid) - + self.bbcovar = cv2d.reshape([self.n_bpws * self.ncross_modes, + self.n_bpws * self.ncross_modes]) self.invcov = np.linalg.inv(self.bbcovar) + return - + def load_cmb(self): """ Loads the CMB BB spectrum as defined in the config file. @@ -227,6 +290,7 @@ def load_cmb(self): mask = (self.cmb_ells <= self.bpw_l.max()) & (self.cmb_ells > 1) self.cmb_ells = self.cmb_ells[mask] + # TODO: this is a patch nell = len(self.cmb_ells) self.cmb_tens = np.zeros([self.npol, self.npol, nell]) self.cmb_lens = np.zeros([self.npol, self.npol, nell]) @@ -241,7 +305,7 @@ def load_cmb(self): self.cmb_scal[ind, ind] = cmb_lensingfile[:, 2][mask] return - def integrate_seds(self, params): + def integrate_seds(self, params): #NB different from NERSC fg_scaling = np.zeros([self.fg_model.n_components, self.nfreqs]) rot_matrices = [] @@ -266,10 +330,10 @@ def sed(nu): return fg_scaling.T, rot_matrices - def evaluate_power_spectra(self, params): + def evaluate_power_spectra(self, params): #NB different from NERSC fg_pspectra = np.zeros([self.fg_model.n_components, self.fg_model.n_components, - self.npol, self.npol, self.n_ell]) + self.npol, self.npol, self.n_ell]) #NB different from NERSC # Fill diagonal for i_c, c_name in enumerate(self.fg_model.component_names): @@ -282,7 +346,7 @@ def evaluate_power_spectra(self, params): for k in clfunc.params] fg_pspectra[i_c, i_c, ip1, ip2, :] = clfunc.eval(self.bpw_l, *pspec_params) * self.dl2cl - # Off diagonals + # Off diagonals #NB different from NERSC for i_c1, c_name1 in enumerate(self.fg_model.component_names): for c_name2, epsname in self.fg_model.components[c_name1]['names_x_dict'].items(): i_c2 = self.fg_model.component_order[c_name2] @@ -293,25 +357,47 @@ def evaluate_power_spectra(self, params): return fg_pspectra - def model(self, params): + def model(self, params): #NB different from NERSC (dimensions) """ Defines the total model and integrates over the bandpasses and windows. """ + # [npol,npol,nell] cmb_cell = (params['r_tensor'] * self.cmb_tens + \ params['A_lens'] * self.cmb_lens + \ - self.cmb_scal) * self.dl2cl # [npol,npol,nell] - fg_scaling, rot_m = self.integrate_seds(params) # [nfreq, ncomp], [ncomp,nfreq,[matrix]] - fg_cell = self.evaluate_power_spectra(params) # [ncomp,ncomp,npol,npol,nell] - - cls_array_fg = np.zeros([self.nfreqs,self.nfreqs,self.n_ell,self.npol,self.npol]) - fg_cell = np.transpose(fg_cell, axes = [0,1,4,2,3]) # [ncomp,ncomp,nell,npol,npol] - cmb_cell = np.transpose(cmb_cell, axes = [2,0,1]) # [nell,npol,npol] + self.cmb_scal) * self.dl2cl + # [nell,npol,npol] + cmb_cell = np.transpose(cmb_cell, axes = [2,0,1]) + if self.config['cmb_model'].get('use_birefringence'): + bi_angle = np.radians(params['birefringence']) + c = np.cos(2*bi_angle) + s = np.sin(2*bi_angle) + bmat = np.array([[c, s], + [-s, c]]) + cmb_cell = rotate_cells_mat(bmat, bmat, cmb_cell) + + # [nfreq, ncomp], [ncomp,nfreq,[matrix]] + fg_scaling, rot_m = self.integrate_seds(params) #NB different from NERSC (dimensions) + # [ncomp,ncomp,npol,npol,nell] + fg_cell = self.evaluate_power_spectra(params) #NB different from NERSC (dimensions) + + # Add all components scaled in frequency (and HWP-rotated if needed) + # [nfreq, nfreq, nell, npol, npol] + cls_array_fg = np.zeros([self.nfreqs, self.nfreqs, + self.n_ell, self.npol, self.npol]) + # [ncomp,ncomp,nell,npol,npol] + fg_cell = np.transpose(fg_cell, axes = [0,1,4,2,3]) #NB different from NERSC + + #NB different from NERSC no CMB scaling + + # Loop over frequencies for f1 in range(self.nfreqs): + # Note that we only need to fill in half of the frequencies for f2 in range(f1,self.nfreqs): - cls=cmb_cell.copy() - + cls=cmb_cell.copy() #NB different from NERSC + + # Loop over component pairs for c1 in range(self.fg_model.n_components): - mat1=rot_m[c1][f1] + mat1=rot_m[c1][f1] #NB different from NERSC a1=fg_scaling[f1,c1] for c2 in range(self.fg_model.n_components): mat2=rot_m[c2][f2] @@ -320,7 +406,70 @@ def model(self, params): cls += clrot*a1*a2 cls_array_fg[f1,f2]=cls - cls_array_list = np.zeros([self.n_bpws, self.nfreqs, self.npol, self.nfreqs, self.npol]) + # Add moment terms if needed + # Start by assuming that there is no cross-correlation between components and betas + if self.config['fg_model'].get('use_moments'): + # TODO: moments work with: + # - B-only + # - No polarization angle business + # - Only power-law beta power spectra at l_pivot=80 + + # Evaluate 1st/2nd order SED derivatives. + # [nfreq, ncomp] + fg_scaling_d1 = self.integrate_seds_der(params, order=1) + fg_scaling_d2 = self.integrate_seds_der(params, order=2) + + # Compute 1x1 for each component + # Compute 0x2 for each component (essentially this is sigma_beta) + # Evaluate beta power spectra. + lmax_mom = self.get_moments_lmax() + # [ncomp, nell, npol, npol] + cls_11 = np.zeros([self.fg_model.n_components, self.n_ell, + self.npol, self.npol]) + # [ncomp, nell, npol, npol] + cls_02 = np.zeros([self.fg_model.n_components, self.n_ell, + self.npol, self.npol]) + + # component_1: dust, component_2: sync + for i_c, c_name in enumerate(self.fg_model.component_names): + comp = self.fg_model.components[c_name] + gamma = params[comp['names_moments_dict']['gamma_beta']] + amp = params[comp['names_moments_dict']['amp_beta']] * 1E-6 + cl_betas = self.bcls(lmax=lmax_mom, gamma=gamma, amp=amp) + #cl_cc = fg_cell[i_c, :] #NERSC version + cl_cc = fg_cell[i_c, i_c, :] #NB different from NERSC + cls_1x1 = self.evaluate_1x1(params, lmax=lmax_mom, + cls_cc=cl_cc, + cls_bb=cl_betas) + cls_11[i_c, :lmax_mom, :, :] = cls_1x1 + + cls_0x2 = self.evaluate_0x2(params, lmax=lmax_mom, + cls_cc=cl_cc, + cls_bb=cl_betas) + cls_02[i_c, :lmax_mom, :, :] = cls_0x2 + + # Add components scaled in frequency + for f1 in range(self.nfreqs): + # Note that we only need to fill in half of the frequencies + for f2 in range(f1, self.nfreqs): + cls = np.zeros([self.n_ell, self.npol, self.npol]) + for c1 in range(self.fg_model.n_components): + cls += (fg_scaling_d1[f1, c1] * fg_scaling_d1[f2, c1] * + cls_11[c1]) + #NERSC version + #cls += 0.5 * (fg_scaling_d2[f1, c1] * + # fg_scaling[c1, c1, f2, f2] + + # fg_scaling_d2[f2, c1] * + # fg_scaling[c1, c1, f1, f1]) * cls_02[c1] + cls += 0.5 * (fg_scaling_d2[f1, c1] * + fg_scaling[f2, c1] + + fg_scaling_d2[f2, c1] * + fg_scaling[f1, c1]) * cls_02[c1] + cls_array_fg[f1, f2] += cls + + # Window convolution + cls_array_list = np.zeros([self.n_bpws, self.nfreqs, + self.npol, self.nfreqs, self.npol]) for f1 in range(self.nfreqs): for p1 in range(self.npol): m1 = f1*self.npol+p1 @@ -333,10 +482,12 @@ def model(self, params): cls_array_list[:, f1, p1, f2, p2] = clband if m1!=m2: cls_array_list[:, f2, p2, f1, p1] = clband - + + # Polarization angle rotation for f1 in range(self.nfreqs): for f2 in range(self.nfreqs): - cls_array_list[:,f1,:,f2,:] = rotate_cells(self.bpss[f2], self.bpss[f1], + cls_array_list[:,f1,:,f2,:] = rotate_cells(self.bpss[f2], + self.bpss[f1], cls_array_list[:,f1,:,f2,:], params) cls_array_list = cls_array_list.reshape([self.n_bpws, self.nmaps, self.nmaps]) @@ -344,6 +495,67 @@ def model(self, params): cls_array_list = np.einsum("aj,bk,ijk->iab", self.R, self.R, cls_array_list) return cls_array_list + def bcls(self, lmax, gamma, amp): + """ + Model beta spectra as power laws. + """ + ls = np.arange(lmax) + bcls = np.zeros(len(ls)) + bcls[2:] = (ls[2:] / 80.)**gamma + return bcls*amp + + def integrate_seds_der(self, params, order=1): + """ + Define the first order derivative of the SED + """ + fg_scaling_der = np.zeros([self.fg_model.n_components, + self.nfreqs]) + + for i_c, c_name in enumerate(self.fg_model.component_names): + comp = self.fg_model.components[c_name] + units = comp['cmb_n0_norm'] + sed_params = [params[comp['names_sed_dict'][k]] + for k in comp['sed'].params] + + # Set SED function with scaling beta + def sed_der(nu): + nu0 = params[comp['names_sed_dict']['nu0']] + x = np.log(nu / nu0) + # This is only valid for spectral indices + return x**order * comp['sed'].eval(nu, *sed_params) + + for tn in range(self.nfreqs): + sed_b = self.bpss[tn].convolve_sed(sed_der, params)[0] + fg_scaling_der[i_c, tn] = sed_b * units + + return fg_scaling_der.T + + def evaluate_1x1(self, params, lmax, cls_cc, cls_bb): + """ + Evaluate the 1x1 moment for auto-spectra + """ + ls = np.arange(lmax) + v_left = (2*ls+1)[:, None, None] * cls_cc[:lmax, :, :] + v_right = (2*ls+1) * cls_bb[:lmax] + + mat = self.big_w3j + v_left = np.transpose(v_left, axes=[1, 0, 2]) + moment1x1 = np.dot(np.dot(mat, v_right), v_left) / (4*np.pi) + return moment1x1 + + def evaluate_0x2(self, params, lmax, cls_cc, cls_bb, + gamma=None): + """ + Evaluate the 0x2 moment for auto-spectra + Assume power law for beta + """ + ls = np.arange(lmax) + if gamma is not None: + prefac = amp * (2 * sp.zeta(-gamma-1) + sp.zeta(-gamma) - 3) / (4 * np.pi * 80**gamma) + else: + prefac = np.sum( (2 * ls + 1) * cls_bb) / (4*np.pi) + return cls_cc[:lmax] * prefac + def chi_sq_dx(self, params): """ Chi^2 likelihood. @@ -358,10 +570,7 @@ def prepare_h_and_l(self): return def h_and_l_dx(self, params): - """ - Hamimeche and Lewis likelihood. - Taken from Cobaya written by H, L and Torrado - See: https://github.com/CobayaSampler/cobaya/blob/master/cobaya/likelihoods/_cmblikes_prototype/_cmblikes_prototype.py + """ Hamimeche and Lewis likelihood. Taken from Cobaya written by H, L and Torrado See: https://github.com/CobayaSampler/cobaya/blob/master/cobaya/likelihoods/_cmblikes_prototype/_cmblikes_prototype.py """ model_cls = self.model(params) dx_vec = [] @@ -451,6 +660,8 @@ def emcee_sampler(self): end = time.time() return sampler, end-start + #NB different from NERSC: Add Polychord + def minimizer(self): """ Find maximum likelihood @@ -460,19 +671,26 @@ def chi2(par): c2=-2*self.lnprob(par) return c2 res=minimize(chi2, self.params.p0, method="Powell") - return res.x + + if len(self.params.p0)<2: + x = np.array([res.x]) + else: + x = res.x + return x def fisher(self): """ Evaluate Fisher matrix """ import numdifftools as nd + def lnprobd(p): l = self.lnprob(p) if l == -np.inf: l = -1E100 return l fisher = - nd.Hessian(lnprobd)(self.params.p0) + #NB different from NERSC: self.params.p0 given by chi^2 minim. return fisher def singlepoint(self): @@ -493,6 +711,7 @@ def timing(self, n_eval=300): end = time.time() return end-start, (end-start)/n_eval + def run(self): from shutil import copyfile self.setup_compsep() diff --git a/bbpower/fg_model.py b/bbpower/fg_model.py index 1fae6db..3fc7ef2 100644 --- a/bbpower/fg_model.py +++ b/bbpower/fg_model.py @@ -1,35 +1,52 @@ -import numpy as np import fgbuster.component_model as fgc import bbpower.fgcls as fgl class FGModel: """ - FGModel loads the foreground models and prepares the unit conversions to K_CMB units. - This creates a class that has an components attribute. The components attribute is a dictionary - of foreground models. Each foreground model is also a dictionary containing the SED function, - SED parameters, SED nu0, CMB nu0 normalization, and the foreground power spectrum parameters. + FGModel loads the foreground models and prepares the unit conversions + to K_CMB units. This creates a class that has an components attribute. + The components attribute is a dictionary of foreground models. Each + foreground model is also a dictionary containing the SED function, + SED parameters, SED nu0, CMB nu0 normalization, and the foreground + power spectrum parameters. """ def __init__(self, config): self.load_foregrounds(config) - return + return + + def component_iterator(self, config): + for key, component in config['fg_model'].items(): + if key.startswith('component_'): + yield key, component def load_foregrounds(self, config): - self.component_names=[] + self.component_names = [] self.components = {} self.component_order = {} i_comp = 0 - for key, component in config['fg_model'].items(): + for key, component in self.component_iterator(config): comp = {} - comp['names_x_dict']={} + decorr = component.get('decorr') + comp['decorr'] = False + if decorr: + comp['decorr'] = True + comp['decorr_param_names'] = {} + for k, l in decorr.items(): + comp['decorr_param_names'][l[0]] = k + + comp['names_x_dict'] = {} d_x = component.get('cross') if d_x: for pn, par in d_x.items(): if par[0] not in config['fg_model'].keys(): - raise KeyError("Component %s is not a valid component" % (par[0]) + + raise KeyError("Component %s " % (par[0]) + + "is not a valid component" + "to correlate %s with" % key) + if par[0] == key: + raise KeyError("%s is cross correlated with itself." % par[0]) comp['names_x_dict'][par[0]] = pn # Loop through SED parameters. @@ -43,17 +60,18 @@ def load_foregrounds(self, config): comp['names_sed_dict'][l[0]] = k # nu0 - if l[0]=='nu0': - if l[1]!='fixed': - raise ValueError("You can't vary reference frequencies!") + if l[0] == 'nu0': + if l[1] != 'fixed': + raise ValueError("You can't vary reference" + " frequencies!") nu0 = l[2][0] # SED parameter - if l[1]=='fixed': + if l[1] == 'fixed': val = l[2][0] else: val = None - params_fgc[l[0]]=val + params_fgc[l[0]] = val # Set units normalization if nu0 is not None: @@ -68,51 +86,52 @@ def load_foregrounds(self, config): # Same thing for C_ell parameters comp['names_cl_dict'] = {} params_fgl = {} - for k,d in component['cl_parameters'].items(): - p1,p2=k + for k, d in component['cl_parameters'].items(): + p1, p2 = k # Add parameters only if we're using both polarization channels - if (p1 in config['pol_channels']) and (p2 in config['pol_channels']): - comp['names_cl_dict'][k]={} - params_fgl[k]={} - for n,l in d.items(): - comp['names_cl_dict'][k][l[0]]=n - if l[0]=='ell0': - if l[1]!='fixed': - raise ValueError("You can't vary reference scales!") - if l[1]=='fixed': + if ((p1 in config['pol_channels']) and + (p2 in config['pol_channels'])): + comp['names_cl_dict'][k] = {} + params_fgl[k] = {} + for n, l in d.items(): + comp['names_cl_dict'][k][l[0]] = n + if l[0] == 'ell0': + if l[1] != 'fixed': + raise ValueError("You can't vary " + "reference scales!") + if l[1] == 'fixed': val = l[2][0] else: val = None - params_fgl[k][l[0]]=val - - ###Added------------------------------------------------------------------------ - ### Moment parameters - ##comp['names_moments_dict'] = {} - ##d = component.get('moments') - ##if d and config['fg_model'].get('use_moments'): - ## comp['moments_pameters'] = component['moments'] - ## for k, l in component['moments'].items(): - ## comp['names_moments_dict'][l[0]] = k - ###------------------------------------------------------------------------ + params_fgl[k][l[0]] = val + + # Moment parameters + comp['names_moments_dict'] = {} + d = component.get('moments') + if d and config['fg_model'].get('use_moments'): + comp['moments_pameters'] = component['moments'] + for k, l in component['moments'].items(): + comp['names_moments_dict'][l[0]] = k # Set Cl functions comp['cl'] = {} - for k,c in component['cl'].items(): - p1,p2=k + for k, c in component['cl'].items(): + p1, p2 = k # Add parameters only if we're using both polarization channels - if (p1 in config['pol_channels']) and (p2 in config['pol_channels']): + if ((p1 in config['pol_channels']) and + (p2 in config['pol_channels'])): cl_fnc = get_function(fgl, c) comp['cl'][k] = cl_fnc(**(params_fgl[k])) self.components[key] = comp self.component_names.append(key) self.component_order[key] = i_comp i_comp += 1 - self.n_components=len(self.component_names) + self.n_components = len(self.component_names) return -def get_function(mod,sed_name): +def get_function(mod, sed_name): try: - return getattr(mod,sed_name) + return getattr(mod, sed_name) except AttributeError: raise KeyError("Function named %s cannot be found" % (sed_name)) diff --git a/bbpower/param_manager.py b/bbpower/param_manager.py index bca10f3..67da208 100644 --- a/bbpower/param_manager.py +++ b/bbpower/param_manager.py @@ -6,7 +6,7 @@ class ParameterManager(object): def _add_parameter(self, p_name, p): # If fixed parameter, just add its name and value if p[1] == 'fixed': - self.p_fixed.append((p_name, p[2][0])) + self.p_fixed.append((p_name, float(p[2][0]))) return # Then move on # Otherwise it's free @@ -17,10 +17,10 @@ def _add_parameter(self, p_name, p): self.p_free_names.append(p_name) self.p_free_priors.append(p) # Add fiducial value to initial vector - if p[1] == 'tophat': - p0 = p[2][1] - elif p[1] == 'Gaussian': - p0 = p[2][0] + if np.char.lower(p[1]) == 'tophat': + p0 = float(p[2][1]) + elif np.char.lower(p[1]) == 'gaussian': + p0 = float(p[2][0]) else: raise ValueError("Unknown prior type %s" % p[1]) self.p0.append(p0) @@ -30,7 +30,14 @@ def _add_parameters(self, params): p = params[p_name] self._add_parameter(p_name, p) - def __init__(self,config): + def get_component_names(self, config): + comps = [] + for c_name in config['fg_model'].keys(): + if c_name.startswith('component_'): + comps.append(c_name) + return sorted(comps) + + def __init__(self, config): self.p_free_names = [] self.p_free_priors = [] self.p_fixed = [] @@ -42,20 +49,27 @@ def __init__(self,config): self._add_parameters(d['params']) # Loop through FG components - for c_name in sorted(config['fg_model'].keys()): + comp_names = self.get_component_names(config) + for c_name in comp_names: c = config['fg_model'][c_name] - for tag in ['sed_parameters', 'cross']: + for tag in ['sed_parameters', 'cross', 'decorr']: d = c.get(tag) if d: self._add_parameters(d) dc = c.get('cl_parameters') if dc: # Power spectra - for cl_name,d in dc.items(): - p1,p2=cl_name - # Add parameters only if we're using both polarization channels - if (p1 in config['pol_channels']) and (p2 in config['pol_channels']): + for cl_name, d in dc.items(): + p1, p2 = cl_name + # Add parameters only if we're using both + # polarization channels + if ((p1 in config['pol_channels']) and + (p2 in config['pol_channels'])): self._add_parameters(d) + dm = c.get('moments') + if dm and config['fg_model'].get('use_moments'): # Moments + self._add_parameters(dm) + # Loop through different systematics if 'systematics' in config.keys(): cnf_sys = config['systematics'] @@ -78,9 +92,9 @@ def build_params(self, par): def lnprior(self, par): lnp = 0 for p, pr in zip(par, self.p_free_priors): - if pr[1] == 'Gaussian': #Gaussian prior + if np.char.lower(pr[1]) == 'gaussian': # Gaussian prior lnp += -0.5 * ((p - pr[2][0])/pr[2][1])**2 - else: #Only other option is top-hat + else: # Only other option is top-hat if not(float(pr[2][0]) <= p <= float(pr[2][2])): return -np.inf return lnp diff --git a/bbpower/power_specter.py b/bbpower/power_specter.py index f97713e..a583f3f 100644 --- a/bbpower/power_specter.py +++ b/bbpower/power_specter.py @@ -1,48 +1,56 @@ from bbpipe import PipelineStage -from .types import FitsFile,TextFile,DummyFile +from .types import FitsFile, TextFile, DummyFile import sacc import numpy as np import healpy as hp import pymaster as nmt import os + class BBPowerSpecter(PipelineStage): """ Template for a power spectrum stage """ - name="BBPowerSpecter" - inputs=[('splits_list',TextFile),('masks_apodized',FitsFile),('bandpasses_list',TextFile), - ('sims_list',TextFile),('beams_list',TextFile)] - outputs=[('cells_all_splits',FitsFile),('cells_all_sims',TextFile),('mcm',DummyFile)] - config_options={'bpw_edges':None, - 'purify_B':True, - 'n_iter':3} + name = "BBPowerSpecter" + inputs = [('splits_list', TextFile), + ('masks_apodized', FitsFile), + ('bandpasses_list', TextFile), + ('sims_list', TextFile), + ('beams_list', TextFile)] + outputs = [('cells_all_splits', FitsFile), + ('cells_all_sims', TextFile), + ('mcm', DummyFile)] + config_options = {'bpw_edges': None, + 'purify_B': True, + 'n_iter': 3} def init_params(self): self.nside = self.config['nside'] self.npix = hp.nside2npix(self.nside) self.prefix_mcm = self.get_output('mcm')[:-4] - def read_beams(self,nbeams): + def read_beams(self, nbeams): from scipy.interpolate import interp1d beam_fnames = [] - with open(self.get_input('beams_list'),'r') as f: + with open(self.get_input('beams_list'), 'r') as f: for fname in f: beam_fnames.append(fname.strip()) # Check that there are enough beams - if len(beam_fnames)!=nbeams: - raise ValueError("Couldn't find enough beams %d != %d" (len(beam_fnames),nbeams)) + if len(beam_fnames) != nbeams: + raise ValueError("Couldn't find enough beams: " + "%d != %d" (len(beam_fnames), nbeams)) self.larr_all = np.arange(3*self.nside) - self.beams={} - for i_f,f in enumerate(beam_fnames): - li,bi=np.loadtxt(f,unpack=True) - bb=interp1d(li,bi,fill_value=0,bounds_error=False)(self.larr_all) - if li[0]!=0: - bb[:int(li[0])]=bi[0] - self.beams['band%d' % (i_f+1)]=bb + self.beams = {} + for i_f, f in enumerate(beam_fnames): + li, bi = np.loadtxt(f, unpack=True) + bb = interp1d(li, bi, fill_value=0, + bounds_error=False)(self.larr_all) + if li[0] != 0: + bb[:int(li[0])] = bi[0] + self.beams['band%d' % (i_f+1)] = bb def compute_cells_from_splits(self, splits_list): # Generate fields @@ -50,88 +58,93 @@ def compute_cells_from_splits(self, splits_list): fields = {} for b in range(self.n_bpss): for s in range(self.nsplits): - name = self.get_map_label(b,s) + name = self.get_map_label(b, s) print(" "+name) fname = splits_list[s] if not os.path.isfile(fname): # See if it's gzipped fname = fname + '.gz' if not os.path.isfile(fname): - raise ValueError("Can't find file ",splits_list[s]) - mp_q,mp_u=hp.read_map(fname, field=[2*b,2*b+1], verbose=False) - fields[name] = self.get_field(b,[mp_q,mp_u]) + raise ValueError("Can't find file ", splits_list[s]) + mp_q, mp_u = hp.read_map(fname, field=[2*b, 2*b+1], + verbose=False) + fields[name] = self.get_field(b, [mp_q, mp_u]) # Iterate over field pairs print(" Computing cross-spectra") cells = {} - for b1,b2,s1,s2,l1,l2 in self.get_cell_iterator(): - wsp = self.workspaces[self.get_workspace_label(b1,b2)] - if cells.get(l1) is None: # Create sub-dictionary if it doesn't exist - cells[l1]={} + for b1, b2, s1, s2, l1, l2 in self.get_cell_iterator(): + wsp = self.workspaces[self.get_workspace_label(b1, b2)] + # Create sub-dictionary if it doesn't exist + if cells.get(l1) is None: + cells[l1] = {} f1 = fields[l1] f2 = fields[l2] # Compute power spectrum print(" "+l1+" "+l2) - cells[l1][l2] = wsp.decouple_cell(nmt.compute_coupled_cell(f1,f2)) + cells[l1][l2] = wsp.decouple_cell(nmt.compute_coupled_cell(f1, f2)) return cells def read_bandpasses(self): bpss_fnames = [] - with open(self.get_input('bandpasses_list'),'r') as f: + with open(self.get_input('bandpasses_list'), 'r') as f: for fname in f: bpss_fnames.append(fname.strip()) self.n_bpss = len(bpss_fnames) - self.bpss={} - for i_f,f in enumerate(bpss_fnames): - nu,bnu=np.loadtxt(f,unpack=True) - dnu=np.zeros_like(nu) - dnu[1:]=np.diff(nu) - dnu[0]=dnu[1] - self.bpss['band%d' % (i_f+1)]={'nu':nu, 'dnu':dnu, 'bnu':bnu} - - def read_masks(self,nbands): - self.masks=[] + self.bpss = {} + for i_f, f in enumerate(bpss_fnames): + nu, bnu = np.loadtxt(f, unpack=True) + dnu = np.zeros_like(nu) + dnu[1:] = np.diff(nu) + dnu[0] = dnu[1] + self.bpss['band%d' % (i_f+1)] = {'nu': nu, + 'dnu': dnu, + 'bnu': bnu} + + def read_masks(self, nbands): + self.masks = [] for i in range(nbands): - m=hp.read_map(self.get_input('masks_apodized'), - verbose=False) - self.masks.append(hp.ud_grade(m,nside_out=self.nside)) + m = hp.read_map(self.get_input('masks_apodized'), + verbose=False) + self.masks.append(hp.ud_grade(m, nside_out=self.nside)) def get_bandpowers(self): # If it's a file containing the bandpower edges - if isinstance(self.config['bpw_edges'],str): + if isinstance(self.config['bpw_edges'], str): # Custom spacing edges = np.loadtxt(self.config['bpw_edges']).astype(int) - bpws = np.zeros(3*self.nside,dtype=int)-1 + bpws = np.zeros(3*self.nside, dtype=int)-1 weights = np.ones(3*self.nside) - for ibpw,(l0,lf) in enumerate(zip(edges[:-1],edges[1:])): - if lf<3*self.nside: - bpws[l0:lf]=ibpw + for ibpw, (l0, lf) in enumerate(zip(edges[:-1], edges[1:])): + if lf < 3*self.nside: + bpws[l0:lf] = ibpw # Add more equi-spaced bandpowers up to the end of the band - if edges[-1]<3*self.nside: + if edges[-1] < 3*self.nside: dell = edges[-1]-edges[-2] l0 = edges[-1] - while l0+dell<3*self.nside: - ibpw+=1 - bpws[l0:l0+dell]=ibpw - l0+=dell + while l0+dell < 3*self.nside: + ibpw += 1 + bpws[l0:l0+dell] = ibpw + l0 += dell is_dell = False if self.config.get('compute_dell'): is_dell = True - self.bins=nmt.NmtBin(self.nside, - bpws=bpws, - ells=self.larr_all, - weights=weights, - is_Dell=is_dell) - else: # otherwise it could be a constant integer interval - self.bins=nmt.NmtBin(self.nside,nlb=int(self.config['bpw_edges'])) - - def get_fname_workspace(self,band1,band2): - b1=min(band1,band2) - b2=max(band1,band2) + self.bins = nmt.NmtBin(self.nside, + bpws=bpws, + ells=self.larr_all, + weights=weights, + is_Dell=is_dell) + else: # otherwise it could be a constant integer interval + self.bins = nmt.NmtBin(self.nside, + nlb=int(self.config['bpw_edges'])) + + def get_fname_workspace(self, band1, band2): + b1 = min(band1, band2) + b2 = max(band1, band2) return self.prefix_mcm+"_%d_%d.fits" % (b1+1, b2+1) - def get_field(self,band,mps): + def get_field(self, band, mps): f = nmt.NmtField(self.masks[band], mps, beam=self.beams['band%d' % (band+1)], @@ -139,32 +152,33 @@ def get_field(self,band,mps): n_iter=self.config['n_iter']) return f - def compute_workspace(self,band1,band2): - b1=min(band1,band2) - b2=max(band1,band2) + def compute_workspace(self, band1, band2): + b1 = min(band1, band2) + b2 = max(band1, band2) w = nmt.NmtWorkspace() - fname = self.get_fname_workspace(b1,b2) + fname = self.get_fname_workspace(b1, b2) # If file exists, just read it if os.path.isfile(fname): print("Reading %d %d" % (b1, b2)) w.read_from(fname) else: print("Computing %d %d" % (b1, b2)) - mdum = np.zeros([2,self.npix]) - f1=self.get_field(b1,mdum) - f2=self.get_field(b2,mdum) - w.compute_coupling_matrix(f1,f2,self.bins,n_iter=self.config['n_iter']) + mdum = np.zeros([2, self.npix]) + f1 = self.get_field(b1, mdum) + f2 = self.get_field(b2, mdum) + w.compute_coupling_matrix(f1, f2, self.bins, + n_iter=self.config['n_iter']) w.write_to(fname) return w - def get_map_label(self,band,split): - return 'band%d_split%d' % (band+1, split+1) + def get_map_label(self, band, split): + return 'band%d_split%d' % (band+1, split+1) - def get_workspace_label(self,band1,band2): - b1=min(band1,band2) - b2=max(band1,band2) + def get_workspace_label(self, band1, band2): + b1 = min(band1, band2) + b2 = max(band1, band2) return 'b%d_b%d' % (b1+1, b2+1) def compute_workspaces(self): @@ -172,24 +186,24 @@ def compute_workspaces(self): # Assumption is that mask is different across bands, # but the same across polarization channels and splits. print("Estimating mode-coupling matrices") - self.workspaces={} + self.workspaces = {} for i1 in range(self.n_bpss): - for i2 in range(i1,self.n_bpss): - name=self.get_workspace_label(i1,i2) - self.workspaces[name] = self.compute_workspace(i1,i2) + for i2 in range(i1, self.n_bpss): + name = self.get_workspace_label(i1, i2) + self.workspaces[name] = self.compute_workspace(i1, i2) def get_cell_iterator(self): for b1 in range(self.n_bpss): - for b2 in range(b1,self.n_bpss): + for b2 in range(b1, self.n_bpss): for s1 in range(self.nsplits): - l1=self.get_map_label(b1,s1) - if b1==b2: - splits_range=range(s1,self.nsplits) + l1 = self.get_map_label(b1, s1) + if b1 == b2: + splits_range = range(s1, self.nsplits) else: - splits_range=range(self.nsplits) + splits_range = range(self.nsplits) for s2 in splits_range: - l2=self.get_map_label(b2,s2) - yield(b1,b2,s1,s2,l1,l2) + l2 = self.get_map_label(b2, s2) + yield(b1, b2, s1, s2, l1, l2) def get_sacc_tracers(self): sacc_t = [] @@ -204,19 +218,23 @@ def get_sacc_tracers(self): bandpass_extra={'dnu': bpss['dnu']}) sacc_t.append(T) return sacc_t - + def get_sacc_windows(self): windows_wsp = {} for b1 in range(self.n_bpss): - for b2 in range(b1,self.n_bpss): + for b2 in range(b1, self.n_bpss): name = self.get_workspace_label(b1, b2) - windows_wsp[name]={} + windows_wsp[name] = {} wsp = self.workspaces[name] bpw_win = wsp.get_bandpower_windows() - windows_wsp[name]['EE'] = sacc.BandpowerWindow(self.larr_all, bpw_win[0, :, 0, :].T) - windows_wsp[name]['EB'] = sacc.BandpowerWindow(self.larr_all, bpw_win[1, :, 1, :].T) - windows_wsp[name]['BE'] = sacc.BandpowerWindow(self.larr_all, bpw_win[2, :, 2, :].T) - windows_wsp[name]['BB'] = sacc.BandpowerWindow(self.larr_all, bpw_win[3, :, 3, :].T) + windows_wsp[name]['EE'] = sacc.BandpowerWindow(self.larr_all, + bpw_win[0, :, 0, :].T) + windows_wsp[name]['EB'] = sacc.BandpowerWindow(self.larr_all, + bpw_win[1, :, 1, :].T) + windows_wsp[name]['BE'] = sacc.BandpowerWindow(self.larr_all, + bpw_win[2, :, 2, :].T) + windows_wsp[name]['BB'] = sacc.BandpowerWindow(self.larr_all, + bpw_win[3, :, 3, :].T) return windows_wsp def save_cell_to_file(self, cell, tracers, fname, with_windows=False): @@ -229,15 +247,15 @@ def save_cell_to_file(self, cell, tracers, fname, with_windows=False): # Add each power spectrum l_eff = self.bins.get_effective_ells() - for b1,b2,s1,s2,l1,l2 in self.get_cell_iterator(): - add_BE = not ((b1==b2) and (s1==s2)) + for b1, b2, s1, s2, l1, l2 in self.get_cell_iterator(): + add_BE = not ((b1 == b2) and (s1 == s2)) if with_windows: wname = self.get_workspace_label(b1, b2) s.add_ell_cl('cl_ee', l1, l2, l_eff, cell[l1][l2][0], window=self.win[wname]['EE']) # EE s.add_ell_cl('cl_eb', l1, l2, l_eff, cell[l1][l2][1], window=self.win[wname]['EB']) # EB - if add_BE: #Only add B1E2 if 1!=2 + if add_BE: # Only add B1E2 if 1!=2 s.add_ell_cl('cl_be', l1, l2, l_eff, cell[l1][l2][2], window=self.win[wname]['BE']) # BE s.add_ell_cl('cl_bb', l1, l2, l_eff, cell[l1][l2][3], @@ -245,14 +263,14 @@ def save_cell_to_file(self, cell, tracers, fname, with_windows=False): else: s.add_ell_cl('cl_ee', l1, l2, l_eff, cell[l1][l2][0]) # EE s.add_ell_cl('cl_eb', l1, l2, l_eff, cell[l1][l2][1]) # EB - if add_BE: #Only add B1E2 if 1!=2 + if add_BE: # Only add B1E2 if 1!=2 s.add_ell_cl('cl_be', l1, l2, l_eff, cell[l1][l2][2]) # BE s.add_ell_cl('cl_bb', l1, l2, l_eff, cell[l1][l2][3]) # EE print("Saving to "+fname) s = s.save_fits(fname, overwrite=True) - def run(self) : + def run(self): self.init_params() # Read bandpasses @@ -275,7 +293,7 @@ def run(self) : # Compile list of splits splits = [] - with open(self.get_input('splits_list'),'r') as f: + with open(self.get_input('splits_list'), 'r') as f: for fname in f: splits.append(fname.strip()) self.nsplits = len(splits) @@ -298,34 +316,48 @@ def run(self) : with_windows=True) # Iterate over simulations sims = [] - with open(self.get_input('sims_list'),'r') as f: + with open(self.get_input('sims_list'), 'r') as f: for dname in f: sims.append(dname.strip()) # Write all output file names into a text file - fo=open(self.get_output('cells_all_sims'),'w') - prefix_out=self.get_output('cells_all_splits')[:-5] - for isim,d in enumerate(sims): - fname=prefix_out + "_sim%d.fits" % isim - fo.write(fname+"\n") - fo.close() - - for isim,d in enumerate(sims): - fname=prefix_out + "_sim%d.fits" % isim - if os.path.isfile(fname): - print("found " + fname) - continue - print("%d-th / %d simulation" % (isim+1, len(sims))) - # Compute list of splits - sim_splits = [d+'/obs_split%dof%d.fits' % (i+1, self.nsplits) - for i in range(self.nsplits)] - # Compute all possible cross-power spectra - cell_sim=self.compute_cells_from_splits(sim_splits) - # Save output - fname=prefix_out + "_sim%d.fits" % isim - self.save_cell_to_file(cell_sim, - self.tracers, - fname, with_windows=False) + # if txt file is not existent + tname = self.get_output('cells_all_sims') + prefix_out = self.get_output('cells_all_splits')[:-5] + if os.path.isfile(tname): + print("found " + tname) + else: + fo = open(self.get_output('cells_all_sims'), 'w') + for isim, d in enumerate(sims): + fname = prefix_out + "_sim%d.fits" % isim + fo.write(fname+"\n") + fo.close() + + # Compile list of cells from txt (added SA) + ccls_all = [] + with open(tname, 'r') as t: + for ttname in t: + ccls_all.append(ttname) + for isim, d in enumerate(sims): + # changed SA + ffname = (ccls_all[isim]).strip() #stripping white spaces + if os.path.isfile(ffname): + print("found " + ffname) #don't re-compute + else: + print("%d-th / %d simulation" % (isim+1, len(sims))) + # Compute list of splits + sim_splits = [d+'/obs_split%dof%d.fits' % (i+1, self.nsplits) + for i in range(self.nsplits)] + #sim_splits = [d+'/SO_SAT_obs_map_split_%dof%d.fits' % (i+1, self.nsplits) + # for i in range(self.nsplits)] + # Compute all possible cross-power spectra + cell_sim = self.compute_cells_from_splits(sim_splits) + # Save output + fname = prefix_out + "_sim%d.fits" % isim + self.save_cell_to_file(cell_sim, + self.tracers, + fname, with_windows=False) + if __name__ == '__main__': cls = PipelineStage.main() diff --git a/bbpower/power_summarizer.py b/bbpower/power_summarizer.py index 0fa249a..fa40c20 100644 --- a/bbpower/power_summarizer.py +++ b/bbpower/power_summarizer.py @@ -15,6 +15,7 @@ class BBPowerSummarizer(PipelineStage): 'data_covar_type':'block_diagonal', 'data_covar_diag_order': 3} + def get_covariance_from_samples(self,v,s,covar_type='dense', off_diagonal_cut=0): """ diff --git a/plotting/labels.py b/plotting/labels.py index 59e9915..65814d5 100644 --- a/plotting/labels.py +++ b/plotting/labels.py @@ -12,7 +12,8 @@ 'shift_10':'\Delta\\nu_{10}', 'shift_11':'\Delta\\nu_{11}', 'shift_12':'\Delta\\nu_{12}', 'angle_1':'\phi_1', 'angle_2':'\phi_2', 'angle_3':'\phi_3', 'angle_4':'\phi_4', 'angle_5':'\phi_5', 'angle_6':'\phi_6', 'dphi1_1':'\Delta\phi_1', 'dphi1_2':'\Delta\phi_2', - 'dphi1_3':'\Delta\phi_3', 'dphi1_4':'\Delta\phi_4', 'dphi1_5':'\Delta\phi_5', 'dphi1_6':'\Delta\phi_6'} + 'dphi1_3':'\Delta\phi_3', 'dphi1_4':'\Delta\phi_4', 'dphi1_5':'\Delta\phi_5', 'dphi1_6':'\Delta\phi_6', + 'amp_d_beta': 'B_d', 'gamma_d_beta': '\\gamma_{d}', 'amp_s_beta': 'B_s', 'gamma_s_beta': '\\gamma_{s}'} allranges={'$r$': [-1.0, 1.0], @@ -34,5 +35,8 @@ '$A^{EB}_s$': [-100.0, 100.0], '$\\alpha^{EB}_s$': [-3.0, 1.0], '$\\theta$': [-30., 30.], - '$\\Delta_s$': [0.9, 1.1]} - + '$\\Delta_s$': [0.9, 1.1], + '$B_d$':[0., 10], + '$\\gamma_{d}$':[-6, -2], + '$B_s$':[0, 10], + '$\\gamma_{s}$':[-6, -2]} From 0a215f5151c1865932bb2e3ea6e17c308a17dc90 Mon Sep 17 00:00:00 2001 From: susannaaz Date: Thu, 25 Aug 2022 15:07:36 +0100 Subject: [PATCH 16/20] hybrid import betas more general --- bbpower/compsep.py | 41 ++++++++++++++++++++++++++++++++-------- bbpower/param_manager.py | 1 + 2 files changed, 34 insertions(+), 8 deletions(-) diff --git a/bbpower/compsep.py b/bbpower/compsep.py index 8c78e00..3850635 100644 --- a/bbpower/compsep.py +++ b/bbpower/compsep.py @@ -42,14 +42,13 @@ def setup_compsep(self): self.hybridparams = np.load(self.config['resid_seds']) hyb_params = self.hybridparams['params_cent'] hyb_sigma = self.hybridparams['params_sigma'] + pnames = np.array(self.params.p_free_names) + ind_bs = int(np.where(pnames == 'beta_s')[0]) + ind_bd = int(np.where(pnames == 'beta_d')[0]) # Beta_s: [centre, sigma] - self.params.p_free_priors[6][2] = [hyb_params[0], hyb_sigma[0]] + self.params.p_free_priors[ind_bs][2] = [hyb_params[0], hyb_sigma[0]] # Beta_d: [centre, sigma] - self.params.p_free_priors[2][2] = [hyb_params[1], hyb_sigma[1]] - if hyb_sigma[0]==np.nan: - hyb_sigma[0]=0.1 - if hyb_sigma[1]==np.nan: - hyb_sigma[1]=0.1 + self.params.p_free_priors[ind_bd][2] = [hyb_params[1], hyb_sigma[1]] if self.use_handl: self.prepare_h_and_l() return @@ -644,8 +643,11 @@ def emcee_sampler(self): if not found_file: backend.reset(nwalkers,ndim) - p0 = self.minimizer() - pos = [p0 + 1.e-3*np.random.randn(ndim) for i in range(nwalkers)] + #p0 = self.minimizer() + #pos = [p0 + 1.e-3*np.random.randn(ndim) for i in range(nwalkers)] + #nsteps_use = n_iters + pos = [self.params.p0 + 1.e-3*np.random.randn(ndim) + for i in range(nwalkers)] nsteps_use = n_iters else: print("Restarting from previous run") @@ -676,6 +678,29 @@ def chi2(par): x = np.array([res.x]) else: x = res.x + + ############################################################# + ### Plot Cls + ##bbdata_cls = self.bbdata # n_bpws, nfreqs, nfreqs + ##names=self.params.p_free_names + ##parameters = dict(zip(list(names),res.x)) + ###print(parameters) + ##parameters['nu0_d'] = 220. + ##parameters['nu0_s'] = 40. + ##model_cls = self.model(parameters) # n_bpws, nfreqs, nfreqs + ## + ### Select the upper triangle + ##for i in range(self.nfreqs): + ## for j in range(i, self.nfreqs): + ## import matplotlib.pyplot as plt + ## plt.figure() + ## plt.plot(self.ell_b, bbdata_cls[:,i,j], label='data') + ## plt.plot(self.ell_b, model_cls[:,i,j], label='theory') + ## plt.legend() + ## plt.savefig(f'03_s000_model_vs_data_{i}_{j}.png',bbox_inches='tight') + ##exit() + ############################################################# + return x def fisher(self): diff --git a/bbpower/param_manager.py b/bbpower/param_manager.py index 67da208..ef9819f 100644 --- a/bbpower/param_manager.py +++ b/bbpower/param_manager.py @@ -92,6 +92,7 @@ def build_params(self, par): def lnprior(self, par): lnp = 0 for p, pr in zip(par, self.p_free_priors): + #if not p[1] == 'fixed': if np.char.lower(pr[1]) == 'gaussian': # Gaussian prior lnp += -0.5 * ((p - pr[2][0])/pr[2][1])**2 else: # Only other option is top-hat From e6086f9d7081061b33e31853e74c8052d30f8f85 Mon Sep 17 00:00:00 2001 From: susannaaz Date: Thu, 25 Aug 2022 16:03:43 +0100 Subject: [PATCH 17/20] Read existing covariance if available --- bbpower/power_summarizer.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/bbpower/power_summarizer.py b/bbpower/power_summarizer.py index fa40c20..f50ae95 100644 --- a/bbpower/power_summarizer.py +++ b/bbpower/power_summarizer.py @@ -38,6 +38,13 @@ def get_covariance_from_samples(self,v,s,covar_type='dense', cuts -= np.diag(np.ones(self.n_bpws-i),k=-i) cov = cov.reshape([nblocks, self.n_bpws, nblocks, self.n_bpws]) cov = (cov * cuts[None, :, None, :]).reshape([nd, nd]) + else: #import pre-computed covariance from sims added SA + import os + if not os.path.isfile(covar_type+'cells_coadded.fits'): + raise ValueError("Can't find file ", covar_type+'cells_coadded.fits') + else: + scv = sacc.Sacc().load_fits(covar_type+'cells_coadded.fits') + cov = scv.covariance.covmat s.add_covariance(cov) def init_params(self): From b3bede68e55d946f0ebfaaa03474864baea860b2 Mon Sep 17 00:00:00 2001 From: susannaaz Date: Thu, 25 Aug 2022 16:51:53 +0100 Subject: [PATCH 18/20] Do not recompute covariance if pre-existing --- bbpower/power_summarizer.py | 82 +++++++++++++++++++++---------------- 1 file changed, 46 insertions(+), 36 deletions(-) diff --git a/bbpower/power_summarizer.py b/bbpower/power_summarizer.py index f50ae95..8a5dbb6 100644 --- a/bbpower/power_summarizer.py +++ b/bbpower/power_summarizer.py @@ -38,12 +38,13 @@ def get_covariance_from_samples(self,v,s,covar_type='dense', cuts -= np.diag(np.ones(self.n_bpws-i),k=-i) cov = cov.reshape([nblocks, self.n_bpws, nblocks, self.n_bpws]) cov = (cov * cuts[None, :, None, :]).reshape([nd, nd]) - else: #import pre-computed covariance from sims added SA + else: #import pre-computed covariance from sims import os - if not os.path.isfile(covar_type+'cells_coadded.fits'): - raise ValueError("Can't find file ", covar_type+'cells_coadded.fits') + if not os.path.isfile(covar_type+'cells_coadded.sacc'): + raise ValueError("Can't find file ", covar_type+'cells_coadded.sacc') else: - scv = sacc.Sacc().load_fits(covar_type+'cells_coadded.fits') + print(f'Importing covariance from {covar_type}') + scv = sacc.Sacc().load_fits(covar_type+'cells_coadded.sacc') cov = scv.covariance.covmat s.add_covariance(cov) @@ -350,38 +351,47 @@ def run(self): # Read data file, coadd and compute nulls print("Reading data") - summ = self.parse_splits_sacc_file(self.s_splits, get_saccs=True, with_windows=True) - - # Read simulations - print("Reading simulations") - sim_cd_t=np.zeros([self.nsims,len(summ['spectra'][0])]) - sim_cd_x=np.zeros([self.nsims,len(summ['spectra'][1])]) - sim_cd_n=np.zeros([self.nsims,len(summ['spectra'][2])]) - sim_null=np.zeros([self.nsims,len(summ['spectra'][3])]) - for i, fn in enumerate(self.fname_sims): - print(fn) - s=sacc.Sacc.load_fits(fn) - sb=self.parse_splits_sacc_file(s) - sim_cd_t[i,:]=sb['spectra'][0] - sim_cd_x[i,:]=sb['spectra'][1] - sim_cd_n[i,:]=sb['spectra'][2] - sim_null[i,:]=sb['spectra'][3] - - # Compute covariance - print("Covariances") - cov_cd_t=self.get_covariance_from_samples(sim_cd_t, summ['saccs'][0], - covar_type=self.config['data_covar_type'], - off_diagonal_cut=self.config['data_covar_diag_order']) - cov_cd_x=self.get_covariance_from_samples(sim_cd_x, summ['saccs'][1], - covar_type=self.config['data_covar_type'], - off_diagonal_cut=self.config['data_covar_diag_order']) - cov_cd_n=self.get_covariance_from_samples(sim_cd_n, summ['saccs'][2], - covar_type=self.config['data_covar_type'], - off_diagonal_cut=self.config['data_covar_diag_order']) - # There are so many nulls that we'll probably run out of memory - cov_null=self.get_covariance_from_samples(sim_null, summ['saccs'][3], - covar_type=self.config['nulls_covar_type'], - off_diagonal_cut=self.config['nulls_covar_diag_order']) + summ = self.parse_splits_sacc_file(self.s_splits, + get_saccs=True, + with_windows=True) + + # Re-compute covariance only if it is not pre-computed + # i.e. if no directory is specified in covar_type in config + if self.config['data_covar_type']=='dense' or self.config['data_covar_type']=='diagonal' or self.config['data_covar_type']=='block_diagonal': + # Read simulations + print("Reading simulations") + sim_cd_t = np.zeros([self.nsims, len(summ['spectra'][0])]) + sim_cd_x = np.zeros([self.nsims, len(summ['spectra'][1])]) + sim_cd_n = np.zeros([self.nsims, len(summ['spectra'][2])]) + sim_null = np.zeros([self.nsims, len(summ['spectra'][3])]) + for i, fn in enumerate(self.fname_sims): + #print(fn) + s = sacc.Sacc.load_fits(fn) + sb = self.parse_splits_sacc_file(s) + sim_cd_t[i, :] = sb['spectra'][0] + sim_cd_x[i, :] = sb['spectra'][1] + sim_cd_n[i, :] = sb['spectra'][2] + sim_null[i, :] = sb['spectra'][3] + + # Compute covariance + print("Covariances") + dctyp = self.config['data_covar_type'] + dcord = self.config['data_covar_diag_order'] + self.get_covariance_from_samples(sim_cd_t, summ['saccs'][0], + covar_type=dctyp, + off_diagonal_cut=dcord) + self.get_covariance_from_samples(sim_cd_x, summ['saccs'][1], + covar_type=dctyp, + off_diagonal_cut=dcord) + self.get_covariance_from_samples(sim_cd_n, summ['saccs'][2], + covar_type=dctyp, + off_diagonal_cut=dcord) + # There are so many nulls that we'll probably run out of memory + nctyp = self.config['nulls_covar_type'] + ncord = self.config['nulls_covar_diag_order'] + self.get_covariance_from_samples(sim_null, summ['saccs'][3], + covar_type=nctyp, + off_diagonal_cut=ncord) # Save data print("Writing output") From 3783da470183145cac4f708c85241091a8b580ac Mon Sep 17 00:00:00 2001 From: susannaaz Date: Tue, 23 May 2023 14:16:57 +0100 Subject: [PATCH 19/20] old changes --- autoqbaseline.sh | 8 ++--- autoqresiduals.sh | 16 +++++----- bbpower/compsep.py | 60 +++++++++++++++++++++++++++++++++++-- bbpower/power_summarizer.py | 1 - 4 files changed, 70 insertions(+), 15 deletions(-) diff --git a/autoqbaseline.sh b/autoqbaseline.sh index a65ed02..d6f2246 100755 --- a/autoqbaseline.sh +++ b/autoqbaseline.sh @@ -21,16 +21,16 @@ ## addqueue -q berg -c "2 hours" -m 1 -s -n 1x8 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config ## done ##done -ellmin=30 #30 #2 #10 +ellmin=2 #30 #2 #10 #simt=d1s1_maskpysm_Bonly #simt=d1s1_maskpysm_Bonly_r0.01 #simt=d1s1_maskpysm_Bonly_r0.01_whitenoiONLY simt=d1s1_maskpysm_Bonly_whitenoiONLY -fsky=0.3 #0.6 #0.8 +fsky=0.8 #0.6 #0.8 #mdir="/mnt/extraspace/susanna/BBHybrid/baseline_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}/" #mdir="/mnt/extraspace/susanna/BBHybrid/baseline_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_r0.01/" #mdir="/mnt/extraspace/susanna/BBHybrid/baseline_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_r0.01_whitenoiONLY/" -mdir="/mnt/extraspace/susanna/BBHybrid/baseline_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_whitenoiONLY/" +mdir="/mnt/extraspace/susanna/BBHybrid/baseline_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_whitenoiONLY_HandLlik/" cf=$mdir"config.yml" for j in {0000..0020} #0000 do @@ -53,7 +53,7 @@ do chain=$output"param_chains.npz" configcopy=$output"config_copy.yml" - addqueue -s -q cmb -c 'baseline fsky0.3 1 hour' -m 4 -s -n 1x12 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config + addqueue -s -q berg -c 'baseline_fsky0.8_ellm2 2.5hr' -m 4 -s -n 1x12 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_coadded_cov=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config done diff --git a/autoqresiduals.sh b/autoqresiduals.sh index be1d22b..a2dda2a 100755 --- a/autoqresiduals.sh +++ b/autoqresiduals.sh @@ -1,12 +1,12 @@ #!/bin/bash -ellmin=30 #2 #10 #30 +ellmin=2 #2 #10 #30 #simt=d1s1_maskpysm_Bonly #simt=d1s1_maskpysm_Bonly_r0.01 -#simt=d1s1_maskpysm_Bonly_r0.01_whitenoiONLY -simt=d1s1_maskpysm_Bonly_whitenoiONLY +simt=d1s1_maskpysm_Bonly_r0.01_whitenoiONLY +#simt=d1s1_maskpysm_Bonly_whitenoiONLY -fsky=0.3 #0.3 0.6 #0.8 +fsky=0.8 #0.3 0.6 #0.8 zz=masked_ #full_ #mdir="/mnt/extraspace/susanna/BBHybrid/hybrid_masked_outs/mask_planck/realistic/d1s1/fsky${fsky}/" @@ -14,10 +14,12 @@ zz=masked_ #full_ #mdir="/mnt/extraspace/susanna/BBHybrid/hybrid_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}/" #mdir="/mnt/extraspace/susanna/BBHybrid/hybrid_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_r0.01/" #mdir="/mnt/extraspace/susanna/BBHybrid/hybrid_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_r0.01_whitenoiONLY/" -mdir="/mnt/extraspace/susanna/BBHybrid/hybrid_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_whitenoiONLY/" +#mdir="/mnt/extraspace/susanna/BBHybrid/hybrid_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_whitenoiONLY/" +#mdir="/mnt/extraspace/susanna/BBHybrid/hybrid_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_whitenoiONLY_HandLlik/" +mdir="/mnt/extraspace/susanna/BBHybrid/hybrid_masked_outs/mask_pysm/realistic/d1s1/fsky${fsky}_ellmin${ellmin}_r0.01_whitenoiONLY_HandLlik/" cf=$mdir"config.yml" -for j in {0000..0020} +for j in {0000..0020} #0000 do output=$mdir"output$j/" echo $output @@ -44,6 +46,6 @@ do chain=$output"param_chains.npz" configcopy=$output"config_copy.yml" - addqueue -s -q cmb -c 'res fsky0.3 1 hour' -m 4 -s -n 1x12 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config + addqueue -s -q cmb -c 'res_fsky0.8_ellm10_h&l 2.5 hour' -m 4 -s -n 1x12 /usr/bin/python3 -m bbpower BBCompSep --cells_coadded=$coadd --cells_coadded_cov=$coadd --cells_noise=$noise --cells_fiducial=$fid --param_chains=$chain --config_copy=$configcopy --config=$config done diff --git a/bbpower/compsep.py b/bbpower/compsep.py index 3850635..4db286b 100644 --- a/bbpower/compsep.py +++ b/bbpower/compsep.py @@ -272,6 +272,15 @@ def parse_sacc_file(self): if self.use_handl: self.bbnoise = self.vector_to_matrix(v2d_noi) self.bbfiducial = self.vector_to_matrix(v2d_fid) + if self.config.get("diff"): + from scipy.linalg import eig, inv, pinv + self.hybridparams = np.load(self.config['resid_seds']) + Qmat = self.hybridparams['Q'] + eval, evec = eig(Qmat) + mask = eval > 0.1 #take away rows equivalent to eigenvalue=0 + self.R = np.linalg.inv(evec)[mask, :] + self.bbnoise = np.einsum("aj,bk,ijk->iab", self.R, self.R, self.bbnoise) + self.bbfiducial = np.einsum("aj,bk,ijk->iab", self.R, self.R, self.bbfiducial) self.bbcovar = cv2d.reshape([self.n_bpws * self.ncross_modes, self.n_bpws * self.ncross_modes]) self.invcov = np.linalg.inv(self.bbcovar) @@ -565,11 +574,16 @@ def chi_sq_dx(self, params): def prepare_h_and_l(self): fiducial_noise = self.bbfiducial + self.bbnoise self.Cfl_sqrt = np.array([sqrtm(f) for f in fiducial_noise]) + print(self.bbdata.shape) + print(self.bbnoise.shape) self.observed_cls = self.bbdata + self.bbnoise return def h_and_l_dx(self, params): - """ Hamimeche and Lewis likelihood. Taken from Cobaya written by H, L and Torrado See: https://github.com/CobayaSampler/cobaya/blob/master/cobaya/likelihoods/_cmblikes_prototype/_cmblikes_prototype.py + """ + Hamimeche and Lewis likelihood. + Taken from Cobaya written by H, L and Torrado + See: https://github.com/CobayaSampler/cobaya/blob/master/cobaya/likelihoods/_cmblikes_prototype/_cmblikes_prototype.py """ model_cls = self.model(params) dx_vec = [] @@ -664,7 +678,7 @@ def emcee_sampler(self): #NB different from NERSC: Add Polychord - def minimizer(self): + def minimizer(self, save_spectra=False): """ Find maximum likelihood """ @@ -679,6 +693,44 @@ def chi2(par): else: x = res.x + if save_spectra: + bbdata_cls = self.bbdata # n_bpws, nfreqs, nfreqs + names=self.params.p_free_names + parameters = dict(zip(list(names),res.x)) + #print(parameters) + parameters['nu0_d'] = 220. + parameters['nu0_s'] = 40. + model_cls = self.model(parameters) # n_bpws, nfreqs, nfreqs + bbcov = self.bbcovar + l = self.ell_b + lmax = 300 + msk = l Date: Tue, 23 May 2023 14:21:02 +0100 Subject: [PATCH 20/20] bash scripts + config files examples --- run_pipeline/bbcompsep_parallel.sh | 37 ++++++ run_pipeline/bbpowerspecter_parallel.sh | 54 ++++++++ run_pipeline/bbpowersummarizer_parallel.sh | 37 ++++++ run_pipeline/config.yml | 137 +++++++++++++++++++++ run_pipeline/config_hybrid.yml | 131 ++++++++++++++++++++ run_pipeline/prepare_directory_parallel.sh | 30 +++++ run_pipeline/run_stage_parallel.sh | 36 ++++++ 7 files changed, 462 insertions(+) create mode 100755 run_pipeline/bbcompsep_parallel.sh create mode 100755 run_pipeline/bbpowerspecter_parallel.sh create mode 100755 run_pipeline/bbpowersummarizer_parallel.sh create mode 100644 run_pipeline/config.yml create mode 100644 run_pipeline/config_hybrid.yml create mode 100755 run_pipeline/prepare_directory_parallel.sh create mode 100755 run_pipeline/run_stage_parallel.sh diff --git a/run_pipeline/bbcompsep_parallel.sh b/run_pipeline/bbcompsep_parallel.sh new file mode 100755 index 0000000..0ea5d67 --- /dev/null +++ b/run_pipeline/bbcompsep_parallel.sh @@ -0,0 +1,37 @@ +#!/bin/bash + +# Reads in seeds to sample +IFS=$'\r\n' GLOBIGNORE='*' command eval 'seeds=($(cat ${dir}/seeds.txt))' + +# Runs pipeline 'nseeds' times in parallel. +for k in $( eval echo {0..$( expr $nseeds / $(( $SLURM_NTASKS + 1 )) )} ) +do + tmp1=$(( $SLURM_PROCID + $k*$SLURM_NTASKS )) + tmp2=$(( $tmp1 < $nseeds ? $tmp1 : $nseeds )) # numerical min + export data_seed=${seeds[$tmp2]} + + export simdir=${dir}${data_seed}/ + export output=$simdir + export config=$output"config.yml" + cp $dir"config.yml" $output"config.yml" + + sed -i "s|ODIR|${dir}|" $config + sed -i "s/0X/${data_seed}/" $config + sed -i "s/0Y/${data_seed}/" $config + + # Only run seed if param_chains file is not there + [ -f ${output}/param_chains.npz ] || + echo "data seed ${data_seed}" + + # Runs BBCompSep + python3 -m bbpower BBCompSep \ + --cells_coadded=${output}cells_coadded.sacc \ + --cells_noise=${output}cells_noise.sacc \ + --cells_fiducial=${simdir}cls_fid_residual_masked_${data_seed}.fits \ + --cells_coadded_cov=${dirall}/cells_coadded.fits \ #this is the precomputed file that contains the covariance + --param_chains=${output}param_chains.npz \ + --config_copy=${output}config_copy.yml \ + --config=$config + + echo "data seed ${data_seed} done at task ${SLURM_PROCID}" +done diff --git a/run_pipeline/bbpowerspecter_parallel.sh b/run_pipeline/bbpowerspecter_parallel.sh new file mode 100755 index 0000000..b4ed51a --- /dev/null +++ b/run_pipeline/bbpowerspecter_parallel.sh @@ -0,0 +1,54 @@ +#!/bin/bash + +# Reads in seeds to sample +IFS=$'\r\n' GLOBIGNORE='*' command eval 'seeds=($(cat ${mdir}/seeds.txt))' + +# Runs pipeline 'nseeds' times in parallel. +for k in $( eval echo {0..$( expr $nseeds / $(( $SLURM_NTASKS + 1 )) )} ) +do + tmp1=$(( $SLURM_PROCID + $k*$SLURM_NTASKS )) + tmp2=$(( $tmp1 < $nseeds ? $tmp1 : $nseeds )) # numerical min + export data_seed=${seeds[$tmp2]} + + echo $data_seed + + export simdir=${mdir}${data_seed}/ + export output=$simdir + + export config=$mdir"config_hybrid.yml" + + # Compile txt file with all Cls in + rm -rf ${output}cells_all_sims.txt + echo ${output}"cells_all_splits_sim0.fits">> ${output}cells_all_sims.txt + for nsim in {1..499} + do + echo ${dirall}"cells_all_splits_sim${nsim}.fits" >> ${output}cells_all_sims.txt + done + + # Creating full simulations list" + rm -rf ${output}list_sims.txt + echo $simdir >> ${output}list_sims.txt + for nsim in {001..499} + do + sim=/dir/to/sim/0${nsim}/ + echo $sim >> ${output}list_sims.txt + done + + # Only run seed if param_chains file is not there + [ -f ${output}/param_chains.npz ] || + echo "data seed ${data_seed}" + + # Runs BBPowerSpecter + python3 -m bbpower BBPowerSpecter \ + --splits_list=${simdir}splits_list.txt \ + --masks_apodized=./test_hybrid/masks_SAT_ns512.fits \ + --bandpasses_list=./examples/data/bpass_list.txt \ + --sims_list=${output}list_sims.txt \ + --beams_list=./examples/data/beams_list.txt \ + --cells_all_splits=${output}cells_all_splits.sacc \ + --cells_all_sims=${output}cells_all_sims.txt \ + --mcm=${output}mcm.dum \ + --config=$config + + echo "data seed ${data_seed} done at task ${SLURM_PROCID}" +done diff --git a/run_pipeline/bbpowersummarizer_parallel.sh b/run_pipeline/bbpowersummarizer_parallel.sh new file mode 100755 index 0000000..f4f405b --- /dev/null +++ b/run_pipeline/bbpowersummarizer_parallel.sh @@ -0,0 +1,37 @@ +#!/bin/bash + +# Reads in seeds to sample +IFS=$'\r\n' GLOBIGNORE='*' command eval 'seeds=($(cat ${dir}/seeds.txt))' + +# Runs pipeline 'nseeds' times in parallel. +for k in $( eval echo {0..$( expr $nseeds / $(( $SLURM_NTASKS + 1 )) )} ) +do + tmp1=$(( $SLURM_PROCID + $k*$SLURM_NTASKS )) + tmp2=$(( $tmp1 < $nseeds ? $tmp1 : $nseeds )) # numerical min + export data_seed=${seeds[$tmp2]} + + echo $data_seed + + export simdir=${dir}${data_seed}/ + export output=$simdir + export covdir=$dirall + cp $dir"config.yml" $output"config.yml" + export config=$output"config_hybrid.yml" + + #sed -i -r "s/ZZ/block_diagonal/" $config #Calculate covariance + sed -i -r "s|ZZ|$covdir|" $config #Use pre-computed covariance + + # Runs BBPowerSummarizer + python3 -m bbpower BBPowerSummarizer \ + --splits_list=${simdir}splits_list.txt \ + --bandpasses_list=./examples/data/bpass_list.txt \ + --cells_all_splits=${output}cells_all_splits.sacc \ + --cells_all_sims=${output}cells_all_sims.txt \ + --cells_coadded_total=${output}cells_coadded_total.sacc \ + --cells_coadded=${output}cells_coadded.sacc \ + --cells_noise=${output}cells_noise.sacc \ + --cells_null=${output}cells_null.sacc \ + --config=${config} + + echo "data seed ${data_seed} done at task ${SLURM_PROCID}" +done diff --git a/run_pipeline/config.yml b/run_pipeline/config.yml new file mode 100644 index 0000000..eb9a259 --- /dev/null +++ b/run_pipeline/config.yml @@ -0,0 +1,137 @@ +# These parameters are accessible to all stages +global: + nside: 512 + compute_dell: True + +BBPowerSpecter: + # Bandpower definition + bpw_edges: "/global/homes/s/susannaz/BBPower/examples/data/bpw_edges.txt" + purify_B: True + n_iter : 3 + # Activates parallelization with mpi4py + parallel: False + +BBPowerSummarizer: + # Covariance types + nulls_covar_type: "diagonal" + nulls_covar_diag_order: 0 + data_covar_type: "block_diagonal" + data_covar_diag_order: 0 + +BBCompSep: + # Sampler type (choose 'emcee', 'polychord', 'maximum_likelihood', + # 'single_point' or 'timing') + sampler: 'emcee' + # If you chose polychord: + nlive: 50 + nrepeat: 50 + # If you chose emcee: + nwalkers: 128 + n_iters: 2000 + # Likelihood type (choose 'chi2' or 'h&l') + likelihood_type: 'chi2' + # What is the starting point? + r_init: 1.e-3 + # Which polarization channels do you want to include? + pol_channels: ['B'] + # Scale cuts (will apply to all frequencies) + l_min: 30 + l_max: 300 + diff: False + + # CMB model + cmb_model: + # Template power spectrum. Should contained the lensed power spectra + # with r=0 and r=1 respectively. + cmb_templates: ['/global/homes/s/susannaz/BBPower/examples/data/camb_lens_nobb_nico.dat', + '/global/homes/s/susannaz/BBPower/examples/data/camb_lens_r1_nico.dat'] + # Free parameters + params: + # tensor-to-scalar ratio + # See below for the meaning of the different elements in the list. + r_tensor: ['r_tensor', 'tophat', [-0.1, 0.00, 0.1]] + # Lensing amplitude + A_lens: ['A_lens', 'tophat', [0.00, 1.0, 2.00]] + + # Foreground model + fg_model: + use_moments: False + # Add one section per component. They should be called `component_X`, + # starting with X=1 + component_1: + # Name for this component + name: Dust + # Type of SED. Should be one of the classes stored in fgbuster.components + # https://github.com/fgbuster/fgbuster/blob/master/fgbuster/component_model.py + sed: Dust + # Type of power spectra for all possible polarization channel combinations. + # Any combinations not added here will be assumed to be zero. + # The names should be one of the classes in bbpower/fgcls.py + cl: + EE: ClPowerLaw + BB: ClPowerLaw + # Parameters of the SED + sed_parameters: + # The key can be anything you want, but each parameter in the model + # must have a different name. + # The first item in the list is the name of the parameter used by fgbuster + # The second item is the type of prior. The last item are the numbers + # necessary to define the prior. They should be: + # - Gaussian: [mean,sigma] + # - tophat: [lower edge, start value, upper edge] + # - fixed: [parameter value] + # nu0-type parameters can only be fixed. + beta_d: ['beta_d', 'Gaussian', [1.54, 0.2]] + temp_d: ['temp', 'fixed', [20.]] + nu0_d: ['nu0', 'fixed', [353.]] + #nu0_d: ['nu0', 'fixed', [220.]] + cl_parameters: + # Same for power spectrum parameters + # (broken down by polarization channel combinations) + EE: + #amp_d_ee: ['amp', 'tophat', [0., 10., "inf"]] + amp_d_ee: ['amp', 'tophat', [40., 56., 80.]] + alpha_d_ee: ['alpha', 'tophat', [-1., -0.32, 0.]] + l0_d_ee: ['ell0', 'fixed', [80.]] + BB: + #amp_d_bb: ['amp', 'tophat', [0., 0.5, "inf"]] + amp_d_bb: ['amp', 'tophat', [0., 10., "inf"]] + #alpha_d_bb: ['alpha', 'tophat', [-1., -0.16, 0.5]] + alpha_d_bb: ['alpha', 'tophat', [-1, -0.42, 0.]] + #alpha_d_bb: ['alpha', 'tophat', [-1, -0.2, 0.]] + l0_d_bb: ['ell0', 'fixed', [80.]] + # If this component should be correlated with any other, list them here + cross: + # In this case the list should contain: + # [component name, prior type, prior parameters] + #epsilon_ds: ['component_2', 'fixed', [0.]] + epsilon_ds: ['component_2', 'tophat', [-0.001, 0., 0.001]] + + component_2: + name: Synchrotron + sed: Synchrotron + cl: + EE: ClPowerLaw + BB: ClPowerLaw + sed_parameters: + beta_s: ['beta_pl', 'Gaussian', [-3.0, 0.3]] + nu0_s: ['nu0', 'fixed', [23.]] + #nu0_s: ['nu0', 'fixed', [40.]] + cl_parameters: + EE: + #amp_s_ee: ['amp', 'tophat', [0., 4., "inf"]] + amp_s_ee: ['amp', 'tophat', [5., 9., 15.]] + alpha_s_ee: ['alpha', 'tophat', [-2., -0.7, 0.]] + l0_s_ee: ['ell0', 'fixed', [80.]] + BB: + amp_s_bb: ['amp', 'tophat', [0., 0.05, "inf"]] + alpha_s_bb: ['alpha', 'tophat', [-2., -0.93, 0.]] + #alpha_s_bb: ['alpha', 'tophat', [-2., -1.5, 0.]] + l0_s_bb: ['ell0', 'fixed', [80.]] + +BBPlotter: + lmax_plot: 300 + plot_coadded_total: False + plot_noise: False + plot_nulls: False + plot_likelihood: True diff --git a/run_pipeline/config_hybrid.yml b/run_pipeline/config_hybrid.yml new file mode 100644 index 0000000..8c5c00d --- /dev/null +++ b/run_pipeline/config_hybrid.yml @@ -0,0 +1,131 @@ +# These parameters are accessible to all stages +global: + nside: 512 + compute_dell: True + +BBPowerSpecter: + # Bandpower definition + bpw_edges: "/global/homes/s/susannaz/BBPower/examples/data/bpw_edges.txt" + purify_B: True + n_iter : 3 + # Activates parallelization with mpi4py + parallel: False + +BBPowerSummarizer: + # Covariance types + nulls_covar_type: "diagonal" + nulls_covar_diag_order: 0 + data_covar_type: "ZZ" + data_covar_diag_order: 0 + +BBCompSep: + # Sampler type (choose 'emcee', 'polychord', 'maximum_likelihood', + # 'single_point' or 'timing') + sampler: 'emcee' + # If you chose polychord: + nlive: 50 + nrepeat: 50 + # If you chose emcee: + nwalkers: 128 + n_iters: 2000 + # Likelihood type (choose 'chi2' or 'h&l') + likelihood_type: 'chi2' + # What is the starting point? + r_init: 1.e-3 + # Which polarization channels do you want to include? + pol_channels: ['B'] + # Scale cuts (will apply to all frequencies) + l_min: 30 + l_max: 300 + diff: True + + # CMB model + cmb_model: + # Template power spectrum. Should contained the lensed power spectra + # with r=0 and r=1 respectively. + cmb_templates: ['/global/homes/s/susannaz/BBPower/examples/data/camb_lens_nobb_nico.dat', + '/global/homes/s/susannaz/BBPower/examples/data/camb_lens_r1_nico.dat'] + # Free parameters + params: + # tensor-to-scalar ratio + # See below for the meaning of the different elements in the list. + r_tensor: ['r_tensor', 'tophat', [-0.1, 0.00, 0.1]] + # Lensing amplitude + A_lens: ['A_lens', 'tophat', [0.00, 1.0, 2.00]] + + # If diff: True + resid_seds: ODIR/0X/masked_hybrid_params_0Y.npz + + # Foreground model + fg_model: + use_moments: False + # Add one section per component. They should be called `component_X`, + # starting with X=1 + component_1: + # Name for this component + name: Dust + # Type of SED. Should be one of the classes stored in fgbuster.components + # https://github.com/fgbuster/fgbuster/blob/master/fgbuster/component_model.py + sed: Dust + # Type of power spectra for all possible polarization channel combinations. + # Any combinations not added here will be assumed to be zero. + # The names should be one of the classes in bbpower/fgcls.py + cl: + EE: ClPowerLaw + BB: ClPowerLaw + # Parameters of the SED + sed_parameters: + # The key can be anything you want, but each parameter in the model + # must have a different name. + # The first item in the list is the name of the parameter used by fgbuster + # The second item is the type of prior. The last item are the numbers + # necessary to define the prior. They should be: + # - Gaussian: [mean,sigma] + # - tophat: [lower edge, start value, upper edge] + # - fixed: [parameter value] + # nu0-type parameters can only be fixed. + beta_d: ['beta_d', 'Gaussian', [1.54, 0.2]] + temp_d: ['temp', 'fixed', [20.]] + nu0_d: ['nu0', 'fixed', [353.]] + cl_parameters: + # Same for power spectrum parameters + # (broken down by polarization channel combinations) + EE: + amp_d_ee: ['amp', 'tophat', [40., 56., 80.]] + alpha_d_ee: ['alpha', 'tophat', [-1., -0.32, 0.]] + l0_d_ee: ['ell0', 'fixed', [80.]] + BB: + amp_d_bb: ['amp', 'tophat', [0., 10., "inf"]] + alpha_d_bb: ['alpha', 'tophat', [-1, -0.42, 0.]] + l0_d_bb: ['ell0', 'fixed', [80.]] + # If this component should be correlated with any other, list them here + cross: + # In this case the list should contain: + # [component name, prior type, prior parameters] + epsilon_ds: ['component_2', 'tophat', [-0.001, 0., 0.001]] + + component_2: + name: Synchrotron + sed: Synchrotron + cl: + EE: ClPowerLaw + BB: ClPowerLaw + sed_parameters: + beta_s: ['beta_pl', 'Gaussian', [-3.0, 0.3]] + nu0_s: ['nu0', 'fixed', [23.]] + cl_parameters: + EE: + amp_s_ee: ['amp', 'tophat', [5., 9., 15.]] + alpha_s_ee: ['alpha', 'tophat', [-2., -0.7, 0.]] + l0_s_ee: ['ell0', 'fixed', [80.]] + BB: + amp_s_bb: ['amp', 'tophat', [0., 0.05, "inf"]] + alpha_s_bb: ['alpha', 'tophat', [-2., -0.93, 0.]] + l0_s_bb: ['ell0', 'fixed', [80.]] + +BBPlotter: + lmax_plot: 300 + plot_coadded_total: False + plot_noise: False + plot_nulls: False + plot_likelihood: True diff --git a/run_pipeline/prepare_directory_parallel.sh b/run_pipeline/prepare_directory_parallel.sh new file mode 100755 index 0000000..dfa04a0 --- /dev/null +++ b/run_pipeline/prepare_directory_parallel.sh @@ -0,0 +1,30 @@ +#!/bin/bash + +# Reads in seeds to sample +IFS=$'\r\n' GLOBIGNORE='*' command eval 'seeds=($(cat ${dir}/seeds.txt))' + +# Runs pipeline 'nseeds' times in parallel. +for k in $( eval echo {0..$( expr $nseeds / $(( $SLURM_NTASKS + 1 )) )} ) +do + tmp1=$(( $SLURM_PROCID + $k*$SLURM_NTASKS )) + tmp2=$(( $tmp1 < $nseeds ? $tmp1 : $nseeds )) # numerical min + export data_seed=${seeds[$tmp2]} + + export simdir=${dir}${data_seed}/ + + # Creating simulated residual directory if not existent + mkdir -p $simdir + cp ${simdir}masked_residualmaps_${data_seed}_split0.fits ${simdir}SO_SAT_obs_map_split_1of4.fits + cp ${simdir}masked_residualmaps_${data_seed}_split1.fits ${simdir}SO_SAT_obs_map_split_2of4.fits + cp ${simdir}masked_residualmaps_${data_seed}_split2.fits ${simdir}SO_SAT_obs_map_split_3of4.fits + cp ${simdir}masked_residualmaps_${data_seed}_split3.fits ${simdir}SO_SAT_obs_map_split_4of4.fits + + # Creating splits list + rm -rf ${simdir}splits_list.txt + echo ${simdir}SO_SAT_obs_map_split_1of4.fits >> ${simdir}splits_list.txt + echo ${simdir}SO_SAT_obs_map_split_2of4.fits >> ${simdir}splits_list.txt + echo ${simdir}SO_SAT_obs_map_split_3of4.fits >> ${simdir}splits_list.txt + echo ${simdir}SO_SAT_obs_map_split_4of4.fits >> ${simdir}splits_list.txt + + echo "prepared directory ${data_seed} at task ${SLURM_PROCID}" +done diff --git a/run_pipeline/run_stage_parallel.sh b/run_pipeline/run_stage_parallel.sh new file mode 100755 index 0000000..a4c484a --- /dev/null +++ b/run_pipeline/run_stage_parallel.sh @@ -0,0 +1,36 @@ +#!/bin/bash + +export HDF5_USE_FILE_LOCKING=FALSE + +## Utils +export idir="/input/dir/" +export dir="/output/dir/" +export dirall="/dir/containing/covariance/" +mkdir -p ${dir} + +nside=512 + +export OMP_NUM_THREADS=1 + +## Config +export cf=/dir/config_file.yml + +# Set data seeds +seeds=( $(seq -f "%04g" 0000 0005 ) ) #run for 5 simulations +printf "%s\n" "${seeds[@]}" > "${dir}/seeds.txt" +export nseeds=${#seeds[@]} + +# Set number of CPUs +ncpus=$(( $SLURM_CPUS_ON_NODE * $SLURM_JOB_NUM_NODES )) +ntasks=$(( $nseeds < $ncpus ? $nseeds : $ncpus )) +echo "Create ${ntasks} tasks" + +# Run each stage +source activate testenv # This loads custom programming environment +# +cp $cf $dir +# Uncomment for stage to be run in queue +srun "--ntasks=${ntasks}" "/global/homes/s/susannaz/BBPower/test_hybrid/prepare_directory_parallel.sh" +#srun "--ntasks=${ntasks}" "/global/homes/s/susannaz/BBPower/test_hybrid/bbpowerspecter_parallel.sh" +#srun "--ntasks=${ntasks}" "/global/homes/s/susannaz/BBPower/test_hybrid/bbpowersummarizer_parallel.sh" +#srun "--ntasks=${ntasks}" "/global/homes/s/susannaz/BBPower/test_hybrid/bbcompsep_parallel.sh"