Commit 4ee3536d authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

llp1 binning debug

parent 67a66070
......@@ -71,7 +71,7 @@ def lmax2nside(lmax, level=2):
return nside
def alm2cov (almtab):
def alm2cov(almtab):
""" Compute second order statistics of spherical harmonic coefficients.
Parameters
......@@ -89,10 +89,10 @@ def alm2cov (almtab):
lmax = int(round((sqrt(1+8*nalm)-3) / 2) )
ncross = nin*(nin+1)/2
cov = zeros((nin, nin, lmax+1))
for i in range(nin):
for j in range(i, nin):
cov[i,j,:] = healpy.alm2cl(almtab[i,:], almtab[j,:])
cov[i,j,:] = hp.alm2cl(almtab[i,:], almtab[j,:])
cov[j,i,:] = cov[i,j,:]
#covtab = zeros((ncross, lmax+1))
......@@ -101,8 +101,7 @@ def alm2cov (almtab):
#for ell in range(lmax+1):
# cov[:,:,ell] = utils.vec2mat(covtab[:,ell])
nmode = bin.get_nmode(bin.uniformbin (0, lmax, 1))
nmode = bin.get_nmode(bin.uniformbin (0, lmax, 1))#, btype)
return cov, nmode
......
......@@ -20,32 +20,35 @@ spectral windows computations.
from numpy import *
import numpy as np
import spherelib_bin
def get_binMatrix(bins, btype='flat'):
def get_binMatrix(bins, btype='flat', lmax=None):
""" Compute binning matrix
"""
lmin = int(bins[ 0, 0])
lmax = int(bins[-1,-1])
lmin = 0
if lmax is None:
lmax = int(bins[-1,-1])
Q = shape(bins)[0]
Mbins = zeros ((Q, lmax+1-lmin))
nmodes = get_nmode(bins, btype)
for q in range(Q):
ellmin = int(bins[q,0])
ellmax = int(bins[q,1])
rg = arange(ellmin,ellmax+1).astype(int32)
rg = arange(ellmin,ellmax+1).astype(int64)
if btype=='flat':
Mbins[q,rg-lmin] = 2*rg+1
##squeeze(lnmodes[rg]/(1.0*sum(lnmodes[rg])))
elif btype=='llp1':
Mbins[q,rg-lmin] = rg*(rg+1)
Mbins[q,rg-lmin] = rg*(rg+1)*(2*rg+1)
nrm = sum(Mbins[q,rg-lmin])
print nrm
Mbins[q,rg-lmin] /= 1.0*nrm
return Mbins
def bin_stat (covmat, bin, nmode=None):
def bin_stat(covmat, bin, nmode=None):
""" Perform the averaging of spectral covariance matrices over flat bins.
Matrices are weighted by their number of modes which is supposed to equal to 2.ell+1.
......@@ -64,27 +67,39 @@ def bin_stat (covmat, bin, nmode=None):
ndet = covmat.shape[0]
bincov = zeros((ndet, ndet, nbin))
lsize = zeros((nbin,1))
for q in range(nbin):
lmin = bin[q,0]
lmax = bin[q,1]
rg = arange(lmin,lmax+1).astype(int32)
rg = arange(lmin,lmax+1).astype(int32)
if nmode is None:
nmodes = 2*rg+1
else:
nmodes = nmode[rg]
tmpcov = covmat[:,:,rg].reshape(ndet*ndet, len(rg))
tmpcov = tmpcov * (ones((ndet*ndet,1))*nmodes)
nmodes = nmodes[rg]
tmpcov = covmat[:,:,rg].reshape(ndet*ndet, len(rg))
tmpcov = tmpcov*(ones((ndet*ndet,1))*nmodes)
tmpcov = sum(tmpcov, axis=1).reshape(ndet, ndet)
ll = sum(nmodes,axis=0)
bincov[:,:,q] = tmpcov * (1.0/(1.0*ll))
ll = sum(nmodes,axis=0)
bincov[:,:,q] = tmpcov * (1.0/(1.0*ll))
lsize[q] = ll
return [bincov,lsize]
def bin_stat_llp1(covmat, Mbin):
"""
"""
nbin = Mbin.shape[0]
ndet = covmat.shape[0]
bincov = zeros((ndet, ndet, nbin))
lsize = zeros((nbin,1))
for q in range(nbin):
tmpcov = covmat.reshape(ndet*ndet, -1)
tmpcov = tmpcov*(ones((ndet*ndet,1))*Mbin[q,:])
tmpcov = sum(tmpcov, axis=1).reshape(ndet, ndet)
bincov[:,:,q] = tmpcov
return bincov
def uniformbin (lstart, lstop, deltal):
""" Return description of flat bins of constant length.
......@@ -173,7 +188,7 @@ def get_lmean (bin):
return lmean
def get_nmode (bin):
def get_nmode(Mbin, btype='flat'):
""" Return the number of modes of flat bins.
Number of modes is supposed to equal to \f$\2ell+1\f$ for each mulitpole.
......@@ -186,11 +201,18 @@ def get_nmode (bin):
-------
array-like, shape (nbin, 1), number of modes of each bin.
"""
nbin = len(bin)
nmode = zeros((nbin, 1))
for q in range(nbin):
rg = arange(bin[q,0] , bin[q,1]+1)
nmode[q] = sum (2 * rg +1, axis=0 )
nbin = Mbin.shape[0]
nmode = np.zeros((nbin, 1))
if btype=='flat':
for q in range(nbin):
rg = np.arange(Mbin[q,0],Mbin[q,1]+1)
nmode[q] = np.sum(2*rg+1, axis=0)
else:
for q in range(nbin):
nmode[q] = np.sum(Mbin[q,:]**2)**2
rg = np.arange(Mbin.shape[1])
nmode[q] /= 1.*np.sum(Mbin[q,:]**4/(2*rg+1))
return nmode
......
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