Merge branch 'doctest' into 'master'

Apply PEP8 recommandation + add initial doctest in simulation.py

See merge request !2

TestEsti.py

-from pylab import *
-import astropy.io.fits as fits
-import xqml
-import timeit
-import numpy as np
-import healpy as hp
-import matplotlib.pyplot as plt
-
-
-nside = 8
-lmax = 3*nside-1
-Slmax = 3*nside-1
-deltal = 1
-clth = hp.read_cl( 'planck_base_planck_2015_TTlowP.fits')
-lth = arange(2,lmax+1)
-
-ellbins = arange(2,lmax+2,deltal)
-ellbins[-1] = lmax+1
-
-P, Q, ell, ellval = xqml.GetBinningMatrix(ellbins, lmax)
-nbins=len(ellbins)-1
-
-
-#create mask
-t,p=hp.pix2ang( nside, range(hp.nside2npix(nside)))
-mask = np.ones(hp.nside2npix(nside),bool)
-mask[abs(90-rad2deg(t)) < 10] = False
-npix = sum(mask)
-
-fwhm = 0.5
-bl=hp.gauss_beam(deg2rad(fwhm), lmax=Slmax+1)
-
-allStoke, der, ind = xqml.getstokes(polar=True,temp=False,EBTB=False)
-nder=len(der)
-
-muKarcmin = 0.1
-
-
-pixvar = 2*xqml.muKarcmin2var(muKarcmin, nside)
-varmap = ones((2*npix))*pixvar
-NoiseVar = np.diag(varmap)
-
-cmb = hp.synfast( clth, nside, fwhm=deg2rad(fwhm), pixwin=True, new=True, verbose=False)
-noise = (randn(len(varmap)) * varmap**.5).reshape(2,-1)
-dm = cmb[1:,mask] + noise
-
-
-###################################### Compute ds_dcb ######################################
-ip=arange(hp.nside2npix(nside))
-ipok=ip[mask]
-
-Pl, S = xqml.compute_ds_dcb(ellbins,nside,ipok,bl,clth,Slmax,polar=True,temp=False, EBTB=False, pixwining=True, timing=True, MC=False)
-Pl = Pl.reshape((nder)*(np.shape(Pl)[1]), 2*npix, 2*npix)
-
-
-
-
-
-
-###################################### Compute spectra ######################################
-esti = xqml.xQML(mask,ellbins, clth, Pl=Pl, fwhm=fwhm)
-esti.construct_esti( NoiseVar, NoiseVar)
-cl = esti.get_spectra( dm, dm)
-V = esti.get_covariance()
-
-
-###################################### Construct MC ######################################
-allcl = []
-esti = xqml.xQML(mask,ellbins, clth, Pl=Pl, fwhm=fwhm)
-esti.construct_esti( NoiseVar, NoiseVar)
-start = timeit.default_timer()
-for n in np.arange(100):
-	xqml.progress_bar( n, 100,timeit.default_timer()-start)
-	cmb = hp.synfast( clth, nside, fwhm=deg2rad(fwhm), pixwin=True, new=True, verbose=False)
-	dm = cmb[1:,mask] + (randn(2*npix)*sqrt(varmap)).reshape( 2, npix)
-	allcl.append(esti.get_spectra( dm, dm))
-
-
-figure()
-subplot( 2,1,1)
-plot( lth, clth.transpose()[lth,1:3], 'k')
-plot( ellval, mean( allcl,0).transpose(), 'r')
-plot( ellval, mean( allcl,0).transpose() + std( allcl,0).transpose(), 'r--')
-plot( ellval, mean( allcl,0).transpose() - std( allcl,0).transpose(), 'r--')
-semilogy()
-subplot( 2,1,2)
-cosmic = sqrt( 2./(2*lth+1))/mean(mask) * clth[1:3,lth]
-plot( lth, cosmic.transpose(), 'k')
-plot( ellval, std( allcl,0).transpose(), 'r')
-plot( ellval, sqrt(diag(V)).reshape(nder,-1).transpose(), 'b')
-semilogy()
+"""
+Test script for xQML
 
+Author: Vanneste
+"""
+from __future__ import division
 
+import timeit
 
+import numpy as np
+import healpy as hp
+from pylab import *
+import astropy.io.fits as fits
 
-show()
+from xqml import xQML
+from libcov import compute_ds_dcb
+from xqml_utils import progress_bar
+from simulation import getstokes, muKarcmin2var, GetBinningMatrix
+
+if __name__ == "__main__":
+    nside = 8
+    lmax = 3 * nside - 1
+    Slmax = 3 * nside - 1
+    deltal = 1
+    clth = hp.read_cl('planck_base_planck_2015_TTlowP.fits')
+    lth = arange(2, lmax+1)
+
+    ellbins = arange(2, lmax + 2, deltal)
+    ellbins[-1] = lmax + 1
+
+    P, Q, ell, ellval = GetBinningMatrix(ellbins, lmax)
+    nbins = len(ellbins) - 1
+
+    # Create mask
+    t, p = hp.pix2ang(nside, range(hp.nside2npix(nside)))
+    mask = np.ones(hp.nside2npix(nside), bool)
+    mask[abs(90 - rad2deg(t)) < 10] = False
+    npix = sum(mask)
+
+    fwhm = 0.5
+    bl = hp.gauss_beam(deg2rad(fwhm), lmax=Slmax + 1)
+
+    allStoke, der, ind = getstokes(polar=True, temp=False, EBTB=False)
+    nder = len(der)
+
+    muKarcmin = 0.1
+
+    pixvar = 2 * muKarcmin2var(muKarcmin, nside)
+    varmap = ones((2 * npix)) * pixvar
+    NoiseVar = np.diag(varmap)
+
+    cmb = hp.synfast(
+        clth, nside, fwhm=deg2rad(fwhm), pixwin=True, new=True, verbose=False)
+    noise = (randn(len(varmap)) * varmap**0.5).reshape(2, -1)
+    dm = cmb[1:, mask] + noise
+
+    ############### Compute ds_dcb ###############
+    ip = arange(hp.nside2npix(nside))
+    ipok = ip[mask]
+
+    Pl, S = compute_ds_dcb(
+        ellbins, nside, ipok, bl, clth, Slmax,
+        polar=True, temp=False, EBTB=False,
+        pixwining=True, timing=True, MC=False)
+    Pl = Pl.reshape((nder)*(np.shape(Pl)[1]), 2 * npix, 2 * npix)
+
+    ############### Compute spectra ###############
+    esti = xQML(mask, ellbins, clth, Pl=Pl, fwhm=fwhm)
+    esti.construct_esti(NoiseVar, NoiseVar)
+    cl = esti.get_spectra(dm, dm)
+    V = esti.get_covariance()
+
+    ############### Construct MC ###############
+    allcl = []
+    esti = xQML(mask, ellbins, clth, Pl=Pl, fwhm=fwhm)
+    esti.construct_esti(NoiseVar, NoiseVar)
+    start = timeit.default_timer()
+    for n in np.arange(100):
+        progress_bar(n, 100, timeit.default_timer() - start)
+        cmb = hp.synfast(
+            clth, nside, fwhm=deg2rad(fwhm),
+            pixwin=True, new=True, verbose=False)
+        dm = cmb[1:, mask] + (randn(2 * npix) * sqrt(varmap)).reshape(2, npix)
+        allcl.append(esti.get_spectra(dm, dm))
+
+    figure()
+    subplot(2, 1, 1)
+    plot(lth, clth.transpose()[lth, 1: 3], 'k')
+    plot(ellval, mean(allcl, 0).transpose(), 'r')
+    plot(ellval, mean(allcl, 0).transpose() + std(allcl, 0).transpose(), 'r--')
+    plot(ellval, mean(allcl, 0).transpose() - std(allcl, 0).transpose(), 'r--')
+    semilogy()
+    subplot(2, 1, 2)
+    cosmic = sqrt(2./(2 * lth + 1)) / mean(mask) * clth[1: 3, lth]
+    plot(lth, cosmic.transpose(), 'k')
+    plot(ellval, std(allcl, 0).transpose(), 'r')
+    plot(ellval, sqrt(diag(V)).reshape(nder, -1).transpose(), 'b')
+    semilogy()
+    show()

