Commit e04935bc authored by Matthieu Tristram's avatar Matthieu Tristram

Update symmetrization for cov

Merge branch 'master' of gitlab.in2p3.fr:tristram/Xpol

Conflicts:
	xpol/xpol.py
parents 2f04d329 29a33ab0
Pipeline #99466 passed with stage
in 2 minutes and 13 seconds
......@@ -4,7 +4,7 @@ import healpy as hp
from numpy import *
import matplotlib.pylab as plt
nside=128
nside=256
lmax=3*nside-1
tagnames = ['TT','EE','BB','TE','TB','EB']
......@@ -30,6 +30,7 @@ xp = xpol.Xpol( mask, binning, verbose=True)
#compute spectra from map dT
pcl,cl = xp.get_spectra( dT)
#plot result
plt.figure()
plt.subplot( 2,1,1)
......@@ -43,39 +44,42 @@ plt.axhline( 1, color='k')
#small MonteCarlo
nsimu = 100
xp = xpol.Xpol( mask, binning)
allcl = []
for i in range(100):
for i in range(nsimu):
print(i)
dT = hp.synfast( clth, nside, new=True, lmax=3*nside-1, pixwin=True, verbose=False)
pcl,cl = xp.get_spectra( dT, pixwin=True)
allcl.append(cl)
#figure
cb = binning.bin_spectra(clth)
cb = list(cb)+list([cb[0]*0.,cb[0]*0.])
## cosmic = sqrt( 2./(2*l+1))/mean(mask) * cb
## cosmic[3] = sqrt( 1./(2*l+1)*(cb[0]*cb[1]+cb[3]*cb[3]))/xpol.fsky(mask)
cosmic = xpol.sample_variance( lb, cb, fsky=xpol.fsky(mask))
tag=1
plt.figure()
for tag in range(6):
plt.subplot( 3,6,0*6+tag+1)
plt.plot( lb, lb*(lb+1)/2/pi*cb[tag], 'k')
# plt.semilogx()
plt.plot( lb, lb*(lb+1)/2/pi*mean(allcl,0)[tag], 'r')
plt.plot( lb, lb*(lb+1)/2/pi*(mean(allcl,0)-std(allcl,0))[tag], 'r--')
plt.plot( lb, lb*(lb+1)/2/pi*(mean(allcl,0)+std(allcl,0))[tag], 'r--')
plt.plot( lb, lb*(lb+1)/2/pi*cb[tag]*1e12, 'k')
plt.plot( lb, lb*(lb+1)/2/pi*mean(allcl,0)[tag]*1e12, 'r')
plt.plot( lb, lb*(lb+1)/2/pi*(mean(allcl,0)-std(allcl,0))[tag]*1e12, 'r--')
plt.plot( lb, lb*(lb+1)/2/pi*(mean(allcl,0)+std(allcl,0))[tag]*1e12, 'r--')
plt.title( tagnames[tag])
plt.ylabel( "$\ell(\ell+1)/2\pi\,C_\ell$")
if tag == 0: plt.ylabel( "$\ell(\ell+1)/2\pi\,C_\ell\,[\mu K^2]$")
plt.subplot( 3,6,1*6+tag+1)
plt.plot( lb, mean(allcl,0)[tag]/cb[tag], 'r')
plt.axhline( 1, color='k')
plt.plot( lb, (mean(allcl,0)[tag]-cb[tag])/(std(allcl,0)[tag]/sqrt(nsimu)), 'r')
plt.ylim( -4,4)
plt.axhline( 0, color='k')
if tag == 0: plt.ylabel( "$R_\ell [\sigma]$")
plt.subplot( 3,6,2*6+tag+1)
plt.plot( lb, lb*(lb+1)/2/pi*std(allcl,0)[tag], 'r')
plt.plot( lb, lb*(lb+1)/2/pi*cosmic[tag], 'k')
plt.ylabel( "$\ell(\ell+1)/2\pi\,\sigma_\ell$")
plt.plot( lb, lb*(lb+1)/2/pi*std(allcl,0)[tag]*1e12, 'r')
plt.plot( lb, lb*(lb+1)/2/pi*cosmic[tag]*1e12, 'k')
if tag == 0: plt.ylabel( "$\ell(\ell+1)/2\pi\,\sigma_\ell\,[\mu K^2]$")
......
......@@ -4,12 +4,12 @@ HEALPIX=/sps/planck/Users/tristram/Soft/Xpol/software/libs/Healpix_3.70
CFITSIO=/sps/planck/Users/tristram/Soft/Xpol/software/libs/
S2HAT=/sps/planck/Users/tristram/Soft/Xpol/software/libs/s2hat
export INC_DIR = -I${S2HAT_DIR} -I${HEALPIX_DIR}/include -I${CFITSIO_DIR}/include
export INC_DIR = -I${MKLROOT}/include -I${S2HAT_DIR} -I${HEALPIX_DIR}/include -I${CFITSIO_DIR}/include
export LIB_DIR = \
-L${S2HAT_DIR} -ls2hat_hlpx_g \
-Wl,-R${HEALPIX_DIR}/lib -L${HEALPIX_DIR}/lib -lchealpix -lhealpix -lsharp -lm \
-Wl,-R${CFITSIO_DIR} -L${CFITSIO_DIR} -lcfitsio -lm \
-Wl,-R${MKLROOT} -L${MKLROOT} -lmkl_blacs_openmpi_lp64 -lmkl_gf_lp64 -lmkl_core -lmkl_gnu_thread -lmkl_def -liomp5 -ldl -lpthread -lm
-Wl,-R${MKLROOT}/lib/intel64 -L${MKLROOT}/lib/intel64 -lmkl_blacs_openmpi_lp64 -lmkl_gf_lp64 -lmkl_core -lmkl_gnu_thread -lmkl_def -ldl -lpthread -lm
export LIB =
......@@ -19,7 +19,8 @@ export F77 = mpif77
export CFLAG = -O3 -DMKL -DHEALPIXDATA=\"${HEALPIX_DIR}/data\" -Wall -Dx86_64 -DOPEN_MPI
export FFLAG = -O3 -fno-second-underscore
export LFLAG = -O3 -lmpi -lmpi_mpifh -lgfortran -lgomp -fopenmp #-lifcore -lintlc -lifport -limf -lirc -ldl -lsvml
export LFLAG = -O3 -lmpi -lmpi_mpifh -lgfortran -lgomp -fopenmp
#-liomp5 -lifcore -lintlc -lifport -limf -lirc -ldl -lsvml
################################################################################################
......
......@@ -811,13 +811,14 @@ void compute_pcl_variance( int nbin, int *bintab,
double BDl = (specBD[b1]+specBD[b2])/2.;
double ADl = (specAD[b1]+specAD[b2])/2.;
double BCl = (specBC[b1]+specBC[b2])/2.;
local_corr[b1*nbin+b2] = ( mll3[b1*nbin+b2]*ACl*BDl + mll4[b1*nbin+b2]*ADl*BCl ) / nu_l;
local_corr[b1*nbin+b2] = ( mll3[b1*nbin+b2] * sqrt( ACl*BDl) + mll4[b1*nbin+b2] * sqrt( ADl*BCl) ) / nu_l;
//old symmetrization
// local_corr[b1*nbin+b2] = (
// mll3[b1*nbin+b2] * sqrt( fabs(specAC[b1]*specAC[b2] * specBD[b1]*specBD[b2])) +
// mll4[b1*nbin+b2] * sqrt( fabs(specAD[b1]*specAD[b2] * specBC[b1]*specBC[b2]))
// ) / nu_l;
}
}
MPI_Reduce( local_corr, corr, nbin*nbin, MPI_DOUBLE, MPI_SUM, root, comm);
......
......@@ -24,7 +24,7 @@ def sample_variance( l, cl, fsky=1.):
if len(cl) > 4:
cosmic[4] = np.sqrt( 1./(2*l+1)*(cl[0]*cl[2]+cl[4]*cl[4]))/fsky
if len(cl) > 5:
cosmic[5] = np.sqrt( 1./(2*l+1)*(cl[1]*cl[2]+cl[5]*cl[5]))/fsky
cosmic[5] = np.sqrt( 1./(2*l+1)*(cl[1]*cl[2]+cl[5]*cl[5]))/fsky
return cosmic
......@@ -122,7 +122,7 @@ class Xpol(object):
biased, unbiased = xpol.get_spectra(map1, map2)
"""
def __init__(self, mask, bins, mask2=None, polar=True, verbose=False):
def __init__(self, mask, bins, mask2=None, polar=True, verbose=False, wlmax=None):
"""
Parameters
----------
......@@ -132,17 +132,13 @@ class Xpol(object):
"""
mask = np.asarray(mask)
lmin = int(bins.lmin)
lmax = int(bins.lmax)
nside = hp.npix2nside(len(mask))
if lmax > 3*nside-1:
raise ValueError('Input lmax is too high for resolution.')
self.lmin = int(bins.lmin)
self.lmax = int(bins.lmax)
self.nside = hp.npix2nside(len(mask))
self.verbose = verbose
self.nside = nside
self.polar = polar
self.lmin = lmin
self.lmax = lmax
if self.lmax > 3*self.nside-1:
raise ValueError('Input lmax is too high for resolution.')
if self.verbose: print( "Compute bin operators")
self._p, self._q = bins._bin_operators()
......@@ -151,7 +147,8 @@ class Xpol(object):
#compute cross-spectrum of weight mask
if self.verbose: print( "alm2map mask")
wlmax = 3*nside-1
if wlmax is None:
wlmax = self.lmax
if mask2 is None:
wl = hp.anafast(mask,lmax=wlmax)
self.mask = self.mask2 = mask
......@@ -166,7 +163,6 @@ class Xpol(object):
if self.verbose: print( "Inverse Mll")
self.mll_binned_inv = self._inv_mll(mll_binned)
# self.mll_binned_inv = np.linalg.inv(mll_binned)
def get_spectra(self, map1, map2=None, bell=None, bell1=None, bell2=None, pixwin=True):
"""
......@@ -199,7 +195,6 @@ class Xpol(object):
unbiased : float array of shape (9, nbins)
The Xpol's (cross-) power spectra for TT,EE,BB,TE,TB,EB,ET,BT,BE.
The corresponding l values are given by `xpol.lbin`.
"""
#Map1
......@@ -213,8 +208,9 @@ class Xpol(object):
map1 = map1.T
self._removeUndef(map1)
map1[0] = hp.remove_dipole( map1[0], verbose=False)
alms1 = hp.map2alm(map1 * self.mask, pol=self.polar, lmax=self.lmax)
#Map2
if map2 is None:
alms2 = alms1
......@@ -224,6 +220,7 @@ class Xpol(object):
if map2.shape[-1] == 3:
map2 = map2.T
self._removeUndef(map2)
map2[0] = hp.remove_dipole( map2[0], verbose=False)
alms2 = hp.map2alm( map2 * self.mask2, pol=self.polar, lmax=self.lmax)
#alm2cl
......@@ -294,9 +291,11 @@ class Xpol(object):
return np.asarray( [TT,EE,BB,TE,TB,EB,ET,BT,BE])
def _replaceUndefWith0(self,mymap):
for i in range(len(mymap)):
if( float(mymap[i]) == hp.UNSEEN):
mymap[i] = 0.
badpix = np.isclose( mymap,hp.UNSEEN)
mymap[badpix] = 0.
# for i in range(len(mymap)):
# if np.isclose( mymap[i],hp.UNSEEN):
# mymap[i] = 0.
def _removeUndef(self,mymap):
if self.polar:
......
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