Commit 888e25c1 authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

doctest

parent fe422e59
......@@ -27,350 +27,261 @@ import spherelib_bin
def get_binMatrix(bins, btype='flat', lmax=None):
""" Compute binning matrix
>>> get_binMatrix(np.array([(1,2),(3,4)]))
array([[ 0. , 0.375 , 0.625 , 0. , 0. ],
[ 0. , 0. , 0. , 0.4375, 0.5625]])
>>> get_binMatrix(np.array([(1,2),(3,4)]), btype="llp1")
array([[ 0. , 0.16666667, 0.83333333, 0. , 0. ],
[ 0. , 0. , 0. , 0.31818182, 0.68181818]])
"""
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)
lmax = int(bins[-1,-1])
Q = np.shape(bins)[0]
Mbins = np.zeros((Q, lmax+1))
for q in range(Q):
ellmin = int(bins[q,0])
ellmax = int(bins[q,1])
rg = arange(ellmin,ellmax+1).astype(int64)
rg = np.arange(ellmin,ellmax+1).astype(int64)
if btype=='flat':
Mbins[q,rg-lmin] = 2*rg+1
Mbins[q,rg] = 2*rg+1
elif btype=='llp1':
Mbins[q,rg-lmin] = rg*(rg+1)*(2*rg+1)
nrm = sum(Mbins[q,rg-lmin])
Mbins[q,rg-lmin] /= 1.0*nrm
Mbins[q,rg] = rg*(rg+1)*(2*rg+1)
else:
return "unknown btype (should be flat or llp1)"
nrm = np.sum(Mbins[q,rg])
Mbins[q,rg] /= 1.0*nrm
return Mbins
def bin_stat(covmat, bin, nmode=None):
def bin_covariance(covmat, bins, 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.
Parameters
----------
covmat : array-like, shape (m,m,lmax+1). Concatenation of lmax+1 symmetric matrix of size m.
bin : array-like, shape (nbin, 2). Two column vectors defining limits (lmin and lmax) of each bin.
nmode: array-like, shape (lmax+1, 1). Number of modes used to build covariance matrices.
Returns
-------
array-like, shape (m,m,nbin). Binned matrices.
"""
nbin = bin.shape[0]
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)
if nmode is None:
nmodes = 2*rg+1
else:
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))
lsize[q] = ll
return [bincov,lsize]
def bin_stat_llp1(covmat, Mbin):
Matrices are weighted by their number of modes ie 2.ell+1.
or using the binning matrix encoded weights
>>> hR,w = bin_covariance(np.ones((2,2,5)), np.array([(1,2),(3,4)]))
>>> print hR[:,:,0]
[[ 1. 1.]
[ 1. 1.]]
>>> print w
[ 8. 16.]
"""
"""
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
nbin = bins.shape[0]
ndet = covmat.shape[0]
bincov = np.zeros((ndet, ndet, nbin))
lsize = np.zeros((nbin))
if bins.shape[1]==2:
for q in range(nbin):
lmin = bins[q,0]
lmax = bins[q,1]
rg = np.arange(lmin,lmax+1).astype(int32)
nmodes = 2*rg+1 if nmode is None else nmode[rg]
tmpcov = covmat[:,:,rg].reshape(ndet*ndet, len(rg))
tmpcov = tmpcov*(np.ones((ndet*ndet,1))*nmodes)
tmpcov = np.sum(tmpcov, axis=1).reshape(ndet, ndet)
ll = np.sum(nmodes,axis=0)
bincov[:,:,q] = tmpcov*(1.0/(1.0*ll))
lsize[q] = ll
return [bincov,lsize]
else: # use binning matrix
for q in range(nbin):
tmpcov = covmat.reshape(ndet*ndet, -1)
tmpcov = tmpcov*(np.ones((ndet*ndet,1))*bins[q,:])
tmpcov = np.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.
Parameters
----------
lstart : scalar, first multipole value.
lstop : scalar, last multipole value.
deltal : scalar, length of bins.
Returns
-------
array-like, shape (nbin, 2). Two column vectors defining limits (lmin and lmax) of each bin.
>>> uniformbin(0,10,2)
array([[ 0., 1.],
[ 2., 3.],
[ 4., 5.],
[ 6., 7.],
[ 8., 10.]])
"""
L = lstop - lstart + 1
nbin = ceil(L / deltal).astype(int32)
nbin = np.ceil(L / deltal).astype(int32)
if nbin==0:
nbin = 1
bin = zeros((nbin, 2))
bin[0,0] = lstart
bins = np.zeros((nbin, 2))
bins[0,0] = lstart
for q in range(1,nbin):
bin[q , 0] = bin[q-1, 0] + deltal
bin[q-1, 1] = bin[q , 0] - 1
bin[-1,1] = lstop
return bin
def wg2bin (lstart, lstop, mind=1):
""" Return description of flat bins as defined by WG2.
Parameters
----------
lstart : scalar, first multipole value.
lstop : scalar, last multipole value.
mind : scalar, minimum bin length.
Returns
-------
array-like, shape (nbin, 2). Two column vectors defining limits (lmin and lmax) of each bin.
"""
wg2_lmin = 0
wg2_lmax = 4000
wg2_ltrans = array([wg2_lmin,10,30,150,420,1200,2500,4000])
wg2_deltal = array([1,2,5,10,20,30,40])
wg2_deltal[wg2_deltal<mind] = mind
## cut ltrans between lmin and lmax
l1 = where(wg2_ltrans>=lstart)[0]
l2 = where(wg2_ltrans<=lstop)[0]
idx = intersect1d(l1, l2)
ltrans = unique1d(hstack((hstack(([max(lstart-1,0)], wg2_ltrans[idx])), [lstop])))
## cut deltal between lmin and lmax
if (len(idx)!=0):
deltal = hstack((array([wg2_deltal[idx[0]-1]]), wg2_deltal[idx]))
deltal = deltal[0:len(ltrans)-1]
else:
deltal = wg2_deltal[where(wg2_ltrans<lstop)[0]][-1]
## build uniform bins
nbin = size(deltal)
if nbin==1:
deltal = [deltal]
bins = zeros((0,2))
for t in range(nbin):
b = uniformbin(ltrans[t]+1,ltrans[t+1],deltal[t])
bins = vstack((bins,b))
bins[q,0] = bins[q-1,0] + deltal
bins[q-1,1] = bins[q,0] - 1
bins[-1,1] = lstop
return bins
def get_lmean (bin):
def get_lmean(bins):
""" Return mean multipole values of flat bins.
Parameters
----------
bin : array-like, shape (nbin, 2). Two column vectors defining limits (lmin and lmax) of each bin.
Returns
-------
array-like, shape (nbin, 1), mean multipole values of each bin.
>>> get_lmean(uniformbin(0,10,2))
array([ 0.5, 2.5, 4.5, 6.5, 9. ])
"""
nbin = len(bin)
lmean = zeros((nbin, 1))
for q in range(nbin):
lmean[q] = (bin[q,0] + bin[q,1] ) / 2
return lmean
nbin = bins.shape[0]
return np.array([(bins[q,0] + bins[q,1] ) / 2. for q in range(nbin)])
def get_nmode(Mbin, btype='flat'):
""" Return the number of modes of flat bins.
def get_nmode(bins, btype='flat'):
""" Return the number of modes.
Number of modes is supposed to equal to \f$\2ell+1\f$ for each mulitpole.
Parameters
----------
bin : array-like, shape (nbin, 2). Two column vectors defining limits (lmin and lmax) of each bin.
Returns
-------
array-like, shape (nbin, 1), number of modes of each bin.
>>> get_nmode(uniformbin(0,10,2))
array([ 4., 12., 20., 28., 57.])
>>> get_nmode(get_binMatrix(uniformbin(0,10,2)))
array([ 4., 12., 20., 28., 57.])
>>> get_nmode(get_binMatrix(uniformbin(0,10,2), btype='llp1'), btype='llp1')
array([ 3. , 8.70140745, 16.71118644, 25.16864036, 49.13735764])
"""
nbin = Mbin.shape[0]
nmode = np.zeros((nbin, 1))
nbin = bins.shape[0]
nmode = np.zeros((nbin))
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:
if bins.shape[1]==2:
rg = np.arange(bins[q,0],bins[q,1]+1)
nmode[q] = np.sum(2*rg+1, axis=0)
else:
nmode[q] = np.sum(bins[q,:])**2
rg = np.arange(bins.shape[1])
nmode[q] /= 1.*np.sum(bins[q,:]**2/(2*rg+1))
elif btype=='llp1':
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))
nmode[q] = np.sum(bins[q,:]**2)**2
rg = np.arange(bins.shape[1])
nmode[q] /= 1.*np.sum(bins[q,:]**4/(2*rg+1))
else:
return 'unknown btype (should be flat or llp1)'
return nmode
def gaussbeam(fwhm_arcmin, lmax):
"""Return the transfer window function of a Gaussian beam up to lmax.
Parameters
----------
fwhm_arcmin : scalar. Full width at half maximum in arcminute.
lmax : maximum multipole value.
Returns
-------
array-like, shape (lmax+1, 1), beam transfert function.
>>> gaussbeam(30, 100)[0:5]
array([ 1. , 0.99998627, 0.9999588 , 0.9999176 , 0.99986268])
"""
fwhm_rad = (array(fwhm_arcmin)*2*pi/60.) /360.
fwhm_rad = reshape(fwhm_rad,(1,size(fwhm_rad)))
l = arange(lmax+1)
l = reshape(l,(size(l),1))
beam = sqrt(exp(-l*(l+1)/8./log(2) * fwhm_rad**2))
fwhm_rad = (np.array(fwhm_arcmin)*2*np.pi/60.) /360.
l = np.arange(lmax+1)
beam = np.sqrt(np.exp(-l*(l+1)/8./np.log(2) * fwhm_rad**2))
return beam
def binavg(x, binsize):
"""Bin a vector into bins of size binsize.
>>> binavg(np.arange(0,10), 2)
array([ 0., 2., 4., 6., 8.])
"""
bin = arange(0,size(x), binsize)
nbin = size(bin)
res = zeros(nbin)
for b in arange(nbin -1):
res[b] = mean(x[bin[b]:bin[b+1]-1])
res[nbin-1] = mean(x[bin[b]:])
bins = np.arange(0,size(x), binsize)
nbin = np.size(bins)
res = np.zeros(nbin)
for b in np.arange(nbin-1):
res[b] = np.mean(x[bins[b]:bins[b+1]-1])
res[nbin-1] = np.mean(x[bins[nbin-1]:-1])
return res
def smooth_cl (cl, dloverl=0.1):
""" Return smoothed power spectrum.
def smooth_cl(cl, dloverl=0.1, width=None, n=0):
"""Return smoothed power spectrum.
Power spectrum is averaged over window of size dloverl (in %)
Power spectrum is either averaged over a sliding window of size
dloverl (in %), or a fixed window (controlled by the width
parameter) in a iterative way (n being the number of iteration.
Parameters
----------
cl : array-like, shape(lmax+1, 1). Power spectrum.
dloverl : averaging window size.
Returns
-------
array-like, shape(lmax+1, 1). Smoothed power spectrum.
>>> smooth_cl(np.array([0,0,1,1,2,2]))
array([ 0. , 0.5, 1. , 1.5, 2. , 2. ])
>>> smooth_cl(np.array([0,0,1,1,2,2]), width=2, n=3)
array([ 0.33796296, 0.4845679 , 0.80864198, 1.19135802, 1.5154321 ,
1.66203704])
"""
lmax = size(cl) - 1
ell = arange(0,lmax+1)
clsmooth = zeros(shape(cl));
lmax = len(cl) - 1
ell = np.arange(0,lmax+1)
clsmooth = np.zeros(np.shape(cl));
for l in range(lmax+1):
llo = max(int(round(ell[l]*(1-dloverl/2.0))), 0)
lhi = min(int(round(ell[l]*(1+dloverl/2.0))), lmax)
if llo==lhi and lhi!=lmax+1:
lhi=llo+1
clsmooth[l] = mean(cl[llo:lhi+1])
return clsmooth
if width is None:
llo = max(int(np.round(ell[l]*(1-dloverl/2.0))), 0)
lhi = min(int(np.round(ell[l]*(1+dloverl/2.0))), lmax)
else:
llo = max(ell[l]-width/2, 0)
lhi = min(ell[l]+width/2, lmax)
if llo==lhi:
lhi=llo+1
clsmooth[l] = np.mean(cl[llo:lhi+1])
if n>0:
return smooth_cl(clsmooth, dloverl=dloverl, width=width, n=n-1)
else:
return clsmooth
def bin_cl (cl, bin):
def bin_cl (cl, bins):
""" Return binned power spectrum.
Power spectrum is averaged over flat bins.
Parameters
----------
cl : array-like, shape(lmax+1, 1). Power spectrum.
bin : array-like, shape (nbin, 2). Two column vectors defining limits (lmin and lmax) of each bin.
Returns
-------
array-like, shape(nbin, 1). Binned power spectrum.
>>> bin_cl(np.arange(10), uniformbin(0,9,2))
array([ 0.5, 2.5, 4.5, 6.5, 8.5])
"""
nbin = shape(bin)[0]
cq = zeros((nbin, 1))
nbin = np.shape(bins)[0]
cq = np.zeros((nbin))
for q in range(nbin):
lmin = bin[q,0]
lmax = bin[q,1]
rg = arange(lmin,lmax+1).astype(int32)
lsize = float(size(rg))
cq[q] = sum(cl[rg])/lsize
return squeeze(cq)
lmin = bins[q,0]
lmax = bins[q,1]
rg = np.arange(lmin,lmax+1).astype(int32)
lsize = np.float(np.size(rg))
cq[q] = np.sum(cl[rg])/lsize
return cq
def unbin_cl (cq, bin):
def unbin_cl (cq, bins):
""" Return unbinned power spectrum.
Supposed that power spectrum has been averaged over flat bins.
Parameters
----------
cq : array-like, shape(nbin, 1). Binned power spectrum.
bin : array-like, shape (nbin, 2). Two column vectors defining limits (lmin and lmax) of each bin.
Returns
-------
array-like, shape(lmax+1, 1). Unbinned power spectrum.
>>> b = uniformbin(0,9,2)
>>> unbin_cl(bin_cl(np.arange(10), b), b)
array([ 0.5, 0.5, 2.5, 2.5, 4.5, 4.5, 6.5, 6.5, 8.5, 8.5])
"""
nbin = shape(bin)[0]
lmax = bin[nbin-1,1]
cl = zeros((lmax+1, 1))
nbin = np.shape(bins)[0]
lmax = bins[nbin-1,1]
cl = np.zeros((lmax+1))
for q in range(nbin):
lmin = bin[q,0]
lmax = bin[q,1]
rg = arange(lmin,lmax+1).astype(int32)
cl[rg] = cq[q]
return squeeze(cl)
lmin = bins[q,0]
lmax = bins[q,1]
rg = np.arange(lmin,lmax+1).astype(int32)
cl[rg] = cq[q]
return cl
def make_prolate (bands, rule=1, thetas=[]):
""" Build prolate window functions.
Parameters
----------
bands: list of integers, windows limits
rule: integer, prolate rule between 1 (acos ( 1 - solid /(2*pi))), 2( 1 / sqrt(lmean)) ,3(1/lmean>1/(2*delta) ? 1/lmean : 1/(2*delta))
thetas: list of float: caps opening values (rule set to 0 in this case).
Returns
-------
array-like (lmax+1, nwin)
bands is a list of integers defining windows limits.
rule can be 1 for (acos ( 1 - solid /(2*pi))),
2 for ( 1 / sqrt(lmean)) ,
3 for (1/lmean>1/(2*delta) ? 1/lmean : 1/(2*delta))
thetas is a list of caps opening values (rule set to 0 in this case).
>>> make_prolate([0,5,10])
array([ 1. , 0.99241115, 0.97735349, 0.95506465, 0.92589514,
0.89030127, 0.84883603, 0.80213815, 0.75091939, 0.69595057,
0.63804657])
"""
nwin = len(bands)-2
if thetas:
rule = 0
thetas = array(thetas, dtype=float32)
thetas = np.array(thetas, dtype=float32)
else:
thetas = ones(nwin, dtype=float32)
bands = array(bands, dtype=int32)
thetas = np.ones(nwin, dtype=float32)
bands = np.array(bands, dtype=int32)
lmax = bands[-1]
bins = zeros((nwin, lmax+1), dtype=float64)
spherelib_bin._make_prolate (bands, bins, thetas, rule)
return bins.transpose()
bins = np.zeros((nwin, lmax+1), dtype=float64)
spherelib_bin._make_prolate(bands, bins, thetas, rule)
return np.squeeze(bins.T)
def make_spline (bands, order=3):
""" Build spline window functions.
Parameters
----------
bands: list of integers, windows limits
order: integer, spline order, default is 3
Returns
-------
array-like (lmax+1, nwin)
>>> make_spline([0,5,10])
array([ 1. , 1. , 1. , 1. , 1. ,
1. , 0.98318666, 0.84273839, 0.53832332, 0.1826034 , 0. ])
"""
order = 3
lmax = bands[-1]
nwin = len(bands)-2
bands = array(bands, dtype=int32)
bins = zeros((nwin, lmax+1), dtype=float64)
spherelib_bin._make_spline (bands, bins, order)
return bins.transpose()
return np.squeeze(bins.T)
def norm_cl (cl, type="sphere"):
""" Return re-normalized spectra.
......@@ -394,3 +305,7 @@ def norm_cl (cl, type="sphere"):
raise Exception ("norm_cl: norm must be one of the following: \"sphere\", \"line\", \"max\" ");
return cl*( norm**(-0.5) )
if __name__ == "__main__":
import doctest
doctest.testmod()
......@@ -72,7 +72,7 @@ def mask_correction (cov, mask, corr="fsky", nmode=None, mll=None):
elif corr=="mll":
lmax = shape(cov)[2]-1
if mll is None:
mll = map.mask2mll (mask, lmax)
mll = map.mask2mll(mask, lmax)
fsky = sum(mll[0,:])
nmap = shape(cov)[0]
for i in range(nmap):
......
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