Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Support passing None for HWP angles to pointing matrix calculation. #345

Merged
merged 2 commits into from
May 6, 2020
Merged

Support passing None for HWP angles to pointing matrix calculation. #345

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9268672
Support passing None for HWP angles to pointing matrix calculation.
tskisner May 6, 2020
1ed3503
Add check for hwp buffer size
tskisner May 6, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 5 additions & 1 deletion src/libtoast/src/toast_tod_pointing.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,11 @@ void toast::pointing_matrix_healpix(toast::HealpixPixels const & hpix,
// FIXME: Switch back to fast version after unit tests improved.
toast::vatan2(n, by, bx, detang.data());

if (hwpang != NULL) {
if (hwpang == NULL) {
for (size_t i = 0; i < n; ++i) {
detang[i] *= 2.0;
}
} else {
for (size_t i = 0; i < n; ++i) {
detang[i] += 2.0 * hwpang[i];
detang[i] *= 2.0;
Expand Down
25 changes: 18 additions & 7 deletions src/toast/_libtoast_tod_pointing.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,15 +9,13 @@
void init_tod_pointing(py::module & m) {
m.def("pointing_matrix_healpix",
[](toast::HealpixPixels const & hpix, bool nest, double eps, double cal,
std::string const & mode, py::buffer pdata, py::buffer hwpang,
std::string const & mode, py::buffer pdata, py::object hwpang,
py::buffer flags, py::buffer pixels, py::buffer weights) {
pybuffer_check_1D <double> (pdata);
pybuffer_check_1D <double> (hwpang);
pybuffer_check_1D <uint8_t> (flags);
pybuffer_check_1D <double> (weights);
pybuffer_check_1D <int64_t> (pixels);
py::buffer_info info_pdata = pdata.request();
py::buffer_info info_hwpang = hwpang.request();
py::buffer_info info_flags = flags.request();
py::buffer_info info_pixels = pixels.request();
py::buffer_info info_weights = weights.request();
Expand All @@ -26,8 +24,7 @@ void init_tod_pointing(py::module & m) {
if (mode.compare("IQU") == 0) {
nw = (size_t)(info_weights.size / 3);
}
if ((info_hwpang.size != n) ||
(info_flags.size != n) ||
if ((info_flags.size != n) ||
(info_pixels.size != n) || (nw != n)) {
auto log = toast::Logger::get();
std::ostringstream o;
Expand All @@ -36,16 +33,30 @@ void init_tod_pointing(py::module & m) {
throw std::runtime_error(o.str().c_str());
}
double * rawpdata = reinterpret_cast <double *> (info_pdata.ptr);
double * rawhwpang = reinterpret_cast <double *> (info_hwpang.ptr);
uint8_t * rawflags = reinterpret_cast <uint8_t *> (info_flags.ptr);
double * rawweights = reinterpret_cast <double *> (info_weights.ptr);
int64_t * rawpixels = reinterpret_cast <int64_t *> (info_pixels.ptr);
double * rawhwpang = NULL;
if (!hwpang.is_none()) {
auto hwpbuf = py::cast <py::buffer> (hwpang);
pybuffer_check_1D <double> (hwpbuf);
py::buffer_info info_hwpang = hwpbuf.request();
if (info_hwpang.size != n) {
auto log = toast::Logger::get();
std::ostringstream o;
o << "HWP buffer size is not consistent.";
log.error(o.str().c_str());
throw std::runtime_error(o.str().c_str());
}
rawhwpang = reinterpret_cast <double *> (info_hwpang.ptr);
tskisner marked this conversation as resolved.
Show resolved Hide resolved
}
toast::pointing_matrix_healpix(hpix, nest, eps, cal, mode, n,
rawpdata, rawhwpang, rawflags, rawpixels,
rawweights);
return;
}, py::arg("hpix"), py::arg("nest"), py::arg("eps"), py::arg("cal"),
py::arg("mode"), py::arg("pdata"), py::arg("hwpang"), py::arg("flags"),
py::arg("mode"), py::arg("pdata"), py::arg("hwpang").none(true),
py::arg("flags"),
py::arg("pixels"), py::arg(
"weights"), R"(
Compute the healpix pixel indices and weights for one detector.
Expand Down
89 changes: 89 additions & 0 deletions src/toast/tests/ops_pmat.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -165,6 +165,95 @@ def test_pointing_matrix_healpix(self):
self.assertFalse(failed)
return

def test_pointing_matrix_healpix_hwp(self):
nside = 64
hpix = HealpixPixels(64)
nest = True
psivec = np.radians([-180, -135, -90, -45, 0, 45, 90, 135, 180])
nsamp = len(psivec)
eps = 0.0
cal = 1.0
mode = "IQU"
nnz = 3
flags = np.zeros(nsamp, dtype=np.uint8)
pix = 49103
theta, phi = hp.pix2ang(nside, pix, nest=nest)
xaxis, yaxis, zaxis = np.eye(3)
thetarot = qa.rotation(yaxis, theta)
phirot = qa.rotation(zaxis, phi)
pixrot = qa.mult(phirot, thetarot)
quats = []
for psi in psivec:
psirot = qa.rotation(zaxis, psi)
quats.append(qa.mult(pixrot, psirot))
quats = np.vstack(quats)

# First with HWP angle == 0.0
hwpang = np.zeros(nsamp)
pixels_zero = np.zeros(nsamp, dtype=np.int64)
weights_zero = np.zeros([nsamp, nnz], dtype=np.float64)
pointing_matrix_healpix(
hpix,
nest,
eps,
cal,
mode,
quats.reshape(-1),
hwpang,
flags,
pixels_zero,
weights_zero.reshape(-1),
)

# Now passing hwpang == None
pixels_none = np.zeros(nsamp, dtype=np.int64)
weights_none = np.zeros([nsamp, nnz], dtype=np.float64)
pointing_matrix_healpix(
hpix,
nest,
eps,
cal,
mode,
quats.reshape(-1),
None,
flags,
pixels_none,
weights_none.reshape(-1),
)
# print("")
# for i in range(nsamp):
# print(
# "HWP zero: {} {} | {} {} {}".format(
# psivec[i],
# pixels_zero[i],
# weights_zero[i][0],
# weights_zero[i][1],
# weights_zero[i][2],
# )
# )
# print(
# " none: {} {} | {} {} {}".format(
# psivec[i],
# pixels_none[i],
# weights_none[i][0],
# weights_none[i][1],
# weights_none[i][2],
# )
# )
failed = False
if not np.all(np.equal(pixels_zero, pixels_none)):
print("HWP pixels do not agree {} != {}".format(pixels_zero, pixels_none))
failed = True

if not np.allclose(weights_zero, weights_none):
print(
"HWP weights do not agree {} != {}".format(weights_zero, weights_none)
)
failed = True

self.assertFalse(failed)
return

def test_hpix_simple(self):
rank = 0
if self.comm is not None:
Expand Down
10 changes: 5 additions & 5 deletions src/toast/todmap/pointing.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ class OpPointingHpix(Operator):
the total response is:

.. math::
d = cal \\left[\\frac{(1+eps)}{2} I + \\frac{(1-eps)}{2} \\left[Q \\cos{4(a+w)} + U \\sin{4(a+w)}\\right]\\right]
d = cal \\left[\\frac{(1+eps)}{2} I + \\frac{(1-eps)}{2} \\left[Q \\cos{2a+4w} + U \\sin{2a+4w}\\right]\\right]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@giuspugl Is this definition of the pointing weights consistent with what Davide used for the LB simulations?


Args:
pixels (str): write pixels to the cache with name <pixels>_<detector>.
Expand Down Expand Up @@ -164,9 +164,7 @@ def exec(self, data):
try:
hwpang = tod.local_hwp_angle()
except:
pass
if hwpang is None:
hwpang = np.zeros(nsamp, dtype=np.float64)
hwpang = None

# read the common flags and apply bitmask

Expand Down Expand Up @@ -227,7 +225,9 @@ def exec(self, data):
# Use cached version
detp = pdata[bslice, :]

hslice = hwpang[bslice].reshape(-1)
hslice = None
if hwpang is not None:
hslice = hwpang[bslice].reshape(-1)
fslice = common[bslice].reshape(-1)

pointing_matrix_healpix(
Expand Down