Skip to content

GitLab

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

Remove Symmarray

parent f2b1113e
Hide whitespace changes
Inline Side-by-side
Showing with 65 additions and 56 deletions
TestEsti.py
View file @ e72cd538
......@@ -68,22 +68,23 @@ lth = arange(2, lmax+1)
##############################
# Create mask
mask = hp.ud_grade(hp.read_map("../Masks/mask_DX12d_galpolco_30pc_ns2048.fits"), nside)
mask[mask > 0.5] = True
mask[mask < 0.5] = False
mask = np.array(mask, bool)
# t, p = hp.pix2ang(nside, range(hp.nside2npix(nside)))
# mask = np.ones(hp.nside2npix(nside), bool)
t, p = hp.pix2ang(nside, range(hp.nside2npix(nside)))
mask = np.ones(hp.nside2npix(nside), bool)
# import random
# random.shuffle(mask)
#mask[abs(90 - rad2deg(t)) < 20] = False
mask[(90 - rad2deg(t)) < glat] = False
if exp == "Big":
mask[abs(90 - rad2deg(t)) < glat] = False
elif exp == "Small":
mask[(90 - rad2deg(t)) < glat] = False
fsky = np.mean(mask)
print("fsky=%.2g %%" % (100*fsky))
npix = sum(mask)
print("fsky=%.2g %% (npix=%d)" % (100*fsky,npix))
toGB = 1024. * 1024. * 1024.
emem = 8.*(npix*2*npix*2) * ( len(lth)*2 ) / toGB
print("mem=%.2g Gb" % emem)
##############################
......@@ -105,16 +106,16 @@ noise = (randn(len(varmap)) * varmap**0.5).reshape(nstoke, -1)
# ############## 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)
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))
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)
s2 = timeit.default_timer()
print( "construct covariance: %d sec" % (s2-s1))
# ############## Construct MC ###############
......
src/libcov.c
View file @ e72cd538
#include "libcov.h"
// void CrossWindowFunction( int nl, int npix, double* El, double* Pl, double* Wll)
// {
// memset( Wll, 0., nl*nl * sizeof(double));
// #pragma omp parallel default(none) shared(nl, El, Pl, Wll)
// {
// #pragma omp for schedule(dynamic)
// for( int l1=0; l1<nl; l1++) {
// for( int l2=0; l2<nl; l2++) {
// for( int p=0; p<npix*npix; p++)
// Wll[l1*nl+l2] += El[l1][
// }
// }
// }
// }
void build_dSdC( int nside, int nstokes, int npix, int nbin, long *ellbins, long *ipix, double *bl, double* dSdC)
{
const int lmax=ellbins[nbin]-1;
......
xqml/estimators.py
View file @ e72cd538
......@@ -7,8 +7,6 @@ from __future__ import division
import numpy as np
from .xqml_utils import symarray
def Pl(ds_dcb):
"""
......@@ -138,16 +136,8 @@ def El(invCAA, invCBB, Pl):
"""
<<<<<<< HEAD
El = np.asarray([np.dot(np.dot(invCAA, symarray(P)), invCBB) for P in Pl])
El = np.asarray([np.dot(np.dot(invCAA, P), invCBB) for P in 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.array([np.linalg.multi_dot((invCAA, Pl[l], invCBB)) for l in lrange]).reshape((lmax, npix, npix))
>>>>>>> 4fd1b451ce149e52585aba6331adb20ead8bcfb5
return El
......@@ -180,9 +170,7 @@ def CrossWindowFunction(El, Pl):
nl = len(El)
# pas de transpose car symm
Wll = np.asarray(
[np.sum(E * symarray(P)) for E in El for P in Pl]
).reshape(nl,nl)
Wll = np.asarray( [np.sum(E * P) for E in El for P in Pl] ).reshape(nl,nl)
return Wll
......@@ -220,7 +208,7 @@ def CrossWindowFunctionLong(invCAA, invCBB, Pl):
lrange = np.arange((lmax))
# Pas de transpose car symm
Wll = np.asarray(
[np.sum(np.dot(np.dot(invCAA, symarray(Pi)), invCBB) * symarray(Pj)) for Pi in Pl for Pj in Pl]
[np.sum(np.dot(np.dot(invCAA, Pi), invCBB) * Pj) for Pi in Pl for Pj in Pl]
).reshape(lmax, lmax)
return Wll
......
xqml/libcov.py
View file @ e72cd538
......@@ -13,13 +13,13 @@ import math
from scipy import special
from .simulation import extrapolpixwin
from .xqml_utils import getstokes, progress_bar, symarray, GetBinningMatrix
from .xqml_utils import getstokes, progress_bar, 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, SymCompress=False):
pixwin=False, timing=False, MC=0, Sonly=False, openMP=True):
"""
Compute the Pl = dS/dCl matrices.
......@@ -117,17 +117,14 @@ def compute_ds_dcb( ellbins, nside, ipok, bl, clth, Slmax, spec,
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))
if SymCompress:
Pl = np.array([symarray(P).packed for P in Pl])
#print( "Total time (npix=%d): %.1f sec" % (len(ipok),timeit.default_timer()-start))
return Pl, S
def SignalCovMatrix(Pl, model, SymCompress=False):
def SignalCovMatrix(Pl, model):
"""
Compute correlation matrix S = sum_l Pl*Cl
......@@ -137,7 +134,7 @@ def SignalCovMatrix(Pl, model, SymCompress=False):
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)
return np.sum(Pl * model[:, None, None], 0)
......
xqml/xqml.py
View file @ e72cd538
......@@ -35,7 +35,7 @@ class xQML(object):
""" Main class to handle the spectrum estimation """
def __init__( self, mask, bins, clth, NA=None, NB=None, lmax=None, Pl=None,
S=None, fwhm=0., bell=None, spec=None, temp=False, polar=True, corr=False,
pixwin=True, SymCompress=False):
pixwin=True):
"""
Parameters
----------
......@@ -108,7 +108,7 @@ class xQML(object):
self.Pl, self.S = compute_ds_dcb(
self.ellbins, self.nside, self.ipok,
self.bl, clth, self.Slmax,
self.spec, pixwin=self.pixwin, timing=True, MC=False, openMP=True, SymCompress=SymCompress)
self.spec, pixwin=self.pixwin, timing=True, MC=False, openMP=True)
else:
self.Pl = Pl
if S is None:
......@@ -121,12 +121,12 @@ class xQML(object):
if NA is not None:
self.construct_esti(NA=NA, NB=NB)
def compute_dSdC( self, clth, lmax=None, timing=True, MC=False, openMP=True, SymCompress=True):
def compute_dSdC( self, clth, lmax=None, timing=True, MC=False, openMP=True):
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.spec, pixwin=self.pixwin, timing=timing, MC=MC, openMP=openMP, SymCompress=SymCompress)
self.spec, pixwin=self.pixwin, timing=timing, MC=MC, openMP=openMP)
return( self.Pl, self.S)
def construct_esti(self, NA, NB=None):
......@@ -154,21 +154,25 @@ class xQML(object):
invCb = pd_inv(self.S + self.NB)
# Compute E using Eq...
self.E = El(invCa, invCb, self.Pl)
self.El = El(invCa, invCb, self.Pl)
if not self.cross:
self.bias = biasQuadEstimator(self.NA, self.E)
# if not self.cross:
# self.bias = biasQuadEstimator(self.NA, self.El)
self.bias = biasQuadEstimator(self.NA, self.El)
# Finally compute invW by inverting...
# s0 = timeit.default_timer()
print("Compute invW")
self.invW = linalg.inv(CrossWindowFunction(self.E, self.Pl))
self.invW = linalg.inv(CrossWindowFunction(self.El, self.Pl))
# s1 = timeit.default_timer()
# self.invW = linalg.inv(CrossWindowFunctionLong(invCa, invCb, self.Pl))
# s2 = timeit.default_timer()
# print( "CrossWindowFunction: %d sec" % (s1-s0))
# print( "CrossWindowFunctionLong: %d sec" % (s2-s1))
#Clean
del(self.Pl)
del(invCa)
del(invCb)
def get_spectra(self, mapA, mapB=None):
"""
......@@ -195,12 +199,10 @@ class xQML(object):
cond_sizeB = np.size(mapB) == self.nstokes * self.npix
dB = mapB if cond_sizeB else mapB[self.istokes][:,self.mask]
yl = yQuadEstimator(dA.ravel(), dB.ravel(), self.E)
yl = yQuadEstimator(dA.ravel(), dB.ravel(), self.El)
cl = ClQuadEstimator(self.invW, yl)
else:
if self.bias is None:
self.bias = biasQuadEstimator(self.NA, self.E)
yl = yQuadEstimator(dA.ravel(), dA.ravel(), self.E) - self.bias
yl = yQuadEstimator(dA.ravel(), dA.ravel(), self.El) - self.bias
cl = ClQuadEstimator(self.invW, yl)
# Return the reshaped set of cls
......@@ -220,9 +222,9 @@ class xQML(object):
# # Do Gll' = S^-1.El.S^-1.El'
cross = self.cross if cross is None else cross
if cross:
G = CrossGisherMatrix(self.E, self.S)
G = CrossGisherMatrix(self.El, self.S)
else:
G = CrossGisherMatrix(self.E, self.S + self.NA)
G = CrossGisherMatrix(self.El, self.S + self.NA)
# # Do V = W^-1.G.W^-1 + W^-1
V = CovAB(self.invW, G)
......
xqml/xqml_utils.py
View file @ e72cd538
......@@ -383,7 +383,7 @@ class symarray(np.ndarray):
m = int(m)
else:
raise ValueError('One dimensional input length must be a triangular number.')
A = np.zeros((m, m))
for i in range(m):
A[i, i:] = l[:(m-i)]
......@@ -391,8 +391,8 @@ class symarray(np.ndarray):
l = l[(m-i):]
obj = np.asarray(A).view(cls)
obj.packed = p
elif len(obj.shape) == 2:
if obj.shape[0] != obj.shape[1]:
raise ValueError('Two dimensional input must be a square matrix.')
......@@ -400,12 +400,12 @@ class symarray(np.ndarray):
for i in range(obj.shape[0]):
packed_out.append(obj[i, i:])
obj.packed = np.concatenate(packed_out)
else:
raise ValueError('Input array must be 1 or 2 dimensional.')
return obj
def __array_finalize__(self, obj):
if obj is None: return
self.packed = getattr(obj, 'packed', None)
......
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