Skip to content

Commit

Permalink
Merge branch 'master' into gp/fix/act_flags
Browse files Browse the repository at this point in the history
  • Loading branch information
@iparask
iparask committed Dec 12, 2024
2 parents d4d6659 + e334326 commit 98449b7
Show file tree
Hide file tree
Showing 26 changed files with 2,350 additions and 968 deletions.
7 changes: 3 additions & 4 deletions docs/preprocess.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,14 +51,12 @@ configuration files and specific manifest databases.

.. autofunction:: sotodlib.site_pipeline.preprocess_tod.preprocess_tod

.. autofunction:: sotodlib.site_pipeline.preprocess_tod.load_preprocess_det_select
.. autofunction:: sotodlib.preprocess.preprocess_util.load_preprocess_det_select

.. autofunction:: sotodlib.site_pipeline.preprocess_tod.load_preprocess_tod
.. autofunction:: sotodlib.preprocess.preprocess_util.load_and_preprocess

.. autofunction:: sotodlib.site_pipeline.preprocess_obs.preprocess_obs

.. autofunction:: sotodlib.site_pipeline.preprocess_obs.load_preprocess_obs

Example TOD Pipeline Configuration File
---------------------------------------

Expand Down Expand Up @@ -254,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
50 changes: 48 additions & 2 deletions docs/site_pipeline.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -850,8 +850,6 @@ Command line arguments
:module: sotodlib.site_pipeline.make_ml_map
:func: get_parser



Default Mapmaker Values
```````````````````````
The following code block contains the hard-coded default values for non-
Expand Down Expand Up @@ -933,6 +931,54 @@ Example of a config file:
verbose: True
quiet: False
make-atomic-filterbin-map
-------------------------

This script will create atomic maps (maps of individual observations by wafer and
frequency, and associated splits). These maps are HWP-demodulated and filtered
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.

Configuration yaml file
````````````````````````

The mapmaker is configured by supplying a yaml file with ``--config_file``.

.. 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
timestreams. A typical configuration file could look like this:

.. code-block:: yaml
context: /global/cfs/projectdirs/sobs/metadata/satp1/contexts/use_this_local.yaml
# Use a pixell area file for rectangular pixel maps or use an nside value for Healpix maps.
# Only use one of these options
area: band_car_fejer1_5arcmin.fits
#nside: 512
# A query can be a file with a list of obs, or an obsdb query
query: obs_list.txt
#query: "subtype == 'cmb' and timestamp >= 1708743600 and timestamp < 1713672000"
odir: output_directory
preprocess_config: preprocess_config.yaml
# Limit the number of obs, map a specific wafer or band
#ntod: 3
#wafer: ws0
#freq: f090
# Plataform to map
site: so_sat1
# Path to housekeeping data (this is used for extracting pwv)
hk_data_path: /global/cfs/cdirs/sobs/data/site/hk/
QDS Monitor
===========
Expand Down
114 changes: 71 additions & 43 deletions sotodlib/coords/helpers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -321,38 +321,39 @@ def get_footprint(tod, wcs_kernel, dets=None, timestamps=None, boresight=None,
fp1 = so3g.proj.FocalPlane.from_xieta(fake_dets, xieta1[0], xieta1[1])

asm = so3g.proj.Assembly.attach(sight, fp1)
output = np.zeros((len(fake_dets), n_samp, 4))
planar = np.zeros((len(fake_dets), n_samp, 4))
proj = so3g.proj.Projectionist.for_geom((1,1), wcs_kernel)
if rot:
# Works whether rot is a quat or a vector of them.
asm.Q = rot * asm.Q
proj.get_planar(asm, output=output)

output2 = output*0
proj.get_coords(asm, output=output2)
proj.get_planar(asm, output=planar)

# Get the pixel extrema in the form [{xmin,ymin},{xmax,ymax}]
delts = wcs_kernel.wcs.cdelt * DEG
planar = output[:,:,:2]
ranges = utils.minmax(planar/delts,(0,1))
# These are in units of pixel *offsets* from crval. crval
# might not correspond to a pixel center, though. So the
# thing that should be integer-valued to preserve pixel compatibility
# is crpix + ranges, not just ranges. Let's add crpix to transform this
# into offsets from the bottom-left pixel to make it easier to reason
# about integers
ranges += wcs_kernel.wcs.crpix
del output

# Start a new WCS and set the lower left corner.
ranges = utils.minmax(planar[:,:,:2]/delts,(0,1))
del planar

# These planar ranges are in units of pixels away from the
# reference point. The reference point is not necessarily on a
# pixel center -- that depends on whether crpix is an integer. So
# transform ranges by +(crpix - 1), making them relative to the
# bottom left pixel of wcs_kernel. Note the -1 here accounts for
# FITS numbering of pixels starting at (1, 1).
ranges += (wcs_kernel.wcs.crpix - 1)

# Round ranges, which in pixel offset units, to nearest integer.
corners = utils.nint(ranges)

# Start a new WCS. Adjust crpix to put our footprint in the
# bottom-left pixel.
w = wcs_kernel.deepcopy()
corners = utils.nint(ranges)
w.wcs.crpix -= corners[0]
shape = tuple(corners[1]-corners[0]+1)[::-1]
shape = tuple(corners[1] - corners[0] + 1)[::-1]
return (shape, w)


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 @@ -365,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 @@ -373,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 @@ -384,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 @@ -476,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 @@ -564,23 +592,23 @@ 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)
# 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 98449b7

Please sign in to comment.