Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/master' into depth1_mapmaker
Browse files Browse the repository at this point in the history
  • Loading branch information
chervias committed Jan 19, 2025
2 parents f8387db + 2347d2d commit 703a3b4
Show file tree
Hide file tree
Showing 60 changed files with 5,644 additions and 1,852 deletions.
6 changes: 2 additions & 4 deletions .github/workflows/test.yml
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,6 @@ jobs:
fail-fast: false
matrix:
include:
- python: "3.8"
pyshort: "38"
- python: "3.9"
pyshort: "39"
- python: "3.11"
Expand Down Expand Up @@ -62,13 +60,13 @@ jobs:
export OMP_NUM_THREADS=2
export OPENBLAS_NUM_THREADS=2
export MPI_DISABLE=1
python3 setup.py test
python3 -m unittest discover
- name: Run MPI Tests
run: |
export OMP_NUM_THREADS=1
export OPENBLAS_NUM_THREADS=1
mpirun -np 2 python3 setup.py test
mpirun -n 2 python3 -m unittest discover
# FIXME: Re-enable after testing this procedure on a local
# apple machine.
Expand Down
3 changes: 3 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -105,3 +105,6 @@ venv.bak/

# pychram
.idea

# vscode
.vscode
11 changes: 11 additions & 0 deletions docs/axisman.rst
Original file line number Diff line number Diff line change
Expand Up @@ -144,6 +144,17 @@ The output of the ``wrap`` cal should be::
Note the boresight entry is marked with a ``*``, indicating that it's
an AxisManager rather than a numpy array.

Data access under an AxisManager is done based on field names. For example::

>>> print(dset.boresight.az)
[0. 0. 0. ... 0. 0. 0.]

Advanced data access is possible by a path like syntax. This is especially useful when
data access is dynamic and the field name is not known in advance. For example::

>>> print(dset["boresight.az"])
[0. 0. 0. ... 0. 0. 0.]

To slice this object, use the restrict() method. First, let's
restrict in the 'dets' axis. Since it's an Axis of type LabelAxis,
the restriction selector must be a list of strings::
Expand Down
1 change: 1 addition & 0 deletions docs/preprocess.rst
Original file line number Diff line number Diff line change
Expand Up @@ -252,6 +252,7 @@ Flagging and Products
.. autoclass:: sotodlib.preprocess.processes.FlagTurnarounds
.. autoclass:: sotodlib.preprocess.processes.DarkDets
.. autoclass:: sotodlib.preprocess.processes.SourceFlags
.. autoclass:: sotodlib.preprocess.processes.GetStats

HWP Related
:::::::::::
Expand Down
13 changes: 5 additions & 8 deletions docs/site_pipeline.rst
Original file line number Diff line number Diff line change
Expand Up @@ -940,16 +940,13 @@ and binned. Every atomic map consist of a ``weights``, ``wmap`` (weighted map),
and ``hits`` map, as well as an information file that is used for adding the map
to an atomic map database.

Command line arguments
``````````````````````
Configuration yaml file
````````````````````````

.. argparse::
:module: sotodlib.site_pipeline.make_atomic_filterbin_map
:func: get_parser
:prog: make-atomic-filterbin-map
The mapmaker is configured by supplying a yaml file with ``--config_file``.

Config file format
``````````````````
.. autoclass:: sotodlib.site_pipeline.make_atomic_filterbin_map.Cfg
:members:

The only mandatory parameters are ``context`` for a context file and ``preprocess_config``,
a preprocess database configuration file that will tell the script how to process the
Expand Down
11 changes: 5 additions & 6 deletions setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,7 @@

import unittest

from setuptools import find_packages, setup, Extension
from setuptools.command.test import test as TestCommand
from setuptools import find_packages, setup, Extension, Command

import versioneer

