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
ecb32a95
Commit
ecb32a95
authored
Jan 21, 2020
by
Matthieu Tristram
Browse files
Update computation for EB, TE, TB
parent
d6698174
Changes
7
Hide whitespace changes
Inline
Side-by-side
Showing
7 changed files
with
166 additions
and
134 deletions
+166
-134
src/_libcov.pyx
src/_libcov.pyx
+3
-1
src/libcov.c
src/libcov.c
+105
-50
src/libcov.h
src/libcov.h
+3
-2
xqml/estimators.py
xqml/estimators.py
+2
-2
xqml/libcov.py
xqml/libcov.py
+29
-23
xqml/xqml.py
xqml/xqml.py
+11
-15
xqml/xqml_utils.py
xqml/xqml_utils.py
+13
-41
No files found.
src/_libcov.pyx
View file @
ecb32a95
...
...
@@ -11,17 +11,19 @@ np.import_array()
cdef
extern
from
"libcov.h"
:
void
build_dSdC
(
int
nside
,
int
nstokes
,
int
npix
,
int
inl
,
long
*
ellbins
,
long
*
ipix
,
double
*
bl
,
double
*
dSdC
)
void
build_dSdC
(
int
nside
,
int
nstokes
,
int
npix
,
int
inl
,
long
*
ispec
,
long
*
ellbins
,
long
*
ipix
,
double
*
bl
,
double
*
dSdC
)
void
dlss
(
double
X
,
int
s1
,
int
s2
,
int
lmax
,
double
*
d
)
int
polrotangle
(
double
*
ri
,
double
*
rj
,
double
*
cos2a
,
double
*
sin2a
)
def
dSdC
(
nside
,
nstokes
,
np
.
ndarray
[
long
,
ndim
=
1
,
mode
=
"c"
]
ispec
not
None
,
np
.
ndarray
[
long
,
ndim
=
1
,
mode
=
"c"
]
ellbins
not
None
,
np
.
ndarray
[
long
,
ndim
=
1
,
mode
=
"c"
]
ipix
not
None
,
np
.
ndarray
[
double
,
ndim
=
1
,
mode
=
"c"
]
bl
not
None
,
np
.
ndarray
[
double
,
ndim
=
1
,
mode
=
"c"
]
dSdC
not
None
):
build_dSdC
(
nside
,
nstokes
,
ipix
.
shape
[
0
],
ellbins
.
shape
[
0
]
-
1
,
<
long
*>
np
.
PyArray_DATA
(
ispec
),
<
long
*>
np
.
PyArray_DATA
(
ellbins
),
<
long
*>
np
.
PyArray_DATA
(
ipix
),
<
double
*>
np
.
PyArray_DATA
(
bl
),
...
...
src/libcov.c
View file @
ecb32a95
...
...
@@ -21,24 +21,26 @@
// }
void
build_dSdC
(
int
nside
,
int
nstokes
,
int
npix
,
int
nbin
,
long
*
ellbins
,
long
*
ipix
,
double
*
bl
,
double
*
dSdC
)
void
build_dSdC
(
int
nside
,
int
nstokes
,
int
npix
,
int
nbin
,
long
*
ispec
,
long
*
ellbins
,
long
*
ipix
,
double
*
bl
,
double
*
dSdC
)
{
const
int
lmax
=
ellbins
[
nbin
]
-
1
;
const
int
lmax1
=
lmax
+
1
;
const
int
ns
=
3
,
nspecall
=
4
;
const
int
nspec
=
nstokes2nspec
(
nstokes
);
const
int
ns
=
3
,
nspecall
=
6
;
const
int
nspec
=
ispec2nspec
(
ispec
);
// const int nspec = nstokes2nspec(nstokes);
const
int64_t
npixtot
=
npix
*
nstokes
;
//
fprintf( stdout, "lmax=%d\n", lmax);
//
fprintf( stdout, "npix=%d\n", npix);
//
fprintf( stdout, "nstokes=%d\n", nstokes);
//
fprintf( stdout, "nspec=%d\n", nspec);
fprintf
(
stdout
,
"lmax=%d
\n
"
,
lmax
);
fprintf
(
stdout
,
"npix=%d
\n
"
,
npix
);
fprintf
(
stdout
,
"nstokes=%d
\n
"
,
nstokes
);
fprintf
(
stdout
,
"nspec=%d
\n
"
,
nspec
);
int64_t
ntot
=
nspec
*
nbin
*
npixtot
*
npixtot
;
//
fprintf( stdout, "memset (%d,%d,%d,%d = %lld)...\n", nspec,nbin,nstokes,npix,ntot);
fprintf
(
stdout
,
"memset (%d,%d,%d,%d = %lld)...
\n
"
,
nspec
,
nbin
,
nstokes
,
npix
,
ntot
);
memset
(
dSdC
,
0
.,
ntot
*
sizeof
(
double
));
fprintf
(
stdout
,
"ispec = (%d,%d,%d,%d,%d,%d)
\n
"
,
ispec
[
0
],
ispec
[
1
],
ispec
[
2
],
ispec
[
3
],
ispec
[
4
],
ispec
[
5
]);
#pragma omp parallel default(none) shared(stdout,nbin,nside,npix,nstokes,dSdC,ipix,bl,ellbins)
#pragma omp parallel default(none) shared(stdout,nbin,nside,npix,nstokes,dSdC,ipix,bl,ellbins
,ispec
)
{
int
s
=
0
;
double
vr
[
3
],
vc
[
3
];
...
...
@@ -49,6 +51,9 @@ void build_dSdC( int nside, int nstokes, int npix, int nbin, long *ellbins, long
dSdCpix
[
il
]
=
(
double
*
)
calloc
(
ns
*
ns
,
sizeof
(
double
));
if
(
dSdCpix
[
il
]
==
NULL
)
EXIT_INFO
(
-
1
,
"Problem allocation dSdCpix (l=%d)...
\n
"
,
il
);
}
int
sI
=
0
,
sQ
=
1
,
sU
=
2
;
if
(
nstokes
==
2
)
{
sQ
=
0
;
sU
=
1
;
}
#pragma omp for schedule(dynamic)
/* loop on local pix to build (3x3) blocks */
...
...
@@ -58,42 +63,60 @@ void build_dSdC( int nside, int nstokes, int npix, int nbin, long *ellbins, long
for
(
int
rpix
=
0
;
rpix
<
npix
;
rpix
++
)
{
pix2vec_ring
(
nside
,
ipix
[
rpix
],
vr
);
QML_compute_dSdC
(
vr
,
vc
,
lmax
,
nstokes
,
dSdCpix
);
QML_compute_dSdC
(
vr
,
vc
,
lmax
,
ispec
,
dSdCpix
);
for
(
int
ib
=
0
;
ib
<
nbin
;
ib
++
)
{
// int l=ellbins[ib];
for
(
int
l
=
ellbins
[
ib
];
l
<=
ellbins
[
ib
+
1
]
-
1
;
l
++
)
{
s
=
0
;
if
(
nstokes
!
=
2
)
{
dSdC
(
s
*
nbin
+
ib
,
0
*
npix
+
cpix
,
0
*
npix
+
rpix
)
+=
dSdCpix
[
0
*
lmax1
+
l
][
0
*
ns
+
0
]
*
bl
[
0
*
lmax1
+
l
]
*
bl
[
0
*
lmax1
+
l
];
//TT on II
if
(
ispec
[
0
]
=
=
1
)
{
dSdC
(
s
*
nbin
+
ib
,
sI
*
npix
+
cpix
,
sI
*
npix
+
rpix
)
+=
dSdCpix
[
0
*
lmax1
+
l
][
0
*
ns
+
0
]
*
bl
[
0
*
lmax1
+
l
]
*
bl
[
0
*
lmax1
+
l
];
//TT on II
s
++
;
}
//EE-BB
if
(
nstokes
>
1
)
{
int
s1
=
s
;
int
s2
=
s
+
1
;
dSdC
(
s
*
nbin
+
ib
,
s1
*
npix
+
cpix
,
s1
*
npix
+
rpix
)
+=
dSdCpix
[
1
*
lmax1
+
l
][
1
*
ns
+
1
]
*
bl
[
1
*
lmax1
+
l
]
*
bl
[
1
*
lmax1
+
l
];
//EE on QQ
dSdC
(
s
*
nbin
+
ib
,
s2
*
npix
+
cpix
,
s1
*
npix
+
rpix
)
+=
dSdCpix
[
1
*
lmax1
+
l
][
2
*
ns
+
1
]
*
bl
[
1
*
lmax1
+
l
]
*
bl
[
1
*
lmax1
+
l
];
//EE on QU
dSdC
(
s
*
nbin
+
ib
,
s1
*
npix
+
cpix
,
s2
*
npix
+
rpix
)
+=
dSdCpix
[
1
*
lmax1
+
l
][
1
*
ns
+
2
]
*
bl
[
1
*
lmax1
+
l
]
*
bl
[
1
*
lmax1
+
l
];
//EE on UQ
dSdC
(
s
*
nbin
+
ib
,
s2
*
npix
+
cpix
,
s2
*
npix
+
rpix
)
+=
dSdCpix
[
1
*
lmax1
+
l
][
2
*
ns
+
2
]
*
bl
[
1
*
lmax1
+
l
]
*
bl
[
1
*
lmax1
+
l
];
//EE on UU
if
(
ispec
[
1
]
==
1
||
ispec
[
2
]
==
1
)
{
dSdC
(
s
*
nbin
+
ib
,
sQ
*
npix
+
cpix
,
sQ
*
npix
+
rpix
)
+=
dSdCpix
[
1
*
lmax1
+
l
][
1
*
ns
+
1
]
*
bl
[
1
*
lmax1
+
l
]
*
bl
[
1
*
lmax1
+
l
];
//EE on QQ
dSdC
(
s
*
nbin
+
ib
,
sU
*
npix
+
cpix
,
sQ
*
npix
+
rpix
)
+=
dSdCpix
[
1
*
lmax1
+
l
][
2
*
ns
+
1
]
*
bl
[
1
*
lmax1
+
l
]
*
bl
[
1
*
lmax1
+
l
];
//EE on QU
dSdC
(
s
*
nbin
+
ib
,
sQ
*
npix
+
cpix
,
sU
*
npix
+
rpix
)
+=
dSdCpix
[
1
*
lmax1
+
l
][
1
*
ns
+
2
]
*
bl
[
1
*
lmax1
+
l
]
*
bl
[
1
*
lmax1
+
l
];
//EE on UQ
dSdC
(
s
*
nbin
+
ib
,
sU
*
npix
+
cpix
,
sU
*
npix
+
rpix
)
+=
dSdCpix
[
1
*
lmax1
+
l
][
2
*
ns
+
2
]
*
bl
[
1
*
lmax1
+
l
]
*
bl
[
1
*
lmax1
+
l
];
//EE on UU
s
++
;
dSdC
(
s
*
nbin
+
ib
,
s
1
*
npix
+
cpix
,
s
1
*
npix
+
rpix
)
+=
dSdCpix
[
2
*
lmax1
+
l
][
1
*
ns
+
1
]
*
bl
[
2
*
lmax1
+
l
]
*
bl
[
2
*
lmax1
+
l
];
//BB on QQ
dSdC
(
s
*
nbin
+
ib
,
s
2
*
npix
+
cpix
,
s
1
*
npix
+
rpix
)
+=
dSdCpix
[
2
*
lmax1
+
l
][
2
*
ns
+
1
]
*
bl
[
2
*
lmax1
+
l
]
*
bl
[
2
*
lmax1
+
l
];
//BB on QU
dSdC
(
s
*
nbin
+
ib
,
s
1
*
npix
+
cpix
,
s
2
*
npix
+
rpix
)
+=
dSdCpix
[
2
*
lmax1
+
l
][
1
*
ns
+
2
]
*
bl
[
2
*
lmax1
+
l
]
*
bl
[
2
*
lmax1
+
l
];
//BB on UQ
dSdC
(
s
*
nbin
+
ib
,
s
2
*
npix
+
cpix
,
s
2
*
npix
+
rpix
)
+=
dSdCpix
[
2
*
lmax1
+
l
][
2
*
ns
+
2
]
*
bl
[
2
*
lmax1
+
l
]
*
bl
[
2
*
lmax1
+
l
];
//BB on UU
dSdC
(
s
*
nbin
+
ib
,
s
Q
*
npix
+
cpix
,
s
Q
*
npix
+
rpix
)
+=
dSdCpix
[
2
*
lmax1
+
l
][
1
*
ns
+
1
]
*
bl
[
2
*
lmax1
+
l
]
*
bl
[
2
*
lmax1
+
l
];
//BB on QQ
dSdC
(
s
*
nbin
+
ib
,
s
U
*
npix
+
cpix
,
s
Q
*
npix
+
rpix
)
+=
dSdCpix
[
2
*
lmax1
+
l
][
2
*
ns
+
1
]
*
bl
[
2
*
lmax1
+
l
]
*
bl
[
2
*
lmax1
+
l
];
//BB on QU
dSdC
(
s
*
nbin
+
ib
,
s
Q
*
npix
+
cpix
,
s
U
*
npix
+
rpix
)
+=
dSdCpix
[
2
*
lmax1
+
l
][
1
*
ns
+
2
]
*
bl
[
2
*
lmax1
+
l
]
*
bl
[
2
*
lmax1
+
l
];
//BB on UQ
dSdC
(
s
*
nbin
+
ib
,
s
U
*
npix
+
cpix
,
s
U
*
npix
+
rpix
)
+=
dSdCpix
[
2
*
lmax1
+
l
][
2
*
ns
+
2
]
*
bl
[
2
*
lmax1
+
l
]
*
bl
[
2
*
lmax1
+
l
];
//BB on UU
s
++
;
}
//TE
if
(
nstokes
==
3
)
{
dSdC
(
s
*
nbin
+
ib
,
1
*
npix
+
cpix
,
0
*
npix
+
rpix
)
+=
dSdCpix
[
3
*
lmax1
+
l
][
1
*
ns
+
0
]
*
bl
[
3
*
lmax1
+
l
]
*
bl
[
3
*
lmax1
+
l
];
//TE on IQ
dSdC
(
s
*
nbin
+
ib
,
2
*
npix
+
cpix
,
0
*
npix
+
rpix
)
+=
dSdCpix
[
3
*
lmax1
+
l
][
2
*
ns
+
0
]
*
bl
[
3
*
lmax1
+
l
]
*
bl
[
3
*
lmax1
+
l
];
//TE on IU
dSdC
(
s
*
nbin
+
ib
,
0
*
npix
+
cpix
,
1
*
npix
+
rpix
)
+=
dSdCpix
[
3
*
lmax1
+
l
][
0
*
ns
+
1
]
*
bl
[
3
*
lmax1
+
l
]
*
bl
[
3
*
lmax1
+
l
];
//TE on QI
dSdC
(
s
*
nbin
+
ib
,
0
*
npix
+
cpix
,
2
*
npix
+
rpix
)
+=
dSdCpix
[
3
*
lmax1
+
l
][
0
*
ns
+
2
]
*
bl
[
3
*
lmax1
+
l
]
*
bl
[
3
*
lmax1
+
l
];
//TE on UI
if
(
ispec
[
3
]
==
1
)
{
dSdC
(
s
*
nbin
+
ib
,
sI
*
npix
+
cpix
,
sQ
*
npix
+
rpix
)
+=
dSdCpix
[
3
*
lmax1
+
l
][
0
*
ns
+
1
]
*
bl
[
3
*
lmax1
+
l
]
*
bl
[
3
*
lmax1
+
l
];
//TE on IQ
dSdC
(
s
*
nbin
+
ib
,
sI
*
npix
+
cpix
,
sU
*
npix
+
rpix
)
+=
dSdCpix
[
3
*
lmax1
+
l
][
0
*
ns
+
2
]
*
bl
[
3
*
lmax1
+
l
]
*
bl
[
3
*
lmax1
+
l
];
//TE on IU
dSdC
(
s
*
nbin
+
ib
,
sQ
*
npix
+
cpix
,
sI
*
npix
+
rpix
)
+=
dSdCpix
[
3
*
lmax1
+
l
][
1
*
ns
+
0
]
*
bl
[
3
*
lmax1
+
l
]
*
bl
[
3
*
lmax1
+
l
];
//TE on QI
dSdC
(
s
*
nbin
+
ib
,
sU
*
npix
+
cpix
,
sI
*
npix
+
rpix
)
+=
dSdCpix
[
3
*
lmax1
+
l
][
2
*
ns
+
0
]
*
bl
[
3
*
lmax1
+
l
]
*
bl
[
3
*
lmax1
+
l
];
//TE on UI
s
++
;
}
//TB
if
(
ispec
[
4
]
==
1
)
{
dSdC
(
s
*
nbin
+
ib
,
sI
*
npix
+
cpix
,
sQ
*
npix
+
rpix
)
+=
dSdCpix
[
4
*
lmax1
+
l
][
0
*
ns
+
1
]
*
bl
[
3
*
lmax1
+
l
]
*
bl
[
3
*
lmax1
+
l
];
//TB on IQ
dSdC
(
s
*
nbin
+
ib
,
sI
*
npix
+
cpix
,
sU
*
npix
+
rpix
)
+=
dSdCpix
[
4
*
lmax1
+
l
][
0
*
ns
+
2
]
*
bl
[
3
*
lmax1
+
l
]
*
bl
[
3
*
lmax1
+
l
];
//TB on IU
dSdC
(
s
*
nbin
+
ib
,
sQ
*
npix
+
cpix
,
sI
*
npix
+
rpix
)
+=
dSdCpix
[
4
*
lmax1
+
l
][
1
*
ns
+
0
]
*
bl
[
3
*
lmax1
+
l
]
*
bl
[
3
*
lmax1
+
l
];
//TB on QI
dSdC
(
s
*
nbin
+
ib
,
sU
*
npix
+
cpix
,
sI
*
npix
+
rpix
)
+=
dSdCpix
[
4
*
lmax1
+
l
][
2
*
ns
+
0
]
*
bl
[
3
*
lmax1
+
l
]
*
bl
[
3
*
lmax1
+
l
];
//TB on UI
s
++
;
}
//EB
if
(
ispec
[
5
]
==
1
)
{
dSdC
(
s
*
nbin
+
ib
,
sQ
*
npix
+
cpix
,
sQ
*
npix
+
rpix
)
+=
dSdCpix
[
5
*
lmax1
+
l
][
1
*
ns
+
1
]
*
bl
[
1
*
lmax1
+
l
]
*
bl
[
2
*
lmax1
+
l
];
//EB on QQ
dSdC
(
s
*
nbin
+
ib
,
sQ
*
npix
+
cpix
,
sU
*
npix
+
rpix
)
+=
dSdCpix
[
5
*
lmax1
+
l
][
1
*
ns
+
2
]
*
bl
[
1
*
lmax1
+
l
]
*
bl
[
2
*
lmax1
+
l
];
//EB on QU
dSdC
(
s
*
nbin
+
ib
,
sU
*
npix
+
cpix
,
sQ
*
npix
+
rpix
)
+=
dSdCpix
[
5
*
lmax1
+
l
][
2
*
ns
+
1
]
*
bl
[
1
*
lmax1
+
l
]
*
bl
[
2
*
lmax1
+
l
];
//EB on UQ
dSdC
(
s
*
nbin
+
ib
,
sU
*
npix
+
cpix
,
sU
*
npix
+
rpix
)
+=
dSdCpix
[
5
*
lmax1
+
l
][
2
*
ns
+
2
]
*
bl
[
1
*
lmax1
+
l
]
*
bl
[
2
*
lmax1
+
l
];
//EB on UU
s
++
;
}
}
/* end loop l */
}
/* end loop bins */
...
...
@@ -116,23 +139,24 @@ void build_dSdC( int nside, int nstokes, int npix, int nbin, long *ellbins, long
/* compute 3x3 matrix correlation matrix for couple (ipix,jpix) */
void
QML_compute_dSdC
(
double
*
vr
,
double
*
vc
,
int
lmax
,
int
nstokes
,
double
**
dSdCpix
)
void
QML_compute_dSdC
(
double
*
vr
,
double
*
vc
,
int
lmax
,
long
*
ispec
,
double
**
dSdCpix
)
{
double
*
pl
=
NULL
,
*
d20
=
NULL
,
*
d2p2
=
NULL
,
*
d2m2
=
NULL
;
double
cos2aij
,
sin2aij
,
cos2aji
,
sin2aji
;
double
norm
,
P02
,
Q22
,
R22
;
double
cos_chi
;
int
l
,
ns
=
3
;
int
spin0
=
0
,
spin2
=
0
,
spin02
=
0
;
/* alloc vectors */
if
(
nstokes
!=
2
)
pl
=
(
double
*
)
calloc
(
(
lmax
+
1
)
,
s
izeof
(
double
))
;
if
(
nstokes
>
1
)
{
d2p2
=
(
double
*
)
calloc
(
(
lmax
+
1
),
sizeof
(
double
));
d2m2
=
(
double
*
)
calloc
(
(
lmax
+
1
),
sizeof
(
double
));
}
if
(
nstokes
>
2
)
d20
=
(
double
*
)
calloc
(
(
lmax
+
1
),
sizeof
(
double
));
if
(
ispec
[
0
]
==
1
)
spin0
=
1
;
if
(
ispec
[
1
]
==
1
||
ispec
[
2
]
==
1
||
ispec
[
5
]
==
1
)
s
pin2
=
1
;
if
(
ispec
[
3
]
==
1
||
ispec
[
4
]
==
1
)
spin02
=
1
;
if
(
spin0
)
pl
=
(
double
*
)
calloc
(
(
lmax
+
1
),
sizeof
(
double
));
if
(
spin2
)
d2p2
=
(
double
*
)
calloc
(
(
lmax
+
1
),
sizeof
(
double
));
if
(
spin2
)
d2m2
=
(
double
*
)
calloc
(
(
lmax
+
1
),
sizeof
(
double
));
if
(
spin02
)
d20
=
(
double
*
)
calloc
(
(
lmax
+
1
),
sizeof
(
double
));
/* generate d_ss'(l) */
cos_chi
=
QML_scal_prod
(
vc
,
vr
);
...
...
@@ -141,23 +165,23 @@ void QML_compute_dSdC( double *vr, double *vc, int lmax, int nstokes, double **d
if
(
cos_chi
<
-
1
)
cos_chi
=
-
1
.;
/* legendre( cos_chi, 0, lmax, pl); */
if
(
nstokes
!=
2
)
dlss
(
cos_chi
,
0
,
0
,
lmax
,
pl
);
if
(
nstokes
>
1
)
dlss
(
cos_chi
,
2
,
2
,
lmax
,
d2p2
);
if
(
nstokes
>
1
)
dlss
(
cos_chi
,
2
,
-
2
,
lmax
,
d2m2
);
if
(
nstokes
>
2
)
dlss
(
cos_chi
,
2
,
0
,
lmax
,
d20
);
if
(
spin0
)
dlss
(
cos_chi
,
0
,
0
,
lmax
,
pl
);
if
(
spin2
)
dlss
(
cos_chi
,
2
,
2
,
lmax
,
d2p2
);
if
(
spin2
)
dlss
(
cos_chi
,
2
,
-
2
,
lmax
,
d2m2
);
if
(
spin0
2
)
dlss
(
cos_chi
,
2
,
0
,
lmax
,
d20
);
/* generate rotation angles */
if
(
nstokes
>
1
)
polrotangle
(
vr
,
vc
,
&
cos2aij
,
&
sin2aij
);
if
(
nstokes
>
1
)
polrotangle
(
vc
,
vr
,
&
cos2aji
,
&
sin2aji
);
if
(
spin2
||
spin02
)
polrotangle
(
vr
,
vc
,
&
cos2aij
,
&
sin2aij
);
if
(
spin2
||
spin02
)
polrotangle
(
vc
,
vr
,
&
cos2aji
,
&
sin2aji
);
/* loop on l */
for
(
l
=
2
;
l
<=
lmax
;
l
++
)
{
norm
=
(
double
)(
2
*
l
+
1
)
/
(
4
.
*
M_PI
);
if
(
nstokes
!
=
2
)
if
(
ispec
[
0
]
=
=
1
)
dSdCpix
[
0
*
(
lmax
+
1
)
+
l
][
0
*
ns
+
0
]
=
norm
*
pl
[
l
]
;
//TT on II
if
(
nstokes
>
1
)
{
if
(
ispec
[
1
]
==
1
||
ispec
[
2
]
==
1
)
{
Q22
=
norm
*
(
d2p2
[
l
]
+
d2m2
[
l
]
)
/
2
.;
R22
=
norm
*
(
d2p2
[
l
]
-
d2m2
[
l
]
)
/
2
.;
...
...
@@ -172,7 +196,7 @@ void QML_compute_dSdC( double *vr, double *vc, int lmax, int nstokes, double **d
dSdCpix
[
2
*
(
lmax
+
1
)
+
l
][
2
*
ns
+
2
]
=
(
sin2aij
*
sin2aji
*
R22
+
cos2aij
*
cos2aji
*
Q22
);
//BB on UU
}
if
(
nstokes
>
2
)
{
if
(
ispec
[
3
]
==
1
)
{
P02
=
-
d20
[
l
];
dSdCpix
[
3
*
(
lmax
+
1
)
+
l
][
1
*
ns
+
0
]
=
P02
*
cos2aji
;
//TE on IQ
dSdCpix
[
3
*
(
lmax
+
1
)
+
l
][
2
*
ns
+
0
]
=
-
P02
*
sin2aji
;
//TE on IU
...
...
@@ -180,12 +204,29 @@ void QML_compute_dSdC( double *vr, double *vc, int lmax, int nstokes, double **d
dSdCpix
[
3
*
(
lmax
+
1
)
+
l
][
0
*
ns
+
2
]
=
-
P02
*
sin2aij
;
//TE on UI
}
if
(
ispec
[
4
]
==
1
)
{
P02
=
-
d20
[
l
];
dSdCpix
[
4
*
(
lmax
+
1
)
+
l
][
1
*
ns
+
0
]
=
P02
*
sin2aji
;
//TB on IQ
dSdCpix
[
4
*
(
lmax
+
1
)
+
l
][
2
*
ns
+
0
]
=
P02
*
cos2aji
;
//TB on IU
dSdCpix
[
4
*
(
lmax
+
1
)
+
l
][
0
*
ns
+
1
]
=
P02
*
sin2aij
;
//TB on QI
dSdCpix
[
4
*
(
lmax
+
1
)
+
l
][
0
*
ns
+
2
]
=
P02
*
cos2aij
;
//TB on UI
}
if
(
ispec
[
5
]
==
1
)
{
Q22
=
norm
*
(
d2p2
[
l
]
+
d2m2
[
l
]
)
/
2
.;
R22
=
norm
*
(
d2p2
[
l
]
-
d2m2
[
l
]
)
/
2
.;
dSdCpix
[
5
*
(
lmax
+
1
)
+
l
][
1
*
ns
+
1
]
=
(
cos2aij
*
sin2aji
*
Q22
-
sin2aij
*
cos2aji
*
R22
+
sin2aij
*
cos2aji
*
Q22
-
cos2aij
*
sin2aji
*
R22
);
//EB on QQ
dSdCpix
[
5
*
(
lmax
+
1
)
+
l
][
1
*
ns
+
2
]
=
(
cos2aij
*
cos2aji
*
Q22
+
sin2aij
*
sin2aji
*
R22
-
sin2aij
*
sin2aji
*
Q22
-
cos2aij
*
cos2aji
*
R22
);
//EB on QU
dSdCpix
[
5
*
(
lmax
+
1
)
+
l
][
2
*
ns
+
1
]
=
(
cos2aij
*
cos2aji
*
Q22
+
sin2aij
*
sin2aji
*
R22
-
sin2aij
*
sin2aji
*
Q22
-
cos2aij
*
cos2aji
*
R22
);
//EB on UQ
dSdCpix
[
5
*
(
lmax
+
1
)
+
l
][
2
*
ns
+
2
]
=
(
-
cos2aij
*
sin2aji
*
Q22
+
sin2aij
*
cos2aji
*
R22
-
sin2aij
*
cos2aji
*
Q22
+
cos2aij
*
sin2aji
*
R22
);
//EB on UU
}
}
if
(
nstokes
!=
2
)
free
(
pl
);
if
(
nstokes
>
1
)
free
(
d2p2
);
if
(
nstokes
>
1
)
free
(
d2m2
);
if
(
nstokes
>
2
)
free
(
d20
);
if
(
spin0
)
free
(
pl
);
if
(
spin2
)
free
(
d2p2
);
if
(
spin2
)
free
(
d2m2
);
if
(
spin0
2
)
free
(
d20
);
}
...
...
@@ -338,6 +379,20 @@ double fact( int n)
/*************************************************************************/
int
ispec2nspec
(
long
*
ispec
)
{
int
nspec
=
0
;
for
(
int
i
=
0
;
i
<
6
;
i
++
)
if
(
ispec
[
i
])
nspec
++
;
//force EE and BB
if
(
ispec
[
1
]
==
1
)
ispec
[
2
]
==
1
;
if
(
ispec
[
2
]
==
1
)
ispec
[
1
]
==
1
;
return
(
nspec
);
}
int
nstokes2nspec
(
int
nstokes
)
{
int
nspec
=
0
;
...
...
src/libcov.h
View file @
ecb32a95
...
...
@@ -12,7 +12,7 @@
#define dSdC(L,I,J) dSdC[ (L)*(npixtot*npixtot) + (I)*npixtot + (J)]
#define halfpi M_PI/2.
void
build_dSdC
(
int
nside
,
int
nstokes
,
int
npix
,
int
inl
,
long
*
ellbins
,
long
*
ipix
,
double
*
bl
,
double
*
dSdC
);
void
build_dSdC
(
int
nside
,
int
nstokes
,
int
npix
,
int
inl
,
long
*
ispec
,
long
*
ellbins
,
long
*
ipix
,
double
*
bl
,
double
*
dSdC
);
/**
Compute the Pl = dS/dCl matrices.
...
...
@@ -39,7 +39,7 @@ void build_dSdC( int nside, int nstokes, int npix, int inl, long *ellbins, long
*/
void
QML_compute_dSdC
(
double
*
vr
,
double
*
vc
,
int
lmax
,
int
nstokes
,
double
**
dSdCpix
);
void
QML_compute_dSdC
(
double
*
vr
,
double
*
vc
,
int
lmax
,
long
*
ispec
,
double
**
dSdCpix
);
int
polrotangle
(
double
*
ri
,
double
*
rj
,
double
*
cos2a
,
double
*
sin2a
);
void
QML_cross_prod
(
double
*
v
,
double
*
w
,
double
*
z
);
double
QML_scal_prod
(
double
*
v
,
double
*
w
);
...
...
@@ -49,6 +49,7 @@ double fact( int n);
int
define_lrange
(
int
lmax
,
int
*
ils
);
int
nstokes2nspec
(
int
nstokes
);
int
ispec2nspec
(
long
*
ispec
);
void
pix2ang_ring
(
long
nside
,
long
ipix
,
double
*
theta
,
double
*
phi
);
void
pix2vec_ring
(
long
nside
,
long
ipix
,
double
*
vec
);
...
...
xqml/estimators.py
View file @
ecb32a95
...
...
@@ -168,10 +168,10 @@ def CrossWindowFunction(El, Pl):
[134 478 822]]
"""
nl
=
len
(
El
)
# pas de transpose car symm
Wll
=
np
.
asarray
(
[
np
.
sum
(
E
*
P
)
for
E
in
El
for
P
in
Pl
]
).
reshape
(
nl
,
nl
)
return
Wll
...
...
xqml/libcov.py
View file @
ecb32a95
...
...
@@ -81,34 +81,38 @@ def compute_ds_dcb( ellbins, nside, ipok, bl, clth, Slmax, spec,
start
=
timeit
.
default_timer
()
temp
=
"TT"
in
spec
clth
=
np
.
asarray
(
clth
)
if
len
(
clth
)
==
4
:
clth
=
np
.
concatenate
((
clth
,
clth
[
0
:
2
]
*
0.
))
temp
=
"TT"
in
spec
polar
=
"EE"
in
spec
or
"BB"
in
spec
corr
=
"TE"
in
spec
or
"TB"
in
spec
or
"EB"
in
spec
corr
=
"TE"
in
spec
or
"TB"
in
spec
or
"EB"
in
spec
if
Sonly
:
if
MC
:
S
=
S_bins_MC
(
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
MC
,
polar
=
polar
,
temp
=
temp
,
corr
=
corr
,
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
MC
,
spec
,
pixwin
=
pixwin
,
timing
=
timing
)
else
:
S
=
compute_S
(
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
polar
=
polar
,
temp
=
temp
,
corr
=
corr
,
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
spec
,
pixwin
=
pixwin
,
timing
=
timing
)
return
S
if
MC
:
Pl
,
S
=
covth_bins_MC
(
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
MC
,
polar
=
polar
,
temp
=
temp
,
corr
=
corr
,
pixwin
=
pixwin
,
timing
=
timing
)
spec
,
pixwin
=
pixwin
,
timing
=
timing
)
elif
openMP
:
fpixwin
=
extrapolpixwin
(
nside
,
Slmax
,
pixwin
)
bell
=
np
.
array
([
bl
*
fpixwin
]
*
4
)[:
Slmax
+
1
].
ravel
()
stokes
,
spec
,
istokes
,
ispecs
=
getstokes
(
spec
)
ispec
=
np
.
zeros
(
6
,
int
)
ispec
[
ispecs
]
=
1
nbins
=
(
len
(
ellbins
)
-
1
)
*
len
(
spec
)
npix
=
len
(
ipok
)
*
len
(
istokes
)
Pl
=
np
.
ndarray
(
nbins
*
npix
**
2
)
clibcov
.
dSdC
(
nside
,
len
(
istokes
),
ellbins
,
ipok
,
bell
,
Pl
)
clibcov
.
dSdC
(
nside
,
len
(
istokes
),
ispec
,
ellbins
,
ipok
,
bell
,
Pl
)
Pl
=
Pl
.
reshape
(
nbins
,
npix
,
npix
)
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
())
...
...
@@ -117,8 +121,9 @@ 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
timing
:
print
(
"Total time (npix=%d): %.1f sec"
%
(
len
(
ipok
),
timeit
.
default_timer
()
-
start
))
return
Pl
,
S
...
...
@@ -235,7 +240,7 @@ def compute_PlS(
for
i
in
np
.
arange
(
npi
):
if
timing
:
progress_bar
(
i
,
npi
,
-
0.5
*
(
start
-
timeit
.
default_timer
())
)
progress_bar
(
i
,
npi
)
for
j
in
np
.
arange
(
i
,
npi
):
cos_chi
=
allcosang
[
i
,
j
]
if
temp
:
...
...
@@ -482,7 +487,7 @@ def compute_S(
start
=
timeit
.
default_timer
()
for
i
in
np
.
arange
(
npi
):
if
timing
:
progress_bar
(
i
,
npi
,
-
0.25
*
(
start
-
timeit
.
default_timer
())
)
progress_bar
(
i
,
npi
)
for
j
in
np
.
arange
(
i
,
npi
):
if
temp
:
pl
=
norm
*
pl0
(
allcosang
[
i
,
j
],
Slmax
)[
2
:]
...
...
@@ -573,7 +578,7 @@ def compute_S(
def
covth_bins_MC
(
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
nsimu
,
polar
=
False
,
temp
=
False
,
corr
=
False
,
pixwin
=
False
,
timing
=
False
):
spec
,
pixwin
=
False
,
timing
=
False
):
"""
Can be particularly slow on sl7.
To be enhanced and extended to TT and correlations.
...
...
@@ -616,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,
...
polar=True, temp=False, corr=False
, pixwin=True, timing=False)
...
spec
, pixwin=True, timing=False)
>>> print(round(np.sum(Pl),5), round(np.sum(S),5))
(12.68135, 406990.12056)
...
...
@@ -625,10 +630,12 @@ 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,
...
polar=False, temp=True, corr=False
, pixwin=True, timing=False)
...
spec
, pixwin=True, timing=False)
>>> print(round(np.sum(Pl),5), round(np.sum(S),5))
(42.35592, 7414.01784)
"""
polar
=
'EE'
in
spec
or
'BB'
in
spec
if
nsimu
==
1
:
nsimu
=
(
12
*
nside
**
2
)
*
10
*
(
int
(
polar
)
+
1
)
...
...
@@ -639,7 +646,7 @@ def covth_bins_MC(
maxell
=
np
.
array
(
ellbins
[
1
:
nbins
+
1
])
-
1
ellval
=
(
minell
+
maxell
)
*
0.5
Stokes
,
spec
,
istokes
,
ispecs
=
getstokes
(
polar
=
polar
,
temp
=
temp
,
corr
=
corr
)
Stokes
,
spec
,
istokes
,
ispecs
=
getstokes
(
spec
)
nspec
=
len
(
spec
)
ll
=
np
.
arange
(
Slmax
+
1
)
fpixwin
=
extrapolpixwin
(
nside
,
Slmax
,
pixwin
)
...
...
@@ -667,8 +674,7 @@ def covth_bins_MC(
start
=
timeit
.
default_timer
()
for
l
in
np
.
arange
((
nspec
*
nbins
)):
if
timing
:
progress_bar
(
l
,
nspec
*
(
nbins
),
-
(
start
-
timeit
.
default_timer
()))
progress_bar
(
l
,
nspec
*
(
nbins
))
data
=
[
np
.
array
(
...
...
@@ -699,8 +705,7 @@ def covth_bins_MC(
Pl
=
np
.
zeros
(((
nbins
),
npix
,
npix
))
for
l
in
np
.
arange
((
nbins
)):
if
timing
:
progress_bar
(
l
,
nspec
*
(
nbins
),
-
(
start
-
timeit
.
default_timer
()))
progress_bar
(
l
,
nspec
*
(
nbins
))
Pl
[
l
]
=
np
.
cov
(
np
.
array
([
hp
.
synfast
(
...
...
@@ -730,7 +735,7 @@ def covth_bins_MC(
def
S_bins_MC
(
ellbins
,
nside
,
ipok
,
allcosang
,
bl
,
clth
,
Slmax
,
nsimu
,
polar
=
True
,
temp
=
True
,
corr
=
False
,
pixwin
=
False
,
timing
=
False
):
spec
,
pixwin
=
False
,
timing
=
False
):
"""
Can be particularly slow on sl7 !
To be enhanced and extended to TT and correlations
...
...
@@ -768,6 +773,7 @@ def S_bins_MC(
Pixel signal covariance matrix S
"""
polar
=
'EE'
in
spec
or
'BB'
in
spec
if
nsimu
==
1
:
nsimu
=
(
12
*
nside
**
2
)
*
10
*
(
int
(
polar
)
+
1
)
lmax
=
ellbins
[
-
1
]
...
...
@@ -776,7 +782,7 @@ def S_bins_MC(
minell
=
np
.
array
(
ellbins
[
0
:
nbins
])
maxell
=
np
.
array
(
ellbins
[
1
:
nbins
+
1
])
-
1
ellval
=
(
minell
+
maxell
)
*
0.5
Stokes
,
spec
,
ind
=
getstokes
(
polar
=
polar
,
temp
=
temp
,
corr
=
corr
)
Stokes
,
spec
,
ind
=
getstokes
(
spec
)
nspec
=
len
(
spec
)
ll
=
np
.
arange
(
Slmax
+
2
)
fpixwin
=
extrapolpixwin
(
nside
,
Slmax
,
pixwin
)
...
...
xqml/xqml.py
View file @
ecb32a95
...
...
@@ -34,8 +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
=
None
,
temp
=
False
,
polar
=
True
,
corr
=
False
,
pixwin
=
True
):
S
=
None
,
fwhm
=
0.
,
bell
=
None
,
spec
=
[
'EE'
,
'BB'
],
pixwin
=
True
):
"""
Parameters
----------
...
...
@@ -53,15 +52,8 @@ class xQML(object):
FWHM of the experiment beam
bell : ndarray, optional
beam transfer function (priority over fwhm)
polar : boolean, optional
Compute the polarisation part E and B. Default: True
temp : boolean, optional
Compute the temperature part T. Default: False
corr : boolean, optional
If True, compute TE, TB, EB spectra. Default: False
pixwin : boolean, optional
If True, applies pixel window function to spectra. Default: True
"""
self
.
bias
=
None
self
.
cross
=
NB
is
not
None
...
...
@@ -94,11 +86,15 @@ class xQML(object):
# For example that would be good to assert that the user
# set at least polar or temp to True.
self
.
stokes
,
self
.
spec
,
self
.
istokes
,
self
.
ispecs
=
getstokes
(
spec
,
temp
,
polar
,
corr
)
self
.
stokes
,
self
.
spec
,
self
.
istokes
,
self
.
ispecs
=
getstokes
(
spec
)
self
.
nstokes
=
len
(
self
.
stokes
)
self
.
nspec
=
len
(
self
.
spec
)
self
.
pixwin
=
pixwin
clth
=
np
.
asarray
(
clth
)
if
len
(
clth
)
==
4
:
clth
=
np
.
concatenate
((
clth
,
clth
[
0
:
2
]
*
0.
))
# If Pl is given by the user, just load it, and then compute the signal
# covariance using the fiducial model.
# Otherwise compute Pl and S from the arguments.
...
...
@@ -146,19 +142,19 @@ class xQML(object):
self
.
cross
=
NB
is
not
None
self
.
NA
=
NA
self
.
NB
=
NB
if
self
.
cross
else
NA
# 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...
self
.
El
=
El
(
invCa
,
invCb
,
self
.
Pl
)
# Finally compute invW by inverting...
self
.
invW
=
linalg
.
inv
(
CrossWindowFunction
(
self
.
El
,
self
.
Pl
))
# Compute bias for auto
# if not self.cross:
# self.bias = biasQuadEstimator(self.NA, self.El)
...
...
xqml/xqml_utils.py
View file @
ecb32a95
...
...
@@ -18,7 +18,7 @@ def pd_inv(a):
return
linalg
.
solve
(
a
,
I
,
sym_pos
=
True
,
overwrite_b
=
True
)
def
getstokes
(
spec
=
None
,
temp
=
False
,
polar
=
False
,
corr
=
False
):
def
getstokes
(
spec
):
"""
Get the Stokes parameters number and name(s)
...
...
@@ -26,12 +26,6 @@ def getstokes(spec=None, temp=False, polar=False, corr=False):
----------
spec : bool
If True, get Stokes parameters for polar (default: True)
polar : bool
If True, get Stokes parameters for polar (default: True)
temp : bool
If True, get Stokes parameters for temperature (default: False)
corr : bool
If True, get Stokes parameters for EB and TB (default: False)
Returns
----------
...
...
@@ -44,50 +38,28 @@ def getstokes(spec=None, temp=False, polar=False, corr=False):
Example
----------
>>> getstokes(
polar=True, temp=False, corr=False
)
>>> getstokes(
['EE','BB']
)
(['Q', 'U'], ['EE', 'BB'], [1, 2])
>>> getstokes(
polar=True, temp=True, corr=False
)
>>> getstokes(
['TT','EE','BB','TE']
)
(['I', 'Q', 'U'], ['TT', 'EE', 'BB', 'TE'], [0, 1, 2, 3])