estimators.py

+"""
+Set of routines to ...
+
+Comment: what is Expd()? It is defined nowhere in the code...
+
+Author: Vanneste
+"""
+from __future__ import division
+
+import numpy as np
+
+def Pl(ds_dcb):
+    """
+    Reshape ds_dcbin Pl
+
+    Parameters
+    ----------
+    ds_dcb : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    nnpix = np.shape(ds_dcb)[-1]
+    return np.copy(ds_dcb).reshape(2 * (np.shape(ds_dcb)[1]), nnpix, nnpix)
+
+def CorrelationMatrix(Clth, Pl, ellbins, polar=True, temp=False, EBTB=False):
+    """
+    Compute correlation matrix S = sum_l Pl*Cl
+
+    Parameters
+    ----------
+    Clth : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    if EBTB:
+        xx = ['TT', 'EE', 'BB', 'TE', 'EB', 'TB']
+        ind = [0, 1, 2, 3, 4, 5]
+    else:
+        xx = ['TT', 'EE', 'BB', 'TE']
+        ind = [0, 1, 2, 3]
+
+    if not temp:
+        allStoke = ['Q', 'U']
+        if EBTB:
+            xx = ['EE', 'BB', 'EB']
+            ind = [1, 2, 5]
+        else:
+            xx = ['EE', 'BB']
+            ind = [1, 2]
+    if not polar:
+        allStoke = ['I']
+        xx = ['TT']
+        ind = [0]
+
+    clth = Clth[ind][:, 2: int(ellbins[-1])].flatten()
+    return np.sum(Pl * clth[:, None, None], 0)
+
+def El(invCAA, invCBB, Pl, Bll=None, expend=False):
+    """
+    Compute El = invCAA.Pl.invCBB
+
+    Parameters
+    ----------
+    invCAA : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    # Tegmark B-matrix useless so far)
+    if Bll is None:
+        Bll = np.diagflat((np.ones(len(Pl))))
+
+    lmax = len(Pl) * 2**int(expend)
+    lrange = np.arange(lmax)
+    npix = len(invCAA)
+    # triangular shape ds_dcb
+    if expend:
+        El = np.array(
+            [np.dot(np.dot(invCAA, Expd(Pl, l)), invCBB) for l in lrange]
+            ).reshape((lmax, npix, npix))
+    else:
+        El = np.array(
+            [np.dot(np.dot(invCAA, Pl[l]), invCBB) for l in lrange]
+            ).reshape((lmax, npix, npix))
+    return El
+
+def ElLong(invCAA, invCBB, Pl, Bll=None, expend=False):
+    """
+    Compute El = invCAA.Pl.invCBB
+
+    Parameters
+    ----------
+    invCAA : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    # Tegmark B-matrix useless so far)
+    if Bll is None:
+        Bll = np.diagflat((np.ones(len(Pl))))
+
+    lmax = len(Pl) * 2**int(expend)
+    lrange = np.arange(lmax)
+    npix = len(invCAA)
+    # triangular shape ds_dcb
+    if expend:
+        for l in lrange:
+            Pl[l] = np.dot(
+                np.dot(invCAA, Expd(Pl, l)), invCBB).reshape((npix, npix))
+    else:
+        for l in lrange:
+            Pl[l] = np.dot(np.dot(invCAA, Pl[l]), invCBB).reshape((npix, npix))
+
+def CrossFisherMatrix(El, Pl, expend=False):
+    """
+    Compute cross (or auto) fisher matrix
+    Fll = Trace[invCAA.Pl.invCBB.Pl] = Trace[El.Pl]
+
+    Parameters
+    ----------
+    El : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    lmax = len(Pl) * 2**int(expend)
+    lrange = np.arange(lmax)
+    if expend:
+        # pas de transpose car symm
+        FAB = np.array(
+            [np.sum(El[il] * Expd(Pl, jl)) for il in lrange for jl in lrange]
+            ).reshape(lmax, lmax)
+    else:
+        # pas de transpose car symm
+        FAB = np.array(
+            [np.sum(El[il] * Pl[jl]) for il in lrange for jl in lrange]
+            ).reshape(lmax, lmax)
+    return FAB
+
+def CrossWindowFunction(El, Pl):
+    """
+    Compute Tegmark cross (or auto) window matrix
+    Wll = Trace[invCAA.Pl.invCBB.Pl] = Trace[El.Pl]
+    Equivalent to Fisher matrix Fll when Tegmark B-matrix = 1
+
+    Parameters
+    ----------
+    El : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    lmax = len(El)
+    lrange = np.arange((lmax))
+    # pas de transpose car symm
+    Wll = np.array(
+        [np.sum(El[il] * Pl[jl]) for il in lrange for jl in lrange]
+        ).reshape(lmax, lmax)
+    return Wll
+
+def CrossWindowFunctionLong(invCAA, invCBB, Pl):
+    """
+    Compute Tegmark cross (or auto) window matrix
+    Wll = Trace[invCAA.Pl.invCBB.Pl] = Trace[El.Pl]
+    Equivalent to Fisher matrix Fll when Tegmark B-matrix = 1
+
+    Parameters
+    ----------
+    invCAA : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    lmax = len(Pl)
+    lrange = np.arange((lmax))
+    # Pas de transpose car symm
+    Wll = np.array(
+        [np.sum(np.dot(np.dot(invCAA, Pl[il]), invCBB) * Pl[jl])
+         for il in lrange for jl in lrange]).reshape(lmax, lmax)
+    return Wll
+
+def CrossMisherMatrix(El, CAA, CBB):
+    """
+    Compute matrix MAB = Trace[El.CAA.El.CBB] (see paper for definition)
+
+    Parameters
+    ----------
+    El : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    lmax = len(El)
+    lrange = np.arange(lmax)
+    El_CAA = np.array([np.dot(CAA, El[il]) for il in lrange])
+    El_CBB = np.array([np.dot(CBB, El[il]) for il in lrange])
+    MAB = np.array(
+        [np.sum(El_CAA[il] * El_CBB[jl].T) for il in lrange for jl in lrange]
+        ).reshape(lmax, lmax)
+    return MAB
+
+def CrossGisherMatrix(El, CAB):
+    """
+    Compute matrix GAB = Trace[El.CAB.El.CAB] (see paper for definition)
+
+    Parameters
+    ----------
+    El : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    lmax = len(El)
+    lrange = np.arange(lmax)
+    El_CAB = np.array([np.dot(CAB, El[il]) for il in lrange])
+    GAB = np.array(
+        [np.sum(El_CAB[il] * El_CAB[jl].T) for il in lrange for jl in lrange]
+        ).reshape(lmax, lmax)
+    return GAB
+
+def CrossGisherMatrixLong(El, CAB):
+    """
+    Compute matrix GAB = Trace[El.CAB.El.CAB] (see paper for definition)
+
+    Parameters
+    ----------
+    El : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    lmax = len(El)
+    lrange = np.arange(lmax)
+    GAB = np.array(
+        [np.sum(np.dot(CAB, El[il]) * np.dot(CAB, El[jl]).T)
+         for il in lrange for jl in lrange]).reshape(lmax, lmax)
+    return GAB
+
+def yQuadEstimator(dA, dB, El):
+    """
+    Compute pre-estimator 'y' such that Cl = Fll^-1 . yl
+
+    Parameters
+    ----------
+    dA : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    npix = len(dA)
+    lrange = np.arange((len(El)))
+    y = np.array([dA.dot(El[l]).dot(dB) for l in lrange])
+    return y
+
+def ClQuadEstimator(invW, y):
+    """
+    Compute estimator 'Cl' such that Cl = Fll^-1 . yl
+
+    Parameters
+    ----------
+    invW : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    Cl = np.dot(invW, y)
+    return Cl
+
+def biasQuadEstimator(NoiseN, El):
+    """
+    Compute bias term bl such that Cl = Fll^-1 . ( yl + bias)
+
+    Parameters
+    ----------
+    NoiseN : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    lrange = np.arange((len(El)))
+    return np.array([np.sum(NoiseN * El[l]) for l in lrange])
+
+def blEstimatorFlat(NoiseN, El):
+    """
+    Compute bias term bl such that Cl = Fll^-1 . ( yl + bl)
+    Not to be confonded with beam function bl(fwhm)
+
+    Parameters
+    ----------
+    NoiseN : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    lrange = np.arange((len(El)))
+
+    return np.array([np.sum(NoiseN * np.diag(El[l])) for l in lrange])
+
+
+def CovAB(invWll, GAB):
+    """
+    Compute analytical covariance matrix Cov(Cl, Cl_prime)
+
+    Parameters
+    ----------
+    invWll : ???
+        ???
+
+    Returns
+    ----------
+    ???
+    """
+    covAB = np.dot(np.dot(invWll, GAB), invWll.T) + invWll
+    return covAB