Expand Down Expand Up @@ -43,7 +42,7 @@
setup_opts["url"] = "https://github.com/simonsobs/sotodlib"
setup_opts["packages"] = find_packages(where=".", exclude="tests")
setup_opts["license"] = "MIT"
setup_opts["python_requires"] = ">=3.8.0"
setup_opts["python_requires"] = ">=3.9.0"
setup_opts["package_data"] = {
"sotodlib": [
"toast/ops/data/*"
Expand Down Expand Up @@ -84,16 +83,16 @@

# Class to run unit tests

class SOTestCommand(TestCommand):
class SOTestCommand(Command):

def __init__(self, *args, **kwargs):
super(SOTestCommand, self).__init__(*args, **kwargs)

def initialize_options(self):
TestCommand.initialize_options(self)
Command.initialize_options(self)

def finalize_options(self):
TestCommand.finalize_options(self)
Command.finalize_options(self)
self.test_suite = True

def mpi_world(self):
Expand Down
10 changes: 7 additions & 3 deletions sotodlib/coords/demod.py
Original file line number Diff line number Diff line change
Expand Up @@ -182,20 +182,24 @@ def from_map(tod, signal_map, cuts=None, flip_gamma=True, wrap=False, modulated=
return signal_sim


def rotate_demodQU(tod, update_focal_plane=True):
def rotate_demodQU(tod, sign=1, offset=0, update_focal_plane=True):
"""
Apply detectors' polarization angle calibration to the HWP demodulated Q and U timestreams to
place all detectors' Q and U timestreams in a common telescope frame. This updates tod.demodQ
and tod.demodU, in place.
Args:
tod : an axisManager object
update_focal_plane (bool, optional): Whether to set focal_plane.gamma angles to zero,
update_focal_plane (bool, optional): Whether to set focal_plane.gamma angles to zero,
consistent with new coordinate reference. Make this true for polarization mapmaking
using make_map.
offset : float, optional
The rotation angle in degrees to apply (default is 0).
sign : int, optional
A sign factor to control the direction of the rotation (default is +1).
"""
demodC = ((tod.demodQ + 1j*tod.demodU).T * np.exp(-2j*tod.focal_plane.gamma)).T
demodC = ((tod.demodQ + 1j*tod.demodU).T * np.exp( sign*(-2j*tod.focal_plane.gamma + 1j*np.deg2rad(offset)) )).T
tod.demodQ = demodC.real
tod.demodU = demodC.imag
del demodC
Expand Down
77 changes: 53 additions & 24 deletions sotodlib/coords/helpers.py
Original file line number Diff line number Diff line change
Expand Up @@ -353,7 +353,7 @@ def get_footprint(tod, wcs_kernel, dets=None, timestamps=None, boresight=None,


def get_focal_plane_cover(tod=None, count=0, focal_plane=None,
xieta=None):
xieta=None, det_weights=None):
"""Process a bunch of detector positions into a center and radius such
that a circle with that center and radius contains all the
detectors. Also return detector positions, arranged approximately
Expand All @@ -366,6 +366,9 @@ def get_focal_plane_cover(tod=None, count=0, focal_plane=None,
not passed in
xieta (array): (2, n) array (or similar) of xi and eta
detector positions.
det_weights (array): If provided, must be same length as the xi
and eta vectors. Only dets with non-zero value for det_weights
will be included in the evaluation of the cover.
Returns:
xieta0: array[2] with array center, (xi0, eta0).
Expand All @@ -374,9 +377,22 @@ def get_focal_plane_cover(tod=None, count=0, focal_plane=None,
coords of the circular convex hull.
Notes:
If count=0, an empty list is returned for xietas. Otherwise,
If count=0, an empty list is returned for xietas. Otherwise,
count must be at least 3 so that the shape is not degenerate.
Any xi, eta that are not finite (e.g. nan or inf) are excluded
from the computation. If no detectors remain after the combined
finiteness and det_weights cuts, a ValueError is raised.
Note that ``det_weights`` can be int, float, or bool
type. Sometimes you might want to only include optical dets in
the result; e.g.::
..., det_weights=(aman.det_info.wafer.type=="OPTC"), ...
In degenerate cases (all dets are in exactly the same place), a
radius of zero may be returned.
"""
if xieta is None:
if focal_plane is None:
Expand All @@ -385,6 +401,17 @@ def get_focal_plane_cover(tod=None, count=0, focal_plane=None,
eta = focal_plane.eta
else:
xi, eta = xieta[:2]
mask = np.isfinite(xi) * np.isfinite(eta)
if det_weights is not None:
mask *= det_weights.astype(bool)

if not np.any(mask):
raise ValueError('All provided (xi, eta) coords are excluded; '
'cannot estimate a focal plane cover.')

# Restrict to only dets under consideration.
xi, eta = xi[mask], eta[mask]

qs = so3g.proj.quat.rotation_xieta(xi, eta)

# Starting guess for center
Expand Down Expand Up @@ -477,23 +504,23 @@ def as_geom(item):
s0 = corner_b - corner_a
return tuple(map(int, s0)), w0

def _confirm_wcs(*maps):
"""Insist that all arguments either have the same .wcs, or do not have
a wcs. Each argument should be either an ndmap (with a .wcs
attribute) or an ndarray (without a .wcs attribute).
def _confirm_wcs(*wcss):
"""Insist that all arguments are either the same wcs, or None.
Raises a ValueError if more than one argument has a .wcs attribute
and they do not all agree. Returns either the first .wcs
attribute encountered, or None if there aren't any.
Raises a ValueError if more than one argument is a wcs, and not
all wcs agree. Returns either the first valid wcs, or None if
there aren't any.
"""
wcs_to_use = None
for i, m in enumerate(maps):
if hasattr(m, 'wcs'):
if wcs_to_use is None:
wcs_to_use = m.wcs
elif not wcsutils.equal(wcs_to_use, m.wcs):
raise ValueError('The wcs from %ith item is discordant with prior ones.' % i)
for i, wcs in enumerate(wcss):
if wcs is None:
continue
if wcs_to_use is None:
wcs_to_use = wcs
elif not wcsutils.equal(wcs_to_use, wcs):
raise ValueError(
f'The wcs from {i}th item ({wcs}) is discordant with prior ones ({wcs_to_use})')
return wcs_to_use

def _invert_weights_map(weights, eigentol=1e-6, kill_partials=True,
Expand Down Expand Up @@ -565,23 +592,25 @@ def _invert_weights_map(weights, eigentol=1e-6, kill_partials=True,
# Reshape the output to match what was passed in.
return iw.transpose(1,2,0).reshape(weights.shape)

def _apply_inverse_weights_map(inverse_weights, target):
def _apply_inverse_weights_map(inverse_weights, target, out=None):
"""Apply a map of matrices to a map of vectors.
Assumes inverse_weights.shape = (a, b, ...) and target.shape =
(b, ...); the result has shape (a, ...).
"""
# master had:
#iw = inverse_weights.transpose((2,3,0,1))
#m = target.transpose((1,2,0)).reshape(
# target.shape[1], target.shape[2], target.shape[0], 1)
#m1 = np.matmul(iw, m)
#return m1.transpose(2,3,0,1).reshape(target.shape)
if out is None:
out = np.empty(inverse_weights.shape[1:],
dtype=target.dtype)
if isinstance(target, enmap.ndmap):
out = enmap.ndmap(out, target.wcs)
# Recall matmul(a, b) operates on the last two axes of (a, b). So
# move axes, and create a second one in target; re-order at end.
iw = np.moveaxis(inverse_weights, (0,1), (-2,-1))
t = np.moveaxis(target[:,None], (0,1), (-2,-1))
m = np.matmul(iw, t)
return np.moveaxis(m, (-2,-1), (0,1))[:,0]
out_moved = np.moveaxis(out[:,None], (0,1), (-2,-1))
np.matmul(iw, t, out=out_moved)
return out

class ScalarLastQuat(np.ndarray):
"""Wrapper class for numpy arrays carrying quaternions with the ijk1
Expand Down
Loading

0 comments on commit 703a3b4

Please sign in to comment.