Commit b3c88e99 authored by Matthieu Tristram's avatar Matthieu Tristram
Browse files

Try to reduce memory

Use symarray class
Clean TestEsti.py
parent 9e452df4
#!/usr/bin/env python
"""
Test script for xQML
......@@ -6,131 +7,146 @@ 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
import timeit
import sys
import xqml
import libcov
import simulation
from xqml_utils import progress_bar
from simulation import muKarcmin2var, GetBinningMatrix
from simulation import extrapolpixwin
from xqml.xqml_utils import progress_bar, getstokes
from xqml.simulation import muKarcmin2var
from xqml.simulation import extrapolpixwin
ion()
show()
# if __name__ == "__main__":
# # Inputs
nside = 64
exp = "Big"
if len(sys.argv) > 1:
if sys.argv[1].lower()[0] == "s":
exp = "Small"
if exp == "Big":
nside = 8
dell = 1
glat = 10
elif exp == "Small":
nside = 64
dell = 10
glat = 80
else:
print( "Need a patch !")
#lmax = nside
lmax = 2 * nside - 1
Slmax = 2 * nside - 1
dell = 10
nsimu = 1000
clth = np.array(hp.read_cl('planck_base_planck_2015_TTlowP.fits'))
Clthshape = zeros(((6,)+shape(clth)[1:]))
Clthshape[:4] = clth
clth = Clthshape
nsimu = 100
MODELFILE = 'planck_base_planck_2015_TTlowP.fits'
Slmax = lmax
# # spectra correlations level
EB = 0.0
clth[4] = EB*sqrt(clth[1]*clth[2])
TB = 0.0
clth[5] = TB*sqrt(clth[0]*clth[2])
# provide list of specs to be computed, and/or options
lth = arange(2, lmax+1)
spec = None # ['EB', 'TE', 'TB']
temp = False
polar = True
corr = False
spec = ['EE','BB'] #'EB', 'TE', 'TB']
pixwin = True
# ellbins = np.append(2, arange(50, lmax + 2, dell))
ellbins = np.arange(2, lmax + 2, dell)
ellbins[-1] = lmax + 1
ellbins[-1] = lmax+1
muKarcmin = 1.0
fwhm = 0.5
##############################
#input model
clth = np.array(hp.read_cl(MODELFILE))
Clthshape = zeros(((6,)+shape(clth)[1:]))
Clthshape[:4] = clth
clth = Clthshape
lth = arange(2, 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)
# # Large scale mask
# mask[abs(90 - rad2deg(t)) < 30] = False
#mask[abs(90 - rad2deg(t)) < 20] = False
mask[(90 - rad2deg(t)) < glat] = False
# # Small scale mask (do not forget to change dell)
mask[(90 - rad2deg(t)) < 78] = False
fsky = np.mean(mask)
print("fsky=%.2g %%" % (100*fsky))
npix = sum(mask)
##############################
fwhm = 0.5
bl = hp.gauss_beam(deg2rad(fwhm), lmax=Slmax + 1)
stokes, spec, istokes, ispecs = simulation.getstokes(
spec=spec, polar=polar, temp=temp, corr=corr)
stokes, spec, istokes, ispecs = getstokes( spec=spec)
print(stokes, spec, istokes, ispecs)
nspec = len(spec)
nstoke = len(stokes)
# ############## Compute ds_dcb ###############
ip = arange(hp.nside2npix(nside))
ipok = ip[mask]
Pl, S = libcov.compute_ds_dcb(ellbins, nside, ipok, bl, clth, Slmax, spec=spec,
pixwin=pixwin, timing=True, MC=False)
# ############## Initialise xqml class ###############
muKarcmin = 1.0
# ############## Generate Noise ###############
pixvar = muKarcmin2var(muKarcmin, nside)
varmap = ones((nstoke * npix)) * pixvar
NoiseVar = np.diag(varmap)
noise = (randn(len(varmap)) * varmap**0.5).reshape(nstoke, -1)
esti = xqml.xQML(mask, ellbins, clth, NA=NoiseVar, NB=NoiseVar, Pl=Pl,
S=S, fwhm=fwhm, spec=spec, temp=temp, polar=polar, corr=corr)
# esti.construct_esti(NoiseVar, NoiseVar)
V = esti.get_covariance(cross=True)
# ############## Initialise xqml class ###############
start = timeit.default_timer()
esti = xqml.xQML(mask, ellbins, clth, NA=NoiseVar, NB=NoiseVar, lmax=lmax, fwhm=fwhm, spec=spec, SymCompress=True)
s1 = timeit.default_timer()
print( "Init: %d sec" % (s1-start))
s2 = timeit.default_timer()
print( "construct_esti: %d sec" % (s2-s1))
ellval = esti.lbin()
# ############## Compute Analytical variance ###############
V = esti.get_covariance(cross=True )
Va = esti.get_covariance(cross=False)
# ############## Construct MC ###############
# ############## Construct MC ###############
allcla = []
allcl = []
# allcmb = []
esti.construct_esti(NoiseVar, NoiseVar)
fpixw = extrapolpixwin(nside, lmax+2, pixwin=pixwin)
start = timeit.default_timer()
t = []
bl = hp.gauss_beam(deg2rad(fwhm), lmax=Slmax)
fpixw = extrapolpixwin(nside, Slmax, pixwin=pixwin)
for n in np.arange(nsimu):
# progress_bar(n, nsimu, timeit.default_timer() - start)
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))
pixwin=False, lmax=Slmax, fwhm=0.0, new=True, verbose=False))
cmbm = cmb[istokes][:, mask]
# allcmb.append(cmbm)
dmA = cmbm + (randn(nstoke * npix) * sqrt(varmap)).reshape(nstoke, npix)
dmB = cmbm + (randn(nstoke * npix) * sqrt(varmap)).reshape(nstoke, npix)
s1 = timeit.default_timer()
allcl.append(esti.get_spectra(dmA, dmB))
t.append( timeit.default_timer() - s1)
allcla.append(esti.get_spectra(dmA))
# myS = cov(np.array(allcmb).reshape(nsimu, -1), rowvar=False)
print( "get_spectra: %d sec" % mean(t))
hcl = mean(allcl, 0)
scl = std(allcl, 0)
hcla = mean(allcla, 0)
scla = std(allcla, 0)
# ############## Plot results ###############
figure(figsize=[12, 8])
clf()
Delta = (ellbins[1:] - ellbins[:-1])/2.
subplot(3, 2, 1)
title("Cross")
plot(lth, (lth*(lth+1)/2./np.pi)[:, None]*clth[ispecs][:, lth].T, '--k')
[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])
for s in np.arange(nspec)]
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])
semilogy()
ylabel(r"$D_\ell$")
legend(loc=4, frameon=True)
......@@ -138,59 +154,51 @@ legend(loc=4, frameon=True)
subplot(3, 2, 2)
title("Auto")
plot(lth,(lth*(lth+1)/2./np.pi)[:, None]*clth[ispecs][:, lth].T, '--k')
[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]) for s in np.arange(nspec)]
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])
semilogy()
subplot(3, 2, 3)
# cosmic = sqrt(2./(2 * lth + 1)) / mean(mask) * clth[ispecs][:, lth]
# plot(lth, cosmic.transpose(), '--k')
[plot(ellval, scl[s], '--', color='C%i' % s, label=r"$\sigma^{%s}_{\rm MC}$" %
spec[s]) for s in np.arange(nspec)]
[plot(ellval, sqrt(diag(V)).reshape(nspec, -1)[s], 'o', color='C%i' % ispecs[s])
for s in np.arange(nspec)]
for s in np.arange(nspec):
plot(ellval, scl[s], '--', color='C%i' % ispecs[s], label=r"$\sigma^{%s}_{\rm MC}$" % spec[s])
plot(ellval, sqrt(diag(V)).reshape(nspec, -1)[s], 'o', color='C%i' % ispecs[s])
ylabel(r"$\sigma(C_\ell)$")
semilogy()
# legend(loc=4, frameon=True)
subplot(3, 2, 4)
[plot(ellval, scla[s], ':', color='C%i' % s, label=r"$\sigma^{%s}_{\rm MC}$" %
spec[s]) for s in np.arange(nspec)]
[plot(ellval, sqrt(diag(Va)).reshape(nspec, -1)[s], 'o', color='C%i' % ispecs[s])
for s in np.arange(nspec)]
for s in np.arange(nspec):
plot(ellval, scla[s], ':', color='C%i' % ispecs[s], label=r"$\sigma^{%s}_{\rm MC}$" % spec[s])
plot(ellval, sqrt(diag(Va)).reshape(nspec, -1)[s], 'o', color='C%i' % ispecs[s])
semilogy()
subplot(3, 2, 5)
[plot(ellval, (hcl[s]-P.dot(clth[ispecs[s]][lth]))/(scl[s]/sqrt(nsimu)), '--o',
color='C%i' % ispecs[s]) for s in np.arange(nspec)]
for s in np.arange(nspec):
plot(ellval, (hcl[s]-esti.BinSpectra(clth)[s])/(scl[s]/sqrt(nsimu)), '--o', color='C%i' % ispecs[s])
ylabel(r"$R[C_\ell]$")
xlabel(r"$\ell$")
ylim(-3, 3)
grid()
subplot(3, 2, 6)
[plot(ellval, (hcla[s]-P.dot(clth[ispecs[s]][lth]))/(scla[s]/sqrt(nsimu)), '--o',
color='C%i' % ispecs[s]) for s in np.arange(nspec)]
for s in np.arange(nspec):
plot(ellval, (hcla[s]-esti.BinSpectra(clth)[s])/(scla[s]/sqrt(nsimu)), '--o', color='C%i' % ispecs[s])
xlabel(r"$\ell$")
ylim(-3, 3)
grid()
show()
savefig("../Plots/Git/Nside%i_dell%i_fsky%.3g_spec%s.png" %
(nside, dell, fsky, "".join(spec)))
if __name__ == "__main__":
"""
Run the doctest using
## if __name__ == "__main__":
## """
## Run the doctest using
python simulation.py
## python simulation.py
If the tests are OK, the script should exit gracefuly, otherwise the
failure(s) will be printed out.
"""
import doctest
if np.__version__ >= "1.14.0":
np.set_printoptions(legacy="1.13")
doctest.testmod()
## If the tests are OK, the script should exit gracefuly, otherwise the
## failure(s) will be printed out.
## """
## import doctest
## if np.__version__ >= "1.14.0":
## np.set_printoptions(legacy="1.13")
## doctest.testmod()
......@@ -7,6 +7,8 @@ from __future__ import division
import numpy as np
from .xqml_utils import symarray
def Pl(ds_dcb):
"""
......@@ -136,12 +138,8 @@ def El(invCAA, invCBB, Pl):
"""
lmax = len(Pl)
lrange = np.arange(lmax)
npix = len(invCAA)
El = np.array(
[np.dot(np.dot(invCAA, Pl[l]), invCBB) for l in lrange]
).reshape((lmax, npix, npix))
El = np.asarray([np.dot(np.dot(invCAA, symarray(P)), invCBB) for P in Pl])
return El
......@@ -149,14 +147,14 @@ def CrossWindowFunction(El, Pl):
"""
Compute mode-mixing matrix (Tegmark's window matrix)
Wll = Trace[invCAA.Pl.invCBB.Pl] = Trace[El.Pl]
Parameters
----------
El : ndarray of floats
Quadratic parameter matrices such that yl = dA.El.dB.T
Pl : ndarray of floats
Rescaled normalize Legendre polynomials dS/dCl
Returns
----------
Wll : 2D square matrix array of floats
......@@ -171,12 +169,13 @@ def CrossWindowFunction(El, Pl):
[110 390 670]
[134 478 822]]
"""
lmax = len(El)
lrange = np.arange((lmax))
nl = len(El)
# 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)
Wll = np.asarray(
[np.sum(E * symarray(P)) for E in El for P in Pl]
).reshape(nl,nl)
return Wll
......@@ -212,9 +211,9 @@ def CrossWindowFunctionLong(invCAA, invCBB, Pl):
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)
Wll = np.asarray(
[np.sum(np.dot(np.dot(invCAA, symarray(Pi)), invCBB) * symarray(Pj)) for Pi in Pl for Pj in Pl]
).reshape(lmax, lmax)
return Wll
......@@ -242,12 +241,12 @@ def CrossGisherMatrix(El, CAB):
[ 701 2205 3709]
[1181 3709 6237]]
"""
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)
nl = len(El)
El_CAB = np.asarray([np.dot(CAB, E) for E in El])
GAB = np.asarray(
[np.sum(Ei * Ej.T) for Ei in El_CAB for Ej in El_CAB]
).reshape(nl,nl)
return GAB
......@@ -278,7 +277,7 @@ def CrossGisherMatrixLong(El, CAB):
"""
lmax = len(El)
lrange = np.arange(lmax)
GAB = np.array(
GAB = np.asarray(
[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
......@@ -305,9 +304,7 @@ def yQuadEstimator(dA, dB, El):
>>> print(yQuadEstimator(dA, dB, El))
[1360788 3356628 5352468]
"""
npix = len(dA)
lrange = np.arange((len(El)))
y = np.array([dA.dot(El[l]).dot(dB) for l in lrange])
y = np.asarray([dA.dot(E).dot(dB) for E in El])
return y
......@@ -349,27 +346,8 @@ def biasQuadEstimator(NoiseN, El):
----------
???
"""
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])
return np.asarray([np.sum(NoiseN * E) for E in El])
def CovAB(invWll, GAB):
......@@ -398,39 +376,11 @@ def CovAB(invWll, GAB):
return covAB
# def GetCorr(F):
# """
# Return correlation of covariance matrix
# corr(F) = F[i, j] / sqrt(F[i, i]*F[j, j])
# Parameters
# ----------
# F : 2D square matrix array (n,n)
# Covariance matrix
# Returns
# ----------
# Corr : 2D square matrix array (n,n)
# Correlation matrix from F
# Example
# ----------
# >>> CorrF = GetCorr(np.matrix('1 .5; .5 4'))
# >>> print(CorrF) # doctest: +NORMALIZE_WHITESPACE
# [[ 1. 0.25]
# [ 0.25 1. ]]
# """
# nbins = len(F)
# Corr = np.array(
# [F[i, j] / (F[i, i]*F[j, j])**.5 for i in np.arange(nbins)
# for j in np.arange(nbins)]).reshape(nbins, nbins)
# return Corr
if __name__ == "__main__":
"""
Run the doctest using
python simulation.py
python estimators.py
If the tests are OK, the script should exit gracefuly, otherwise the
failure(s) will be printed out.
......
......@@ -12,14 +12,14 @@ import healpy as hp
import math
from scipy import special
from simulation import extrapolpixwin
from xqml_utils import getstokes, progress_bar
from .simulation import extrapolpixwin
from .xqml_utils import getstokes, progress_bar, symarray, GetBinningMatrix
import _libcov as clibcov
def compute_ds_dcb( ellbins, nside, ipok, bl, clth, Slmax, spec,
pixwin=False, timing=False, MC=0, Sonly=False, openMP=True):
pixwin=False, timing=False, MC=0, Sonly=False, openMP=True, SymCompress=False):
"""
Compute the Pl = dS/dCl matrices.
......@@ -63,12 +63,12 @@ def compute_ds_dcb( ellbins, nside, ipok, bl, clth, Slmax, spec,
Example
----------
>>> Pl, S = compute_ds_dcb(
... np.array([2,4,5,8]), 2, np.array([0,1,4,10,11]), np.arange(10),
... clth=np.arange(60).reshape(6,-1), Slmax=8,
... np.array([2,4,5,10]), 2, np.array([0,1,4,10,11]), np.arange(10),
... clth=np.arange(60).reshape(6,-1), Slmax=9,
... spec=["TT", "EE", "BB", "TE", "EB", "TB"],
... pixwin=True, timing=False, MC=0, Sonly=False)
>>> print(round(np.sum(Pl),5), round(np.sum(S),5))
(1128.61613, 29501.83265)
(1149.18805, 23675.10206)
"""
if Slmax < ellbins[-1]-1:
print("WARNING : Slmax < lmax")
......@@ -102,27 +102,48 @@ def compute_ds_dcb( ellbins, nside, ipok, bl, clth, Slmax, spec,
ellbins, nside, ipok, allcosang, bl, clth, Slmax, MC,
polar=polar, temp=temp, corr=corr, pixwin=pixwin, timing=timing)
elif openMP:
fpixwin = extrapolpixwin(nside, Slmax+1, pixwin)
fpixwin = extrapolpixwin(nside, Slmax, pixwin)
bell = np.array([bl*fpixwin]*4)[:Slmax+1].ravel()
stokes, spec, istokes, ispecs = getstokes(spec)
nbins = (len(ellbins)-1)*len(spec)
npix = len(ipok)*len(istokes)
Pl = np.ndarray( nbins*npix**2)
## print( "size dSdC: %d" % len(Pl))
clibcov.dSdC( nside, len(istokes), ellbins, ipok, bell, Pl)
Pl = Pl.reshape( nbins, npix, npix)
S = SignalCovMatrix(Pl,clth[ispecs][:,2:int(ellbins[-1])].ravel())
P, Q, ell, ellval = GetBinningMatrix(ellbins, Slmax)
S = SignalCovMatrix(Pl,np.array([P.dot(clth[isp,2:int(ellbins[-1])]) for isp in ispecs]).ravel())
else:
Pl, S = compute_PlS(
ellbins, nside, ipok, allcosang, bl, clth, Slmax,
spec=spec, pixwin=pixwin, timing=timing)
print( "Total time (npix=%d): %.1f sec" % (len(ipok),timeit.default_timer()-start))
#print( "Total time (npix=%d): %.1f sec" % (len(ipok),timeit.default_timer()-start))
if SymCompress:
Pl = np.array([symarray(P).packed for P in Pl])
return Pl, S
def SignalCovMatrix(Pl, model, SymCompress=False):
"""
Compute correlation matrix S = sum_l Pl*Cl
Parameters
----------
clth : ndarray of floats
Array containing fiducial CMB spectra (unbinned).
"""
# Return scalar product btw Pl and the fiducial spectra.
return np.sum([symarray(P) for P in Pl] * model[:, None, None], 0)
def compute_PlS(
ellbins, nside, ipok, allcosang, bl, clth, Slmax,
spec, pixwin=True, timing=False):
......@@ -162,20 +183,19 @@ def compute_PlS(
Example
----------
>>> Pl, S = compute_PlS(np.array([2,3,4,7]),4,ipok=np.array([0,1,3]),
>>> Pl, S = compute_PlS(np.array([2,3,4,10]),4,ipok=np.array([0,1,3]),
... allcosang=np.linspace(0,1,15).reshape(3,-1), bl=np.arange(13),
... clth=np.arange(6*13).reshape(6,-1), Slmax=11,
... clth=np.arange(6*13).reshape(6,-1), Slmax=9,
... spec=["TT", "EE", "BB", "TE", "EB", "TB"], pixwin=True, timing=False)
>>> print(round(np.sum(Pl),5), round(np.sum(S),5))
(14.38023, -4125.5419)
(-429.8591, -17502.8982)
>>> Pl, S = compute_PlS(np.array([2,3,4,7]),4,ipok=np.array([0,1,3]),
>>> Pl, S = compute_PlS(np.array([2,3,4,10]),4,ipok=np.array([0,1,3]),
... allcosang=np.linspace(0,1,15).reshape(3,-1), bl=np.arange(13),
... clth=np.arange(6*13).reshape(6,-1), Slmax=11,
... clth=np.arange(6*13).reshape(6,-1), Slmax=9,
... spec=["TT", "EE", "BB", "TE", "EB", "TB"], pixwin=False, timing=False)
>>> print(round(np.sum(Pl),5), round(np.sum(S),5))
(17.44139, 3086.13152)
(-756.35517, -31333.69722)
"""
lmax = ellbins[-1]
......@@ -205,7 +225,7 @@ def compute_PlS(
rpix = np.array(hp.pix2vec(nside, ipok))
ll = np.arange(Slmax+1)
fpixwin = extrapolpixwin(nside, Slmax+1, pixwin)
fpixwin = extrapolpixwin(nside, Slmax, pixwin)
norm = (2*ll[2:]+1)/(4.*np.pi)*(fpixwin[2:]**2)*(bl[2:Slmax+1]**2)
clthn = clth[:, 2: Slmax+1]
masks = []
......@@ -417,20 +437,19 @@ def compute_S(
Example
----------
>>> S = compute_S(np.array([2,3,4,7]),4,ipok=np.array([0,1,3]),
... allcosang=np.linspace(0,1,15).reshape(3,-1), bl=np.arange(13),
... clth=np.arange(6*13).reshape(6,-1), Slmax=11,
>>> S = compute_S(np.array([2,3,4,10]),4,ipok=np.array([0,1,3]),