Commit 2f04d329 authored by Matthieu Tristram's avatar Matthieu Tristram

Update symmetrization

parent affa59c9
......@@ -144,6 +144,12 @@ class Xpol(object):
self.lmin = lmin
self.lmax = lmax
if self.verbose: print( "Compute bin operators")
self._p, self._q = bins._bin_operators()
self.lbin = bins.lbin
self.bin_spectra = bins.bin_spectra
#compute cross-spectrum of weight mask
if self.verbose: print( "alm2map mask")
wlmax = 3*nside-1
if mask2 is None:
......@@ -152,17 +158,11 @@ class Xpol(object):
else:
wl = hp.anafast(mask,map2=mask2,lmax=wlmax)
self.mask = mask
self.mask2 = mask2
self.wl = wl
self.mask2 = mask2
if self.verbose: print( "Compute bin operators")
self._p, self._q = bins._bin_operators()
self.lbin = bins.lbin
self.bin_spectra = bins.bin_spectra
if self.verbose: print( "Compute Mll")
mll_binned = self._get_Mbb()
#compute coupling kernels for covariance
if self.verbose: print( "Compute coupling kernels Mll")
mll_binned = self._get_Mbb( wl)
if self.verbose: print( "Inverse Mll")
self.mll_binned_inv = self._inv_mll(mll_binned)
......@@ -329,11 +329,11 @@ class Xpol(object):
return TT_TT, EE_EE, EE_BB, TE_TE, EB_EB
def _get_Mll_blocks(self):
def _get_Mll_blocks(self, wl):
if self.polar:
TT_TT, EE_EE, EE_BB, TE_TE, EB_EB, ier = flib.xpol.mll_blocks_pol(self.lmax, self.wl)
TT_TT, EE_EE, EE_BB, TE_TE, EB_EB, ier = flib.xpol.mll_blocks_pol(self.lmax, wl)
else:
TT_TT, ier = flib.xpol.mll_blocks(self.lmax, self.wl)
TT_TT, ier = flib.xpol.mll_blocks(self.lmax, wl)
if ier > 0:
msg = ['Either L2 < ABS(M2) or L3 < ABS(M3).',
'Either L2+ABS(M2) or L3+ABS(M3) non-integer.'
......@@ -347,11 +347,11 @@ class Xpol(object):
else:
return TT_TT
def _get_Mbb(self):
def _get_Mbb(self, wl):
if self.polar:
TT_TT, EE_EE, EE_BB, TE_TE, EB_EB = self._get_Mll_blocks()
TT_TT, EE_EE, EE_BB, TE_TE, EB_EB = self._get_Mll_blocks( wl)
else:
TT_TT = self._get_Mll_blocks()
TT_TT = self._get_Mll_blocks( wl)
def func(x):
return np.dot(np.dot(self._p, x), self._q)
......@@ -380,7 +380,7 @@ class Xcov(Xpol):
xcov = Xcov(mask, lmin, lmax, delta_ell)
"""
def __init__(self, mask, bins):
def __init__(self, mask, bins, mask2=None, polar=True, verbose=False):
"""
Parameters
----------
......@@ -388,73 +388,98 @@ class Xcov(Xpol):
Mask defining the region of interest (of value True)
"""
Xpol.__init__(self,mask,bins)
Xpol.__init__( self, mask, bins, mask2=mask2, polar=polar, verbose=verbose)
self.dl = bins.dl
self.wl = hp.anafast(mask**2)
self.TT_TT, self.EE_EE, self.EE_BB, self.TE_TE, self.EB_EB = self._get_Mbb()
#compute coupling kernels for covariance
wlmax = 3*nside-1
if mask2 is None:
wl = hp.anafast(mask**2,lmax=wlmax)
else:
wl = hp.anafast(mask**2,map2=mask2**2,lmax=wlmax)
self.TT_TT, self.EE_EE, self.EE_BB, self.TE_TE, self.EB_EB = self._get_Mbb( wl)
def _get_pcl_cov(self, clAC, clBD, clAD, clBC):
"""
Compute covariance matrix for pseudo-spectra <AB,CD>
tags=[TT,EE,BB,TE,ET]
inputs:
cl: 2d array (ntag,nbin)
output:
out: 2d array (5*nbin,5*nbin)
"""
n = len(self.lbin)
out = np.zeros((5*n, 5*n))
nu_l = (2.*self.lbin+1.) * self.dl
#symmetrization
tags = dict( zip( ["TT","EE","BB","TE","ET"], [0,1,2,3,4]))
symcl = lambda cl,t: np.add.outer(cl[tags[t]],cl[tags[t]])/2.
#TT_TT
out[0*n:1*n, 0*n:1*n] = np.sqrt(abs(np.outer(clAC[0],clAC[0])*np.outer(clBD[0],clBD[0])))*self.TT_TT/nu_l + \
np.sqrt(abs(np.outer(clAD[0],clAD[0])*np.outer(clBC[0],clBC[0])))*self.TT_TT/nu_l
out[0*n:1*n, 0*n:1*n] = symcl(clAC,"TT")*symcl(clBD,"TT")*self.TT_TT/nu_l + \
symcl(clAD,"TT")*symcl(clBC,"TT")*self.TT_TT/nu_l
#TT_EE
out[0*n:1*n, 1*n:2*n] = np.sqrt(abs(np.outer(clAC[3],clAC[3])*np.outer(clBD[3],clBD[3])))*self.TT_TT/nu_l + \
np.sqrt(abs(np.outer(clAD[3],clAD[3])*np.outer(clBC[3],clBC[3])))*self.TT_TT/nu_l
out[0*n:1*n, 1*n:2*n] = symcl(clAC,"TE")*symcl(clBD,"TE")*self.TT_TT/nu_l + \
symcl(clAD,"TE")*symcl(clBC,"TE")*self.TT_TT/nu_l
#TT_BB
out[0*n:1*n, 2*n:3*n] = 0.
#TT_TE
out[0*n:1*n, 3*n:4*n]= np.sqrt(abs(np.outer(clAC[0],clAC[0])*np.outer(clBD[3],clBD[3])))*self.TT_TT/nu_l + \
np.sqrt(abs(np.outer(clAD[3],clAD[3])*np.outer(clBC[0],clBC[0])))*self.TT_TT/nu_l
out[0*n:1*n, 3*n:4*n]= symcl(clAC,"TT")*symcl(clBD,"TE")*self.TT_TT/nu_l + \
symcl(clAD,"TE")*symcl(clBC,"TT")*self.TT_TT/nu_l
#TT_ET
out[0*n:1*n, 4*n:5*n] = np.sqrt(abs(np.outer(clAC[3],clAC[3])*np.outer(clBD[0],clBD[0])))*self.TT_TT/nu_l + \
np.sqrt(abs(np.outer(clAD[0],clAD[0])*np.outer(clBC[3],clBC[3])))*self.TT_TT/nu_l
out[0*n:1*n, 4*n:5*n] = symcl(clAC,"TE")*symcl(clBD,"TT")*self.TT_TT/nu_l + \
symcl(clAD,"TT")*symcl(clBC,"TE")*self.TT_TT/nu_l
#EE_TT
out[1*n:2*n, 0*n:1*n] = np.sqrt(abs(np.outer(clAC[4],clAC[4])*np.outer(clBD[4],clBD[4])))*self.TT_TT/nu_l + \
np.sqrt(abs(np.outer(clAD[4],clAD[4])*np.outer(clBC[4],clBC[4])))*self.TT_TT/nu_l
out[1*n:2*n, 0*n:1*n] = symcl(clAC,"ET")*symcl(clBD,"ET")*self.TT_TT/nu_l + \
symcl(clAD,"ET")*symcl(clBC,"ET")*self.TT_TT/nu_l
#EE_EE
out[1*n:2*n, 1*n:2*n] = np.sqrt(abs(np.outer(clAC[1],clAC[1])*np.outer(clBD[1],clBD[1])))*self.EE_EE/nu_l + \
np.sqrt(abs(np.outer(clAD[1],clAD[1])*np.outer(clBC[1],clBC[1])))*self.EE_EE/nu_l
out[1*n:2*n, 1*n:2*n] = symcl(clAC,"EE")*symcl(clBD,"EE")*self.EE_EE/nu_l + \
symcl(clAD,"EE")*symcl(clBC,"EE")*self.EE_EE/nu_l
#EE_BB
out[1*n:2*n, 2*n:3*n] = ( np.sqrt(abs(np.outer(clAC[1],clAC[1])*np.outer(clBD[1],clBD[1]))) + \
np.sqrt(abs(np.outer(clAC[1],clAC[1])*np.outer(clBD[2],clBD[2]))) + \
np.sqrt(abs(np.outer(clAC[2],clAC[2])*np.outer(clBD[1],clBD[1]))) + \
np.sqrt(abs(np.outer(clAC[2],clAC[2])*np.outer(clBD[2],clBD[2]))) )*self.EE_BB/nu_l + \
( np.sqrt(abs(np.outer(clAD[1],clAD[1])*np.outer(clBC[1],clBC[1]))) + \
np.sqrt(abs(np.outer(clAD[1],clAD[1])*np.outer(clBC[2],clBC[2]))) + \
np.sqrt(abs(np.outer(clAD[2],clAD[2])*np.outer(clBC[1],clBC[1]))) + \
np.sqrt(abs(np.outer(clAD[2],clAD[2])*np.outer(clBC[2],clBC[2]))) )*self.EE_BB/nu_l
out[1*n:2*n, 2*n:3*n] = ( symcl(clAC,"EE")*symcl(clBD,"EE") + \
symcl(clAC,"EE")*symcl(clBD,"BB") + \
symcl(clAC,"BB")*symcl(clBD,"EE") + \
symcl(clAC,"BB")*symcl(clBD,"BB") )*self.EE_BB/nu_l + \
( symcl(clAD,"EE")*symcl(clBC,"EE") + \
symcl(clAD,"EE")*symcl(clBC,"BB") + \
symcl(clAD,"BB")*symcl(clBC,"EE") + \
symcl(clAD,"BB")*symcl(clBC,"BB") )*self.EE_BB/nu_l
#EE_TE
out[1*n:2*n, 3*n:4*n] = np.sqrt(abs(np.outer(clAC[4],clAC[4])*np.outer(clBD[1],clBD[1])))*self.TE_TE/nu_l + \
np.sqrt(abs(np.outer(clAD[1],clAD[1])*np.outer(clBC[4],clBC[4])))*self.TE_TE/nu_l
out[1*n:2*n, 3*n:4*n] = symcl(clAC,"ET")*symcl(clBD,"EE")*self.TE_TE/nu_l + \
symcl(clAD,"EE")*symcl(clBC,"ET")*self.TE_TE/nu_l
#EE_ET
out[1*n:2*n, 4*n:5*n] = np.sqrt(abs(np.outer(clAC[1],clAC[1])*np.outer(clBD[4],clBD[4])))*self.TE_TE/nu_l + \
np.sqrt(abs(np.outer(clAD[4],clAD[4])*np.outer(clBC[1],clBC[1])))*self.TE_TE/nu_l
out[1*n:2*n, 4*n:5*n] = symcl(clAC,"EE")*symcl(clBD,"ET")*self.TE_TE/nu_l + \
symcl(clAD,"ET")*symcl(clBC,"EE")*self.TE_TE/nu_l
#BB_TT
out[2*n:3*n, 0*n:1*n] = 0.
#BB_EE
out[2*n:3*n, 1*n:2*n] = np.sqrt(abs(np.outer(clAC[1],clAC[1])*np.outer(clBD[1],clBD[1])))*self.EE_BB/nu_l + \
np.sqrt(abs(np.outer(clAD[1],clAD[1])*np.outer(clBC[1],clBC[1])))*self.EE_BB/nu_l
out[2*n:3*n, 1*n:2*n] = ( symcl(clAC,"BB")*symcl(clBD,"BB") + \
symcl(clAC,"BB")*symcl(clBD,"EE") + \
symcl(clAC,"EE")*symcl(clBD,"BB") + \
symcl(clAC,"EE")*symcl(clBD,"EE") )*self.EE_BB/nu_l + \
( symcl(clAD,"BB")*symcl(clBC,"BB") + \
symcl(clAD,"BB")*symcl(clBC,"EE") + \
symcl(clAD,"EE")*symcl(clBC,"BB") + \
symcl(clAD,"EE")*symcl(clBC,"EE") )*self.EE_BB/nu_l
#BB_BB
out[2*n:3*n, 2*n:3*n] = np.sqrt(abs(np.outer(clAC[2],clAC[2])*np.outer(clBD[2],clBD[2])))*self.EE_EE/nu_l + \
np.sqrt(abs(np.outer(clAD[2],clAD[2])*np.outer(clBC[2],clBC[2])))*self.EE_EE/nu_l
out[2*n:3*n, 2*n:3*n] = symcl(clAC,"BB")*symcl(clBD,"BB")*self.EE_EE/nu_l + \
symcl(clAD,"BB")*symcl(clBC,"BB")*self.EE_EE/nu_l
#BB_TE
out[2*n:3*n, 3*n:4*n] = 0.
......@@ -463,48 +488,59 @@ class Xcov(Xpol):
out[2*n:3*n, 4*n:5*n] = 0.
#TE_TT
out[3*n:4*n, 0*n:1*n] = np.sqrt(abs(np.outer(clAC[0],clAC[0])*np.outer(clBD[4],clBD[4])))*self.TT_TT/nu_l + \
np.sqrt(abs(np.outer(clAD[0],clAD[0])*np.outer(clBC[4],clBC[4])))*self.TT_TT/nu_l
out[3*n:4*n, 0*n:1*n] = symcl(clAC,"TT")*symcl(clBD,"ET")*self.TT_TT/nu_l + \
symcl(clAD,"TT")*symcl(clBC,"ET")*self.TT_TT/nu_l
#TE_EE
out[3*n:4*n, 1*n:2*n] = np.sqrt(abs(np.outer(clAC[3],clAC[3])*np.outer(clBD[1],clBD[1])))*self.TE_TE/nu_l + \
np.sqrt(abs(np.outer(clAD[3],clAD[3])*np.outer(clBC[1],clBC[1])))*self.TE_TE/nu_l
out[3*n:4*n, 1*n:2*n] = symcl(clAC,"TE")*symcl(clBD,"EE")*self.TE_TE/nu_l + \
symcl(clAD,"TE")*symcl(clBC,"EE")*self.TE_TE/nu_l
#TE_BB
out[3*n:4*n, 2*n:3*n] = 0.
#TE_TE
out[3*n:4*n, 3*n:4*n] = np.sqrt(abs(np.outer(clAC[0],clAC[0])*np.outer(clBD[1],clBD[1])))*self.TE_TE/nu_l + \
np.sqrt(abs(np.outer(clAD[3],clAD[3])*np.outer(clBC[4],clBC[4])))*self.TT_TT/nu_l
out[3*n:4*n, 3*n:4*n] = symcl(clAC,"TT")*symcl(clBD,"EE")*self.TE_TE/nu_l + \
symcl(clAD,"TE")*symcl(clBC,"ET")*self.TT_TT/nu_l
#TE_ET
out[3*n:4*n, 4*n:5*n] = np.sqrt(abs(np.outer(clAC[3],clAC[3])*np.outer(clBD[4],clBD[4])))*self.TT_TT/nu_l + \
np.sqrt(abs(np.outer(clAD[0],clAD[0])*np.outer(clBC[1],clBC[1])))*self.TE_TE/nu_l
out[3*n:4*n, 4*n:5*n] = symcl(clAC,"TE")*symcl(clBD,"ET")*self.TT_TT/nu_l + \
symcl(clAD,"TT")*symcl(clBC,"EE")*self.TE_TE/nu_l
#ET_TT
out[4*n:5*n, 0*n:1*n] = np.sqrt(abs(np.outer(clAC[4],clAC[4])*np.outer(clBD[0],clBD[0])))*self.TT_TT/nu_l + \
np.sqrt(abs(np.outer(clAD[4],clAD[4])*np.outer(clBC[0],clBC[0])))*self.TT_TT/nu_l
out[4*n:5*n, 0*n:1*n] = symcl(clAC,"ET")*symcl(clBD,"TT")*self.TT_TT/nu_l + \
symcl(clAD,"ET")*symcl(clBC,"TT")*self.TT_TT/nu_l
#ET_EE
out[4*n:5*n, 1*n:2*n] = np.sqrt(abs(np.outer(clAC[1],clAC[1])*np.outer(clBD[3],clBD[3])))*self.TE_TE/nu_l + \
np.sqrt(abs(np.outer(clAD[1],clAD[1])*np.outer(clBC[3],clBC[3])))*self.TE_TE/nu_l
out[4*n:5*n, 1*n:2*n] = symcl(clAC,"EE")*symcl(clBD,"TE")*self.TE_TE/nu_l + \
symcl(clAD,"EE")*symcl(clBC,"TE")*self.TE_TE/nu_l
#ET_BB
out[4*n:5*n, 2*n:3*n] = 0.
#ET_TE
out[4*n:5*n, 3*n:4*n] = np.sqrt(abs(np.outer(clAC[4],clAC[4])*np.outer(clBD[3],clBD[3])))*self.TT_TT/nu_l + \
np.sqrt(abs(np.outer(clAD[1],clAD[1])*np.outer(clBC[0],clBC[0])))*self.TE_TE/nu_l
out[4*n:5*n, 3*n:4*n] = symcl(clAC,"ET")*symcl(clBD,"TE")*self.TT_TT/nu_l + \
symcl(clAD,"EE")*symcl(clBC,"TT")*self.TE_TE/nu_l
#ET_ET
out[4*n:5*n, 4*n:5*n] = np.sqrt(abs(np.outer(clAC[1],clAC[1])*np.outer(clBD[0],clBD[0])))*self.TE_TE/nu_l + \
np.sqrt(abs(np.outer(clAD[4],clAD[4])*np.outer(clBC[3],clBC[3])))*self.TT_TT/nu_l
out[4*n:5*n, 4*n:5*n] = symcl(clAC,"EE")*symcl(clBD,"TT")*self.TE_TE/nu_l + \
symcl(clAD,"ET")*symcl(clBC,"TE")*self.TT_TT/nu_l
return out
def get_clcov_blocks(self, clAC, clBD, clAD, clBC, nsm=2):
"""
Compute analytical covariance matrix for cross spectra: <AB,CD>
tags=[TT,EE,BB,TE,ET]
inputs:
cl: 2d array (ntag,nbin)
output:
out: 2d array (5*nbin,5*nbin)
"""
#use smooth data as signal model for covariance
sclAC = self._smooth_cl( clAC, nsm=nsm)
sclBD = self._smooth_cl( clBD, nsm=nsm)
sclAD = self._smooth_cl( clAD, nsm=nsm)
......@@ -592,13 +628,13 @@ class Xpol_Wrapper(object):
self._verbose = verbose
mask = np.asarray(mask)
self._MLLFILE = "%s/mll_%d" % (self._tmpdir,self._runnb)
self._MLLFILE = "{:s}/mll_{:d}".format(self._tmpdir,self._runnb)
self._MASKFILE = []
self._MASKFILE.append( "%s/mask_%d_0.fits" % (self._tmpdir,self._runnb))
self._MASKFILE.append( "{:s}/mask_{:d}_0.fits".format(self._tmpdir,self._runnb))
hp.write_map( self._MASKFILE[0], mask)
if mask2 is not None:
self._MASKFILE.append( "%s/mask_%d_1.fits" % (self._tmpdir,self._runnb))
self._MASKFILE.append( "{:s}/mask_{:d}_1.fits".format(self._tmpdir,self._runnb))
hp.write_map( self._MASKFILE[1], mask2)
self._compute_mll()
......@@ -607,14 +643,14 @@ class Xpol_Wrapper(object):
raise ValueError('Mll not found.')
def __del__(self):
self._clean( "mask_%d_0.fits" % self._runnb)
self._clean( "mask_%d_1.fits" % self._runnb)
self._clean( "mll_%d_*.fits" % self._runnb)
self._clean( "pseudo_%d_*.fits" % self._runnb)
self._clean( "cross_%d_*.fits" % self._runnb)
self._clean( "mask_{:d}_0.fits".format(self._runnb))
self._clean( "mask_{:d}_1.fits".format(self._runnb))
self._clean( "mll_{:d}_*.fits".format(self._runnb))
self._clean( "pseudo_{:d}_*.fits".format(self._runnb))
self._clean( "cross_{:d}_*.fits".format(self._runnb))
def _clean(self, tmpfile):
os.system( "rm -f %s/%s" % (self._tmpdir,tmpfile))
os.system( "rm -f {}/{}".format(self._tmpdir,tmpfile))
def _check_exe( self, exe, path=None):
bindir = None
......@@ -634,27 +670,27 @@ class Xpol_Wrapper(object):
def _compute_mll(self,verbose=True):
# lmax = min( [2*self.lmax,3*self.nside-1])
f = open( "%s/xpol_create_mll_%d.par" % (self._tmpdir,self._runnb), "w")
f.write( "nside = %d\n" % self.nside )
f.write( "nstokes = %d\n" % self.nstokes )
f.write( "lmax = %d\n" % self.lmax )
f.write( "nmaps = %d\n" % 2)
f.write( "weightI1 = %s\n" % self._MASKFILE[0])
f.write( "weightP1 = %s\n" % self._MASKFILE[0])
f = open( "{:s}/xpol_create_mll_{:d}.par".format(self._tmpdir,self._runnb), "w")
f.write( "nside = {}\n".format(self.nside ))
f.write( "nstokes = {}\n".format(self.nstokes ))
f.write( "lmax = {}\n".format(self.lmax ))
f.write( "nmaps = {}\n".format( 2))
f.write( "weightI1 = {}\n".format(self._MASKFILE[0]))
f.write( "weightP1 = {}\n".format(self._MASKFILE[0]))
if len(self._MASKFILE) == 2:
f.write( "weightI2 = %s\n" % self._MASKFILE[1])
f.write( "weightP2 = %s\n" % self._MASKFILE[1])
f.write( "mlloutfile = %s\n" % self._MLLFILE )
f.write( "weightI2 = {}\n".format(self._MASKFILE[1]))
f.write( "weightP2 = {}\n".format(self._MASKFILE[1]))
f.write( "mlloutfile = {}\n".format(self._MLLFILE ))
f.close()
str_exe = "%s -n %d %s/xpol_create_mll %s/xpol_create_mll_%d.par" % (self.mpirun,self.nprocs,self.bindir,self._tmpdir,self._runnb)
str_exe = "{} -n {:d} {}/xpol_create_mll {}/xpol_create_mll_{:d}.par".format(self.mpirun,self.nprocs,self.bindir,self._tmpdir,self._runnb)
if not self._verbose:
str_exe += " >& %s/xpol_create_mll_%d.out" % (self._tmpdir,self._runnb)
str_exe += " >& {}/xpol_create_mll_{:d}.out".format(self._tmpdir,self._runnb)
err = os.system( str_exe)
if err:
raiseError( "Error create mll")
self._clean( "xpol_create_mll_%d.par" % self._runnb)
self._clean( "xpol_create_mll_%d.out" % self._runnb)
self._clean( "xpol_create_mll_{:d}.par".format(self._runnb))
self._clean( "xpol_create_mll_{:d}.out".format(self._runnb))
def _compute_spectra( self, map1, map2, bell1, bell2):
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment