# Bayesian signal reconstruction with Gaussian Random Fields (Wiener Filtering) - CMB example¶

Florent Leclercq,
Institute of Cosmology and Gravitation, University of Portsmouth,
[email protected]

In [1]:
import numpy as np
import scipy.linalg
import os.path
import matplotlib.pyplot as plt
from classy import Class
import healpy as hp
%matplotlib inline


## Setup Cls using CLASS¶

In [2]:
# Define cosmology (what is not specified will be set to CLASS default parameters)
paramsLCDM = {
'h':0.702,
'n_s':0.9619,
'Omega_b':0.045,
'Omega_cdm':0.272-0.045,
'output':'tCl lCl',
'l_max_scalars': 1000,
}
# Create an instance of the CLASS wrapper
cosmoLCDM = Class()
# Set the parameters to the cosmological code
cosmoLCDM.set(paramsLCDM)
# Run the whole code.
cosmoLCDM.compute()

In [3]:
# Compute the power spectrum for the current set of cosmological parameters
res=cosmoLCDM.raw_cl(1000)
ell=res['ell']
Cl=res['tt']

In [4]:
# Clean CLASS
cosmoLCDM.struct_cleanup()

In [5]:
# Define a function of l from the arrays
from scipy.interpolate import InterpolatedUnivariateSpline
Clfunc = InterpolatedUnivariateSpline(ell, Cl, k=2)

In [6]:
# Plot l*(l+1)*Cl/(2*pi)
plt.xlabel("$\ell$",size=18)
plt.ylabel("$\ell(\ell+1) C_\ell$",size=18)
plt.loglog(ell,Cl*(ell*(ell+1))/(2*np.pi),label='CLASS output')
plt.loglog(np.arange(0,999,2),Clfunc(np.arange(0,999,2))*(np.arange(0,999,2)*(np.arange(0,999,2)+1))/(2*np.pi)
,label='Spline function')
plt.legend(loc='best')
plt.show()

In [7]:
# Setup notations
nside=32
npix=hp.pixelfunc.nside2npix(nside)
lmax=min(len(Cl)-1, 3*nside-1)
nalm=hp.Alm.getsize(lmax)


Above $\mathrm{nside}=64$ the covariance matrix can't be stored...

## Setup signal covariance matrix¶

In [8]:
# Setup covariance matrix in harmonic space
covarFourier=np.zeros(nalm,dtype=complex)
for l in xrange(lmax+1):
for m in xrange(lmax+1):
idx=hp.Alm.getidx(lmax,l,m)
covarFourier[idx]=Clfunc(l)

In [9]:
# Plot diagonal of the harmonic-space covariance matrix
plt.xlabel("$a_{\ell m}$",size=18)
plt.ylabel("$C_{\ell}$",size=18)
plt.loglog(np.arange(nalm),covarFourier.real)
plt.show()


Each row of the covariance matrix in pixel space is given by: $C_i = \mathcal{F}^{-1} \left[ C_{a_{\ell m}a_{\ell m}}a_{\ell m}(i) \right]$ where $a_{\ell m}(i) = \mathcal{F} \left[ e(i) \right]$.
$e(i) = [\delta^{i,j}]_{0 \leq j < N_\mathrm{pix}}^\mathrm{T}$ a the unit vector of the pixel space basis.
$\mathcal{F}$ and $\mathcal{F}^{-1}$ are transformations from harmonic space to pixel space: \begin{eqnarray} m(n) & = & \sum_{\ell m} a_{\ell m} Y_{\ell m}(n)\\ a_{\ell m} & = & \sum_n m(n) Y_{\ell m}^*(n). \end{eqnarray}

In [11]:
%%capture
if os.path.isfile("data/sqrtsignalcovarPix.npy"):
# Load precomputed sqrt covariance matrix in pixel space
else:
# Setup covariance matrix in pixel space
covarFouriermat=np.diagflat(covarFourier)
signalcovarPix=np.zeros((npix,npix))
for i in xrange(npix):
ei=np.zeros(npix)
ei[i]=1.
alm=hp.sphtfunc.map2alm(ei)
Calm=covarFouriermat.dot(alm)
signalcovarPix[i]=hp.sphtfunc.alm2map(Calm,nside)
signalcovarPix=signalcovarPix.T
sqrtsignalcovarPix=scipy.linalg.sqrtm(signalcovarPix)
del covarFouriermat, signalcovarPix
np.save("data/sqrtsignalcovarPix",signalcovarPix)


The signal covariance matrix is diagonal in harmonic space, but dense in pixel space...

