Quadratic Estimator Module
Quadratic Estimator for Cosmic Birefringence Module
This module implements quadratic estimator (QE) methods for reconstructing cosmic birefringence rotation angle fields from CMB polarization observations.
The module provides: - Wiener filtering of E and B modes using C⁻¹ filters - Quadratic estimator reconstruction of rotation angle alms - N0 and RDN0 bias estimation - Mean field subtraction - Likelihood analysis for cosmic birefringence amplitude
Classes
- FilterEB
Implements C⁻¹ Wiener filtering for E and B mode polarization using the curvedsky library. Handles component-separated maps with noise and beam deconvolution.
- QE
Quadratic estimator class for cosmic birefringence reconstruction. Computes rotation angle power spectra with various bias correction methods (analytical N0, realization-dependent N0, mean field).
Simulation Index Structure (sim_config is required):
From sim_config={‘set1’: 400, ‘reuse_last’: 100} and nsims_mf=100: - stat_index: 0-299 (set1 - nsims_mf) - Statistics and OCL computation - mf_index: 300-399 (last nsims_mf from set1) - Mean field simulations - vary_index: 300-399 (last reuse_last from set1) - Varying alpha (CMB mode=’vary’) - const_index: 400-499 (reuse_last sims) - Constant alpha (CMB mode=’constant’) - null_index: 500-599 (reuse_last sims) - Null alpha (CMB mode=’null’)
Note: When nsims_mf equals reuse_last, mf_index and vary_index are identical.
Bias Estimation: - MCN0(‘stat’): Uses stat_index (0-299) - MCN0(‘vary’): Uses vary_index (300-399) - MCN0(‘const’): Uses const_index (400-499) - N1 = MCN0(‘const’) - MCN0(‘vary’) - Nlens: Uses null_index (500-599) and MCN0(‘vary’)
- AcbLikelihood
Likelihood analysis class for constraining the cosmic birefringence amplitude A_CB from reconstructed rotation angle power spectra.
Algorithm
The quadratic estimator exploits EB correlations induced by cosmic birefringence:
Apply C⁻¹ Wiener filter to observed E and B maps
Compute quadratic combination: qlm ∝ Ẽlm B̃*lm
Apply normalization from Fisher information
Estimate and subtract biases (N0, mean field)
Compare to theoretical prediction for A_CB
Requirements
Requires the curvedsky library for efficient curved-sky operations. Install from: https://github.com/toshiyan/cmblensplus
Example
from cobi.quest import FilterEB, QE from cobi.simulation import LATsky, Mask
# Set up sky and mask sky = LATsky(libdir, nside=512) mask = Mask(libdir, nside=512, mask_str=’LAT’)
# Initialize Wiener filter filter_eb = FilterEB(
sky=sky, mask=mask, lmax=3000, fwhm=2.0, # arcmin sht_backend=’ducc0’
)
# Create quadratic estimator qe = QE(
filter=filter_eb, lmin=30, lmax=300, recon_lmax=100, norm_nsim=100
)
# Reconstruct rotation angle for simulation 0 qlm = qe.qlm(idx=0)
# Get bias-corrected power spectrum cl_aa = qe.qcl(idx=0, rm_bias=True, rdn0=True)
# Run likelihood analysis likelihood = AcbLikelihood(qe, lmin=2, lmax=50) samples = likelihood.samples(nwalkers=100, nsamples=2000)
Notes
The curvedsky library must be installed to use this module. MPI parallelization is supported for RDN0 computation.
- class FilterEB(sky, mask, lmax, fwhm=2, sht_backend='healpy')[source]
Bases:
object- property Bl
- convolved_EB(idx, split=0)[source]
convolve the component separated map with the beam
- Parameters:
idx (int :
indexofthe simulation)
- NL(idx, split=0)[source]
array manipulation of noise spectra obtained by ILC weight for the filtering process
- QU(idx, split=0)[source]
deconvolve the beam from the QU map
- Parameters:
idx (int :
indexofthe simulation)
- class QE(filter, lmin, lmax, recon_lmax, nsims_mf=100, nlb=10, lmax_bin=1024)[source]
Bases:
object- property cl_aa
- RDN0_mpi(idx)[source]
MPI-parallel version of RDN0(…) using mpi4py.
Parallelizes the 100 Monte-Carlo iterations; safe for any number of ranks. Only rank 0 touches the on-disk cache.
- N0_sim(idx, which='vary', debug=False)[source]
Calculate the N0 bias from the Reconstructed potential using filtered Fields with different CMB fields. If E modes is from ith simulation then B modes is from (i+1)th simulation
idx: int : index of the N0 which: str : ‘stat’, ‘vary’, or ‘const’ to select which index range to use debug: bool : if True, print the indices used for computation and return None
Index ranges and wrapping: - ‘stat’: stat_index, wraps within stat range - ‘vary’: vary_index, wraps within vary range - ‘const’: const_index, wraps within const range
Requires sim_config to be set, or manually set the corresponding index array.
- MCN0(which='vary')[source]
Monte Carlo average of N0_sim over specified index range
- which: str‘stat’, ‘vary’, or ‘const’ to select which index range to use
‘stat’ uses stat_index ‘vary’ uses vary_index ‘const’ uses const_index
Requires sim_config to be set, or manually set the corresponding index arrays. Note: vary_index overlaps with mf_index only when nsims_mf equals reuse_last.
- N1(binned=False)[source]
N1 bias: difference between MCN0 for constant and varying alpha N1 = MCN0(‘const’) - MCN0(‘vary’)
- class CrossQE(filter, lmin, lmax, recon_lmax, nsims_mf=100, nlb=10, lmax_bin=1024)[source]
Bases:
object- property cl_aa
- precompute_split_ocl(splits=(1, 2, 3, 4))[source]
Computes split-level total spectra: oclE[s] = ClEE_lensed + <N_EE_split_s> oclB[s] = ClBB_lensed + <N_BB_split_s> Saves dict {s: (oclE, oclB)}.
- precompute_pair_norms(pairs=((1, 2), (3, 4), (1, 3), (2, 4), (1, 4), (2, 3)))[source]
Precomputes normalization arrays for each split pair.
Returns dict: norms[(i,j)] = {“EiBj”: norm_ij, “EjBi”: norm_ji} where norm_ij normalizes Q(E^i, B^j).
- mean_field_sim(idx, splits=(1, 2, 3, 4))[source]
Cross-only lensing power estimate: average over cross-spectra of reconstructions from disjoint split pairs. For 4 splits, uses the 3 disjoint pairings.
- qcl_cross_only(idx, splits=(1, 2, 3, 4))[source]
Cross-only lensing power estimate: average over cross-spectra of reconstructions from disjoint split pairs. For 4 splits, uses the 3 disjoint pairings.
- Nlens(binned=False)[source]
Lensing bias: average qcl on null_index minus MCN0(‘vary’) Nlens = <qcl(null_alpha)> - MCN0(‘vary’)
- N1(binned=False)[source]
N1 bias: difference between MCN0 for constant and varying alpha N1 = MCN0(‘const’) - MCN0(‘vary’)
- qcl_stat(n0=None, n1=False, mf=False, nlens=False, binned=False)[source]
Compute statistics of qcl over stat_sims with various bias corrections.
- Parameters:
- Returns:
cl – Array of shape (nsims, nlmax+1) or (nsims, nbins) if binned
- Return type:
array
- N0_sim(idx, which='vary', debug=False)[source]
Calculate the N0 bias from the Reconstructed potential using filtered Fields with different CMB fields. If E modes is from ith simulation then B modes is from (i+1)th simulation
idx: int : index of the N0 which: str : ‘stat’, ‘vary’, or ‘const’ to select which index range to use debug: bool : if True, print the indices used for computation and return None
Index ranges and wrapping: - ‘stat’: stat_index, wraps within stat range - ‘vary’: vary_index, wraps within vary range - ‘const’: const_index, wraps within const range
Requires sim_config to be set, or manually set the corresponding index array.
- RDN0(idx, navg=100)[source]
Realization-dependent N0 (RDN0) compatible with qcl_cross_only() (symmetric estimator).
Uses the same Eq.(43)-style structure: RDN0 = Cx(ds,ds) + Cx(ds,sd) + Cx(sd,ds) + Cx(sd,sd) - Cx(ss0,ss0) - Cx(ss0,s0s) where Cx is the cross-only spectrum operator (12×34+13×24+14×23)/3.
Returns: array C_L[0..recon_lmax]
- RDN0_mpi(idx, navg=100)[source]
MPI-parallel cross-only RDN0 for SYM estimator (compatible with qcl_cross_only()).
Implements the Eq.(43)-style combination (same algebra as your rot version), but with SYM estimator: RDN0 = Cx(ds,ds)+Cx(ds,sd)+Cx(sd,ds)+Cx(sd,sd) - Cx(ss0,ss0) - Cx(ss0,s0s)
Parallelizes navg Monte-Carlo draws over (s, s0) pairs. Only rank 0 reads/writes the cache.
Filtering Class
- class FilterEB(sky, mask, lmax, fwhm=2, sht_backend='healpy')[source]
-
- property Bl
- convolved_EB(idx, split=0)[source]
convolve the component separated map with the beam
- Parameters:
idx (int :
indexofthe simulation)
- NL(idx, split=0)[source]
array manipulation of noise spectra obtained by ILC weight for the filtering process
- QU(idx, split=0)[source]
deconvolve the beam from the QU map
- Parameters:
idx (int :
indexofthe simulation)
Quadratic Estimator Class
- class QE(filter, lmin, lmax, recon_lmax, nsims_mf=100, nlb=10, lmax_bin=1024)[source]
-
- property cl_aa
- RDN0_mpi(idx)[source]
MPI-parallel version of RDN0(…) using mpi4py.
Parallelizes the 100 Monte-Carlo iterations; safe for any number of ranks. Only rank 0 touches the on-disk cache.
- N0_sim(idx, which='vary', debug=False)[source]
Calculate the N0 bias from the Reconstructed potential using filtered Fields with different CMB fields. If E modes is from ith simulation then B modes is from (i+1)th simulation
idx: int : index of the N0 which: str : ‘stat’, ‘vary’, or ‘const’ to select which index range to use debug: bool : if True, print the indices used for computation and return None
Index ranges and wrapping: - ‘stat’: stat_index, wraps within stat range - ‘vary’: vary_index, wraps within vary range - ‘const’: const_index, wraps within const range
Requires sim_config to be set, or manually set the corresponding index array.
- MCN0(which='vary')[source]
Monte Carlo average of N0_sim over specified index range
- which: str‘stat’, ‘vary’, or ‘const’ to select which index range to use
‘stat’ uses stat_index ‘vary’ uses vary_index ‘const’ uses const_index
Requires sim_config to be set, or manually set the corresponding index arrays. Note: vary_index overlaps with mf_index only when nsims_mf equals reuse_last.
- N1(binned=False)[source]
N1 bias: difference between MCN0 for constant and varying alpha N1 = MCN0(‘const’) - MCN0(‘vary’)