libangles.py

+"""
+Set of routines to ...
+
+Author: Vanneste
+"""
+from __future__ import division
+
+import numpy as np
+
+def polrotangle(ri, rj):
+    """
+    ???
+
+    Parameters
+    ----------
+    ri : ???
+        ???
+    rj : ???
+        ???
+
+    Returns
+    ----------
+    cos2a : 1D array of floats
+        ???
+    sin2a : 1D array of floats
+        ???
+    """
+    z = np.array([0.0, 0.0, 1.0])
+
+    # Compute ri^rj : unit vector for the great circle connecting i and j
+    rij = np.cross(ri, rj)
+    norm = np.sqrt(np.dot(rij, np.transpose(rij)))
+
+    # case where pixels are identical or diametrically opposed on the sky
+    if norm <= 1e-15:
+        cos2a = 1.0
+        sin2a = 0.0
+        return cos2a, sin2a
+    rij = rij / norm
+
+    # Compute z^ri : unit vector for the meridian passing through pixel i
+    ris = np.cross(z, ri)
+    norm = np.sqrt(np.dot(ris, np.transpose(ris)))
+
+    # case where pixels is at the pole
+    if norm <= 1e-15:
+        cos2a = 1.0
+        sin2a = 0.0
+        return cos2a, sin2a
+    ris = ris / norm
+
+    # Now, the angle we want is that
+    # between these two great circles: defined by
+    cosa = np.dot(rij, np.transpose(ris))
+
+    # the sign is more subtle : see tegmark et de oliveira costa 2000 eq. A6
+    rijris = np.cross(rij, ris)
+    sina = np.dot(rijris, np.transpose(ri))
+
+    # so now we have directly cos2a and sin2a
+    cos2a = 2.0 * cosa * cosa - 1.0
+    sin2a = 2.0 * cosa * sina
+
+    return cos2a, sin2a

simulation.py

+"""
+Set of routines to ...
+
+Author: Vanneste
+"""
+from __future__ import division
+
+import numpy as np
+import healpy as hp
+
+from scipy import sparse
+
+def muKarcmin2var(muKarcmin, nside):
+    """
+    Return pixel variance for a given nside and noise level [1e-6 K . arcmin]
+