diff --git a/examples/script_xpol.py b/examples/script_xpol.py index ada1327ee75f028d3bd8f300b876cffc5db5b35c..4e8e5c6027952ee08ff068f58832edf1e18b21bc 100644 --- a/examples/script_xpol.py +++ b/examples/script_xpol.py @@ -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]$") + + diff --git a/makefile_gnu b/makefile_gnu index 64c722f6e4111565a5edb44fb5ed0de2e1a1aa8d..7e5bc291a81b966d65b0ad78a1ab2c4df578f4f5 100644 --- a/makefile_gnu +++ b/makefile_gnu @@ -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 ################################################################################################ diff --git a/src/xpol_fisher_tools.c b/src/xpol_fisher_tools.c index 8c659a368ad5a338abcda83088bccc6d24b22fbd..2a000d44301086959307b5d6ccb62f5c01786b50 100644 --- a/src/xpol_fisher_tools.c +++ b/src/xpol_fisher_tools.c @@ -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); diff --git a/xpol/xpol.py b/xpol/xpol.py index 32c0b0a32174eae0eaaca6424086270bb2263cf5..0aca9d80c866be4d8b9f225b6f6f3194d409bdc9 100644 --- a/xpol/xpol.py +++ b/xpol/xpol.py @@ -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: