Commit 3e62b456 by Matthieu Tristram

 ... ... @@ -5,7 +5,6 @@ Test script for xQML from __future__ import division import numpy as np import healpy as hp from pylab import * ... ... @@ -19,38 +18,36 @@ from xqml.simulation import extrapolpixwin #ion() #show() exp = "Big" patch = "Big" if len(sys.argv) > 1: if sys.argv[1].lower()[0] == "s": exp = "Small" patch = "Small" if exp == "Big": if patch == "Big": nside = 8 dell = 1 glat = 10 elif exp == "Small": fwhm = 1 #deg lmin=2 elif patch == "Small": nside = 64 dell = 10 glat = 80 glat = 70 fwhm = 1 #deg lmin = 2 else: print( "Need a patch !") #lmax = nside lmax = 2 * nside - 1 lmax = 3 * nside - 1 nsimu = 100 MODELFILE = 'planck_base_planck_2015_TTlowP.fits' Slmax = 3*nside-1 # provide list of specs to be computed, and/or options spec = ['EE','BB','EB'] #spec = ['TT','EE','BB','TE']#'EE','BB', 'EB']#, 'TE', 'TB'] pixwin = True ellbins = np.arange(2, lmax + 2, dell) ellbins[-1] = lmax+1 muKarcmin = 1.0 fwhm = 10 muKarcmin = 0.1 ... ... @@ -65,20 +62,19 @@ lth = arange(2, lmax+1) ############################## # Create mask t, p = hp.pix2ang(nside, range(hp.nside2npix(nside))) mask = np.ones(hp.nside2npix(nside), bool) # import random # random.shuffle(mask) if exp == "Big": if patch == "Big": mask[abs(90 - rad2deg(t)) < glat] = False elif exp == "Small": elif patch == "Small": mask[(90 - rad2deg(t)) < glat] = False fsky = np.mean(mask) npix = sum(mask) print("%s patch: fsky=%.2g %% (npix=%d)" % (exp,100*fsky,npix)) print("%s patch: fsky=%.2g %% (npix=%d)" % (patch,100*fsky,npix)) toGB = 1024. * 1024. * 1024. emem = 8.*(npix*2*npix*2) * ( len(lth)*2 ) / toGB print("mem=%.2g Gb" % emem) ... ... @@ -98,34 +94,31 @@ varmap = ones((nstoke * npix)) * pixvar NoiseVar = np.diag(varmap) # ############## Initialise xqml class ############### start = timeit.default_timer() esti = xqml.xQML(mask, ellbins, clth, NA=NoiseVar, NB=NoiseVar, lmax=lmax, fwhm=fwhm, spec=spec) s1 = timeit.default_timer() print( "Init: %d sec" % (s1-start)) ellval = esti.lbin() bins = xqml.Bins.fromdeltal( lmin, lmax, dell) esti = xqml.xQML(mask, bins, clth, NA=NoiseVar, NB=NoiseVar, lmax=lmax, fwhm=fwhm, spec=spec) lb = bins.lbin # ############## Compute Analytical variance ############### #V = esti.get_covariance(cross=True ) #Va = esti.get_covariance(cross=False) s2 = timeit.default_timer() print( "construct covariance: %d sec" % (s2-s1)) #s2 = timeit.default_timer() #print( "construct covariance: %d sec" % (s2-s1)) # ############## Construct MC ############### allcla = [] allcl = [] t = [] bl = hp.gauss_beam(deg2rad(fwhm), lmax=Slmax) fpixw = extrapolpixwin(nside, Slmax, pixwin=pixwin) for n in np.arange(nsimu): bl = hp.gauss_beam(deg2rad(fwhm), lmax=lmax) fpixw = extrapolpixwin(nside, lmax, pixwin=pixwin) for n in range(nsimu): progress_bar(n, nsimu) cmb = np.array(hp.synfast(clth[:, :len(fpixw)]*(fpixw*bl)**2, nside, pixwin=False, lmax=Slmax, fwhm=0.0, new=True, verbose=False)) cmb = np.array(hp.synfast(clth[:, :len(fpixw)]*(fpixw*bl)**2, nside, pixwin=False, lmax=lmax, fwhm=0.0, new=True, verbose=False)) cmbm = cmb[istokes][:, mask] dmA = cmbm + (randn(nstoke * npix) * sqrt(varmap)).reshape(nstoke, npix) dmB = cmbm + (randn(nstoke * npix) * sqrt(varmap)).reshape(nstoke, npix) dmA = cmbm + np.random.randn(nstoke, npix) * sqrt(pixvar) dmB = cmbm + np.random.randn(nstoke, npix) * sqrt(pixvar) s1 = timeit.default_timer() allcl.append(esti.get_spectra(dmA, dmB)) t.append( timeit.default_timer() - s1) ... ... @@ -142,14 +135,13 @@ scla = std(allcla, 0) # ############## Plot results ############### figure(figsize=[12, 8]) clf() Delta = (ellbins[1:] - ellbins[:-1])/2. toDl = lb*(lb+1)/2./np.pi subplot(3, 2, 1) title("Cross") plot(lth, (lth*(lth+1)/2./np.pi)[:, None]*clth[ispecs][:, lth].T, '--k') for s in np.arange(nspec): errorbar(ellval, ellval*(ellval+1)/2./np.pi*hcl[s], yerr=scl[s], xerr=Delta, fmt='o', color='C%i' % ispecs[s], label=r"$%s$" % spec[s]) errorbar(lb, toDl*hcl[s], yerr=toDl*scl[s], xerr=bins.dl/2, fmt='o', color='C%i' % ispecs[s], label=r"$%s$" % spec[s]) semilogy() ylabel(r"$D_\ell$") legend(loc=4, frameon=True) ... ... @@ -158,25 +150,25 @@ subplot(3, 2, 2) title("Auto") plot(lth,(lth*(lth+1)/2./np.pi)[:, None]*clth[ispecs][:, lth].T, '--k') for s in np.arange(nspec): errorbar(ellval, ellval*(ellval+1)/2./np.pi*hcla[s], yerr=scla[s], xerr=Delta, fmt='o', color='C%i' % ispecs[s], label=r"$%s$" % spec[s]) errorbar(lb, toDl*hcla[s], yerr=toDl*scla[s], xerr=bins.dl/2, fmt='o', color='C%i' % ispecs[s], label=r"$%s$" % spec[s]) semilogy() subplot(3, 2, 3) for s in np.arange(nspec): plot(ellval, scl[s], 'o', color='C%i' % ispecs[s], label=r"$\sigma^{%s}_{\rm MC}$" % spec[s]) # plot(ellval, sqrt(diag(V)).reshape(nspec, -1)[s], '-', color='C%i' % ispecs[s]) ylabel(r"$\sigma(C_\ell)$") plot(lb, toDl*scl[s], color='C%i' % ispecs[s], label=r"$\sigma^{%s}_{\rm MC}$" % spec[s]) # plot(lb, sqrt(diag(V)).reshape(nspec, -1)[s], '--', color='C%i' % ispecs[s]) ylabel(r"$\sigma(D_\ell)$") semilogy() subplot(3, 2, 4) for s in np.arange(nspec): plot(ellval, scla[s], 'o', color='C%i' % ispecs[s], label=r"$\sigma^{%s}_{\rm MC}$" % spec[s]) # plot(ellval, sqrt(diag(Va)).reshape(nspec, -1)[s], '-', color='C%i' % ispecs[s]) plot(lb, toDl*scla[s], color='C%i' % ispecs[s], label=r"$\sigma^{%s}_{\rm MC}$" % spec[s]) # plot(lb, sqrt(diag(Va)).reshape(nspec, -1)[s], '--', color='C%i' % ispecs[s]) semilogy() subplot(3, 2, 5) for s in np.arange(nspec): plot(ellval, (hcl[s]-esti.BinSpectra(clth)[s])/(scl[s]/sqrt(nsimu)), '--o', color='C%i' % ispecs[s]) plot(lb, (hcl[s]-bins.bin_spectra(clth)[ispecs[s]])/(scl[s]/sqrt(nsimu)), '--o', color='C%i' % ispecs[s]) ylabel(r"$R[C_\ell]$") xlabel(r"$\ell$") ylim(-3, 3) ... ... @@ -184,7 +176,7 @@ grid() subplot(3, 2, 6) for s in np.arange(nspec): plot(ellval, (hcla[s]-esti.BinSpectra(clth)[s])/(scla[s]/sqrt(nsimu)), '--o', color='C%i' % ispecs[s]) plot(lb, (hcla[s]-bins.bin_spectra(clth)[ispecs[s]])/(scla[s]/sqrt(nsimu)), '--o', color='C%i' % ispecs[s]) xlabel(r"$\ell$") ylim(-3, 3) grid() ... ...
 ... ... @@ -8,12 +8,14 @@ import sys import timeit import string #from scipy import linalg from scipy import linalg import numpy as np import healpy as hp import random as rd from .bins import Bins from .xqml_utils import getstokes, pd_inv, GetBinningMatrix from .estimators import El ... ... @@ -38,8 +40,8 @@ class xQML(object): ---------- mask : 1D array of booleans Mask defining the region of interest (of value True) bins : 1D array of floats (nbin+1) lower multipole for each bin bins : Bins class object Contains informations about bins clth : ndarray of floats Array containing fiducial CMB spectra (unbinned) lmax : int ... ... @@ -47,7 +49,7 @@ class xQML(object): Pl : ndarray or None, optional Normalize Legendre polynomials dS/dCl. Default: None fwhm : float, optional FWHM of the experiment beam FWHM of the experiment beam in degree bell : ndarray, optional beam transfer function (priority over fwhm) pixwin : boolean, optional ... ... @@ -69,12 +71,12 @@ class xQML(object): self.ipok = np.arange(npixtot)[self.mask] self.npix = len(self.ipok) # lower multipole bin self.ellbins = bins # binning (Bins class) self.bins = bins # Maximum multipole based on nside (rule of thumb to avoid aliasing) self.Slmax = np.max(bins)-1 if lmax is None else lmax self.Slmax = bins.lmax if lmax is None else lmax # Beam 2pt function (Gaussian) self.bl = hp.gauss_beam(np.deg2rad(fwhm), lmax=self.Slmax) if bell is not None: ... ... @@ -90,7 +92,7 @@ class xQML(object): if len(clth) == 4: clth = np.concatenate((clth,clth[0:2]*0.)) nbin = len(self.ellbins) nbin = bins.nbins nmem = self.nspec*nbin*(self.nstokes*self.npix)**2 toGb = 1024. * 1024. * 1024. if verbose: ... ... @@ -104,10 +106,16 @@ class xQML(object): # Otherwise compute Pl and S from the arguments. # Ok, but Pl cannot be binned, otherwise S construction is not valid if Pl is None: self.Pl, self.S = compute_ds_dcb( self.ellbins, self.nside, self.ipok, self.bl, clth, self.Slmax, self.spec, pixwin=self.pixwin, verbose=verbose, MC=False, openMP=True) self.Pl, self.S = compute_ds_dcb( self.bins, self.nside, self.ipok, self.bl, clth, self.Slmax, self.spec, pixwin=self.pixwin, verbose=verbose, openMP=True) else: self.Pl = Pl if S is None: ... ... @@ -122,10 +130,10 @@ class xQML(object): if lmax is None: lmax = 2*self.nside-1 #Why ? self.Pl, self.S = compute_ds_dcb( self.ellbins, self.nside, self.ipok, self.bl, clth, lmax, self.Pl, self.S = compute_ds_dcb( self.bins, self.nside, self.ipok, self.bl, clth, lmax, self.spec, pixwin=self.pixwin, verbose=verbose, MC=MC, openMP=openMP) return( self.Pl, self.S) def construct_esti(self, NA, NB=None, verbose=False, thread=False): """ Compute the inverse of the datasets pixel covariance matrices C, ... ... @@ -151,18 +159,18 @@ class xQML(object): # Invert (signalB + noise) matrix invCb = pd_inv(self.S + self.NB) # Compute El = Ca^-1.Pl.Cb^-1 (long) self.El = El(invCa, invCb, self.Pl, thread=thread, verbose=verbose) # Finally compute invW by inverting (longer) self.invW = np.linalg.inv(CrossWindowFunction(self.El, self.Pl, thread=thread, verbose=verbose)) self.invW = linalg.inv(CrossWindowFunction(self.El, self.Pl, thread=thread, verbose=verbose)) # Compute bias for auto # if not self.cross: # self.bias = biasQuadEstimator(self.NA, self.El) self.bias = biasQuadEstimator(self.NA, self.El) if verbose: print( "Construct estimator: %.1f sec" % (timeit.default_timer()-tstart)) ... ... @@ -248,80 +256,3 @@ class xQML(object): P, Q, ell, ellval = GetBinningMatrix(self.ellbins, self.Slmax) return( ellval) ##Not USED YET class Bins(object): """ lmins : list of integers Lower bound of the bins lmaxs : list of integers Upper bound of the bins (not included) """ def __init__( self, lmins, lmaxs): if not(len(lmins) == len(lmaxs)): raise ValueError('Incoherent inputs') cutfirst = np.where( lmaxs>2)[0] self.lmins = lmins[cutfirst] self.lmaxs = lmaxs[cutfirst] if self.lmins[0] < 2: self.lmins[0] = 2 self._derive_ext() @classmethod def fromdeltal( cls, lmin, lmax, delta_ell): nbins = (lmax - lmin + 1) // delta_ell lmins = lmin + np.arange(nbins) * delta_ell lmaxs = lmins + delta_ell return( cls( lmins, lmaxs)) def _derive_ext( self): self.lmin = min(self.lmins) self.lmax = max(self.lmaxs)-1 if self.lmin < 1: raise ValueError('Input lmin is less than 1.') if self.lmax < self.lmin: raise ValueError('Input lmax is less than lmin.') self.nbins = len(self.lmins) self.lbin = (self.lmins + self.lmaxs - 1) / 2 self.dl = (self.lmaxs - self.lmins) def bins(self): return (self.lmins,self.lmaxs) def cut_binning(self, lmin, lmax): sel = np.where( (self.lmins > lmin) & (self.lmaxs <= lmax+1) )[0] self.lmins = self.lmins[sel] self.lmaxs = self.lmaxs[sel] self._derive_ext() def _bin_operators(self): ell2 = np.arange(self.lmax+1) ell2 = ell2 * (ell2 + 1) / (2 * np.pi) p = np.zeros((self.nbins, self.lmax+1)) q = np.zeros((self.lmax+1, self.nbins)) for b, (a, z) in enumerate(zip(self.lmins, self.lmaxs)): p[b, a:z] = ell2[a:z] / (z - a) q[a:z, b] = 1 / ell2[a:z] return p, q def bin_spectra(self, spectra): """ Average spectra in bins specified by lmin, lmax and delta_ell, weighted by l(l+1)/2pi. """ spectra = np.asarray(spectra) minlmax = min([spectra.shape[-1] - 1,self.lmax]) fact_binned = 2 * np.pi / (self.lbin * (self.lbin + 1)) _p, _q = self._bin_operators() return np.dot(spectra[..., :minlmax+1], _p.T[:minlmax+1,...]) * fact_binned