In [12]:
# Show sqrt covariance matrix in pixel space
i=1000
ei=np.zeros(npix)
ei[i]=1.
hp.mollview(ei, coord='G', title='one pixel', unit='mK')
hp.mollview(sqrtsignalcovarPix.real[i], coord='G', title='one row of the sqrt covariance matrix', unit='mK')
plt.show()


## Setup mask and noise covariance¶

In [13]:
%%capture
# Simplify to binary

In [14]:
hp.mollview(tmask_map, coord='G', title='mask', unit='mK')
plt.show()

In [15]:
%%capture
if os.path.isfile("data/invnoisecovar.npy"):
# Load precomputed inverse noise covariance matrix in pixel space
else:
# Setup the noise covariance (infinite in masked regions)
noisepower=1e-12
invnoisecovar=1./noisecovar
np.save("data/invnoisecovar",invnoisecovar)


The noise covariance matrix is assumed diagonal in pixel space.

## Generate mock data¶

In [16]:
%%capture
# Generate the true CMB signal
signal_map, signal_alms = hp.synfast(Cl, nside, alm=True, lmax=lmax, new=True)

In [17]:
hp.mollview(signal_map, coord='G', title='true signal', unit='mK')
plt.show()


Data model: $d=s+n$

In [19]:
# Generate simulated data
normalsim = np.random.normal(0., 1., len(signal_map))
noise_map = normalsim/np.sqrt(invnoisecovar+1e-12)
data_map = signal_map + noise_map
data_map_vis = hp.ma(data_map)

In [20]:
hp.mollview(data_map_vis, coord='G', title='simulated data', unit='mK')
plt.show()


## Setup Wiener filter¶

Covariance of the Wiener Filter: $$\mathrm{Cov}_\mathrm{WF} = (N^{-1}+S^{-1})^{-1} = S^{1/2}(I+S^{1/2}N^{-1}S^{1/2})^{-1}S^{1/2}$$ For applications we only need to keep in memory $\mathrm{Cov}_\mathrm{WF}N^{-1}$ and $\sqrt{\mathrm{Cov}_\mathrm{WF}}$.

In [22]:
%%capture
if os.path.isfile("data/CovWFinvnoisecovarmat.npy") and os.path.isfile("data/sqrtCovWF.npy"):
# Load precomputed covariance of the Wiener filter
del sqrtsignalcovarPix
else:
# Setup covariance of the Wiener filter
M=np.identity(npix)+sqrtsignalcovarPix.dot(invnoisecovarmat).dot(sqrtsignalcovarPix)
CovWF=sqrtsignalcovarPix.dot(np.linalg.inv(M)).dot(sqrtsignalcovarPix)
del sqrtsignalcovarPix, M
# Setup CovWF*N^-1
CovWFinvnoisecovarmat=CovWF.dot(invnoisecovarmat)
np.save("data/CovWFinvnoisecovarmat",CovWFinvnoisecovarmat)
# Setup sqrtCovWF for simulating signals
CovWFsym=(CovWF+CovWF.T)/2
del CovWF
sqrtCovWF=scipy.linalg.sqrtm(CovWFsym)
del CovWFsim
np.save("data/sqrtCovWF",sqrtCovWF)


## Perform signal reconstruction (apply Wiener Filtering)¶

Mean of the Wiener posterior: $$s_\mathrm{WF} = \mathrm{Cov}_\mathrm{WF} N^{-1} d$$

In [24]:
# Wiener filtered map (posterior mean)
sWF=CovWFinvnoisecovarmat.dot(data_map).real

In [25]:
hp.mollview(sWF, coord='G', title='Wiener-filtered data', unit='mK')
plt.show()


## Generate constrained realizations (draw samples from the Wiener posterior)¶

Samples of the Wiener posterior: $$s=s_\mathrm{WF}+\sqrt{C_\mathrm{WF}} \, G(0,1)$$ so that $\left\langle s \right\rangle = s_\mathrm{WF}$ and $\mathrm{Cov}(s) = C_\mathrm{WF}$

In [27]:
# Constrained realizations (samples of the Wiener posterior)
cr1=sqrtCovWF.dot(np.random.normal(0.,1.,npix)).real+sWF
cr2=sqrtCovWF.dot(np.random.normal(0.,1.,npix)).real+sWF
cr3=sqrtCovWF.dot(np.random.normal(0.,1.,npix)).real+sWF

In [28]:
hp.mollview(cr1.real, coord='G', title='One simulated signal', unit='mK')
hp.mollview(cr2.real, coord='G', title='One simulated signal', unit='mK')
hp.mollview(cr3.real, coord='G', title='One simulated signal', unit='mK')
plt.show()