Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Open sidebar
xQML
xQML
Commits
211cac05
Commit
211cac05
authored
Mar 04, 2020
by
Matthieu Tristram
Browse files
Change "timing" to "verbose"
parent
04a55c92
Changes
2
Hide whitespace changes
Inline
Side-by-side
Showing
2 changed files
with
53 additions
and
47 deletions
+53
-47
xqml/libcov.py
xqml/libcov.py
+33
-33
xqml/xqml.py
xqml/xqml.py
+20
-14
No files found.
xqml/libcov.py
View file @
211cac05
...
...
@@ -19,7 +19,7 @@ import _libcov as clibcov
def
compute_ds_dcb
(
ellbins
,
nside
,
ipok
,
bl
,
clth
,
Slmax
,
spec
,
pixwin
=
True
,
timing
=
False
,
MC
=
0
,
Sonly
=
False
,
openMP
=
True
):
pixwin
=
True
,
verbose
=
False
,
MC
=
0
,
Sonly
=
False
,
openMP
=
True
):
"""
Compute the Pl = dS/dCl matrices.
...
...
@@ -43,7 +43,7 @@ def compute_ds_dcb( ellbins, nside, ipok, bl, clth, Slmax, spec,
pixwin : bool
If True, multiplies the beam window function by the pixel
window function. Default: True
timing
: bool
verbose
: bool
If True, displays timmer. Default: False
MC : int
If not 0, computes Pl using Monte-Carlo method from MC simulations.
...
...
@@ -66,7 +66,7 @@ def compute_ds_dcb( ellbins, nside, ipok, bl, clth, Slmax, spec,
... 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)
... pixwin=True,
verbose
=False, MC=0, Sonly=False)
>>> print(round(np.sum(Pl),5), round(np.sum(S),5))
(1149.18805, 23675.10206)
"""
...
...
@@ -79,7 +79,7 @@ def compute_ds_dcb( ellbins, nside, ipok, bl, clth, Slmax, spec,
allcosang
[
allcosang
>
1
]
=
1.0
allcosang
[
allcosang
<
-
1
]
=
-
1.0
start
=
timeit
.
default_timer
()
t
start
=
timeit
.
default_timer
()
clth
=
np
.
asarray
(
clth
)
if
len
(
clth
)
==
4
:
...
...
@@ -92,17 +92,17 @@ def compute_ds_dcb( ellbins, nside, ipok, bl, clth, Slmax, spec,
if
MC
:
S
=
S_bins_MC
(
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
MC
,
spec
,
pixwin
=
pixwin
,
timing
=
timing
)
pixwin
=
pixwin
,
verbose
=
verbose
)
else
:
S
=
compute_S
(
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
spec
,
pixwin
=
pixwin
,
timing
=
timing
)
pixwin
=
pixwin
,
verbose
=
verbose
)
return
S
if
MC
:
Pl
,
S
=
covth_bins_MC
(
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
MC
,
spec
,
pixwin
=
pixwin
,
timing
=
timing
)
spec
,
pixwin
=
pixwin
,
verbose
=
verbose
)
elif
openMP
:
fpixwin
=
extrapolpixwin
(
nside
,
Slmax
,
pixwin
)
bell
=
np
.
array
([
bl
*
fpixwin
]
*
4
)[:
Slmax
+
1
].
ravel
()
...
...
@@ -119,10 +119,10 @@ def compute_ds_dcb( ellbins, nside, ipok, bl, clth, Slmax, spec,
else
:
Pl
,
S
=
compute_PlS
(
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
spec
=
spec
,
pixwin
=
pixwin
,
timing
=
timing
)
spec
=
spec
,
pixwin
=
pixwin
,
verbose
=
verbose
)
if
timing
:
print
(
"
Total time
(npix=%d): %.1f sec"
%
(
len
(
ipok
),
timeit
.
default_timer
()
-
start
))
if
verbose
:
print
(
"
Construct Pl
(npix=%d): %.1f sec"
%
(
len
(
ipok
),
timeit
.
default_timer
()
-
t
start
))
return
Pl
,
S
...
...
@@ -146,7 +146,7 @@ def SignalCovMatrix(Pl, model):
def
compute_PlS
(
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
spec
,
pixwin
=
True
,
timing
=
False
):
spec
,
pixwin
=
True
,
verbose
=
False
):
"""
Computes Legendre polynomes Pl = dS/dCb and signal matrix S.
...
...
@@ -170,7 +170,7 @@ def compute_PlS(
pixwin : bool
If True, multiplies the beam window function by the pixel
window function. Default: True
timing
: bool
verbose
: bool
If True, displays timmer. Default: False
Returns
...
...
@@ -186,14 +186,14 @@ def compute_PlS(
>>> 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=9,
... spec=["TT", "EE", "BB", "TE", "EB", "TB"], pixwin=True,
timing
=False)
... spec=["TT", "EE", "BB", "TE", "EB", "TB"], pixwin=True,
verbose
=False)
>>> print(round(np.sum(Pl),5), round(np.sum(S),5))
(-429.8591, -17502.8982)
>>> 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=9,
... spec=["TT", "EE", "BB", "TE", "EB", "TB"], pixwin=False,
timing
=False)
... spec=["TT", "EE", "BB", "TE", "EB", "TB"], pixwin=False,
verbose
=False)
>>> print(round(np.sum(Pl),5), round(np.sum(S),5))
(-756.35517, -31333.69722)
"""
...
...
@@ -236,10 +236,10 @@ def compute_PlS(
Pl
=
np
.
zeros
((
nspec
*
nbins
,
nsto
*
npi
,
nsto
*
npi
))
S
=
np
.
zeros
((
nsto
*
npi
,
nsto
*
npi
))
start
=
timeit
.
default_timer
()
t
start
=
timeit
.
default_timer
()
for
i
in
np
.
arange
(
npi
):
if
timing
:
if
verbose
:
progress_bar
(
i
,
npi
)
for
j
in
np
.
arange
(
i
,
npi
):
cos_chi
=
allcosang
[
i
,
j
]
...
...
@@ -403,7 +403,7 @@ def compute_PlS(
def
compute_S
(
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
spec
,
pixwin
=
True
,
timing
=
False
):
spec
,
pixwin
=
True
,
verbose
=
False
):
"""
Computes signal matrix S.
...
...
@@ -427,7 +427,7 @@ def compute_S(
pixwin : bool
If True, multiplies the beam window function by the pixel
window function. Default: True
timing
: bool
verbose
: bool
If True, displays timmer. Default: False
Returns
...
...
@@ -440,14 +440,14 @@ def compute_S(
>>> S = compute_S(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(10),
... clth=np.arange(6*13).reshape(6,-1), Slmax=9,
... spec=["TT", "EE", "BB", "TE", "EB", "TB"], pixwin=True,
timing
=False)
... spec=["TT", "EE", "BB", "TE", "EB", "TB"], pixwin=True,
verbose
=False)
>>> print(round(np.sum(S),5))
-17502.8982
>>> S = compute_S(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=9,
... spec=["TT", "EE", "BB", "TE", "EB", "TB"], pixwin=False,
timing
=False)
... spec=["TT", "EE", "BB", "TE", "EB", "TB"], pixwin=False,
verbose
=False)
>>> print(round(np.sum(S),5))
-31333.69722
"""
...
...
@@ -484,9 +484,9 @@ def compute_S(
clthn
=
clth
[:,
2
:
Slmax
+
1
]
S
=
np
.
zeros
((
nsto
*
npi
,
nsto
*
npi
))
start
=
timeit
.
default_timer
()
t
start
=
timeit
.
default_timer
()
for
i
in
np
.
arange
(
npi
):
if
timing
:
if
verbose
:
progress_bar
(
i
,
npi
)
for
j
in
np
.
arange
(
i
,
npi
):
if
temp
:
...
...
@@ -578,7 +578,7 @@ def compute_S(
def
covth_bins_MC
(
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
nsimu
,
spec
,
pixwin
=
False
,
timing
=
False
):
spec
,
pixwin
=
False
,
verbose
=
False
):
"""
Can be particularly slow on sl7.
To be enhanced and extended to TT and correlations.
...
...
@@ -603,7 +603,7 @@ def covth_bins_MC(
pixwin : bool
If True, multiplies the beam window function by the pixel
window function. Default: True
timing
: bool
verbose
: bool
If True, displays timmer. Default: False
Returns
...
...
@@ -621,7 +621,7 @@ def covth_bins_MC(
>>> Pl, S = covth_bins_MC(np.array([2,3,4]), 4, ipok=np.array([0,1,3]),
... allcosang=np.linspace(0,1,13), bl=np.arange(13),
... clth=np.arange(4*13).reshape(4,-1), Slmax=11, nsimu=100,
... spec, pixwin=True,
timing
=False)
... spec, pixwin=True,
verbose
=False)
>>> print(round(np.sum(Pl),5), round(np.sum(S),5))
(12.68135, 406990.12056)
...
...
@@ -630,7 +630,7 @@ def covth_bins_MC(
>>> Pl, S = covth_bins_MC(np.array([2,3,4]), 4, ipok=np.array([0,1,3]),
... allcosang=np.linspace(0,1,13), bl=np.arange(13),
... clth=np.arange(4*13).reshape(4,-1), Slmax=11, nsimu=100,
... spec, pixwin=True,
timing
=False)
... spec, pixwin=True,
verbose
=False)
>>> print(round(np.sum(Pl),5), round(np.sum(S),5))
(42.35592, 7414.01784)
"""
...
...
@@ -655,7 +655,7 @@ def covth_bins_MC(
masks
.
append
((
ll
[:]
>=
minell
[
i
])
&
(
ll
[:]
<=
maxell
[
i
]))
masks
=
np
.
array
(
masks
)
npix
=
len
(
ipok
)
start
=
timeit
.
default_timer
()
t
start
=
timeit
.
default_timer
()
norm
=
bl
[
0
:
Slmax
+
1
]
**
2
*
fpixwin
[
0
:
Slmax
+
1
]
**
2
if
polar
:
...
...
@@ -671,9 +671,9 @@ def covth_bins_MC(
ClthOne
[
l
,
4
]
=
masks
[
l
%
nbins
]
*
norm
Pl
=
np
.
zeros
((
nspec
*
(
nbins
),
2
*
npix
,
2
*
npix
))
start
=
timeit
.
default_timer
()
t
start
=
timeit
.
default_timer
()
for
l
in
np
.
arange
((
nspec
*
nbins
)):
if
timing
:
if
verbose
:
progress_bar
(
l
,
nspec
*
(
nbins
))
data
=
[
...
...
@@ -704,7 +704,7 @@ def covth_bins_MC(
ClthOne
[
l
]
=
masks
[
l
]
*
norm
Pl
=
np
.
zeros
(((
nbins
),
npix
,
npix
))
for
l
in
np
.
arange
((
nbins
)):
if
timing
:
if
verbose
:
progress_bar
(
l
,
nspec
*
(
nbins
))
Pl
[
l
]
=
np
.
cov
(
np
.
array
([
...
...
@@ -735,7 +735,7 @@ def covth_bins_MC(
def
S_bins_MC
(
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
nsimu
,
spec
,
pixwin
=
False
,
timing
=
False
):
spec
,
pixwin
=
False
,
verbose
=
False
):
"""
Can be particularly slow on sl7 !
To be enhanced and extended to TT and correlations
...
...
@@ -764,7 +764,7 @@ def S_bins_MC(
pixwin : bool
If True, multiplies the beam window function by the pixel
window function. Default: True
timing
: bool
verbose
: bool
If True, displays timmer. Default: False
Returns
...
...
@@ -791,7 +791,7 @@ def S_bins_MC(
masks
.
append
((
ll
[:]
>=
minell
[
i
])
&
(
ll
[:]
<=
maxell
[
i
]))
masks
=
np
.
array
(
masks
)
npix
=
len
(
ipok
)
start
=
timeit
.
default_timer
()
t
start
=
timeit
.
default_timer
()
norm
=
bl
[
0
:
Slmax
+
1
]
**
2
*
fpixwin
[
0
:
Slmax
+
1
]
**
2
if
polar
:
...
...
xqml/xqml.py
View file @
211cac05
...
...
@@ -34,7 +34,7 @@ __all__ = ['xQML','Bins']
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
=
[
'EE'
,
'BB'
],
pixwin
=
True
):
S
=
None
,
fwhm
=
0.
,
bell
=
None
,
spec
=
[
'EE'
,
'BB'
],
pixwin
=
True
,
verbose
=
True
):
"""
Parameters
----------
...
...
@@ -95,10 +95,11 @@ class xQML(object):
nbin
=
len
(
self
.
ellbins
)
nmem
=
self
.
nspec
*
nbin
*
(
self
.
nstokes
*
self
.
npix
)
**
2
toGb
=
1024.
*
1024.
*
1024.
print
(
"xQML"
)
print
(
"spec: "
,
spec
)
print
(
"nbin: "
,
nbin
)
print
(
"Memset: %.2f Gb (%d,%d,%d,%d)"
%
(
8.
*
nmem
/
toGb
,
self
.
nspec
,
nbin
,
self
.
nstokes
,
self
.
npix
))
if
verbose
:
print
(
"xQML"
)
print
(
"spec: "
,
spec
)
print
(
"nbin: "
,
nbin
)
print
(
"Memset: %.2f Gb (%d,%d,%d,%d)"
%
(
8.
*
nmem
/
toGb
,
self
.
nspec
,
nbin
,
self
.
nstokes
,
self
.
npix
))
# If Pl is given by the user, just load it, and then compute the signal
# covariance using the fiducial model.
...
...
@@ -108,7 +109,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
=
Tru
e
,
MC
=
False
,
openMP
=
True
)
self
.
spec
,
pixwin
=
self
.
pixwin
,
verbose
=
verbos
e
,
MC
=
False
,
openMP
=
True
)
else
:
self
.
Pl
=
Pl
if
S
is
None
:
...
...
@@ -117,17 +118,17 @@ class xQML(object):
self
.
S
=
S
if
NA
is
not
None
:
self
.
construct_esti
(
NA
=
NA
,
NB
=
NB
)
self
.
construct_esti
(
NA
=
NA
,
NB
=
NB
,
verbose
=
verbose
)
def
compute_dSdC
(
self
,
clth
,
lmax
=
None
,
timing
=
True
,
MC
=
False
,
openMP
=
True
):
def
compute_dSdC
(
self
,
clth
,
lmax
=
None
,
verbose
=
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
)
self
.
spec
,
pixwin
=
self
.
pixwin
,
verbose
=
verbose
,
MC
=
MC
,
openMP
=
openMP
)
return
(
self
.
Pl
,
self
.
S
)
def
construct_esti
(
self
,
NA
,
NB
=
None
):
def
construct_esti
(
self
,
NA
,
NB
=
None
,
verbose
=
False
):
"""
Compute the inverse of the datasets pixel covariance matrices C,
the quadratic matrix parameter E, and inverse of the window
...
...
@@ -144,17 +145,19 @@ class xQML(object):
self
.
cross
=
NB
is
not
None
self
.
NA
=
NA
self
.
NB
=
NB
if
self
.
cross
else
NA
tstart
=
timeit
.
default_timer
()
# Invert (signalA + noise) matrix
invCa
=
pd_inv
(
self
.
S
+
self
.
NA
)
# Invert (signalB + noise) matrix
invCb
=
pd_inv
(
self
.
S
+
self
.
NB
)
# Compute E
using Eq...
# Compute E
l = Ca^-11.Pl.Cb^-1 (long)
self
.
El
=
El
(
invCa
,
invCb
,
self
.
Pl
)
# Finally compute invW by inverting
...
# Finally compute invW by inverting
(long)
self
.
invW
=
linalg
.
inv
(
CrossWindowFunction
(
self
.
El
,
self
.
Pl
))
# Compute bias for auto
...
...
@@ -162,6 +165,9 @@ class xQML(object):
# self.bias = biasQuadEstimator(self.NA, self.El)
self
.
bias
=
biasQuadEstimator
(
self
.
NA
,
self
.
El
)
if
verbose
:
print
(
"Construct estimator: %.1f sec"
%
(
timeit
.
default_timer
()
-
tstart
))
def
get_spectra
(
self
,
mapA
,
mapB
=
None
):
"""
...
...
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment