From 6a10d5bae861ff12a7671c9acbb9adae1f44a81a Mon Sep 17 00:00:00 2001 From: Tomoki Terasaki Date: Tue, 19 Mar 2024 02:05:11 +0000 Subject: [PATCH 01/27] initial commit --- sotodlib/coords/mapbased_pointing.py | 545 ++++++++++++++++++ .../site_pipeline/combine_focal_planes.py | 1 + .../site_pipeline/make_mapbased_pointing.py | 102 ++++ sotodlib/site_pipeline/update_pointing.py | 196 +++++++ 4 files changed, 844 insertions(+) create mode 100644 sotodlib/coords/mapbased_pointing.py create mode 100644 sotodlib/site_pipeline/combine_focal_planes.py create mode 100644 sotodlib/site_pipeline/make_mapbased_pointing.py create mode 100644 sotodlib/site_pipeline/update_pointing.py diff --git a/sotodlib/coords/mapbased_pointing.py b/sotodlib/coords/mapbased_pointing.py new file mode 100644 index 000000000..276c4dd7a --- /dev/null +++ b/sotodlib/coords/mapbased_pointing.py @@ -0,0 +1,545 @@ +import os +import re +from tqdm import tqdm +import numpy as np +from scipy import interpolate +from scipy.optimize import curve_fit + +from sotodlib import core +from sotodlib import coords +from sotodlib.coords import optics +from sotodlib.core import metadata +from sotodlib.io.metadata import write_dataset, read_dataset + +from so3g.proj import quat +from pixell import enmap +import h5py +from scipy.ndimage import maximum_filter + +def get_planet_trajectry(tod, planet, _split=20, return_model=False): + """ + Generate the trajectory of a given planet over a specified time range. + + Parameters: + tod : An axis manager + planet (str): The name of the planet for which to generate the trajectory. + _split (int, optional): Number of points to interpolate the trajectory. Defaults to 20. + return_model (bool, optional): If True, returns interpolation functions of az and el. Defaults to False. + + Returns: + If return_model is True: + tuple: Tuple containing interpolation functions for azimuth and elevation. + If return_model is False: + array: Array of quaternions representing trajectry of the planet at each timestamp. + """ + timestamps_sparse = np.linspace(tod.timestamps[0], tod.timestamps[-1], _split) + + planet_az_sparse = np.zeros_like(timestamps_sparse) + planet_el_sparse = np.zeros_like(timestamps_sparse) + for i, timestamp in enumerate(timestamps_sparse): + az, el, _ = coords.planets.get_source_azel(planet, timestamp) + planet_az_sparse[i] = az + planet_el_sparse[i] = el + planet_az_func = interpolate.interp1d(timestamps_sparse, planet_az_sparse, kind="quadratic", fill_value='extrapolate') + planet_el_func = interpolate.interp1d(timestamps_sparse, planet_el_sparse, kind="quadratic", fill_value='extrapolate') + if return_model: + return planet_az_func, planet_el_func + else: + planet_az = planet_az_func(tod.timestamps) + planet_el = planet_el_func(tod.timestamps) + q_planet = quat.rotation_lonlat(planet_az, planet_el) + return q_planet + +def get_wafer_centered_sight(tod, planet, q_planet=None, q_bs=None, q_wafer=None): + """ + Calculate the sightline vector from the focal plane, centered on the wafer, to a planet. + + Parameters: + tod : An axis manager + planet (str): The name of the planet to calculate the sightline vector. + q_planet (optional): Quaternion representing the trajectry of the planet. + If None, it will be computed using get_planet_trajectory. Defaults to None. + q_bs (optional): Quaternion representing the trajectry of the boresight. + If None, it will be computed using the current boresight angles from tod. Defaults to None. + q_wafer (optional): Quaternion representing the center of wafer to the center of boresight. + If None, it will be computed using the median of the focal plane xi and eta from tod.focal_plane. + Defaults to None. + + Returns: + Sightline vector for the planet trajectry centered on the center of the wafer. + """ + if q_planet is None: + q_planet = get_planet_trajectry(tod, planet) + if q_bs is None: + q_bs = quat.rotation_lonlat(tod.boresight.az, tod.boresight.el) + if q_wafer is None: + q_wafer = quat.rotation_xieta(np.nanmedian(tod.focal_plane.xi), + np.nanmedian(tod.focal_plane.eta)) + + xi_wafer, eta_wafer, _ = quat.decompose_xieta(q_wafer) + q_wafer_f = quat.rotation_xieta(-xi_wafer, eta_wafer) + z_to_x = quat.rotation_lonlat(0, 0) + sight = z_to_x * ~(q_bs * q_wafer_f) * q_planet + return sight + +def get_wafer_xieta(wafer_slot, optics_config_fn, xieta_bs_offset=(0., 0.), + roll_bs_offset=0., tod=None, wrap_to_tod=True,): + """ + Calculate the xi and eta coordinates for a given wafer slot on the focal plane. + + Parameters: + wafer_slot (str): The slot identifier of the wafer. + optics_config_fn (str): File name containing the optics configuration. + xieta_bs_offset (tuple): Offset in xieta coordinates for the focal plane, default is (0., 0.). + roll_bs_offset (float): Boresight roll offset. Default is 0 + tod (TimeOrderedData): TOD object to which focal plane infomation that all detectors have uniform pointing at center of the wafer is wrapped. + wrap_to_tod (bool): If True, wrap the calculated xi and eta coordinates to the Time-Ordered Data (TOD), default is True. + + Returns: + tuple: A tuple containing the calculated xi and eta coordinates for the specified wafer slot. + """ + + optics_config = optics.load_ufm_to_fp_config(optics_config_fn)['SAT'] + wafer_x, wafer_y = optics_config[wafer_slot]['dx'], optics_config[wafer_slot]['dy'] + wafer_r = np.sqrt(wafer_x**2 + wafer_y**2) + wafer_theta = np.arctan2(wafer_y, wafer_x) + + fp_to_sky = optics.sat_to_sky(optics.SAT_X, optics.SAT_LON) + lon = fp_to_sky(wafer_r) + + q1 = quat.rotation_iso(lon, 0) + + q2 = quat.rotation_iso(0, 0, np.pi/2 - wafer_theta - roll_bs_offset) + q3 = quat.rotation_xieta(xieta_bs_offset[0], xieta_bs_offset[1]) + q = q3 * q2 * q1 + + xi_wafer, eta_wafer, _ = quat.decompose_xieta(q) + if wrap_to_tod: + if tod is None: + raise ValueError('tod is not provided.') + if 'focal_plane' in tod._fields.keys(): + tod.move('focal_plane', None) + focal_plane = core.AxisManager(tod.dets) + focal_plane.wrap('xi', np.ones(tod.dets.count, dtype='float32') * xi_wafer, [(0, 'dets')]) + focal_plane.wrap('eta', np.ones(tod.dets.count, dtype='float32') * eta_wafer, [(0, 'dets')]) + focal_plane.wrap('gamma', np.zeros(tod.dets.count, dtype='float32'), [(0, 'dets')]) + tod.wrap('focal_plane', focal_plane) + tod.boresight.roll *= 0. + return xi_wafer, eta_wafer + + +def make_wafer_centered_maps(tod, planet, optics_config_fn, map_hdf, + xieta_bs_offset=(0., 0.), roll_bs_offset=None, + signal='signal', wcs_kernel=None, res=0.3*coords.DEG, cuts=None,): + """ + Generate boresight-centered maps from Time-Ordered Data (TOD) for each individual detector. + + Parameters: + tod : an axismanager object + planet (str): Name of the planet for which the trajectory is calculated. + optics_config_fn (str): File name containing the optics configuration. + map_hdf (str): Path to the HDF5 file where the maps will be saved. + xieta_bs_offset (tuple): Offset in xieta coordinates for the boresight, default is (0., 0.). + roll_bs_offset (float): Offset in roll angle for the boresight, default is None. + signal (str): Name of the signal to be used, default is 'signal'. + wcs_kernel (ndarray): WCS kernel for mapping, default is None. + res (float): Resolution of the map in degrees, default is 0.3 degrees. + cuts (tuple): Cuts to be applied to the map, default is None. + + Returns: + None + """ + + if wcs_kernel is None: + wcs_kernel = coords.get_wcs_kernel('car', 0, 0, res) + + q_planet = get_planet_trajectry(tod, planet) + q_bs = quat.rotation_lonlat(tod.boresight.az, tod.boresight.el) + + if roll_bs_offset is None: + roll_bs_offset = np.mean(tod.boresight.roll) + + # wafer + if np.unique(tod.det_info.wafer_slot).shape[0] > 1: + raise ValueError('tod include detectors from more than one wafer') + wafer_slot = tod.det_info.wafer_slot[0] + xi_wafer, eta_wafer = get_wafer_xieta(wafer_slot=wafer_slot, + xieta_bs_offset=xieta_bs_offset, + roll_bs_offset=roll_bs_offset, + tod=tod, + optics_config_fn=optics_config_fn, + wrap_to_tod=True) + + coords.planets.compute_source_flags(tod, center_on=planet, max_pix=100000000, + wrap=planet, mask={'shape':'circle', 'xyr':[0,0,8]}) + + q_wafer = quat.rotation_xieta(xi_wafer, eta_wafer) + sight = get_wafer_centered_sight(tod, planet, q_planet, q_bs, q_wafer) + + + + xi0 = tod.focal_plane.xi[0] + eta0 = tod.focal_plane.eta[0] + xi_bs_offset, eta_bs_offset = xieta_bs_offset + + tod.focal_plane.xi *= 0. + tod.focal_plane.eta *= 0. + tod.boresight.roll *= 0. + + if cuts is None: + cuts = ~tod.flags[planet] + P = coords.P.for_tod(tod=tod, wcs_kernel=wcs_kernel, comps='T', cuts=cuts, sight=sight, threads=False) + for di, det in enumerate(tqdm(tod.dets.vals)): + det_weights = np.zeros(tod.dets.count, dtype='float32') + det_weights[di] = 1. + mT_weighted = P.to_map(tod=tod, signal=signal, comps='T', det_weights=det_weights) + wT = P.to_weights(tod, signal=signal, comps='T', det_weights=det_weights) + mT = P.remove_weights(signal_map=mT_weighted, weights_map=wT, comps='T')[0] + + enmap.write_hdf(map_hdf, mT, address=det, + extra={'xi0': xi0, 'eta0': eta0, + 'xi_bs_offset': xi_bs_offset, 'eta_bs_offset': eta_bs_offset, 'roll_bs_offset': roll_bs_offset}) + return + +def detect_peak_xieta(mT, filter_size=None): + """ + Detects the peak in a given pixcell map and converts it to ξ and η coordinates. + + Parameters: + - mT (enmap.ndmap): a map object + - filter_size (int, optional): Size of the filter window for peak detection. + If not provided, it's calculated as a fraction + of the minimum dimension of mT. + + Returns: + - xi_peak (float): xi coordinate of the peak. + - eta_peak (float): eta coordinate of the peak. + - ra_peak (float): ra coordinate of the peak. + - dec_peak (float): dec coordinate of the peak. + - peak_i (int): Row index of the peak. + - peak_j (int): Column index of the peak. + """ + if filter_size is None: + filter_size = int(np.min(mT.shape)//10) + local_max = maximum_filter(mT, footprint=np.ones((filter_size, filter_size)), + mode='constant', cval=np.nan) + peak_i, peak_j = np.where(mT == np.nanmax(local_max)) + peak_i = int(np.median(peak_i)) + peak_j = int(np.median(peak_j)) + dec_grid, ra_grid = mT.posmap() + + ra_peak = ra_grid[peak_i][peak_j] + dec_peak = dec_grid[peak_i][peak_j] + xi_peak, eta_peak = _radec2xieta(ra_peak, dec_peak) + return xi_peak, eta_peak, ra_peak, dec_peak, peak_i, peak_j + +def get_center_of_mass(x, y, z, + circle_mask={'x0':0, 'y0':0, 'r_circle':3.0*coords.DEG}, + percentile_mask = {'q': 50}): + """ + Calculates the center of mass of a dataset within specified masks. + + Parameters: + - x (ndarray): Array of x-coordinates. + - y (ndarray): Array of y-coordinates. + - z (ndarray): Array of data values corresponding to the coordinates. + - circle_mask (dict, optional): Parameters defining circular mask. + Should contain keys 'x0', 'y0', and 'r_circle'. + Defaults to a circle centered at (0, 0) with radius 3.0 degrees. + - percentile_mask (dict, optional): Parameters defining percentile mask. + Should contain key 'q' representing the percentile threshold. + Defaults to the 50th percentile. + + Returns: + - x_center (float): x-coordinate of the center of mass. + - y_center (float): y-coordinate of the center of mass. + """ + mask = ~np.isnan(z) + if circle_mask is not None: + x0, y0 = circle_mask['x0'], circle_mask['y0'] + r_circle = circle_mask['r_circle'] + r = np.sqrt((x-x0)**2 + (y-y0)**2) + mask = np.logical_and(mask, rnp.nanpercentile(z[mask], q)) + + _x = x[mask] + _y = y[mask] + _z = z[mask] + + total_mass = np.nansum(_z) + x_center = np.nansum(_x * _z) / total_mass + y_center = np.nansum(_y * _z) / total_mass + return x_center, y_center + +def get_edgemap(mT, edge_avoidance=1*coords.DEG, edge_check='nan'): + """ + Generates an edge map for a given map, marking regions near the edges where data is potentially unreliable. + + Parameters: + - mT (enmap.ndmap): a map object + - edge_avoidance (float, optional): Size of the edge avoidance region, defaults to 1 degree. + - edge_check (str, optional): Method for checking edges. Should be one of {'nan', 'zero'}. + 'nan': Checks for NaN values at edges. + 'zero': Checks for zero values at edges. + Defaults to 'nan'. + + Returns: + - edge_map (enmap.ndmap): 2D boolean array representing the edge map, where True indicates regions near the edges. + """ + if edge_check not in ('nan', 'zero'): + raise ValueError('only `nan` or `zero` is supported') + + edge_map = enmap.zeros(mT.shape, mT.wcs) + edge_margin_size = int(edge_avoidance/np.mean(mT.pixshape())) + + for i, row in enumerate(mT): + if edge_check == 'nan': + nonzero_idxes = np.where(~np.isnan(row))[0] + elif edge_check == 'zero': + nonzero_idxes = np.where(row != 0)[0] + if len(nonzero_idxes>0): + edge_map[i, :nonzero_idxes[0] + edge_margin_size] = True + edge_map[i, nonzero_idxes[-1] - edge_margin_size:] = True + else: + edge_map[i, :] = True + + for j, col in enumerate(mT.T): + if edge_check == 'nan': + nonzero_idxes = np.where(~np.isnan(col))[0] + elif edge_check == 'zero': + nonzero_idxes = np.where(col != 0)[0] + if len(nonzero_idxes>0): + edge_map[:nonzero_idxes[0] + edge_margin_size, j] = True + edge_map[nonzero_idxes[-1] - edge_margin_size:, j] = True + else: + edge_map[:, j] = True + return edge_map + + + +def map_to_xieta(mT, edge_avoidance=1.0*coords.DEG, edge_check='nan', + r_tune_circle=1.0*coords.DEG, q_tune=50, + r_fit_circle=3.0*coords.DEG, beam_sigma_init=0.5*coords.DEG, ): + """ + Derive (xi,eta) coordinate of a peak from a given map and calculates the coefficient of determination (R^2) + as a measure of how well the data fits a Gaussian model around the peak. + + Parameters: + - mT (enmap.ndmap): a map object. + - edge_avoidance (float, optional): Size of the edge avoidance region, defaults to 1 degree. + - edge_check (str, optional): Method for checking edges. Should be one of {'nan', 'zero'}. Defaults to 'nan'. + - r_tune_circle (float, optional): Radius of the circle used for tuning the peak position, specified in radians. Defaults to 1 degree. + - q_tune (int, optional): Percentile threshold used for tuning the peak position. Defaults to 50. + - r_fit_circle (float, optional): Radius of the circle used for fitting the Gaussian model, specified in radians. Defaults to 3 degrees. + - beam_sigma_init (float, optional): Initial guess for the sigma parameter of the Gaussian beam, specified in radians. Defaults to 0.5 degree. + + Returns: + - xi_det (float): ξ coordinate of the detected peak. + - eta_det (float): η coordinate of the detected peak. + - R2_det (float): Coefficient of determination (R^2) indicating the goodness of fit of the data around the peak. + If no valid peak is detected or if fitting fails, returns NaN. + """ + if np.all(np.isnan(mT)): + xi_det, eta_det, R2_det = np.nan, np.nan, np.nan + + else: + xi_peak, eta_peak, ra_peak, dec_peak, peak_i, peak_j = detect_peak_xieta(mT) + if edge_avoidance > 0.: + edge_map = get_edgemap(mT, edge_avoidance=edge_avoidance, edge_check=edge_check) + edge_valid = not edge_map[peak_i, peak_j] + else: + edge_valid = True + + if edge_valid: + dec_flat, ra_flat = mT.posmap() + dec_flat, ra_flat = dec_flat.flatten(), ra_flat.flatten() + xi_flat, eta_flat = _radec2xieta(ra_flat, dec_flat) + + circle_mask = {'x0':xi_peak, 'y0':eta_peak, 'r_circle':r_tune_circle} + percentile_mask = {'q': q_tune} + xi_peak, eta_peak = get_center_of_mass(xi_flat, eta_flat, mT.flatten(), + circle_mask=circle_mask, percentile_mask=percentile_mask) + + # check R2(=coefficient of determination) + r = np.sqrt((xi_flat - xi_peak)**2 + (eta_flat - eta_peak)**2) + z = mT.flatten() + mask_fit = np.logical_and(~np.isnan(z), r R2_threshold + if np.all(~mask): + xi, eta, gamma, R2 = np.nan, np.nan, np.nan, np.nan + else: + if method == 'highest_R2': + idx = np.argmax(val['R2'][mask]) + xi, eta, gamma, R2 = val['xi'][mask][idx], val['eta'][mask][idx], val['gamma'][mask][idx], val['R2'][mask][idx] + elif method == 'mean': + xi, eta, gamma = np.mean(val['xi'][mask]), np.mean(val['eta'][mask]), np.mean(val['gamma'][mask]) + R2 = np.nan + elif method == 'median': + xi, eta, gamma = np.median(val['xi'][mask]), np.median(val['eta'][mask]), np.median(val['gamma'][mask]) + R2 = np.nan + else: + raise ValueError('Not supported method. Supported methods are `highest_R2`, `mean` or `median`') + focal_plane.rows.append((det, band, channel, R2, xi, eta, gamma)) + if save: + if output_dir is None: + output_dir = os.path.join(os.getcwd(), 'combined_pointing_results') + if not os.path.exists(output_dir): + os.makedirs(output_dir) + if save_name is None: + ctimes = np.atleast_1d([]) + wafer_slots = np.atleast_1d([]) + + for file in pointing_result_files: + filename = os.path.basename(file) + match = re.search('\d{10}', filename) + ctime = int(match.group(0) if match else None) + match = re.search('ws\d{1}', filename) + ws = match.group(0) + ctimes = np.append(ctimes, ctime) + wafer_slots = np.append(wafer_slots, ws) + ctimes = ctimes.astype('int') + wafer_slots = np.sort(np.unique(wafer_slots.astype('U3'))) + save_name = f'focal_plane_{ctimes.min()}_{ctimes.max()}_' + ''.join(wafer_slots) + '.hdf' + + write_dataset(focal_plane, os.path.join(output_dir, save_name), 'focal_plane', overwrite=True) + return focal_plane \ No newline at end of file diff --git a/sotodlib/site_pipeline/combine_focal_planes.py b/sotodlib/site_pipeline/combine_focal_planes.py new file mode 100644 index 000000000..b190074fa --- /dev/null +++ b/sotodlib/site_pipeline/combine_focal_planes.py @@ -0,0 +1 @@ +# combine multiple results \ No newline at end of file diff --git a/sotodlib/site_pipeline/make_mapbased_pointing.py b/sotodlib/site_pipeline/make_mapbased_pointing.py new file mode 100644 index 000000000..82feaf79e --- /dev/null +++ b/sotodlib/site_pipeline/make_mapbased_pointing.py @@ -0,0 +1,102 @@ +import os +import numpy as np +import yaml +import argparse + +from sotodlib import core +from sotodlib import coords +from sotodlib import tod_ops +from sotodlib.tod_ops.filters import high_pass_sine2, low_pass_sine2, fourier_filter +from sotodlib.coords import mapbased_pointing as mbp +from sotodlib.site_pipeline import update_pointing as up +from sotodlib.io.metadata import write_dataset + +from sotodlib.site_pipeline import util +logger = util.init_logger(__name__, 'make_mapbased_pointing: ') + +def filter_tod(tod, cutoff_high=0.01, cutoff_low=1.8): + if cutoff_low is not None: + tod.signal = fourier_filter(tod, filt_function=low_pass_sine2(cutoff=cutoff_low),) + if cutoff_high is not None: + tod.signal = fourier_filter(tod, filt_function=high_pass_sine2(cutoff=cutoff_high),) + return + +def tod_process(tod): + tod_ops.detrend_tod(tod) + tod_ops.apodize_cosine(tod, apodize_samps=2000) + filter_tod(tod) + tod.restrict('samps', (tod.samps.offset+2000, tod.samps.offset+tod.samps.count-2000)) + return + +def main(ctx_file, obs_id, wafer_slot, + sso_name, optics_config_fn, + map_dir, mapbased_result_dir, todbased_result_dir, + tune_by_tod=True, R2_threshold=0.3, restrict_dets=False): + + ctx = core.Context(ctx_file) + meta = ctx.get_meta(obs_id) + meta.restrict('dets', meta.dets.vals[meta.det_info.wafer_slot == wafer_slot]) + if restrict_dets: + meta.restrict('dets', meta.dets.vals[:100]) + logger.info('loading data') + tod = ctx.get_obs(meta) + logger.info('tod processing') + tod_process(tod) + + if not os.path.exists(map_dir): + logger.info(f'Make a directory: f{map_dir}') + os.makedirs(map_dir) + + logger.info(f'Making single detector maps') + map_hdf = os.path.join(map_dir, f'{obs_id}_{wafer_slot}.hdf') + mbp.make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf=map_hdf,) + + logger.info(f'Making map-based pointing results') + result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' + focal_plane_rset_mapbased = mbp.get_xieta_from_maps(map_hdf, + save=True, + output_dir=mapbased_result_dir, + filename=result_filename, + force_zero_roll=False, + edge_avoidance=1.0*coords.DEG) + + if tune_by_tod: + focal_plane = core.AxisManager(tod.dets) + focal_plane.wrap('xi', focal_plane_rset_mapbased['xi'], [(0, 'dets')]) + focal_plane.wrap('eta', focal_plane_rset_mapbased['eta'], [(0, 'dets')]) + focal_plane.wrap('gamma', focal_plane_rset_mapbased['gamma'], [(0, 'dets')]) + is_low_R2 = focal_plane_rset_mapbased['R2'] < R2_threshold + focal_plane.xi[is_low_R2] = np.nan + focal_plane.eta[is_low_R2] = np.nan + + tod.focal_plane = focal_plane + tod.flags.move(sso_name, None) + logger.info(f'Making tod-based pointing results') + focal_plane_rset_todbased = up.update_xieta(tod, sso_name, ds_factor=10, + save=True, + result_dir=todbased_result_dir, + filename=result_filename) + return + + + +def get_parser(): + parser = argparse.ArgumentParser(description="Process TOD data and update pointing") + parser.add_argument("ctx_file", type=str, help="Path to the context file") + parser.add_argument("obs_id", type=str, help="Observation ID") + parser.add_argument("wafer_slot", type=int, help="Wafer slot number") + parser.add_argument("sso_name", type=str, help="Name of Solar System Object (SSO)") + parser.add_argument("optics_config_fn", type=str, help="Path to optics configuration file") + parser.add_argument("map_dir", type=str, help="Directory to save map data") + parser.add_argument("mapbased_result_dir", type=str, help="Directory to save map-based result") + parser.add_argument("todbased_result_dir", type=str, help="Directory to save TOD-based result") + parser.add_argument("--tune_by_tod", action="store_true", help="Whether to tune by TOD data") + parser.add_argument("--R2_threshold", type=float, default=0.3, + help="Threshold for R2 value. If R2 of map-domain result is lower than the threshold,\ + the tod-fitting for that detector is skipped.") + parser.add_argument("--restrict_dets", action="store_true", + help="If specified, number of detectors are restricted to 100") + return parser + +if __name__ == '__main__': + util.main_launcher(main, get_parser) diff --git a/sotodlib/site_pipeline/update_pointing.py b/sotodlib/site_pipeline/update_pointing.py new file mode 100644 index 000000000..2fe72094b --- /dev/null +++ b/sotodlib/site_pipeline/update_pointing.py @@ -0,0 +1,196 @@ +import os +import numpy as np +import yaml +import h5py +import matplotlib.pyplot as plt +from tqdm import tqdm +import scipy +from scipy.optimize import minimize +from sotodlib.core import metadata +from sotodlib.io.metadata import write_dataset + +from sotodlib import core +from sotodlib import coords +from sotodlib import tod_ops +import so3g +from so3g.proj import quat +import sotodlib.coords.planets as planets + +from sotodlib.tod_ops import pca +from so3g.proj import Ranges, RangesMatrix +from pixell import enmap, enplot +from sotodlib.tod_ops.filters import high_pass_sine2, low_pass_sine2, fourier_filter + +from sotodlib.site_pipeline import util +logger = util.init_logger(__name__, 'update_pointing: ') + +def _gaussian2d(xi, eta, a, xi0, eta0, fwhm_xi, fwhm_eta, phi): + """Simulate a time stream with an Gaussian beam model + Args + ------ + xi, eta: cordinates in the detector's system + a: float + amplitude of the Gaussian beam model + xi0, eta0: float, float + center position of the Gaussian beam model + fwhm_xi, fwhm_eta, phi: float, float, float + fwhm along the xi, eta axis (rotated) + and the rotation angle (in radians) + + Ouput: + ------ + sim_data: 1d array of float + Time stream at sampling points given by xieta + """ + xi_rot = xi * np.cos(phi) - eta * np.sin(phi) + eta_rot = xi * np.sin(phi) + eta * np.cos(phi) + factor = 2 * np.sqrt(2 * np.log(2)) + xi_coef = -0.5 * (xi_rot - xi0) ** 2 / (fwhm_xi / factor) ** 2 + eta_coef = -0.5 * (eta_rot - eta0) ** 2 / (fwhm_eta / factor) ** 2 + sim_data = a * np.exp(xi_coef + eta_coef) + return sim_data + +def filter_tod(tod, cutoff_high=0.01, cutoff_low=1.8): + if cutoff_low is not None: + tod.signal = fourier_filter(tod, filt_function=low_pass_sine2(cutoff=cutoff_low),) + if cutoff_high is not None: + tod.signal = fourier_filter(tod, filt_function=high_pass_sine2(cutoff=cutoff_high),) + return + +def tod_process(tod): + tod_ops.detrend_tod(tod) + tod_ops.apodize_cosine(tod, apodize_samps=2000) + filter_tod(tod) + tod.restrict('samps', (tod.samps.offset+2000, tod.samps.offset+tod.samps.count-2000)) + return + +def update_xieta(tod, sso_name='moon', ds_factor=10, fwhm = 1.*coords.DEG, + save=False, result_dir=None, filename=None): + """ + Update xieta parameters for each detector by tod fitting of a point source observation + + Parameters: + - tod : an Axismanager object + - sso_name (str): Name of the Solar System Object (SSO). + - ds_factor (int): Downsampling factor for processing TOD. + - fwhm (float): Full width at half maximum of the Gaussian model. + - save (bool): Flag indicating whether to save the updated focal plane data. + - result_dir (str): Directory where the updated data will be saved. + - filename (str): Name of the file to save the updated data. + + Returns: + - focal_plane (ResultSet): ResultSet containing updated xieta parameters for each detector. + """ + mask_ds = slice(None, None, ds_factor) + + fp_isnan = np.isnan(tod.focal_plane.xi) + if np.any(fp_isnan): + tod.focal_plane.xi[fp_isnan] = 0. + tod.focal_plane.eta[fp_isnan] = 0. + tod.focal_plane.gamma[fp_isnan] = 0. + + coords.planets.compute_source_flags(tod, center_on=sso_name, max_pix=100000000, + wrap=sso_name, mask={'shape':'circle', 'xyr':[0,0,3]}) + + summed_flag = np.sum(tod.flags[sso_name].mask()[~fp_isnan], axis=0).astype('bool') + idx_hit = np.where(summed_flag)[0] + idx_first, idx_last = idx_hit[0], idx_hit[-1] + tod.restrict('samps', (tod.samps.offset+idx_first, tod.samps.offset+idx_last)) + csl = so3g.proj.CelestialSightLine.az_el(tod.timestamps[mask_ds], tod.boresight.az[mask_ds], + tod.boresight.el[mask_ds], weather="typical") + q_bore = csl.Q + + ts_ds = tod.timestamps[mask_ds] + sig_ds = tod.signal[:, mask_ds] + source_flags_ds = tod.flags[sso_name].mask()[:, mask_ds] + xieta_dict = {} + for di, det in enumerate(tqdm(tod.dets.vals)): + if fp_isnan[di]: + xieta_dict[det] = {'xi':np.nan, 'eta':np.nan, 'R2':np.nan} + else: + mask_di = source_flags_ds[di] + ts = ts_ds[mask_di] + + xieta_det = np.array([tod.focal_plane.xi[di], tod.focal_plane.eta[di]]) + q_det = so3g.proj.quat.rotation_xieta(xieta_det[0], xieta_det[1]) + d1_unix = np.median(ts) + planet = planets.SlowSource.for_named_source(sso_name, d1_unix * 1.) + ra0, dec0 = planet.pos(d1_unix) + q_obj = so3g.proj.quat.rotation_lonlat(ra0, dec0) + q_total = ~q_det * ~q_bore * q_obj + + xi_src, eta_src, _ = quat.decompose_xieta(q_total) + xieta_src = np.array([xi_src, eta_src]) + xieta_src = xieta_src[:, mask_di] + + + + sig = sig_ds[di][mask_di] + amp = np.ptp(sig) + def fit_func(xi0, eta0): + model_tod = _gaussian2d(xieta_src[0], xieta_src[1], amp, xi0, eta0, fwhm, fwhm, 0) + residual = sig - model_tod + return np.sum(residual ** 2) + + res = minimize(lambda x: fit_func(*x), [0, 0]) + R2 = 1 - res.fun/np.sum((sig - np.mean(sig))**2) + + if np.rad2deg(np.sqrt(np.sum(res.x**2))) > 1.0: + xieta_dict[det] = {'xi':np.nan, 'eta':np.nan, 'R2':np.nan} + else: + xieta_det += res.x + xieta_dict[det] = {'xi':xieta_det[0], 'eta':xieta_det[1], 'R2':R2} + + focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'band', 'channel', 'R2', 'xi', 'eta', 'gamma']) + for det in tod.dets.vals: + band = int(det.split('_')[-2]) + channel = int(det.split('_')[-1]) + focal_plane.rows.append((det, band, channel, xieta_dict[det]['R2'], + xieta_dict[det]['xi'], xieta_dict[det]['eta'], 0.)) + if save: + assert result_dir is not None + assert filename is not None + if not os.path.exists(result_dir): + os.makedirs(result_dir) + write_dataset(focal_plane, + filename=os.path.join(result_dir, filename), + address='focal_plane', + overwrite=True) + return focal_plane + +def main(ctx_file, obs_id, wafer_slot, sso_name, result_dir, + ds_factor=10, fwhm = 1.*coords.DEG, restrict_dets=False): + ctx = core.Context(ctx_file) + meta = ctx.get_meta(obs_id) + meta.restrict('dets', meta.dets.vals[meta.det_info.wafer_slot == wafer_slot]) + if restrict_dets: + meta.restrict('dets', meta.dets.vals[:100]) + + logger.info('loading data') + tod = ctx.get_obs(meta) + logger.info('tod processing') + tod_process(tod) + + if not os.path.exists(result_dir): + logger.info(f'Make a directory: f{result_dir}') + os.makedirs(result_dir) + + result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' + focal_plane_rset = update_xieta(tod=tod, sso_name=sso_name, ds_factor=ds_factor, fwhm=fwhm, + save=True, result_dir=result_dir, filename=result_filename) + return + +def get_parser(): + parser = argparse.ArgumentParser(description="Description of the script.") + parser.add_argument("ctx_file", type=str, help="Path to the context file.") + parser.add_argument("obs_id", type=str, help="Observation ID.") + parser.add_argument("wafer_slot", type=int, help="Wafer slot number.") + parser.add_argument("sso_name", type=str, help="Name of the Solar System Object (SSO).") + parser.add_argument("result_dir", type=str, help="Directory to save the result.") + parser.add_argument("--ds_factor", type=int, default=10, help="Downsampling factor for TOD processing.") + parser.add_argument("--fwhm", type=float, default=1.0, help="Full width at half maximum of the Gaussian model.") + parser.add_argument("--restrict_dets", action="store_true", help="Flag to restrict the number of detectors.") + return parser + +if __name__ == '__main__': + util.main_launcher(main, get_parser) From 0af181a458224aca53ee095a7466101da75bcc9b Mon Sep 17 00:00:00 2001 From: Tomoki Terasaki Date: Tue, 19 Mar 2024 05:36:51 +0000 Subject: [PATCH 02/27] add combine_focal_planes.py --- sotodlib/coords/mapbased_pointing.py | 65 --------- .../site_pipeline/combine_focal_planes.py | 131 +++++++++++++++++- sotodlib/site_pipeline/update_pointing.py | 2 +- 3 files changed, 131 insertions(+), 67 deletions(-) diff --git a/sotodlib/coords/mapbased_pointing.py b/sotodlib/coords/mapbased_pointing.py index 276c4dd7a..a374a621e 100644 --- a/sotodlib/coords/mapbased_pointing.py +++ b/sotodlib/coords/mapbased_pointing.py @@ -478,68 +478,3 @@ def _add_xieta(xieta1, xieta2): xi_add, eta_add, _ = quat.decompose_xieta(q_add) return xi_add, eta_add -def combine_pointings(pointing_result_files, method='mean', R2_threshold=0.3, - save=False, output_dir=None, save_name=None): - combined_dict = {} - for file in pointing_result_files: - rset = read_dataset(file, 'focal_plane') - for row in rset[:]: - if row['dets:readout_id'] not in combined_dict.keys(): - combined_dict[row['dets:readout_id']] = {} - combined_dict[row['dets:readout_id']]['band'] = row['band'] - combined_dict[row['dets:readout_id']]['channel'] = row['channel'] - - combined_dict[row['dets:readout_id']]['R2'] = np.atleast_1d([]) - combined_dict[row['dets:readout_id']]['xi'] = np.atleast_1d([]) - combined_dict[row['dets:readout_id']]['eta'] = np.atleast_1d([]) - combined_dict[row['dets:readout_id']]['gamma'] = np.atleast_1d([]) - - combined_dict[row['dets:readout_id']]['R2'] = np.append(combined_dict[row['dets:readout_id']]['R2'], row['R2']) - combined_dict[row['dets:readout_id']]['xi'] = np.append(combined_dict[row['dets:readout_id']]['xi'], row['xi']) - combined_dict[row['dets:readout_id']]['eta'] = np.append(combined_dict[row['dets:readout_id']]['eta'], row['eta']) - combined_dict[row['dets:readout_id']]['gamma'] = np.append(combined_dict[row['dets:readout_id']]['gamma'], row['gamma']) - - focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'band', 'channel', 'R2', 'xi', 'eta', 'gamma']) - for det, val in combined_dict.items(): - band = int(val['band']) - channel = int(val['channel']) - - mask = val['R2'] > R2_threshold - if np.all(~mask): - xi, eta, gamma, R2 = np.nan, np.nan, np.nan, np.nan - else: - if method == 'highest_R2': - idx = np.argmax(val['R2'][mask]) - xi, eta, gamma, R2 = val['xi'][mask][idx], val['eta'][mask][idx], val['gamma'][mask][idx], val['R2'][mask][idx] - elif method == 'mean': - xi, eta, gamma = np.mean(val['xi'][mask]), np.mean(val['eta'][mask]), np.mean(val['gamma'][mask]) - R2 = np.nan - elif method == 'median': - xi, eta, gamma = np.median(val['xi'][mask]), np.median(val['eta'][mask]), np.median(val['gamma'][mask]) - R2 = np.nan - else: - raise ValueError('Not supported method. Supported methods are `highest_R2`, `mean` or `median`') - focal_plane.rows.append((det, band, channel, R2, xi, eta, gamma)) - if save: - if output_dir is None: - output_dir = os.path.join(os.getcwd(), 'combined_pointing_results') - if not os.path.exists(output_dir): - os.makedirs(output_dir) - if save_name is None: - ctimes = np.atleast_1d([]) - wafer_slots = np.atleast_1d([]) - - for file in pointing_result_files: - filename = os.path.basename(file) - match = re.search('\d{10}', filename) - ctime = int(match.group(0) if match else None) - match = re.search('ws\d{1}', filename) - ws = match.group(0) - ctimes = np.append(ctimes, ctime) - wafer_slots = np.append(wafer_slots, ws) - ctimes = ctimes.astype('int') - wafer_slots = np.sort(np.unique(wafer_slots.astype('U3'))) - save_name = f'focal_plane_{ctimes.min()}_{ctimes.max()}_' + ''.join(wafer_slots) + '.hdf' - - write_dataset(focal_plane, os.path.join(output_dir, save_name), 'focal_plane', overwrite=True) - return focal_plane \ No newline at end of file diff --git a/sotodlib/site_pipeline/combine_focal_planes.py b/sotodlib/site_pipeline/combine_focal_planes.py index b190074fa..b8a1184a4 100644 --- a/sotodlib/site_pipeline/combine_focal_planes.py +++ b/sotodlib/site_pipeline/combine_focal_planes.py @@ -1 +1,130 @@ -# combine multiple results \ No newline at end of file +import os +import re +import glob +import numpy as np + +from sotodlib.core import metadata +from sotodlib.io.metadata import write_dataset, read_dataset + +from sotodlib.site_pipeline import util +logger = util.init_logger(__name__, 'combine_focal_planes: ') + +def combine_pointings(pointing_result_files, method='highest_R2', R2_threshold=0.3, + save=False, output_dir=None, save_name=None): + combined_dict = {} + for file in pointing_result_files: + rset = read_dataset(file, 'focal_plane') + for row in rset[:]: + if row['dets:readout_id'] not in combined_dict.keys(): + combined_dict[row['dets:readout_id']] = {} + combined_dict[row['dets:readout_id']]['band'] = row['band'] + combined_dict[row['dets:readout_id']]['channel'] = row['channel'] + + combined_dict[row['dets:readout_id']]['R2'] = np.atleast_1d([]) + combined_dict[row['dets:readout_id']]['xi'] = np.atleast_1d([]) + combined_dict[row['dets:readout_id']]['eta'] = np.atleast_1d([]) + combined_dict[row['dets:readout_id']]['gamma'] = np.atleast_1d([]) + + combined_dict[row['dets:readout_id']]['R2'] = np.append(combined_dict[row['dets:readout_id']]['R2'], row['R2']) + combined_dict[row['dets:readout_id']]['xi'] = np.append(combined_dict[row['dets:readout_id']]['xi'], row['xi']) + combined_dict[row['dets:readout_id']]['eta'] = np.append(combined_dict[row['dets:readout_id']]['eta'], row['eta']) + combined_dict[row['dets:readout_id']]['gamma'] = np.append(combined_dict[row['dets:readout_id']]['gamma'], row['gamma']) + + focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'band', 'channel', 'R2', 'xi', 'eta', 'gamma']) + for det, val in combined_dict.items(): + band = int(val['band']) + channel = int(val['channel']) + + mask = val['R2'] > R2_threshold + if np.all(~mask): + xi, eta, gamma, R2 = np.nan, np.nan, np.nan, np.nan + else: + if method == 'highest_R2': + idx = np.argmax(val['R2'][mask]) + xi, eta, gamma, R2 = val['xi'][mask][idx], val['eta'][mask][idx], val['gamma'][mask][idx], val['R2'][mask][idx] + elif method == 'mean': + xi, eta, gamma = np.mean(val['xi'][mask]), np.mean(val['eta'][mask]), np.mean(val['gamma'][mask]) + R2 = np.nan + elif method == 'median': + xi, eta, gamma = np.median(val['xi'][mask]), np.median(val['eta'][mask]), np.median(val['gamma'][mask]) + R2 = np.nan + else: + raise ValueError('Not supported method. Supported methods are `highest_R2`, `mean` or `median`') + focal_plane.rows.append((det, band, channel, R2, xi, eta, gamma)) + if save: + if output_dir is None: + output_dir = os.path.join(os.getcwd(), 'combined_pointing_results') + if not os.path.exists(output_dir): + os.makedirs(output_dir) + if save_name is None: + ctimes = np.atleast_1d([]) + wafer_slots = np.atleast_1d([]) + for file in pointing_result_files: + filename = os.path.basename(file) + match = re.search('\d{10}', filename) + ctime = int(match.group(0) if match else None) + match = re.search('ws\d{1}', filename) + ws = match.group(0) + ctimes = np.append(ctimes, ctime) + wafer_slots = np.append(wafer_slots, ws) + ctimes = ctimes.astype('int') + wafer_slots = np.sort(np.unique(wafer_slots.astype('U3'))) + save_name = f'focal_plane_{ctimes.min()}_{ctimes.max()}_' + ''.join(wafer_slots) + '.hdf' + + write_dataset(focal_plane, os.path.join(output_dir, save_name), 'focal_plane', overwrite=True) + return focal_plane + +def combine_onewafer_results(pointing_dir, ws, output_dir, filename=None, + method='highest_R2', R2_threshold=0.3,): + pointing_result_files = glob.glob(os.path.join(pointing_dir, f'focal_plane*{ws}.hdf')) + if filename is None: + filename = f'focal_plane_{ws}_combined.hdf' + _ = combine_pointings(pointing_result_files, save=True, output_dir=output_dir, save_name=filename) + return + +def combine_allwafer_results(pointing_dir, output_dir, filename=None, + method='highest_R2', R2_threshold=0.3,): + pointing_result_files = glob.glob(os.path.join(pointing_dir, 'focal_plane*.hdf')) + if filename is None: + filename = f'focal_plane_combined.hdf' + _ = combine_pointings(pointing_result_files, save=True, output_dir=output_dir, save_name=filename) + return + +def make_detabase(focal_plane_file, db_file,): + scheme = metadata.ManifestScheme().add_data_field('dataset') + db = metadata.ManifestDb(scheme=scheme) + db.add_entry({'dataset': 'focal_plane'}, filename=focal_plane_file) + db.to_file(db_file) + return + +def main(pointing_dir, output_dir=None, method='highest_R2', R2_threshold=0.3,): + if output_dir is None: + output_dir = os.path.join(os.getcwd(), 'combined_results') + + logger.info('Combining each wafer resluts') + wafer_slots = [f'ws{i}' for i in range(7)] + for ws in wafer_slots: + combine_onewafer_results(pointing_dir=pointing_dir, ws=ws, + output_dir=output_dir, filename=None, + method=method, R2_threshold=R2_threshold) + + logger.info('Combining all wafer resluts') + combine_allwafer_results(pointing_dir=pointing_dir, output_dir=output_dir, filename='focal_plane_combined.hdf', + method=method, R2_threshold=R2_threshold) + + logger.info('Making a database') + focal_plane_file = os.path.join(output_dir, 'focal_plane_combined.hdf') + db_file = os.path.join(output_dir, 'focal_plane_combined.sqlite') + make_detabase(focal_plane_file, db_file,) + return + +def get_parser(): + parser = argparse.ArgumentParser(description="Combine multiple result of pointing.") + parser.add_argument('--pointing_dir', type=str, required=True, help='Directory containing pointing result files.') + parser.add_argument('--output_dir', type=str, default=None, help='Directory to save combined results. Default is "combined_results".') + parser.add_argument('--method', type=str, default='highest_R2', choices=['highest_R2', 'mean', 'median'], help='Combination method. Default is "highest_R2".') + parser.add_argument('--R2_threshold', type=float, default=0.3, help='Threshold for R2 value. Default is 0.3.') + return parser + +if __name__ == '__main__': + util.main_launcher(main, get_parser) diff --git a/sotodlib/site_pipeline/update_pointing.py b/sotodlib/site_pipeline/update_pointing.py index 2fe72094b..d3015d688 100644 --- a/sotodlib/site_pipeline/update_pointing.py +++ b/sotodlib/site_pipeline/update_pointing.py @@ -181,7 +181,7 @@ def main(ctx_file, obs_id, wafer_slot, sso_name, result_dir, return def get_parser(): - parser = argparse.ArgumentParser(description="Description of the script.") + parser = argparse.ArgumentParser(description="Get updated result of pointings with tod-based results") parser.add_argument("ctx_file", type=str, help="Path to the context file.") parser.add_argument("obs_id", type=str, help="Observation ID.") parser.add_argument("wafer_slot", type=int, help="Wafer slot number.") From 32e24af07371958171d1134f3bfb4170ef0d6edd Mon Sep 17 00:00:00 2001 From: Tomoki Terasaki Date: Tue, 19 Mar 2024 05:47:57 +0000 Subject: [PATCH 03/27] resolve wrong git operations --- sotodlib/coords/demod.py | 42 +---------------------------------- tests/test_demod_map.py | 48 ---------------------------------------- 2 files changed, 1 insertion(+), 89 deletions(-) diff --git a/sotodlib/coords/demod.py b/sotodlib/coords/demod.py index 709c34462..c6b9dcd66 100644 --- a/sotodlib/coords/demod.py +++ b/sotodlib/coords/demod.py @@ -140,44 +140,4 @@ def make_map(tod, output = {'map': mTQU, 'weighted_map': mTQU_weighted, 'weight': wTQU} - return output - -def from_map(tod, signal_map, cuts=None, flip_gamma=True, wrap=False, modulated=False): - """ - Generate simulated TOD with HWP from a given signal map. - - Args: - tod : an axisManager object - signal_map: pixell.enmap.ndmap containing (Tmap, Qmap, Umap) representing the signal. - cuts (RangesMatrix, optional): Cuts to apply to the data. Default is None. - flip_gamma (bool, optional): Whether to flip detector coordinate. If you use the HWP, keep it `True`. Default is True. - wrap (bool, optional): Whether to wrap the simulated data. Default is False. - modulated (bool, optional): If True, return modulated signal. If False, return the demodulated signal - (`dsT`, `demodQ`, and `demodU`). Default is False. - - Returns: - `modulate==False`: A tuple containing the TOD (np.array) of dsT, demodQ and demodU. - `modulate==True` : The modulated TOD (np.array) - - """ - Tmap, Qmap, Umap = signal_map - - P = coords.P.for_tod(tod=tod, geom=signal_map.geometry, cuts=cuts, - comps='QU', hwp=flip_gamma) - dsT_sim = P.from_map(Tmap, comps='T') - demodQ_sim = P.from_map(enmap.enmap([Qmap, Umap]), comps='QU') - demodU_sim = P.from_map(enmap.enmap([Umap, -Qmap]), comps='QU') - - if modulated is False: - if wrap: - tod.wrap('dsT', dsT_sim, [(0, 'dets'), (1, 'samps')]) - tod.wrap('demodQ', demodQ_sim, [(0, 'dets'), (1, 'samps')]) - tod.wrap('demodU', demodU_sim, [(0, 'dets'), (1, 'samps')]) - return dsT_sim, demodQ_sim, demodU_sim - else: - assert 'hwp_angle' in tod._fields - signal_sim = dsT_sim + demodQ_sim*np.cos(4*tod.hwp_angle) + demodU_sim*np.sin(4*tod.hwp_angle) - if wrap: - tod.wrap('signal', signal_sim, [(0, 'dets'), (1, 'samps')]) - return signal_sim - \ No newline at end of file + return output \ No newline at end of file diff --git a/tests/test_demod_map.py b/tests/test_demod_map.py index 525f35467..36d1b7853 100644 --- a/tests/test_demod_map.py +++ b/tests/test_demod_map.py @@ -130,52 +130,4 @@ def test_10_mod_demod(self): means = [m[s].mean() for m in m0] print(means) assert(abs(means[1] - Q_stream) < TOL) - assert(abs(means[2] - U_stream) < TOL) - - def test_from_map_demodulated(self): - """Test the coords.demod.from_map function of demodulated signal. - - """ - tod = quick_tod(10, 10000) - TOL = 0.0001 - - shape, wcs = enmap.fullsky_geometry(res=0.5*coords.DEG) - signal_map = enmap.zeros((3, *shape), wcs) - T_stream, Q_stream, U_stream = 1., 0.25, 0.01 - signal_map[0] += T_stream - signal_map[1] += Q_stream - signal_map[2] += U_stream - _ = coords.demod.from_map(tod, signal_map, modulated=False, wrap=True) - - results = coords.demod.make_map(tod) - m0 = results['map'] - s = m0[1] != 0 - means = [m[s].mean() for m in m0] - assert(abs(means[1] - Q_stream) < TOL) - assert(abs(means[2] - U_stream) < TOL) - - def test_from_map_modulated(self): - """Test the coords.demod.from_map function of modulated signal. - - """ - tod = quick_tod(10, 10000) - tod.move('signal', None) - TOL = .01 - - shape, wcs = enmap.fullsky_geometry(res=0.5*coords.DEG) - signal_map = enmap.zeros((3, *shape), wcs) - - T_stream, Q_stream, U_stream = 1., 0.25, 0.01 - signal_map[0] += T_stream - signal_map[1] += Q_stream - signal_map[2] += U_stream - - _ = coords.demod.from_map(tod, signal_map, modulated=True, wrap=True) - hwp.demod_tod(tod) - results = coords.demod.make_map(tod) - - m0 = results['map'] - s = m0[1] != 0. - means = [m[s].mean() for m in m0] - assert(abs(means[1] - Q_stream) < TOL) assert(abs(means[2] - U_stream) < TOL) \ No newline at end of file From 76db4c563d54c1ab9a7029bd6498f8ff9396bda1 Mon Sep 17 00:00:00 2001 From: Tomoki Terasaki Date: Tue, 30 Apr 2024 09:02:35 +0000 Subject: [PATCH 04/27] using a new preprocess style in map_based pointing reconstruction --- ...ased_pointing.py => map_based_pointing.py} | 0 .../site_pipeline/make_mapbased_pointing.py | 128 ++++++++++-------- 2 files changed, 68 insertions(+), 60 deletions(-) rename sotodlib/coords/{mapbased_pointing.py => map_based_pointing.py} (100%) diff --git a/sotodlib/coords/mapbased_pointing.py b/sotodlib/coords/map_based_pointing.py similarity index 100% rename from sotodlib/coords/mapbased_pointing.py rename to sotodlib/coords/map_based_pointing.py diff --git a/sotodlib/site_pipeline/make_mapbased_pointing.py b/sotodlib/site_pipeline/make_mapbased_pointing.py index 82feaf79e..c898e015f 100644 --- a/sotodlib/site_pipeline/make_mapbased_pointing.py +++ b/sotodlib/site_pipeline/make_mapbased_pointing.py @@ -7,95 +7,103 @@ from sotodlib import coords from sotodlib import tod_ops from sotodlib.tod_ops.filters import high_pass_sine2, low_pass_sine2, fourier_filter -from sotodlib.coords import mapbased_pointing as mbp +from sotodlib.coords import map_based_pointing as mbp from sotodlib.site_pipeline import update_pointing as up from sotodlib.io.metadata import write_dataset from sotodlib.site_pipeline import util -logger = util.init_logger(__name__, 'make_mapbased_pointing: ') - -def filter_tod(tod, cutoff_high=0.01, cutoff_low=1.8): - if cutoff_low is not None: - tod.signal = fourier_filter(tod, filt_function=low_pass_sine2(cutoff=cutoff_low),) - if cutoff_high is not None: - tod.signal = fourier_filter(tod, filt_function=high_pass_sine2(cutoff=cutoff_high),) - return - -def tod_process(tod): - tod_ops.detrend_tod(tod) - tod_ops.apodize_cosine(tod, apodize_samps=2000) - filter_tod(tod) - tod.restrict('samps', (tod.samps.offset+2000, tod.samps.offset+tod.samps.count-2000)) - return +from sotodlib.preprocess import Pipeline +logger = util.init_logger(__name__, 'make_map_based_pointing: ') + +def main(configs, obs_id, wafer_slot, + sso_name=None, optics_config_fn=None, + single_det_maps_dir=None, map_based_result_dir=None, tod_based_result_dir=None, + tune_by_tod=None, restrict_dets_for_debug=False): + + if type(configs) == str: + configs = yaml.safe_load(open(configs, "r")) + + # Derive parameters from config file + if optics_config_fn is None: + optics_config_fn = configs.get('optics_config_fn') + if single_det_maps_dir is None: + single_det_maps_dir = configs.get('single_det_maps_dir') + if map_based_result_dir is None: + map_based_result_dir = configs.get('map_based_result_dir') + if tod_based_result_dir is None: + tod_based_result_dir = configs.get('tod_based_result_dir') -def main(ctx_file, obs_id, wafer_slot, - sso_name, optics_config_fn, - map_dir, mapbased_result_dir, todbased_result_dir, - tune_by_tod=True, R2_threshold=0.3, restrict_dets=False): + res_deg = configs.get('res_deg') + edge_avoidance_deg = configs.get('edge_avoidance_deg') + tune_by_tod = configs.get('tune_by_tod') + R2_threshold = configs.get('R2_threshold') + ds_factor = configs.get('ds_factor') - ctx = core.Context(ctx_file) - meta = ctx.get_meta(obs_id) - meta.restrict('dets', meta.dets.vals[meta.det_info.wafer_slot == wafer_slot]) - if restrict_dets: - meta.restrict('dets', meta.dets.vals[:100]) + + ctx = core.Context(configs.get('context_file')) + # If sso_name is not specified, get sso name from observation tags + obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] + if sso_name is None: + if 'moon' in obs_tags: + sso_name = 'moon' + elif 'jupiter' in obs_tags: + sso_name = 'jupiter' + else: + raise ValueError('sso_name is not specified') + + # Load data logger.info('loading data') + meta = ctx.get_meta(obs_id, dets={'wafer_slot': wafer_slot}) + if restrict_dets_for_debug is not False: + meta.restrict('dets', meta.dets.vals[:restrict_dets_for_debug]) tod = ctx.get_obs(meta) - logger.info('tod processing') - tod_process(tod) - if not os.path.exists(map_dir): - logger.info(f'Make a directory: f{map_dir}') - os.makedirs(map_dir) + # tod processing + logger.info('tod processing') + pipe = Pipeline(configs["process_pipe"], logger=logger) + proc_aman, success = pipe.run(tod) - logger.info(f'Making single detector maps') - map_hdf = os.path.join(map_dir, f'{obs_id}_{wafer_slot}.hdf') - mbp.make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf=map_hdf,) + # make single detecctor maps + logger.info(f'Making single detector maps') + os.makedirs(single_det_maps_dir, exist_ok=True) + map_hdf = os.path.join(single_det_maps_dir, f'{obs_id}_{wafer_slot}.hdf') + mbp.make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf=map_hdf, res=res_deg*coords.DEG) + # reconstruct pointing from single detector maps logger.info(f'Making map-based pointing results') result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' - focal_plane_rset_mapbased = mbp.get_xieta_from_maps(map_hdf, - save=True, - output_dir=mapbased_result_dir, + focal_plane_rset_map_based = mbp.get_xieta_from_maps(map_hdf, save=True, + output_dir=map_based_result_dir, filename=result_filename, force_zero_roll=False, - edge_avoidance=1.0*coords.DEG) + edge_avoidance = edge_avoidance_deg*coords.DEG) if tune_by_tod: focal_plane = core.AxisManager(tod.dets) - focal_plane.wrap('xi', focal_plane_rset_mapbased['xi'], [(0, 'dets')]) - focal_plane.wrap('eta', focal_plane_rset_mapbased['eta'], [(0, 'dets')]) - focal_plane.wrap('gamma', focal_plane_rset_mapbased['gamma'], [(0, 'dets')]) - is_low_R2 = focal_plane_rset_mapbased['R2'] < R2_threshold + focal_plane.wrap('xi', focal_plane_rset_map_based['xi'], [(0, 'dets')]) + focal_plane.wrap('eta', focal_plane_rset_map_based['eta'], [(0, 'dets')]) + focal_plane.wrap('gamma', focal_plane_rset_map_based['gamma'], [(0, 'dets')]) + is_low_R2 = focal_plane_rset_map_based['R2'] < R2_threshold focal_plane.xi[is_low_R2] = np.nan focal_plane.eta[is_low_R2] = np.nan tod.focal_plane = focal_plane tod.flags.move(sso_name, None) logger.info(f'Making tod-based pointing results') - focal_plane_rset_todbased = up.update_xieta(tod, sso_name, ds_factor=10, - save=True, - result_dir=todbased_result_dir, - filename=result_filename) + focal_plane_rset_tod_based = up.update_xieta(tod, sso_name, ds_factor=ds_factor, save=True, + result_dir=tod_based_result_dir, filename=result_filename) return - - def get_parser(): parser = argparse.ArgumentParser(description="Process TOD data and update pointing") - parser.add_argument("ctx_file", type=str, help="Path to the context file") - parser.add_argument("obs_id", type=str, help="Observation ID") + parser.add_argument("configs", type=str, help="Path to the configuration file") + parser.add_argument("obs_id", type=int, help="Observation ID") parser.add_argument("wafer_slot", type=int, help="Wafer slot number") - parser.add_argument("sso_name", type=str, help="Name of Solar System Object (SSO)") - parser.add_argument("optics_config_fn", type=str, help="Path to optics configuration file") - parser.add_argument("map_dir", type=str, help="Directory to save map data") - parser.add_argument("mapbased_result_dir", type=str, help="Directory to save map-based result") - parser.add_argument("todbased_result_dir", type=str, help="Directory to save TOD-based result") - parser.add_argument("--tune_by_tod", action="store_true", help="Whether to tune by TOD data") - parser.add_argument("--R2_threshold", type=float, default=0.3, - help="Threshold for R2 value. If R2 of map-domain result is lower than the threshold,\ - the tod-fitting for that detector is skipped.") - parser.add_argument("--restrict_dets", action="store_true", - help="If specified, number of detectors are restricted to 100") + parser.add_argument("--sso_name", type=str, default=None, help="Name of solar system object (e.g., 'moon', 'jupiter')") + parser.add_argument("--optics_config_fn", type=str, default=None, help="Path to optics configuration file") + parser.add_argument("--single_det_maps_dir", type=str, default=None, help="Directory to save single detector maps") + parser.add_argument("--map_based_result_dir", type=str, default=None, help="Directory to save map-based pointing results") + parser.add_argument("--tod_based_result_dir", type=str, default=None, help="Directory to save TOD-based pointing results") return parser if __name__ == '__main__': From 20b522cf00b1eec129600d17c604acdeff78c713 Mon Sep 17 00:00:00 2001 From: Tomoki Terasaki Date: Tue, 30 Apr 2024 12:23:26 +0000 Subject: [PATCH 05/27] added except method in flags.reduce --- sotodlib/core/flagman.py | 4 +++- sotodlib/preprocess/processes.py | 7 ++++++- sotodlib/site_pipeline/update_pointing.py | 2 -- 3 files changed, 9 insertions(+), 4 deletions(-) diff --git a/sotodlib/core/flagman.py b/sotodlib/core/flagman.py index 4a513278b..3a64faf50 100644 --- a/sotodlib/core/flagman.py +++ b/sotodlib/core/flagman.py @@ -172,7 +172,7 @@ def reduce(self, flags=None, method='union', wrap=False, new_flag=None, flags: List of flags to collapse together. Uses their names. If flags is None then all flags are reduced method: How to collapse the data. Accepts 'union','intersect', - or function. + 'except', or function. wrap: if True, add reduced flag to self new_flag: name of new flag, required if wrap is True remove_reduced: if True, remove all reduced flags from self @@ -198,6 +198,8 @@ def reduce(self, flags=None, method='union', wrap=False, new_flag=None, op = lambda x, y: x+y elif method == 'intersect': op = lambda x, y: x*y + elif method == 'except': + op = lambda x, y: x*~y else: op = method out = reduce(op, to_reduce) diff --git a/sotodlib/preprocess/processes.py b/sotodlib/preprocess/processes.py index 1e1293344..baae3cae6 100644 --- a/sotodlib/preprocess/processes.py +++ b/sotodlib/preprocess/processes.py @@ -782,7 +782,11 @@ def plot(self, aman, proc_aman, filename): for sso in proc_aman.sso_footprint._assignments.keys(): planet_aman = proc_aman.sso_footprint[sso] plot_sso_footprint(aman, planet_aman, sso, filename=filename.replace('{name}', f'{sso}_sso_footprint'), **self.plot_cfgs) - + +class ComputeSourceFlags(_Preprocess): + name = 'compute_source_flags' + def process(self, aman, proc_aman): + planets.compute_source_flags(aman, **self.process_cfgs) class FourierFilter(_Preprocess): """ @@ -882,3 +886,4 @@ def save(self, proc_aman, rc_aman): _Preprocess.register(SubPolyf) _Preprocess.register(DetBiasFlags) _Preprocess.register(SSOFootprint) +_Preprocess.register(ComputeSourceFlags) diff --git a/sotodlib/site_pipeline/update_pointing.py b/sotodlib/site_pipeline/update_pointing.py index d3015d688..a1cfade81 100644 --- a/sotodlib/site_pipeline/update_pointing.py +++ b/sotodlib/site_pipeline/update_pointing.py @@ -122,8 +122,6 @@ def update_xieta(tod, sso_name='moon', ds_factor=10, fwhm = 1.*coords.DEG, xi_src, eta_src, _ = quat.decompose_xieta(q_total) xieta_src = np.array([xi_src, eta_src]) xieta_src = xieta_src[:, mask_di] - - sig = sig_ds[di][mask_di] amp = np.ptp(sig) From 740f1e30d804936731b1bb32af151fa7b8f9d17b Mon Sep 17 00:00:00 2001 From: Tomoki Terasaki Date: Tue, 30 Apr 2024 12:50:03 +0000 Subject: [PATCH 06/27] added ReduceFlags function in preprocessing --- sotodlib/preprocess/processes.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/sotodlib/preprocess/processes.py b/sotodlib/preprocess/processes.py index baae3cae6..14adb046f 100644 --- a/sotodlib/preprocess/processes.py +++ b/sotodlib/preprocess/processes.py @@ -783,6 +783,11 @@ def plot(self, aman, proc_aman, filename): planet_aman = proc_aman.sso_footprint[sso] plot_sso_footprint(aman, planet_aman, sso, filename=filename.replace('{name}', f'{sso}_sso_footprint'), **self.plot_cfgs) +class ReduceFlags(_Preprocess): + name = 'reduce_flags' + def process(self, aman, proc_aman): + aman.flags.reduce(**self.process_cfgs) + class ComputeSourceFlags(_Preprocess): name = 'compute_source_flags' def process(self, aman, proc_aman): @@ -887,3 +892,4 @@ def save(self, proc_aman, rc_aman): _Preprocess.register(DetBiasFlags) _Preprocess.register(SSOFootprint) _Preprocess.register(ComputeSourceFlags) +_Preprocess.register(ReduceFlags) From 7b403a78b50a681241dd9a1d2cbacfbc328a40c5 Mon Sep 17 00:00:00 2001 From: Tomoki Terasaki Date: Thu, 2 May 2024 12:44:39 +0000 Subject: [PATCH 07/27] new fitting code --- sotodlib/coords/map_based_pointing.py | 47 ++-- .../site_pipeline/make_mapbased_pointing.py | 37 ++- sotodlib/site_pipeline/update_pointing.py | 224 +++++++++++++----- 3 files changed, 220 insertions(+), 88 deletions(-) diff --git a/sotodlib/coords/map_based_pointing.py b/sotodlib/coords/map_based_pointing.py index a374a621e..294a350a2 100644 --- a/sotodlib/coords/map_based_pointing.py +++ b/sotodlib/coords/map_based_pointing.py @@ -108,8 +108,7 @@ def get_wafer_xieta(wafer_slot, optics_config_fn, xieta_bs_offset=(0., 0.), lon = fp_to_sky(wafer_r) q1 = quat.rotation_iso(lon, 0) - - q2 = quat.rotation_iso(0, 0, np.pi/2 - wafer_theta - roll_bs_offset) + q2 = quat.rotation_iso(0, 0, np.pi/2 - wafer_theta + roll_bs_offset) q3 = quat.rotation_xieta(xieta_bs_offset[0], xieta_bs_offset[1]) q = q3 * q2 * q1 @@ -124,36 +123,39 @@ def get_wafer_xieta(wafer_slot, optics_config_fn, xieta_bs_offset=(0., 0.), focal_plane.wrap('eta', np.ones(tod.dets.count, dtype='float32') * eta_wafer, [(0, 'dets')]) focal_plane.wrap('gamma', np.zeros(tod.dets.count, dtype='float32'), [(0, 'dets')]) tod.wrap('focal_plane', focal_plane) + + # set boresight roll to zero + tod.boresight.wrap('roll_original', tod.boresight.roll, [(0, 'samps')]) tod.boresight.roll *= 0. + return xi_wafer, eta_wafer -def make_wafer_centered_maps(tod, planet, optics_config_fn, map_hdf, +def make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf, xieta_bs_offset=(0., 0.), roll_bs_offset=None, - signal='signal', wcs_kernel=None, res=0.3*coords.DEG, cuts=None,): + signal='signal', wafer_mask_deg=8., res_deg=0.3, cuts=None,): """ Generate boresight-centered maps from Time-Ordered Data (TOD) for each individual detector. Parameters: tod : an axismanager object - planet (str): Name of the planet for which the trajectory is calculated. + sso_name (str): Name of the planet for which the trajectory is calculated. optics_config_fn (str): File name containing the optics configuration. map_hdf (str): Path to the HDF5 file where the maps will be saved. xieta_bs_offset (tuple): Offset in xieta coordinates for the boresight, default is (0., 0.). roll_bs_offset (float): Offset in roll angle for the boresight, default is None. signal (str): Name of the signal to be used, default is 'signal'. wcs_kernel (ndarray): WCS kernel for mapping, default is None. - res (float): Resolution of the map in degrees, default is 0.3 degrees. + res_deg (float): Resolution of the map in degrees, default is 0.3 degrees. cuts (tuple): Cuts to be applied to the map, default is None. Returns: None """ - if wcs_kernel is None: - wcs_kernel = coords.get_wcs_kernel('car', 0, 0, res) + - q_planet = get_planet_trajectry(tod, planet) + q_planet = get_planet_trajectry(tod, sso_name) q_bs = quat.rotation_lonlat(tod.boresight.az, tod.boresight.el) if roll_bs_offset is None: @@ -170,25 +172,26 @@ def make_wafer_centered_maps(tod, planet, optics_config_fn, map_hdf, optics_config_fn=optics_config_fn, wrap_to_tod=True) - coords.planets.compute_source_flags(tod, center_on=planet, max_pix=100000000, - wrap=planet, mask={'shape':'circle', 'xyr':[0,0,8]}) + coords.planets.compute_source_flags(tod, center_on=sso_name, max_pix=100000000, + wrap='source', mask={'shape':'circle', 'xyr':[0., 0., wafer_mask_deg]}) - q_wafer = quat.rotation_xieta(xi_wafer, eta_wafer) - sight = get_wafer_centered_sight(tod, planet, q_planet, q_bs, q_wafer) - - + + q_wafer = quat.rotation_xieta(xi_wafer, eta_wafer) + sight = get_wafer_centered_sight(tod, sso_name, q_planet, q_bs, q_wafer) xi0 = tod.focal_plane.xi[0] eta0 = tod.focal_plane.eta[0] - xi_bs_offset, eta_bs_offset = xieta_bs_offset - + xi_bs_offset, eta_bs_offset = xieta_bs_offset tod.focal_plane.xi *= 0. tod.focal_plane.eta *= 0. tod.boresight.roll *= 0. + + box = np.deg2rad([[-wafer_mask_deg, -wafer_mask_deg], [wafer_mask_deg, wafer_mask_deg]]) + geom = enmap.geometry(pos=box, res=res_deg*coords.DEG) if cuts is None: - cuts = ~tod.flags[planet] - P = coords.P.for_tod(tod=tod, wcs_kernel=wcs_kernel, comps='T', cuts=cuts, sight=sight, threads=False) + cuts = ~tod.flags['source'] + P = coords.P.for_tod(tod=tod, geom=geom, comps='T', cuts=cuts, sight=sight, threads=False) for di, det in enumerate(tqdm(tod.dets.vals)): det_weights = np.zeros(tod.dets.count, dtype='float32') det_weights[di] = 1. @@ -337,8 +340,8 @@ def map_to_xieta(mT, edge_avoidance=1.0*coords.DEG, edge_check='nan', - beam_sigma_init (float, optional): Initial guess for the sigma parameter of the Gaussian beam, specified in radians. Defaults to 0.5 degree. Returns: - - xi_det (float): ξ coordinate of the detected peak. - - eta_det (float): η coordinate of the detected peak. + - xi_det (float): xi coordinate of the detected peak. + - eta_det (float): eta coordinate of the detected peak. - R2_det (float): Coefficient of determination (R^2) indicating the goodness of fit of the data around the peak. If no valid peak is detected or if fitting fails, returns NaN. """ @@ -432,7 +435,7 @@ def get_xieta_from_maps(map_hdf_file, q1 = quat.rotation_xieta(xi, eta) q2 = quat.rotation_xieta(xi_bs_offset, eta_bs_offset) - q3 = quat.rotation_iso(0, 0, roll_bs_offset) + q3 = quat.rotation_iso(0, 0, -roll_bs_offset) # q = q3 * ~q2 * q1 xieta = quat.decompose_xieta(q) xi, eta = xieta[0], xieta[1] diff --git a/sotodlib/site_pipeline/make_mapbased_pointing.py b/sotodlib/site_pipeline/make_mapbased_pointing.py index c898e015f..e566dbf21 100644 --- a/sotodlib/site_pipeline/make_mapbased_pointing.py +++ b/sotodlib/site_pipeline/make_mapbased_pointing.py @@ -33,8 +33,12 @@ def main(configs, obs_id, wafer_slot, if tod_based_result_dir is None: tod_based_result_dir = configs.get('tod_based_result_dir') - res_deg = configs.get('res_deg') - edge_avoidance_deg = configs.get('edge_avoidance_deg') + xieta_bs_offset = configs.get('xieta_bs_offset', [0., 0.]) + wafer_mask_deg = configs.get('wafer_mask_deg', 8.) + res_deg = configs.get('res_deg', 0.3) + edge_avoidance_deg = configs.get('edge_avoidance_deg', 0.3) + save_force_zero_roll = configs.get('save_force_zero_roll', True) + tune_by_tod = configs.get('tune_by_tod') R2_threshold = configs.get('R2_threshold') ds_factor = configs.get('ds_factor') @@ -67,12 +71,14 @@ def main(configs, obs_id, wafer_slot, logger.info(f'Making single detector maps') os.makedirs(single_det_maps_dir, exist_ok=True) map_hdf = os.path.join(single_det_maps_dir, f'{obs_id}_{wafer_slot}.hdf') - mbp.make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf=map_hdf, res=res_deg*coords.DEG) + mbp.make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf=map_hdf, + xieta_bs_offset=xieta_bs_offset, + wafer_mask_deg=wafer_mask_deg, res_deg=res_deg) # reconstruct pointing from single detector maps - logger.info(f'Making map-based pointing results') + logger.info(f'Saving map-based pointing results') result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' - focal_plane_rset_map_based = mbp.get_xieta_from_maps(map_hdf, save=True, + fp_rset_map_based = mbp.get_xieta_from_maps(map_hdf, save=True, output_dir=map_based_result_dir, filename=result_filename, force_zero_roll=False, @@ -80,18 +86,27 @@ def main(configs, obs_id, wafer_slot, if tune_by_tod: focal_plane = core.AxisManager(tod.dets) - focal_plane.wrap('xi', focal_plane_rset_map_based['xi'], [(0, 'dets')]) - focal_plane.wrap('eta', focal_plane_rset_map_based['eta'], [(0, 'dets')]) - focal_plane.wrap('gamma', focal_plane_rset_map_based['gamma'], [(0, 'dets')]) - is_low_R2 = focal_plane_rset_map_based['R2'] < R2_threshold + focal_plane.wrap('xi', fp_rset_map_based['xi'], [(0, 'dets')]) + focal_plane.wrap('eta', fp_rset_map_based['eta'], [(0, 'dets')]) + focal_plane.wrap('gamma', fp_rset_map_based['gamma'], [(0, 'dets')]) + is_low_R2 = fp_rset_map_based['R2'] < R2_threshold focal_plane.xi[is_low_R2] = np.nan focal_plane.eta[is_low_R2] = np.nan tod.focal_plane = focal_plane tod.flags.move(sso_name, None) logger.info(f'Making tod-based pointing results') - focal_plane_rset_tod_based = up.update_xieta(tod, sso_name, ds_factor=ds_factor, save=True, - result_dir=tod_based_result_dir, filename=result_filename) + fp_rset_tod_based = up.update_xieta(tod, sso_name, ds_factor=ds_factor, save=True, + result_dir=tod_based_result_dir, filename=result_filename) + + if save_force_zero_roll: + logger.info(f'Saving map-based pointing results (force-zero-roll)') + output_dir = map_based_result_dir + '_force_zero_roll' + fp_rset_map_based_force_zero_roll = mbp.get_xieta_from_maps(map_hdf, save=True, + output_dir=output_dir, + filename=result_filename, + force_zero_roll=True, + edge_avoidance = edge_avoidance_deg*coords.DEG) return def get_parser(): diff --git a/sotodlib/site_pipeline/update_pointing.py b/sotodlib/site_pipeline/update_pointing.py index a1cfade81..65472819c 100644 --- a/sotodlib/site_pipeline/update_pointing.py +++ b/sotodlib/site_pipeline/update_pointing.py @@ -5,9 +5,9 @@ import matplotlib.pyplot as plt from tqdm import tqdm import scipy -from scipy.optimize import minimize +from scipy.optimize import curve_fit from sotodlib.core import metadata -from sotodlib.io.metadata import write_dataset +from sotodlib.io.metadata import read_dataset, write_dataset from sotodlib import core from sotodlib import coords @@ -24,48 +24,110 @@ from sotodlib.site_pipeline import util logger = util.init_logger(__name__, 'update_pointing: ') -def _gaussian2d(xi, eta, a, xi0, eta0, fwhm_xi, fwhm_eta, phi): +def _gaussian2d(xieta, xi0, eta0, fwhm_xi, fwhm_eta, phi, a): """Simulate a time stream with an Gaussian beam model Args ------ xi, eta: cordinates in the detector's system + xi0, eta0: float, float + center position of the Gaussian beam model + fwhm_xi, fwhm_eta, phi: float, float, float + fwhm along the xi, eta axis (rotated) + and the rotation angle (in radians) a: float amplitude of the Gaussian beam model + + Ouput: + ------ + sim_data: 1d array of float + Time stream at sampling points given by xieta + """ + xi, eta = xieta + xi_rot = xi * np.cos(phi) - eta * np.sin(phi) + eta_rot = xi * np.sin(phi) + eta * np.cos(phi) + factor = 2 * np.sqrt(2 * np.log(2)) + xi_coef = -0.5 * (xi_rot - xi0) ** 2 / (fwhm_xi / factor) ** 2 + eta_coef = -0.5 * (eta_rot - eta0) ** 2 / (fwhm_eta / factor) ** 2 + sim_data = a * np.exp(xi_coef + eta_coef) + return sim_data + +def _gaussian2d_nonlin(xieta, xi0, eta0, fwhm_xi, fwhm_eta, phi, a, b2,): + """Simulate a time stream with an Gaussian beam model with non-linear response + Args + ------ + xi, eta: cordinates in the detector's system xi0, eta0: float, float center position of the Gaussian beam model fwhm_xi, fwhm_eta, phi: float, float, float fwhm along the xi, eta axis (rotated) and the rotation angle (in radians) + a: float + amplitude of the Gaussian beam model + b2: float + coefficient of 2nd-order term Ouput: ------ sim_data: 1d array of float Time stream at sampling points given by xieta """ + xi, eta = xieta xi_rot = xi * np.cos(phi) - eta * np.sin(phi) eta_rot = xi * np.sin(phi) + eta * np.cos(phi) factor = 2 * np.sqrt(2 * np.log(2)) xi_coef = -0.5 * (xi_rot - xi0) ** 2 / (fwhm_xi / factor) ** 2 eta_coef = -0.5 * (eta_rot - eta0) ** 2 / (fwhm_eta / factor) ** 2 - sim_data = a * np.exp(xi_coef + eta_coef) + _y = np.exp(xi_coef + eta_coef) + sim_data = a * (_y + b2*_y**2) return sim_data -def filter_tod(tod, cutoff_high=0.01, cutoff_low=1.8): - if cutoff_low is not None: - tod.signal = fourier_filter(tod, filt_function=low_pass_sine2(cutoff=cutoff_low),) - if cutoff_high is not None: - tod.signal = fourier_filter(tod, filt_function=high_pass_sine2(cutoff=cutoff_high),) - return +# def filter_tod(tod, cutoff_high=0.01, cutoff_low=1.8): +# if cutoff_low is not None: +# tod.signal = fourier_filter(tod, filt_function=low_pass_sine2(cutoff=cutoff_low),) +# if cutoff_high is not None: +# tod.signal = fourier_filter(tod, filt_function=high_pass_sine2(cutoff=cutoff_high),) +# return -def tod_process(tod): - tod_ops.detrend_tod(tod) - tod_ops.apodize_cosine(tod, apodize_samps=2000) - filter_tod(tod) - tod.restrict('samps', (tod.samps.offset+2000, tod.samps.offset+tod.samps.count-2000)) - return +# def tod_process(tod): +# tod_ops.detrend_tod(tod) +# tod_ops.apodize_cosine(tod, apodize_samps=2000) +# filter_tod(tod) +# tod.restrict('samps', (tod.samps.offset+2000, tod.samps.offset+tod.samps.count-2000)) +# return + +def wrap_fp_from_hdf(tod, fp_hdf_file, data_set='focal_plane'): + fp_rset = read_dataset(fp_hdf_file, data_set) + tod.restrict('dets', tod.dets.vals[np.in1d(tod.dets.vals, fp_rset['dets:readout_id'])]) + focal_plane = core.AxisManager(tod.dets) + focal_plane.wrap_new('xi', shape=('dets', )) + focal_plane.wrap_new('eta', shape=('dets', )) + focal_plane.wrap_new('gamma', shape=('dets', )) -def update_xieta(tod, sso_name='moon', ds_factor=10, fwhm = 1.*coords.DEG, - save=False, result_dir=None, filename=None): + for di, det in enumerate(tod.dets.vals): + di_rset = np.where(fp_rset['dets:readout_id'] == det)[0][0] + focal_plane.xi[di] = fp_rset['xi'][di_rset] + focal_plane.eta[di] = fp_rset['eta'][di_rset] + focal_plane.gamma[di] = fp_rset['gamma'][di_rset] + + if 'focal_plane' in tod._fields.keys(): + tod.move('focal_plane', None) + tod.wrap('focal_plane', focal_plane) + return + + +def update_xieta(tod, + sso_name='moon', + fp_hdf_file=None, + input_force_zero_roll=False, + pipe=None, + ds_factor=10, + mask_deg=3, + fit_func_name = '_gaussian2d_nonlin', + fwhm_init_deg = 0.5, + error_estimation_method='force_one_redchi2', # rms_from_data + flag_name_rms_calc = 'source', + flag_rms_calc_exclusive = True, + save=False, result_dir=None, filename=None): """ Update xieta parameters for each detector by tod fitting of a point source observation @@ -81,70 +143,122 @@ def update_xieta(tod, sso_name='moon', ds_factor=10, fwhm = 1.*coords.DEG, Returns: - focal_plane (ResultSet): ResultSet containing updated xieta parameters for each detector. """ - mask_ds = slice(None, None, ds_factor) + # if focal_plane result is specified, use the information as a prior + if fp_hdf_file is not None: + wrap_fp_from_hdf(tod, fp_hdf_file) + + # set dets without focal_plane info to have (xi, eta, gamma) = (0, 0, 0), just to avoid error + xieta_isnan = (np.isnan(tod.focal_plane.xi)) | (np.isnan(tod.focal_plane.eta)) + gamma_isnan = np.isnan(tod.focal_plane.gamma) + tod.focal_plane.xi[xieta_isnan] = 0. + tod.focal_plane.eta[xieta_isnan] = 0. + tod.focal_plane.gamma[gamma_isnan] = 0. - fp_isnan = np.isnan(tod.focal_plane.xi) - if np.any(fp_isnan): - tod.focal_plane.xi[fp_isnan] = 0. - tod.focal_plane.eta[fp_isnan] = 0. - tod.focal_plane.gamma[fp_isnan] = 0. + # If input focal_plane is a result with `force_zero_roll`, set the roll to be zero + # Original value is stored to `roll_original` + if input_force_zero_roll: + if 'roll_original' in tod.boresight._fields.keys(): + pass + else: + tod.boresight.wrap('roll_original', tod.boresight.roll, [(0, 'samps')]) + tod.boresight.roll *= 0. - coords.planets.compute_source_flags(tod, center_on=sso_name, max_pix=100000000, - wrap=sso_name, mask={'shape':'circle', 'xyr':[0,0,3]}) + # compute source flags + if 'source' in tod.flags._fields.keys(): + tod.flags.move('source', None) + coords.planets.compute_source_flags(tod, + center_on=sso_name, + max_pix=1e10, + wrap='source', + mask={'shape':'circle', 'xyr':[0.,0.,mask_deg]}) - summed_flag = np.sum(tod.flags[sso_name].mask()[~fp_isnan], axis=0).astype('bool') + # restrict data to duration when at least one detector hit the source + summed_flag = np.sum(tod.flags['source'].mask()[~xieta_isnan], axis=0).astype('bool') idx_hit = np.where(summed_flag)[0] idx_first, idx_last = idx_hit[0], idx_hit[-1] tod.restrict('samps', (tod.samps.offset+idx_first, tod.samps.offset+idx_last)) - csl = so3g.proj.CelestialSightLine.az_el(tod.timestamps[mask_ds], tod.boresight.az[mask_ds], - tod.boresight.el[mask_ds], weather="typical") - q_bore = csl.Q + # run preprocess pipeline if provided + if pipe is not None: + proc_aman, success = pipe.run(tod) + + # get rms of flagged region for later error estimation + if flag_rms_calc_exclusive: + mask_for_rms_calc = tod.flags[flag_name_rms_calc].mask() + else: + mask_for_rms_calc = ~tod.flags[flag_name_rms_calc].mask() + rms = np.ma.std(np.ma.masked_array(tod.signal, mask_for_rms_calc), axis=1).data + tod.wrap('rms', rms, [(0, 'dets')]) + + # use downsampled data for faster fitting + mask_ds = slice(None, None, ds_factor) ts_ds = tod.timestamps[mask_ds] + q_bore = so3g.proj.CelestialSightLine.az_el(ts_ds, tod.boresight.az[mask_ds], + tod.boresight.el[mask_ds], weather="typical").Q + q_bore_roll = quat.rotation_iso(0, 0, np.median(tod.boresight.roll)) sig_ds = tod.signal[:, mask_ds] - source_flags_ds = tod.flags[sso_name].mask()[:, mask_ds] + source_flags_ds = tod.flags['source'].mask()[:, mask_ds] + + # fit each detector data xieta_dict = {} for di, det in enumerate(tqdm(tod.dets.vals)): - if fp_isnan[di]: - xieta_dict[det] = {'xi':np.nan, 'eta':np.nan, 'R2':np.nan} + mask_di = source_flags_ds[di] + if np.any([xieta_isnan[di], np.all(mask_di==False), tod.rms[di]==0.]): + xieta_dict[det] = {'xi': np.nan, 'eta': np.nan, 'xi_err': np.nan, 'eta_err': np.nan, + 'R2': np.nan, 'redchi2': np.nan} else: - mask_di = source_flags_ds[di] ts = ts_ds[mask_di] - + d1_unix = np.median(ts) xieta_det = np.array([tod.focal_plane.xi[di], tod.focal_plane.eta[di]]) + q_det = so3g.proj.quat.rotation_xieta(xieta_det[0], xieta_det[1]) - d1_unix = np.median(ts) planet = planets.SlowSource.for_named_source(sso_name, d1_unix * 1.) ra0, dec0 = planet.pos(d1_unix) q_obj = so3g.proj.quat.rotation_lonlat(ra0, dec0) - q_total = ~q_det * ~q_bore * q_obj + q_total = ~q_det * ~q_bore_roll * ~q_bore * q_obj xi_src, eta_src, _ = quat.decompose_xieta(q_total) xieta_src = np.array([xi_src, eta_src]) xieta_src = xieta_src[:, mask_di] - sig = sig_ds[di][mask_di] - amp = np.ptp(sig) - def fit_func(xi0, eta0): - model_tod = _gaussian2d(xieta_src[0], xieta_src[1], amp, xi0, eta0, fwhm, fwhm, 0) - residual = sig - model_tod - return np.sum(residual ** 2) - res = minimize(lambda x: fit_func(*x), [0, 0]) - R2 = 1 - res.fun/np.sum((sig - np.mean(sig))**2) + ptp_val = np.ptp(np.percentile(sig, [0.1, 99.9])) + + if fit_func_name == '_gaussian2d': + p0 = (0., 0., fwhm_init_deg*coords.DEG, fwhm_init_deg*coords.DEG, 0., ptp_val) + fit_func = _gaussian2d + elif fit_func_name == '_gaussian2d_nonlin': + p0 = (0., 0., fwhm_init_deg*coords.DEG, fwhm_init_deg*coords.DEG, 0., ptp_val, -0.1,) + fit_func = _gaussian2d_nonlin + + popt, pcov = curve_fit(fit_func, xdata=xieta_src, ydata=sig, p0=p0, sigma=tod.rms[di]*np.ones_like(sig), + absolute_sigma=True, maxfev=int(1e5)) - if np.rad2deg(np.sqrt(np.sum(res.x**2))) > 1.0: - xieta_dict[det] = {'xi':np.nan, 'eta':np.nan, 'R2':np.nan} + chi2 = np.sum(((fit_func(xieta_src, *popt) - sig)/tod.rms[di])**2) + redchi2 = chi2 / (np.prod(xieta_src.shape) - popt.shape[0]) + R2 = 1. - np.sum((fit_func(xieta_src, *popt) - sig)**2) / np.sum((sig - sig.mean())**2) + xi_opt, eta_opt = popt[0], popt[1] + + if error_estimation_method == 'rms_from_data': + xi_err, eta_err = np.sqrt(pcov[0,0]), np.sqrt(pcov[1,1]) + elif error_estimation_method == 'force_one_redchi2': + # The error of (xi, eta) is equivalent the case if the error bar of each data point is set + # as the reduced chi-square is equal to unity. + xi_err, eta_err = np.sqrt(pcov[0,0] * redchi2), np.sqrt(pcov[1,1] * redchi2) + redchi2 = 1. else: - xieta_det += res.x - xieta_dict[det] = {'xi':xieta_det[0], 'eta':xieta_det[1], 'R2':R2} - - focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'band', 'channel', 'R2', 'xi', 'eta', 'gamma']) + raise NameError("Unsupported name for 'error_estimation_method'") + + xieta_det += np.array([xi_opt, eta_opt]) + xieta_dict[det] = {'xi': xieta_det[0], 'eta': xieta_det[1], 'xi_err': xi_err, 'eta_err': eta_err, + 'R2': R2, 'redchi2': redchi2} + + focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', 'xi_err', 'eta_err', 'R2', 'redchi2']) for det in tod.dets.vals: - band = int(det.split('_')[-2]) - channel = int(det.split('_')[-1]) - focal_plane.rows.append((det, band, channel, xieta_dict[det]['R2'], - xieta_dict[det]['xi'], xieta_dict[det]['eta'], 0.)) + focal_plane.rows.append((det, xieta_dict[det]['xi'], xieta_dict[det]['eta'], 0., + xieta_dict[det]['xi_err'], xieta_dict[det]['eta_err'], + xieta_dict[det]['R2'], xieta_dict[det]['redchi2'], + )) if save: assert result_dir is not None assert filename is not None From 5101736a1950e3df2a0ccdfaa731271519f9e889 Mon Sep 17 00:00:00 2001 From: Tomoki Terasaki Date: Tue, 7 May 2024 02:27:21 +0000 Subject: [PATCH 08/27] added argument save_normal_roll --- .../site_pipeline/make_mapbased_pointing.py | 121 ++++--- sotodlib/site_pipeline/update_pointing.py | 296 +++++++++++------- 2 files changed, 256 insertions(+), 161 deletions(-) diff --git a/sotodlib/site_pipeline/make_mapbased_pointing.py b/sotodlib/site_pipeline/make_mapbased_pointing.py index e566dbf21..23f15416d 100644 --- a/sotodlib/site_pipeline/make_mapbased_pointing.py +++ b/sotodlib/site_pipeline/make_mapbased_pointing.py @@ -16,7 +16,7 @@ logger = util.init_logger(__name__, 'make_map_based_pointing: ') def main(configs, obs_id, wafer_slot, - sso_name=None, optics_config_fn=None, + sso_name=None, single_det_maps_dir=None, map_based_result_dir=None, tod_based_result_dir=None, tune_by_tod=None, restrict_dets_for_debug=False): @@ -24,35 +24,42 @@ def main(configs, obs_id, wafer_slot, configs = yaml.safe_load(open(configs, "r")) # Derive parameters from config file - if optics_config_fn is None: - optics_config_fn = configs.get('optics_config_fn') + ctx = core.Context(configs.get('context_file')) if single_det_maps_dir is None: single_det_maps_dir = configs.get('single_det_maps_dir') if map_based_result_dir is None: map_based_result_dir = configs.get('map_based_result_dir') if tod_based_result_dir is None: tod_based_result_dir = configs.get('tod_based_result_dir') - + optics_config_fn = configs.get('optics_config_fn') xieta_bs_offset = configs.get('xieta_bs_offset', [0., 0.]) wafer_mask_deg = configs.get('wafer_mask_deg', 8.) res_deg = configs.get('res_deg', 0.3) edge_avoidance_deg = configs.get('edge_avoidance_deg', 0.3) + save_normal_roll = configs.get('save_normal_roll', True) save_force_zero_roll = configs.get('save_force_zero_roll', True) + # parameters for tod tuning tune_by_tod = configs.get('tune_by_tod') - R2_threshold = configs.get('R2_threshold') - ds_factor = configs.get('ds_factor') + if tune_by_tod: + tod_ds_factor = configs.get('tod_ds_factor') + tod_mask_deg = configs.get('tod_mask_deg') + tod_fit_func_name = configs.get('tod_fit_func_name') + tod_max_non_linear_order = configs.get('tod_max_non_linear_order') + tod_fwhm_init_deg = configs.get('tod_fwhm_init_deg') + tod_error_estimation_method = configs.get('tod_error_estimation_method') + tod_flag_name_rms_calc = configs.get('tod_flag_name_rms_calc') + tod_flag_rms_calc_exclusive = configs.get('tod_flag_rms_calc_exclusive') - ctx = core.Context(configs.get('context_file')) # If sso_name is not specified, get sso name from observation tags obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] if sso_name is None: - if 'moon' in obs_tags: - sso_name = 'moon' - elif 'jupiter' in obs_tags: - sso_name = 'jupiter' - else: + known_source_names = ['moon', 'jupiter'] + for _source_name in known_source_names: + if _source_name in obs_tags: + sso_name = _source_name + if _source_name is None: raise ValueError('sso_name is not specified') # Load data @@ -76,49 +83,83 @@ def main(configs, obs_id, wafer_slot, wafer_mask_deg=wafer_mask_deg, res_deg=res_deg) # reconstruct pointing from single detector maps - logger.info(f'Saving map-based pointing results') - result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' - fp_rset_map_based = mbp.get_xieta_from_maps(map_hdf, save=True, - output_dir=map_based_result_dir, - filename=result_filename, - force_zero_roll=False, - edge_avoidance = edge_avoidance_deg*coords.DEG) - - if tune_by_tod: - focal_plane = core.AxisManager(tod.dets) - focal_plane.wrap('xi', fp_rset_map_based['xi'], [(0, 'dets')]) - focal_plane.wrap('eta', fp_rset_map_based['eta'], [(0, 'dets')]) - focal_plane.wrap('gamma', fp_rset_map_based['gamma'], [(0, 'dets')]) - is_low_R2 = fp_rset_map_based['R2'] < R2_threshold - focal_plane.xi[is_low_R2] = np.nan - focal_plane.eta[is_low_R2] = np.nan - - tod.focal_plane = focal_plane - tod.flags.move(sso_name, None) - logger.info(f'Making tod-based pointing results') - fp_rset_tod_based = up.update_xieta(tod, sso_name, ds_factor=ds_factor, save=True, - result_dir=tod_based_result_dir, filename=result_filename) + if save_normal_roll: + logger.info(f'Saving map-based pointing results') + result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' + fp_rset_map_based = mbp.get_xieta_from_maps(map_hdf, save=True, + output_dir=map_based_result_dir, + filename=result_filename, + force_zero_roll=False, + edge_avoidance = edge_avoidance_deg*coords.DEG) + + if tune_by_tod: + logger.info(f'Making tod-based pointing results') + up.wrap_fp_rset(tod, fp_rset_map_based) + fp_rset_tod_based = up.update_xieta( tod, + sso_name=sso_name, + fp_hdf_file=None, + force_zero_roll=False, + pipe=None, + ds_factor=tod_ds_factor, + mask_deg=tod_mask_deg, + fit_func_name = tod_fit_func_name, + max_non_linear_order = tod_max_non_linear_order, + fwhm_init_deg = tod_fwhm_init_deg, + error_estimation_method=tod_error_estimation_method, + flag_name_rms_calc = tod_flag_name_rms_calc, + flag_rms_calc_exclusive = tod_flag_rms_calc_exclusive, + ) + os.makedirs(tod_based_result_dir, exist_ok=True) + write_dataset(fp_rset_tod_based, + filename=os.path.join(tod_based_result_dir, f'focal_plane_{obs_id}_{wafer_slot}.hdf'), + address='focal_plane', + overwrite=True) if save_force_zero_roll: logger.info(f'Saving map-based pointing results (force-zero-roll)') - output_dir = map_based_result_dir + '_force_zero_roll' + map_based_result_dir_force_zero_roll = map_based_result_dir + '_force_zero_roll' fp_rset_map_based_force_zero_roll = mbp.get_xieta_from_maps(map_hdf, save=True, - output_dir=output_dir, + output_dir=map_based_result_dir_force_zero_roll, filename=result_filename, force_zero_roll=True, - edge_avoidance = edge_avoidance_deg*coords.DEG) + edge_avoidance = edge_avoidance_deg*coords.DEG) + if tune_by_tod: + logger.info(f'Making tod-based pointing results (force-zero-roll)') + up.wrap_fp_rset(tod, fp_rset_map_based_force_zero_roll) + tod_based_result_dir_force_zero_roll = tod_based_result_dir + '_force_zero_roll' + fp_rset_tod_based_force_zero_roll = up.update_xieta( tod, + sso_name=sso_name, + fp_hdf_file=None, + force_zero_roll=False, + pipe=None, + ds_factor=tod_ds_factor, + mask_deg=tod_mask_deg, + fit_func_name = tod_fit_func_name, + max_non_linear_order = tod_max_non_linear_order, + fwhm_init_deg = tod_fwhm_init_deg, + error_estimation_method=tod_error_estimation_method, + flag_name_rms_calc = tod_flag_name_rms_calc, + flag_rms_calc_exclusive = tod_flag_rms_calc_exclusive, + ) + os.makedirs(tod_based_result_dir_force_zero_roll, exist_ok=True) + write_dataset(fp_rset_tod_based_force_zero_roll, + filename=os.path.join(tod_based_result_dir_force_zero_roll, f'focal_plane_{obs_id}_{wafer_slot}.hdf'), + address='focal_plane', + overwrite=True) + + return def get_parser(): parser = argparse.ArgumentParser(description="Process TOD data and update pointing") parser.add_argument("configs", type=str, help="Path to the configuration file") - parser.add_argument("obs_id", type=int, help="Observation ID") - parser.add_argument("wafer_slot", type=int, help="Wafer slot number") + parser.add_argument("obs_id", type=str, help="Observation id") + parser.add_argument("wafer_slot", type=str, help="Wafer slot") parser.add_argument("--sso_name", type=str, default=None, help="Name of solar system object (e.g., 'moon', 'jupiter')") - parser.add_argument("--optics_config_fn", type=str, default=None, help="Path to optics configuration file") parser.add_argument("--single_det_maps_dir", type=str, default=None, help="Directory to save single detector maps") parser.add_argument("--map_based_result_dir", type=str, default=None, help="Directory to save map-based pointing results") parser.add_argument("--tod_based_result_dir", type=str, default=None, help="Directory to save TOD-based pointing results") + parser.add_argument("--restrict_dets_for_debug", type=int, default=False) return parser if __name__ == '__main__': diff --git a/sotodlib/site_pipeline/update_pointing.py b/sotodlib/site_pipeline/update_pointing.py index 65472819c..1d107c0e7 100644 --- a/sotodlib/site_pipeline/update_pointing.py +++ b/sotodlib/site_pipeline/update_pointing.py @@ -24,8 +24,8 @@ from sotodlib.site_pipeline import util logger = util.init_logger(__name__, 'update_pointing: ') -def _gaussian2d(xieta, xi0, eta0, fwhm_xi, fwhm_eta, phi, a): - """Simulate a time stream with an Gaussian beam model +def gaussian2d_nonlin(xieta, xi0, eta0, fwhm_xi, fwhm_eta, phi, a, nonlin_coeffs): + """ An Gaussian beam model with non-linear response Args ------ xi, eta: cordinates in the detector's system @@ -36,11 +36,12 @@ def _gaussian2d(xieta, xi0, eta0, fwhm_xi, fwhm_eta, phi, a): and the rotation angle (in radians) a: float amplitude of the Gaussian beam model - + nonlin_coeffs: float + Coefficient of non-linear term normalized by linear term (from 2nd term). + The order is ascending. Ouput: ------ - sim_data: 1d array of float - Time stream at sampling points given by xieta + Model at xieta """ xi, eta = xieta xi_rot = xi * np.cos(phi) - eta * np.sin(phi) @@ -48,55 +49,18 @@ def _gaussian2d(xieta, xi0, eta0, fwhm_xi, fwhm_eta, phi, a): factor = 2 * np.sqrt(2 * np.log(2)) xi_coef = -0.5 * (xi_rot - xi0) ** 2 / (fwhm_xi / factor) ** 2 eta_coef = -0.5 * (eta_rot - eta0) ** 2 / (fwhm_eta / factor) ** 2 - sim_data = a * np.exp(xi_coef + eta_coef) - return sim_data - -def _gaussian2d_nonlin(xieta, xi0, eta0, fwhm_xi, fwhm_eta, phi, a, b2,): - """Simulate a time stream with an Gaussian beam model with non-linear response - Args - ------ - xi, eta: cordinates in the detector's system - xi0, eta0: float, float - center position of the Gaussian beam model - fwhm_xi, fwhm_eta, phi: float, float, float - fwhm along the xi, eta axis (rotated) - and the rotation angle (in radians) - a: float - amplitude of the Gaussian beam model - b2: float - coefficient of 2nd-order term + lin_gauss = np.exp(xi_coef + eta_coef) + polycoeffs_discending = np.hstack([nonlin_coeffs[::-1], [1, 0]]) + return a * np.poly1d(polycoeffs_discending)(lin_gauss) - Ouput: - ------ - sim_data: 1d array of float - Time stream at sampling points given by xieta +def wrapper_gaussian2d_nonlin(xieta, xi0, eta0, fwhm_xi, fwhm_eta, phi, a, *args): """ - xi, eta = xieta - xi_rot = xi * np.cos(phi) - eta * np.sin(phi) - eta_rot = xi * np.sin(phi) + eta * np.cos(phi) - factor = 2 * np.sqrt(2 * np.log(2)) - xi_coef = -0.5 * (xi_rot - xi0) ** 2 / (fwhm_xi / factor) ** 2 - eta_coef = -0.5 * (eta_rot - eta0) ** 2 / (fwhm_eta / factor) ** 2 - _y = np.exp(xi_coef + eta_coef) - sim_data = a * (_y + b2*_y**2) - return sim_data - -# def filter_tod(tod, cutoff_high=0.01, cutoff_low=1.8): -# if cutoff_low is not None: -# tod.signal = fourier_filter(tod, filt_function=low_pass_sine2(cutoff=cutoff_low),) -# if cutoff_high is not None: -# tod.signal = fourier_filter(tod, filt_function=high_pass_sine2(cutoff=cutoff_high),) -# return - -# def tod_process(tod): -# tod_ops.detrend_tod(tod) -# tod_ops.apodize_cosine(tod, apodize_samps=2000) -# filter_tod(tod) -# tod.restrict('samps', (tod.samps.offset+2000, tod.samps.offset+tod.samps.count-2000)) -# return + A wrapper for `gaussian2d_nonlin` + """ + nonlin_coeffs = np.array(args) + return gaussian2d_nonlin(xieta, xi0, eta0, fwhm_xi, fwhm_eta, phi, a, nonlin_coeffs) -def wrap_fp_from_hdf(tod, fp_hdf_file, data_set='focal_plane'): - fp_rset = read_dataset(fp_hdf_file, data_set) +def wrap_fp_rset(tod, fp_rset): tod.restrict('dets', tod.dets.vals[np.in1d(tod.dets.vals, fp_rset['dets:readout_id'])]) focal_plane = core.AxisManager(tod.dets) focal_plane.wrap_new('xi', shape=('dets', )) @@ -114,31 +78,68 @@ def wrap_fp_from_hdf(tod, fp_hdf_file, data_set='focal_plane'): tod.wrap('focal_plane', focal_plane) return +def wrap_fp_from_hdf(tod, fp_hdf_file, data_set='focal_plane'): + fp_rset = read_dataset(fp_hdf_file, data_set) + wrap_fp_rset(tod, fp_rset) + return + def update_xieta(tod, sso_name='moon', fp_hdf_file=None, - input_force_zero_roll=False, + save_dir=None, pipe=None, + force_zero_roll=False, ds_factor=10, mask_deg=3, - fit_func_name = '_gaussian2d_nonlin', + fit_func_name = 'gaussian2d_nonlin', + max_non_linear_order = 1, fwhm_init_deg = 0.5, error_estimation_method='force_one_redchi2', # rms_from_data flag_name_rms_calc = 'source', flag_rms_calc_exclusive = True, - save=False, result_dir=None, filename=None): + save=True, ): """ - Update xieta parameters for each detector by tod fitting of a point source observation + Update xieta parameters for each detector by TOD fitting of a point source observation. Parameters: - - tod : an Axismanager object - - sso_name (str): Name of the Solar System Object (SSO). - - ds_factor (int): Downsampling factor for processing TOD. - - fwhm (float): Full width at half maximum of the Gaussian model. - - save (bool): Flag indicating whether to save the updated focal plane data. - - result_dir (str): Directory where the updated data will be saved. - - filename (str): Name of the file to save the updated data. + - tod : + an Axismanager object + - sso_name (str): + Name of the Solar System Object (SSO). + - fp_hdf_file (str or None): + Path to the HDF file containing focal plane information. Default is None. + If None, tod.focal_plane is used for focal plane information. + - save_dir (str or None): + Directory where the updated data will be saved. Required if save is True. + - force_zero_roll (bool): + Flag indicating whether to force the roll to be zero. Default is False. + If True, input and output focal plane information assumes force_zero_roll condition. + - pipe (Pipeline or None): + Preprocessing pipeline to be applied to the TOD. Default is None, which + do not apply any processing. + - ds_factor (int): + Downsampling factor for fitting. Default is 10. + - mask_deg (float): + Mask radius in degrees for source flagging. Default is 3. + - fit_func_name (str): + Name of the fitting function. Default is 'gaussian2d_nonlin'. 'gaussian2d_nonlin' is only supported. + - max_non_linear_order (int): + Maximum non-linear order for fitting function. Default is 1. If you want to use simple gaussian set it to be 1. + Higher order is for the case that detector response is distorted by non-point-like source or too-strogng source, such as the Moon. + - fwhm_init_deg (float): + Initial guess for full width at half maximum in degrees. Default is 0.5. + - error_estimation_method (str): + Method for error estimation. Default is 'rms_from_data'. 'rms_from_data' and 'force_one_redchi2' are supported. + If 'rms_from_data', errorbar of each data point is set by root-mean-square of the data points flaged by 'flag_name_rms_calc', + and errorbar of xi,eta is set from the fit covariance matrix. If 'force_one_redchi2', the errorbar of (xi,eta) is equivalent the case + if the error bar of each data point is set as the reduced chi-square is equal to unity. + - flag_name_rms_calc (str): + Name of the flag used for RMS calculation. Default is 'source'. + - flag_rms_calc_exclusive (bool): + Flag indicating whether the RMS calculation is exclusive to the flag. Default is True. + - save (bool): + Flag indicating whether to save the updated focal plane data. Default is True. Returns: - focal_plane (ResultSet): ResultSet containing updated xieta parameters for each detector. @@ -156,7 +157,7 @@ def update_xieta(tod, # If input focal_plane is a result with `force_zero_roll`, set the roll to be zero # Original value is stored to `roll_original` - if input_force_zero_roll: + if force_zero_roll: if 'roll_original' in tod.boresight._fields.keys(): pass else: @@ -188,6 +189,8 @@ def update_xieta(tod, else: mask_for_rms_calc = ~tod.flags[flag_name_rms_calc].mask() rms = np.ma.std(np.ma.masked_array(tod.signal, mask_for_rms_calc), axis=1).data + if 'rms' in tod._fields.keys(): + tod.move('rms', None) tod.wrap('rms', rms, [(0, 'dets')]) # use downsampled data for faster fitting @@ -220,38 +223,50 @@ def update_xieta(tod, xi_src, eta_src, _ = quat.decompose_xieta(q_total) xieta_src = np.array([xi_src, eta_src]) xieta_src = xieta_src[:, mask_di] - sig = sig_ds[di][mask_di] - + sig = sig_ds[di][mask_di] ptp_val = np.ptp(np.percentile(sig, [0.1, 99.9])) - if fit_func_name == '_gaussian2d': - p0 = (0., 0., fwhm_init_deg*coords.DEG, fwhm_init_deg*coords.DEG, 0., ptp_val) - fit_func = _gaussian2d - elif fit_func_name == '_gaussian2d_nonlin': - p0 = (0., 0., fwhm_init_deg*coords.DEG, fwhm_init_deg*coords.DEG, 0., ptp_val, -0.1,) - fit_func = _gaussian2d_nonlin - - popt, pcov = curve_fit(fit_func, xdata=xieta_src, ydata=sig, p0=p0, sigma=tod.rms[di]*np.ones_like(sig), - absolute_sigma=True, maxfev=int(1e5)) - - chi2 = np.sum(((fit_func(xieta_src, *popt) - sig)/tod.rms[di])**2) - redchi2 = chi2 / (np.prod(xieta_src.shape) - popt.shape[0]) - R2 = 1. - np.sum((fit_func(xieta_src, *popt) - sig)**2) / np.sum((sig - sig.mean())**2) - xi_opt, eta_opt = popt[0], popt[1] - - if error_estimation_method == 'rms_from_data': - xi_err, eta_err = np.sqrt(pcov[0,0]), np.sqrt(pcov[1,1]) - elif error_estimation_method == 'force_one_redchi2': - # The error of (xi, eta) is equivalent the case if the error bar of each data point is set - # as the reduced chi-square is equal to unity. - xi_err, eta_err = np.sqrt(pcov[0,0] * redchi2), np.sqrt(pcov[1,1] * redchi2) - redchi2 = 1. + if fit_func_name == 'gaussian2d_nonlin': + p0 = np.array([0., 0., fwhm_init_deg*coords.DEG, fwhm_init_deg*coords.DEG, 0., ptp_val]) + bounds = np.array( + [[-np.inf, -np.inf, 0.1*fwhm_init_deg*coords.DEG, 0.1*fwhm_init_deg*coords.DEG, -np.pi, 0.1*ptp_val], + [np.inf, np.inf, 10*fwhm_init_deg*coords.DEG, 10*fwhm_init_deg*coords.DEG, np.pi, 10*ptp_val]] + ) + if max_non_linear_order >= 2: + p0 = np.append(p0, np.zeros(max_non_linear_order-1)) + bounds = np.hstack([bounds, + np.vstack([[-np.inf * np.ones(max_non_linear_order-1), + np.inf * np.ones(max_non_linear_order-1)]]) + ]) + fit_func = wrapper_gaussian2d_nonlin else: - raise NameError("Unsupported name for 'error_estimation_method'") + raise NameError("Unsupported name for 'fit_func_name'") - xieta_det += np.array([xi_opt, eta_opt]) - xieta_dict[det] = {'xi': xieta_det[0], 'eta': xieta_det[1], 'xi_err': xi_err, 'eta_err': eta_err, - 'R2': R2, 'redchi2': redchi2} + try: + popt, pcov = curve_fit(fit_func, xdata=xieta_src, ydata=sig, sigma=tod.rms[di]*np.ones_like(sig), + p0=p0, bounds=bounds, absolute_sigma=True) + + chi2 = np.sum(((fit_func(xieta_src, *popt) - sig)/tod.rms[di])**2) + redchi2 = chi2 / (np.prod(xieta_src.shape) - popt.shape[0]) + R2 = 1. - np.sum((fit_func(xieta_src, *popt) - sig)**2) / np.sum((sig - sig.mean())**2) + xi_opt, eta_opt = popt[0], popt[1] + + if error_estimation_method == 'rms_from_data': + xi_err, eta_err = np.sqrt(pcov[0,0]), np.sqrt(pcov[1,1]) + elif error_estimation_method == 'force_one_redchi2': + # The error of (xi, eta) is equivalent the case if the error bar of each data point is set + # as the reduced chi-square is equal to unity. + xi_err, eta_err = np.sqrt(pcov[0,0] * redchi2), np.sqrt(pcov[1,1] * redchi2) + redchi2 = 1. + else: + raise NameError("Unsupported name for 'error_estimation_method'") + + xieta_det += np.array([xi_opt, eta_opt]) + xieta_dict[det] = {'xi': xieta_det[0], 'eta': xieta_det[1], 'xi_err': xi_err, 'eta_err': eta_err, + 'R2': R2, 'redchi2': redchi2} + except RuntimeError: + xieta_dict[det] = {'xi': np.nan, 'eta': np.nan, 'xi_err': np.nan, 'eta_err': np.nan, + 'R2': np.nan, 'redchi2': np.nan} focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', 'xi_err', 'eta_err', 'R2', 'redchi2']) for det in tod.dets.vals: @@ -259,37 +274,78 @@ def update_xieta(tod, xieta_dict[det]['xi_err'], xieta_dict[det]['eta_err'], xieta_dict[det]['R2'], xieta_dict[det]['redchi2'], )) - if save: - assert result_dir is not None - assert filename is not None - if not os.path.exists(result_dir): - os.makedirs(result_dir) - write_dataset(focal_plane, - filename=os.path.join(result_dir, filename), - address='focal_plane', - overwrite=True) + return focal_plane -def main(ctx_file, obs_id, wafer_slot, sso_name, result_dir, - ds_factor=10, fwhm = 1.*coords.DEG, restrict_dets=False): - ctx = core.Context(ctx_file) - meta = ctx.get_meta(obs_id) - meta.restrict('dets', meta.dets.vals[meta.det_info.wafer_slot == wafer_slot]) - if restrict_dets: - meta.restrict('dets', meta.dets.vals[:100]) - +def main(configs, obs_id, wafer_slot, sso_name=None, + fp_hdf_file=None, save_dir=None, restrict_dets_for_debug=False): + if type(configs) == str: + configs = yaml.safe_load(open(configs, "r")) + + # Derive parameters from config file + ctx = core.Context(configs.get('context_file')) + if fp_hdf_file is None: + fp_hdf_file = configs.get('fp_hdf_file', None) + if save_dir is None: + save_dir = configs.get('save_dir', None) + + # get sso_name if it is not specified + obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] + if sso_name is None: + known_source_names = ['moon', 'jupiter'] + for _source_name in known_source_names: + if _source_name in obs_tags: + sso_name = _source_name + if _source_name is None: + raise ValueError('sso_name is not specified') + + # construct pipeline from configs + pipe = Pipeline(configs["process_pipe"]) + for pipe_component in pipe: + if pipe_component.name == 'compute_source_flags': + pipe_component.process_cfgs['center_on'] = sso_name + + # Other parameters + force_zero_roll = configs.get('force_zero_roll') + ds_factor = configs.get('ds_factor') + mask_deg = configs.get('mask_deg') + fit_func_name = configs.get('fit_func_name') + max_non_linear_order = configs.get('max_non_linear_order') + fwhm_init_deg = configs.get('fwhm_init_deg') + error_estimation_method = configs.get('error_estimation_method') + flag_name_rms_calc = configls.get('flag_name_rms_calc') + flag_rms_calc_exclusive = configls.get('flag_rms_calc_exclusive') + + + # Load data logger.info('loading data') + meta = ctx.get_meta(obs_id, dets={'wafer_slot': wafer_slot}) + if restrict_dets_for_debug is not False: + meta.restrict('dets', meta.dets.vals[:restrict_dets_for_debug]) tod = ctx.get_obs(meta) - logger.info('tod processing') - tod_process(tod) - if not os.path.exists(result_dir): - logger.info(f'Make a directory: f{result_dir}') - os.makedirs(result_dir) + # get pointing + focal_plane_rset = update_xieta( tod, + sso_name=sso_name, + fp_hdf_file=fp_hdf_file, + force_zero_roll=force_zero_roll, + pipe=pipe, + ds_factor=ds_factor, + mask_deg=mask_deg, + fit_func_name = fit_func_name, + max_non_linear_order = max_non_linear_order, + fwhm_init_deg = fwhm_init_deg, + error_estimation_method=error_estimation_method, + flag_name_rms_calc = flag_name_rms_calc, + flag_rms_calc_exclusive = flag_rms_calc_exclusive, + ) - result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' - focal_plane_rset = update_xieta(tod=tod, sso_name=sso_name, ds_factor=ds_factor, fwhm=fwhm, - save=True, result_dir=result_dir, filename=result_filename) + os.makedirs(save_dir, exist_ok=True) + write_dataset(focal_plane_rset, + filename=os.path.join(save_dir, f'focal_plane_{obs_id}_{wafer_slot}.hdf'), + address='focal_plane', + overwrite=True) + return def get_parser(): @@ -297,11 +353,9 @@ def get_parser(): parser.add_argument("ctx_file", type=str, help="Path to the context file.") parser.add_argument("obs_id", type=str, help="Observation ID.") parser.add_argument("wafer_slot", type=int, help="Wafer slot number.") - parser.add_argument("sso_name", type=str, help="Name of the Solar System Object (SSO).") - parser.add_argument("result_dir", type=str, help="Directory to save the result.") - parser.add_argument("--ds_factor", type=int, default=10, help="Downsampling factor for TOD processing.") - parser.add_argument("--fwhm", type=float, default=1.0, help="Full width at half maximum of the Gaussian model.") - parser.add_argument("--restrict_dets", action="store_true", help="Flag to restrict the number of detectors.") + parser.add_argument("sso_name", type=str, default=None, help="Name of the Solar System Object (SSO).") + parser.add_argument("save_dir", type=str, help="Directory to save the result.") + parser.add_argument("restrict_dets_for_debug", action="store_true", help="Flag to restrict the number of detectors.") return parser if __name__ == '__main__': From d708edb7f980f51b2e5c0aab0400e8e14b37e031 Mon Sep 17 00:00:00 2001 From: Tomoki Terasaki Date: Tue, 7 May 2024 09:49:37 +0000 Subject: [PATCH 09/27] tod based fitting is automated --- sotodlib/coords/map_based_pointing.py | 35 ++++++- .../site_pipeline/make_mapbased_pointing.py | 60 ++++++++++-- sotodlib/site_pipeline/update_pointing.py | 97 ++++++++++++++----- 3 files changed, 159 insertions(+), 33 deletions(-) diff --git a/sotodlib/coords/map_based_pointing.py b/sotodlib/coords/map_based_pointing.py index 294a350a2..9f87ae406 100644 --- a/sotodlib/coords/map_based_pointing.py +++ b/sotodlib/coords/map_based_pointing.py @@ -50,7 +50,7 @@ def get_planet_trajectry(tod, planet, _split=20, return_model=False): q_planet = quat.rotation_lonlat(planet_az, planet_el) return q_planet -def get_wafer_centered_sight(tod, planet, q_planet=None, q_bs=None, q_wafer=None): +def get_wafer_centered_sight(tod=None, planet=None, q_planet=None, q_bs=None, q_wafer=None): """ Calculate the sightline vector from the focal plane, centered on the wafer, to a planet. @@ -131,6 +131,34 @@ def get_wafer_xieta(wafer_slot, optics_config_fn, xieta_bs_offset=(0., 0.), return xi_wafer, eta_wafer +def get_rough_hit_time(tod, wafer_slot, sso_name, circle_r_deg=7.,optics_config_fn=None): + """ + Estimate the rough hit time for a axismanager, wafer_slot, and sso_name. + + Parameters: + tod : An AxisManager object + wafer_slot (str): Identifier for the wafer slot. + sso_name (str): Name of the Solar System Object (e.g., 'moon', 'jupiter'). + circle_r_deg (float, optional): Radius in degrees defining the circular region around the wafer center. + Defaults to 7 degrees. + + Returns: + float: Estimated rough hit time within the circular region around the wafer center. + """ + q_bs = quat.rotation_lonlat(tod.boresight.az, tod.boresight.el) + q_planet = get_planet_trajectry(tod, sso_name) + xi_wafer, eta_wafer = get_wafer_xieta(wafer_slot, optics_config_fn=optics_config_fn, + roll_bs_offset=np.median(tod.boresight.roll), wrap_to_tod=False) + q_wafer = quat.rotation_xieta(xi_wafer, eta_wafer) + + q_wafer_centered = get_wafer_centered_sight(q_planet=q_planet, q_bs=q_bs, q_wafer=q_wafer) + x_to_z = ~quat.rotation_lonlat(0, 0) + xi_wafer_centered, eta_wafer_centered, _ = quat.decompose_xieta(x_to_z * q_wafer_centered) + r_wafer_centered = np.sqrt(xi_wafer_centered**2 + eta_wafer_centered**2) + hit_time = (tod.timestamps[-1] - tod.timestamps[0]) * np.mean(np.rad2deg(r_wafer_centered) < circle_r_deg) + return hit_time + + def make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf, xieta_bs_offset=(0., 0.), roll_bs_offset=None, signal='signal', wafer_mask_deg=8., res_deg=0.3, cuts=None,): @@ -151,10 +179,7 @@ def make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf, Returns: None - """ - - - + """ q_planet = get_planet_trajectry(tod, sso_name) q_bs = quat.rotation_lonlat(tod.boresight.az, tod.boresight.el) diff --git a/sotodlib/site_pipeline/make_mapbased_pointing.py b/sotodlib/site_pipeline/make_mapbased_pointing.py index 23f15416d..3adc0a401 100644 --- a/sotodlib/site_pipeline/make_mapbased_pointing.py +++ b/sotodlib/site_pipeline/make_mapbased_pointing.py @@ -14,8 +14,8 @@ from sotodlib.site_pipeline import util from sotodlib.preprocess import Pipeline logger = util.init_logger(__name__, 'make_map_based_pointing: ') - -def main(configs, obs_id, wafer_slot, + +def main_one_wafer(configs, obs_id, wafer_slot, sso_name=None, single_det_maps_dir=None, map_based_result_dir=None, tod_based_result_dir=None, tune_by_tod=None, restrict_dets_for_debug=False): @@ -82,10 +82,11 @@ def main(configs, obs_id, wafer_slot, xieta_bs_offset=xieta_bs_offset, wafer_mask_deg=wafer_mask_deg, res_deg=res_deg) + result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' # reconstruct pointing from single detector maps if save_normal_roll: logger.info(f'Saving map-based pointing results') - result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' + fp_rset_map_based = mbp.get_xieta_from_maps(map_hdf, save=True, output_dir=map_based_result_dir, filename=result_filename, @@ -146,19 +147,66 @@ def main(configs, obs_id, wafer_slot, filename=os.path.join(tod_based_result_dir_force_zero_roll, f'focal_plane_{obs_id}_{wafer_slot}.hdf'), address='focal_plane', overwrite=True) - - return +def main(configs, obs_id, wafer_slots, + sso_name=None, + single_det_maps_dir=None, map_based_result_dir=None, tod_based_result_dir=None, + tune_by_tod=None, hit_time_threshold=1200, hit_circle_r_deg=7.0, + restrict_dets_for_debug=False): + + logger.info('get wafer_slots which hit the source because wafer_slots are not specified') + if wafer_slots is None: + if type(configs) == str: + configs = yaml.safe_load(open(configs, "r")) + + ctx = core.Context(configs.get('context_file')) + optics_config_fn = configs.get('optics_config_fn') + + obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] + if sso_name is None: + known_source_names = ['moon', 'jupiter'] + for _source_name in known_source_names: + if _source_name in obs_tags: + sso_name = _source_name + if _source_name is None: + raise ValueError('sso_name is not specified') + + wafer_slots = [] + tod = ctx.get_obs(obs_id, dets=[]) + for ws in [f'ws{i}' for i in range(7)]: + hit_time = mbp.get_rough_hit_time(tod, wafer_slot=ws, sso_name=sso_name, circle_r_deg=hit_circle_r_deg, + optics_config_fn=optics_config_fn) + logger.info(f'hit_time for {ws} is {hit_time:.1f} [sec]') + if hit_time > hit_time_threshold: + wafer_slots.append(ws) + assert np.all(np.array(wafer_slots, dtype='U2') == 'ws') + + logger.info(f'wafer_slots which pointing calculated: {wafer_slots}') + for wafer_slot in wafer_slots: + main_one_wafer(configs=configs, + obs_id=obs_id, + wafer_slot=wafer_slot, + sso_name=sso_name, + single_det_maps_dir=single_det_maps_dir, + map_based_result_dir=map_based_result_dir, + tod_based_result_dir=tod_based_result_dir, + tune_by_tod=tune_by_tod, + restrict_dets_for_debug=restrict_dets_for_debug) + def get_parser(): parser = argparse.ArgumentParser(description="Process TOD data and update pointing") parser.add_argument("configs", type=str, help="Path to the configuration file") parser.add_argument("obs_id", type=str, help="Observation id") - parser.add_argument("wafer_slot", type=str, help="Wafer slot") + parser.add_argument("--wafer_slots", nargs='*', default=None, help="Wafer slots to be processed") parser.add_argument("--sso_name", type=str, default=None, help="Name of solar system object (e.g., 'moon', 'jupiter')") parser.add_argument("--single_det_maps_dir", type=str, default=None, help="Directory to save single detector maps") parser.add_argument("--map_based_result_dir", type=str, default=None, help="Directory to save map-based pointing results") parser.add_argument("--tod_based_result_dir", type=str, default=None, help="Directory to save TOD-based pointing results") + parser.add_argument("--hit_time_threshold", type=float, default=1200, + help="Minimum hit time. If calculated wafer hit time is smaller than that, pointing calculation for that wafer is skipped") + parser.add_argument("--hit_circle_r_deg", type=float, default=7., + help="circle radius for wafer hit time calculation") parser.add_argument("--restrict_dets_for_debug", type=int, default=False) return parser diff --git a/sotodlib/site_pipeline/update_pointing.py b/sotodlib/site_pipeline/update_pointing.py index 1d107c0e7..9c98c2f69 100644 --- a/sotodlib/site_pipeline/update_pointing.py +++ b/sotodlib/site_pipeline/update_pointing.py @@ -2,26 +2,21 @@ import numpy as np import yaml import h5py -import matplotlib.pyplot as plt +import argparse from tqdm import tqdm import scipy from scipy.optimize import curve_fit from sotodlib.core import metadata from sotodlib.io.metadata import read_dataset, write_dataset - +from sotodlib.coords import map_based_pointing as mbp from sotodlib import core from sotodlib import coords from sotodlib import tod_ops import so3g from so3g.proj import quat import sotodlib.coords.planets as planets - -from sotodlib.tod_ops import pca -from so3g.proj import Ranges, RangesMatrix -from pixell import enmap, enplot -from sotodlib.tod_ops.filters import high_pass_sine2, low_pass_sine2, fourier_filter - from sotodlib.site_pipeline import util +from sotodlib.preprocess import Pipeline logger = util.init_logger(__name__, 'update_pointing: ') def gaussian2d_nonlin(xieta, xi0, eta0, fwhm_xi, fwhm_eta, phi, a, nonlin_coeffs): @@ -109,9 +104,7 @@ def update_xieta(tod, Name of the Solar System Object (SSO). - fp_hdf_file (str or None): Path to the HDF file containing focal plane information. Default is None. - If None, tod.focal_plane is used for focal plane information. - - save_dir (str or None): - Directory where the updated data will be saved. Required if save is True. + If None, tod.focal_plane is used for focal plane information. - force_zero_roll (bool): Flag indicating whether to force the roll to be zero. Default is False. If True, input and output focal plane information assumes force_zero_roll condition. @@ -229,8 +222,8 @@ def update_xieta(tod, if fit_func_name == 'gaussian2d_nonlin': p0 = np.array([0., 0., fwhm_init_deg*coords.DEG, fwhm_init_deg*coords.DEG, 0., ptp_val]) bounds = np.array( - [[-np.inf, -np.inf, 0.1*fwhm_init_deg*coords.DEG, 0.1*fwhm_init_deg*coords.DEG, -np.pi, 0.1*ptp_val], - [np.inf, np.inf, 10*fwhm_init_deg*coords.DEG, 10*fwhm_init_deg*coords.DEG, np.pi, 10*ptp_val]] + [[-np.inf, -np.inf, fwhm_init_deg*coords.DEG/5., fwhm_init_deg*coords.DEG/5., -np.pi, 0.1*ptp_val], + [np.inf, np.inf, fwhm_init_deg*coords.DEG*5, fwhm_init_deg*coords.DEG*5, np.pi, 10*ptp_val]] ) if max_non_linear_order >= 2: p0 = np.append(p0, np.zeros(max_non_linear_order-1)) @@ -277,15 +270,28 @@ def update_xieta(tod, return focal_plane -def main(configs, obs_id, wafer_slot, sso_name=None, - fp_hdf_file=None, save_dir=None, restrict_dets_for_debug=False): +def main_one_wafer(configs, obs_id, wafer_slot, + sso_name=None, fp_hdf_file=None, fp_hdf_dir=None, + save_dir=None, restrict_dets_for_debug=False): if type(configs) == str: configs = yaml.safe_load(open(configs, "r")) # Derive parameters from config file ctx = core.Context(configs.get('context_file')) + + # get prior if fp_hdf_file is None: fp_hdf_file = configs.get('fp_hdf_file', None) + if fp_hdf_dir is None: + fp_hdf_dir = configs.get('fp_hdf_dir', None) + if fp_hdf_file is None: + if fp_hdf_dir is None: + pass + else: + fp_hdf_file = os.path.join(fp_hdf_dir, f'focal_plane_{obs_id}_{wafer_slot}.hdf') + if not os.path.exists(fp_hdf_file): + fp_hdf_file = None + if save_dir is None: save_dir = configs.get('save_dir', None) @@ -313,8 +319,8 @@ def main(configs, obs_id, wafer_slot, sso_name=None, max_non_linear_order = configs.get('max_non_linear_order') fwhm_init_deg = configs.get('fwhm_init_deg') error_estimation_method = configs.get('error_estimation_method') - flag_name_rms_calc = configls.get('flag_name_rms_calc') - flag_rms_calc_exclusive = configls.get('flag_rms_calc_exclusive') + flag_name_rms_calc = configs.get('flag_name_rms_calc') + flag_rms_calc_exclusive = configs.get('flag_rms_calc_exclusive') # Load data @@ -348,14 +354,61 @@ def main(configs, obs_id, wafer_slot, sso_name=None, return +def main(configs, obs_id, wafer_slots=None, + sso_name=None, fp_hdf_file=None, fp_hdf_dir=None, save_dir=None, + hit_time_threshold=1200, hit_circle_r_deg=7.0, + restrict_dets_for_debug=False): + logger.info('get wafer_slots which hit the source because wafer_slots are not specified') + if wafer_slots is None: + if type(configs) == str: + configs = yaml.safe_load(open(configs, "r")) + + ctx = core.Context(configs.get('context_file')) + optics_config_fn = configs.get('optics_config_fn') + + obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] + if sso_name is None: + known_source_names = ['moon', 'jupiter'] + for _source_name in known_source_names: + if _source_name in obs_tags: + sso_name = _source_name + if _source_name is None: + raise ValueError('sso_name is not specified') + + wafer_slots = [] + tod = ctx.get_obs(obs_id, dets=[]) + for ws in [f'ws{i}' for i in range(7)]: + hit_time = mbp.get_rough_hit_time(tod, wafer_slot=ws, sso_name=sso_name, circle_r_deg=hit_circle_r_deg, + optics_config_fn=optics_config_fn) + logger.info(f'hit_time for {ws} is {hit_time:.1f} [sec]') + if hit_time > hit_time_threshold: + wafer_slots.append(ws) + assert np.all(np.array(wafer_slots, dtype='U2') == 'ws') + + logger.info(f'wafer_slots which pointing calculated: {wafer_slots}') + for wafer_slot in wafer_slots: + main_one_wafer(configs, + obs_id, + wafer_slot, + sso_name=sso_name, + fp_hdf_file=fp_hdf_file, + fp_hdf_dir=fp_hdf_dir, + save_dir=save_dir, + restrict_dets_for_debug=restrict_dets_for_debug) + return + + def get_parser(): parser = argparse.ArgumentParser(description="Get updated result of pointings with tod-based results") - parser.add_argument("ctx_file", type=str, help="Path to the context file.") + parser.add_argument("configs", type=str, help="Path to the configuration file") parser.add_argument("obs_id", type=str, help="Observation ID.") - parser.add_argument("wafer_slot", type=int, help="Wafer slot number.") - parser.add_argument("sso_name", type=str, default=None, help="Name of the Solar System Object (SSO).") - parser.add_argument("save_dir", type=str, help="Directory to save the result.") - parser.add_argument("restrict_dets_for_debug", action="store_true", help="Flag to restrict the number of detectors.") + parser.add_argument("--wafer_slots", nargs='*', default=None, help="Wafer slots to be processed") + parser.add_argument("--sso_name", type=str, default=None, help="Name of the Solar System Object (SSO).") + parser.add_argument("--fp_hdf_file", type=str, default=None, help="File path to the focal_plane hdf file used as a prior") + parser.add_argument("--fp_hdf_dir", type=str, default=None, + help="Directory path where focal_plane hdf file of each observation are stored. Used only fp_hdf_file is not specified.") + parser.add_argument("--save_dir", type=str, help="Directory to save the result.") + parser.add_argument("--restrict_dets_for_debug", type=int, help="Flag to restrict the number of detectors.") return parser if __name__ == '__main__': From 34481e7108541aeef800e1a48b8c4c5f444ef3b7 Mon Sep 17 00:00:00 2001 From: tterasaki Date: Thu, 9 May 2024 07:45:20 +0000 Subject: [PATCH 10/27] main function is done --- sotodlib/coords/map_based_pointing.py | 8 +- .../site_pipeline/make_mapbased_pointing.py | 343 +++++++++++------- sotodlib/site_pipeline/update_pointing.py | 288 ++++++++++----- 3 files changed, 414 insertions(+), 225 deletions(-) diff --git a/sotodlib/coords/map_based_pointing.py b/sotodlib/coords/map_based_pointing.py index 9f87ae406..23e814713 100644 --- a/sotodlib/coords/map_based_pointing.py +++ b/sotodlib/coords/map_based_pointing.py @@ -477,12 +477,10 @@ def get_xieta_from_maps(map_hdf_file, filename = 'focal_plane_' + os.path.splitext(os.path.basename(map_hdf_file))[0] + '.hdf' output_file = os.path.join(output_dir, filename) - focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'band', 'channel', 'R2', 'xi', 'eta', 'gamma']) + focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', 'R2']) for det in dets: - band = int(det.split('_')[-2]) - channel = int(det.split('_')[-1]) - focal_plane.rows.append((det, band, channel, xieta_dict[det]['R2'], - xieta_dict[det]['xi'], xieta_dict[det]['eta'], 0.)) + focal_plane.rows.append((det, xieta_dict[det]['xi'], xieta_dict[det]['eta'], 0., + xieta_dict[det]['R2'])) write_dataset(focal_plane, output_file, 'focal_plane', overwrite=True) return focal_plane diff --git a/sotodlib/site_pipeline/make_mapbased_pointing.py b/sotodlib/site_pipeline/make_mapbased_pointing.py index 3adc0a401..faed09a43 100644 --- a/sotodlib/site_pipeline/make_mapbased_pointing.py +++ b/sotodlib/site_pipeline/make_mapbased_pointing.py @@ -2,66 +2,61 @@ import numpy as np import yaml import argparse +import time +import glob from sotodlib import core from sotodlib import coords from sotodlib import tod_ops -from sotodlib.tod_ops.filters import high_pass_sine2, low_pass_sine2, fourier_filter from sotodlib.coords import map_based_pointing as mbp from sotodlib.site_pipeline import update_pointing as up -from sotodlib.io.metadata import write_dataset +from sotodlib.io import metadata +from sotodlib.io.metadata import read_dataset, write_dataset from sotodlib.site_pipeline import util from sotodlib.preprocess import Pipeline logger = util.init_logger(__name__, 'make_map_based_pointing: ') -def main_one_wafer(configs, obs_id, wafer_slot, - sso_name=None, - single_det_maps_dir=None, map_based_result_dir=None, tod_based_result_dir=None, - tune_by_tod=None, restrict_dets_for_debug=False): +def _get_sso_names_from_tags(ctx, obs_id, candidate_names=['moon', 'jupiter']): + obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] + sso_names = [] + for _name in candidate_names: + if _name in obs_tags: + sso_names.append(_name) + if len(sso_names) == 0: + raise NameError('Could not find sso_name from observation tags') + else: + return sso_names +def main_one_wafer(configs, obs_id, wafer_slot, sso_name=None, + restrict_dets_for_debug=False): if type(configs) == str: configs = yaml.safe_load(open(configs, "r")) # Derive parameters from config file + # required parameters ctx = core.Context(configs.get('context_file')) - if single_det_maps_dir is None: - single_det_maps_dir = configs.get('single_det_maps_dir') - if map_based_result_dir is None: - map_based_result_dir = configs.get('map_based_result_dir') - if tod_based_result_dir is None: - tod_based_result_dir = configs.get('tod_based_result_dir') + single_det_maps_dir = configs.get('single_det_maps_dir') + result_dir = configs.get('result_dir') optics_config_fn = configs.get('optics_config_fn') + save_normal_roll = configs.get('save_normal_roll') + save_force_zero_roll = configs.get('save_force_zero_roll') + + # optional parameters xieta_bs_offset = configs.get('xieta_bs_offset', [0., 0.]) wafer_mask_deg = configs.get('wafer_mask_deg', 8.) res_deg = configs.get('res_deg', 0.3) edge_avoidance_deg = configs.get('edge_avoidance_deg', 0.3) - save_normal_roll = configs.get('save_normal_roll', True) - save_force_zero_roll = configs.get('save_force_zero_roll', True) - # parameters for tod tuning - tune_by_tod = configs.get('tune_by_tod') - if tune_by_tod: - tod_ds_factor = configs.get('tod_ds_factor') - tod_mask_deg = configs.get('tod_mask_deg') - tod_fit_func_name = configs.get('tod_fit_func_name') - tod_max_non_linear_order = configs.get('tod_max_non_linear_order') - tod_fwhm_init_deg = configs.get('tod_fwhm_init_deg') - tod_error_estimation_method = configs.get('tod_error_estimation_method') - tod_flag_name_rms_calc = configs.get('tod_flag_name_rms_calc') - tod_flag_rms_calc_exclusive = configs.get('tod_flag_rms_calc_exclusive') - - - # If sso_name is not specified, get sso name from observation tags - obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] if sso_name is None: - known_source_names = ['moon', 'jupiter'] - for _source_name in known_source_names: - if _source_name in obs_tags: - sso_name = _source_name - if _source_name is None: - raise ValueError('sso_name is not specified') - + logger.info('deriving sso_name from observation tag') + obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] + sso_names = _get_sso_names_from_tags(ctx, obs_id) + sso_name = sso_names[0] + if len(sso_names) >= 2: + logger.info(f'sso_names of {sso_names} are found from observation tags.' + + f'Processing only {sso_name}') + # Load data logger.info('loading data') meta = ctx.get_meta(obs_id, dets={'wafer_slot': wafer_slot}) @@ -88,125 +83,203 @@ def main_one_wafer(configs, obs_id, wafer_slot, logger.info(f'Saving map-based pointing results') fp_rset_map_based = mbp.get_xieta_from_maps(map_hdf, save=True, - output_dir=map_based_result_dir, - filename=result_filename, - force_zero_roll=False, - edge_avoidance = edge_avoidance_deg*coords.DEG) - - if tune_by_tod: - logger.info(f'Making tod-based pointing results') - up.wrap_fp_rset(tod, fp_rset_map_based) - fp_rset_tod_based = up.update_xieta( tod, - sso_name=sso_name, - fp_hdf_file=None, - force_zero_roll=False, - pipe=None, - ds_factor=tod_ds_factor, - mask_deg=tod_mask_deg, - fit_func_name = tod_fit_func_name, - max_non_linear_order = tod_max_non_linear_order, - fwhm_init_deg = tod_fwhm_init_deg, - error_estimation_method=tod_error_estimation_method, - flag_name_rms_calc = tod_flag_name_rms_calc, - flag_rms_calc_exclusive = tod_flag_rms_calc_exclusive, - ) - os.makedirs(tod_based_result_dir, exist_ok=True) - write_dataset(fp_rset_tod_based, - filename=os.path.join(tod_based_result_dir, f'focal_plane_{obs_id}_{wafer_slot}.hdf'), - address='focal_plane', - overwrite=True) + output_dir=result_dir, + filename=result_filename, + force_zero_roll=False, + edge_avoidance = edge_avoidance_deg*coords.DEG) if save_force_zero_roll: logger.info(f'Saving map-based pointing results (force-zero-roll)') - map_based_result_dir_force_zero_roll = map_based_result_dir + '_force_zero_roll' + result_dir_force_zero_roll = result_dir + '_force_zero_roll' fp_rset_map_based_force_zero_roll = mbp.get_xieta_from_maps(map_hdf, save=True, - output_dir=map_based_result_dir_force_zero_roll, + output_dir=result_dir_force_zero_roll, filename=result_filename, force_zero_roll=True, edge_avoidance = edge_avoidance_deg*coords.DEG) - if tune_by_tod: - logger.info(f'Making tod-based pointing results (force-zero-roll)') - up.wrap_fp_rset(tod, fp_rset_map_based_force_zero_roll) - tod_based_result_dir_force_zero_roll = tod_based_result_dir + '_force_zero_roll' - fp_rset_tod_based_force_zero_roll = up.update_xieta( tod, - sso_name=sso_name, - fp_hdf_file=None, - force_zero_roll=False, - pipe=None, - ds_factor=tod_ds_factor, - mask_deg=tod_mask_deg, - fit_func_name = tod_fit_func_name, - max_non_linear_order = tod_max_non_linear_order, - fwhm_init_deg = tod_fwhm_init_deg, - error_estimation_method=tod_error_estimation_method, - flag_name_rms_calc = tod_flag_name_rms_calc, - flag_rms_calc_exclusive = tod_flag_rms_calc_exclusive, - ) - os.makedirs(tod_based_result_dir_force_zero_roll, exist_ok=True) - write_dataset(fp_rset_tod_based_force_zero_roll, - filename=os.path.join(tod_based_result_dir_force_zero_roll, f'focal_plane_{obs_id}_{wafer_slot}.hdf'), - address='focal_plane', - overwrite=True) return -def main(configs, obs_id, wafer_slots, - sso_name=None, - single_det_maps_dir=None, map_based_result_dir=None, tod_based_result_dir=None, - tune_by_tod=None, hit_time_threshold=1200, hit_circle_r_deg=7.0, - restrict_dets_for_debug=False): - - logger.info('get wafer_slots which hit the source because wafer_slots are not specified') - if wafer_slots is None: - if type(configs) == str: - configs = yaml.safe_load(open(configs, "r")) - - ctx = core.Context(configs.get('context_file')) - optics_config_fn = configs.get('optics_config_fn') +def main_one_wafer_dummy(configs, obs_id, wafer_slot, restrict_dets_for_debug=False): + if type(configs) == str: + configs = yaml.safe_load(open(configs, "r")) + ctx = core.Context(configs.get('context_file')) + single_det_maps_dir = configs.get('single_det_maps_dir') + result_dir = configs.get('result_dir') + save_normal_roll = configs.get('save_normal_roll', True) + save_force_zero_roll = configs.get('save_force_zero_roll', True) + + meta = ctx.get_meta(obs_id, dets={'wafer_slot': wafer_slot}) + if restrict_dets_for_debug is not False: + meta.restrict('dets', meta.dets.vals[:restrict_dets_for_debug]) + result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' + + fp_rset_dummy_map_based = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', 'R2']) + for det in meta.dets.vals: + fp_rset_dummy_map_based.rows.append((det, np.nan, np.nan, np.nan, np.nan)) + if save_normal_roll: + os.makedirs(result_dir, exist_ok=True) + write_dataset(fp_rset_dummy_map_based, + filename=os.path.join(result_dir, result_filename), + address='focal_plane', + overwrite=True) + + if save_force_zero_roll: + result_dir_force_zero_roll = result_dir + '_force_zero_roll' + os.makedirs(result_dir_force_zero_roll, exist_ok=True) + write_dataset(fp_rset_dummy_map_based, + filename=os.path.join(result_dir_force_zero_roll, result_filename), + address='focal_plane', + overwrite=True) + return + +def combine_pointings(pointing_result_files): + combined_dict = {} + for file in pointing_result_files: + rset = read_dataset(file, 'focal_plane') + for row in rset[:]: + if row['dets:readout_id'] not in combined_dict.keys(): + combined_dict[row['dets:readout_id']] = {} + combined_dict[row['dets:readout_id']]['xi'] = row['xi'] + combined_dict[row['dets:readout_id']]['eta'] = row['eta'] + combined_dict[row['dets:readout_id']]['gamma'] = row['gamma'] + combined_dict[row['dets:readout_id']]['R2'] = row['R2'] + + focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', 'R2']) + + for det, val in combined_dict.items(): + focal_plane.rows.append((det, val['xi'], val['eta'], val['gamma'], val['R2'])) + return focal_plane + +def main_one_obs(configs, obs_id, sso_name=None, + restrict_dets_for_debug=False): + if type(configs) == str: + configs = yaml.safe_load(open(configs, "r")) + ctx = core.Context(configs.get('context_file')) + optics_config_fn = configs.get('optics_config_fn') + + result_dir = configs.get('result_dir') + save_normal_roll = configs.get('save_normal_roll') + save_force_zero_roll = configs.get('save_force_zero_roll') + + hit_time_threshold = configs.get('hit_time_threshold', 600) + hit_circle_r_deg = configs.get('hit_circle_r_deg', 7.0) + + if sso_name is None: + logger.info('deriving sso_name from observation tag') obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] - if sso_name is None: - known_source_names = ['moon', 'jupiter'] - for _source_name in known_source_names: - if _source_name in obs_tags: - sso_name = _source_name - if _source_name is None: - raise ValueError('sso_name is not specified') - - wafer_slots = [] - tod = ctx.get_obs(obs_id, dets=[]) - for ws in [f'ws{i}' for i in range(7)]: - hit_time = mbp.get_rough_hit_time(tod, wafer_slot=ws, sso_name=sso_name, circle_r_deg=hit_circle_r_deg, - optics_config_fn=optics_config_fn) - logger.info(f'hit_time for {ws} is {hit_time:.1f} [sec]') - if hit_time > hit_time_threshold: - wafer_slots.append(ws) - assert np.all(np.array(wafer_slots, dtype='U2') == 'ws') - - logger.info(f'wafer_slots which pointing calculated: {wafer_slots}') - for wafer_slot in wafer_slots: + sso_names = _get_sso_names_from_tags(ctx, obs_id) + sso_name = sso_names[0] + if len(sso_names) >= 2: + logger.info(f'sso_names of {sso_names} are found from observation tags.' + + f'Processing only {sso_name}') + + tod = ctx.get_obs(obs_id, dets=[]) + streamed_wafer_slots = ['ws{}'.format(index) for index, bit in enumerate(obs_id.split('_')[-1]) if bit == '1'] + processed_wafer_slots = [] + skipped_wafer_slots = [] + + for ws in streamed_wafer_slots: + hit_time = mbp.get_rough_hit_time(tod, wafer_slot=ws, sso_name=sso_name, circle_r_deg=hit_circle_r_deg, + optics_config_fn=optics_config_fn) + logger.info(f'hit_time for {ws} is {hit_time:.1f} [sec]') + if hit_time > hit_time_threshold: + processed_wafer_slots.append(ws) + else: + skipped_wafer_slots.append(ws) + + logger.info(f'wafer_slots which pointing calculated: {processed_wafer_slots}') + for wafer_slot in processed_wafer_slots: + logger.info(f'Processing {obs_id}, {wafer_slot}') main_one_wafer(configs=configs, obs_id=obs_id, wafer_slot=wafer_slot, sso_name=sso_name, - single_det_maps_dir=single_det_maps_dir, - map_based_result_dir=map_based_result_dir, - tod_based_result_dir=tod_based_result_dir, - tune_by_tod=tune_by_tod, restrict_dets_for_debug=restrict_dets_for_debug) + + logger.info(f'create dummy hdf for non-hitting wafer: {skipped_wafer_slots}') + for wafer_slot in skipped_wafer_slots: + main_one_wafer_dummy(configs=configs, + obs_id=obs_id, + wafer_slot=wafer_slot, + restrict_dets_for_debug=restrict_dets_for_debug) + + logger.info('making combined result') + if save_normal_roll: + pointing_result_files = glob.glob(os.path.join(result_dir, f'focal_plane_{obs_id}_ws[0-6].hdf')) + fp_rset_full = combine_pointings(pointing_result_files) + fp_rset_full_file = os.path.join(os.path.join(result_dir, f'focal_plane_{obs_id}_all.hdf')) + write_dataset(fp_rset_full, filename=fp_rset_full_file, + address='focal_plane', overwrite=True) + + + if save_force_zero_roll: + result_dir_force_zero_roll = result_dir + '_force_zero_roll' + pointing_result_files = glob.glob(os.path.join(result_dir_force_zero_roll, f'focal_plane_{obs_id}_ws[0-6].hdf')) + fp_rset_full = combine_pointings(pointing_result_files) + fp_rset_full_file = os.path.join(os.path.join(result_dir_force_zero_roll, f'focal_plane_{obs_id}_all.hdf')) + write_dataset(fp_rset_full, filename=fp_rset_full_file, + address='focal_plane', overwrite=True) + + + return + +def main(configs, min_ctime=None, max_ctime=None, update_delay=None, + obs_id=None, wafer_slot=None, sso_name=None, restrict_dets_for_debug=False): + if (min_ctime is None) and (update_delay is not None): + # If min_ctime is provided it will use that.. + # Otherwise it will use update_delay to set min_ctime. + min_ctime = int(time.time()) - update_delay*86400 + + if type(configs) == str: + configs = yaml.safe_load(open(configs, "r")) + ctx = core.Context(configs.get('context_file')) + + if obs_id is None: + query_text = configs.get('query_text', None) + query_tags = configs.get('query_tags', None) + tot_query = "and " + if query_text is not None: + tot_query += f"{query_text} and " + if min_ctime is not None: + tot_query += f"timestamp>={min_ctime} and " + if max_ctime is not None: + tot_query += f"timestamp<={max_ctime} and " + tot_query = tot_query[4:-4] + if tot_query == "": + tot_query = "1" + + logger.info(f'tot_query: {tot_query}') + obs_list= ctx.obsdb.query(tot_query, query_tags) + + for obs in obs_list: + obs_id = obs['obs_id'] + logger.info(f'Processing {obs_id}') + main_one_obs(configs=configs, obs_id=obs_id, + restrict_dets_for_debug=restrict_dets_for_debug) + + elif obs_id is not None: + if wafer_slot is None: + main_one_obs(configs=configs, obs_id=obs_id, sso_name=sso_name, + restrict_dets_for_debug=restrict_dets_for_debug) + else: + main_one_wafer(configs=configs, obs_id=obs_id, wafer_slot=wafer_slot, sso_name=sso_name, + restrict_dets_for_debug=restrict_dets_for_debug) def get_parser(): parser = argparse.ArgumentParser(description="Process TOD data and update pointing") parser.add_argument("configs", type=str, help="Path to the configuration file") - parser.add_argument("obs_id", type=str, help="Observation id") - parser.add_argument("--wafer_slots", nargs='*', default=None, help="Wafer slots to be processed") - parser.add_argument("--sso_name", type=str, default=None, help="Name of solar system object (e.g., 'moon', 'jupiter')") - parser.add_argument("--single_det_maps_dir", type=str, default=None, help="Directory to save single detector maps") - parser.add_argument("--map_based_result_dir", type=str, default=None, help="Directory to save map-based pointing results") - parser.add_argument("--tod_based_result_dir", type=str, default=None, help="Directory to save TOD-based pointing results") - parser.add_argument("--hit_time_threshold", type=float, default=1200, - help="Minimum hit time. If calculated wafer hit time is smaller than that, pointing calculation for that wafer is skipped") - parser.add_argument("--hit_circle_r_deg", type=float, default=7., - help="circle radius for wafer hit time calculation") + parser.add_argument('--min_ctime', type=int, help="Minimum timestamp for the beginning of an observation list") + parser.add_argument('--max_ctime', type=int, help="Maximum timestamp for the beginning of an observation list") + parser.add_argument('--update-delay', type=int, help="Number of days (unit is days) in the past to start observation list.") + parser.add_argument("--obs_id", type=str, + help="Specific observation obs_id to process. If provided, overrides other filtering parameters.") + + parser.add_argument("--wafer_slot", type=str, default=None, + help="Wafer slot to be processed (e.g., 'ws0', 'ws3'). Valid only when obs_id is specified.") + + parser.add_argument("--sso_name", type=str, default=None, + help="Name of solar system object (e.g., 'moon', 'jupiter'). If not specified, get sso_name from observation tags. "\ + + "Valid only when obs_id is specified") parser.add_argument("--restrict_dets_for_debug", type=int, default=False) return parser diff --git a/sotodlib/site_pipeline/update_pointing.py b/sotodlib/site_pipeline/update_pointing.py index 9c98c2f69..7f1c73b31 100644 --- a/sotodlib/site_pipeline/update_pointing.py +++ b/sotodlib/site_pipeline/update_pointing.py @@ -3,8 +3,10 @@ import yaml import h5py import argparse +import time +import glob from tqdm import tqdm -import scipy + from scipy.optimize import curve_fit from sotodlib.core import metadata from sotodlib.io.metadata import read_dataset, write_dataset @@ -19,6 +21,17 @@ from sotodlib.preprocess import Pipeline logger = util.init_logger(__name__, 'update_pointing: ') +def _get_sso_names_from_tags(ctx, obs_id, candidate_names=['moon', 'jupiter']): + obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] + sso_names = [] + for _name in candidate_names: + if _name in obs_tags: + sso_names.append(_name) + if len(sso_names) == 0: + raise NameError('Could not find sso_name from observation tags') + else: + return sso_names + def gaussian2d_nonlin(xieta, xi0, eta0, fwhm_xi, fwhm_eta, phi, a, nonlin_coeffs): """ An Gaussian beam model with non-linear response Args @@ -80,20 +93,19 @@ def wrap_fp_from_hdf(tod, fp_hdf_file, data_set='focal_plane'): def update_xieta(tod, - sso_name='moon', + sso_name=None, fp_hdf_file=None, - save_dir=None, - pipe=None, force_zero_roll=False, + pipe=None, ds_factor=10, mask_deg=3, fit_func_name = 'gaussian2d_nonlin', max_non_linear_order = 1, fwhm_init_deg = 0.5, - error_estimation_method='force_one_redchi2', # rms_from_data + error_estimation_method='force_one_redchi2', flag_name_rms_calc = 'source', flag_rms_calc_exclusive = True, - save=True, ): + ): """ Update xieta parameters for each detector by TOD fitting of a point source observation. @@ -123,7 +135,7 @@ def update_xieta(tod, - fwhm_init_deg (float): Initial guess for full width at half maximum in degrees. Default is 0.5. - error_estimation_method (str): - Method for error estimation. Default is 'rms_from_data'. 'rms_from_data' and 'force_one_redchi2' are supported. + Method for error estimation. Default is 'force_one_redchi2'. 'force_one_redchi2' and 'rms_from_data' are supported. If 'rms_from_data', errorbar of each data point is set by root-mean-square of the data points flaged by 'flag_name_rms_calc', and errorbar of xi,eta is set from the fit covariance matrix. If 'force_one_redchi2', the errorbar of (xi,eta) is equivalent the case if the error bar of each data point is set as the reduced chi-square is equal to unity. @@ -131,8 +143,6 @@ def update_xieta(tod, Name of the flag used for RMS calculation. Default is 'source'. - flag_rms_calc_exclusive (bool): Flag indicating whether the RMS calculation is exclusive to the flag. Default is True. - - save (bool): - Flag indicating whether to save the updated focal plane data. Default is True. Returns: - focal_plane (ResultSet): ResultSet containing updated xieta parameters for each detector. @@ -270,9 +280,8 @@ def update_xieta(tod, return focal_plane -def main_one_wafer(configs, obs_id, wafer_slot, - sso_name=None, fp_hdf_file=None, fp_hdf_dir=None, - save_dir=None, restrict_dets_for_debug=False): +def main_one_wafer(configs, obs_id, wafer_slot, sso_name=None, + restrict_dets_for_debug=False): if type(configs) == str: configs = yaml.safe_load(open(configs, "r")) @@ -280,47 +289,45 @@ def main_one_wafer(configs, obs_id, wafer_slot, ctx = core.Context(configs.get('context_file')) # get prior + fp_hdf_file = configs.get('fp_hdf_file', None) + fp_hdf_dir = configs.get('fp_hdf_dir', None) if fp_hdf_file is None: - fp_hdf_file = configs.get('fp_hdf_file', None) - if fp_hdf_dir is None: - fp_hdf_dir = configs.get('fp_hdf_dir', None) - if fp_hdf_file is None: - if fp_hdf_dir is None: - pass - else: + if fp_hdf_dir is not None: fp_hdf_file = os.path.join(fp_hdf_dir, f'focal_plane_{obs_id}_{wafer_slot}.hdf') if not os.path.exists(fp_hdf_file): fp_hdf_file = None - - if save_dir is None: - save_dir = configs.get('save_dir', None) + + result_dir = configs.get('result_dir') + force_zero_roll = configs.get('force_zero_roll', True) + if force_zero_roll: + result_dir = result_dir + '_force_zero_roll' # get sso_name if it is not specified - obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] if sso_name is None: - known_source_names = ['moon', 'jupiter'] - for _source_name in known_source_names: - if _source_name in obs_tags: - sso_name = _source_name - if _source_name is None: - raise ValueError('sso_name is not specified') + logger.info('deriving sso_name from observation tag') + obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] + sso_names = _get_sso_names_from_tags(ctx, obs_id) + sso_name = sso_names[0] + if len(sso_names) >= 2: + logger.info(f'sso_names of {sso_names} are found from observation tags.' + + f'Processing only {sso_name}') # construct pipeline from configs - pipe = Pipeline(configs["process_pipe"]) + pipe = Pipeline(configs["process_pipe"], logger=logger) for pipe_component in pipe: if pipe_component.name == 'compute_source_flags': pipe_component.process_cfgs['center_on'] = sso_name # Other parameters force_zero_roll = configs.get('force_zero_roll') - ds_factor = configs.get('ds_factor') - mask_deg = configs.get('mask_deg') - fit_func_name = configs.get('fit_func_name') - max_non_linear_order = configs.get('max_non_linear_order') - fwhm_init_deg = configs.get('fwhm_init_deg') - error_estimation_method = configs.get('error_estimation_method') - flag_name_rms_calc = configs.get('flag_name_rms_calc') - flag_rms_calc_exclusive = configs.get('flag_rms_calc_exclusive') + ds_factor = configs.get('ds_factor', 20) + mask_deg = configs.get('mask_deg', 3.0) + fit_func_name = configs.get('fit_func_name', 'gaussian2d_nonlin') + max_non_linear_order = configs.get('max_non_linear_order', 2) + fwhm_init_deg = configs.get('fwhm_init_deg', 0.5) + error_estimation_method = configs.get('error_estimation_method', 'force_one_redchi2') + flag_name_rms_calc = configs.get('flag_name_rms_calc', 'source') + flag_rms_calc_exclusive = configs.get('flag_rms_calc_exclusive', True) # Load data @@ -346,69 +353,180 @@ def main_one_wafer(configs, obs_id, wafer_slot, flag_rms_calc_exclusive = flag_rms_calc_exclusive, ) - os.makedirs(save_dir, exist_ok=True) + os.makedirs(result_dir, exist_ok=True) write_dataset(focal_plane_rset, - filename=os.path.join(save_dir, f'focal_plane_{obs_id}_{wafer_slot}.hdf'), + filename=os.path.join(result_dir, f'focal_plane_{obs_id}_{wafer_slot}.hdf'), address='focal_plane', overwrite=True) return -def main(configs, obs_id, wafer_slots=None, - sso_name=None, fp_hdf_file=None, fp_hdf_dir=None, save_dir=None, - hit_time_threshold=1200, hit_circle_r_deg=7.0, - restrict_dets_for_debug=False): - logger.info('get wafer_slots which hit the source because wafer_slots are not specified') - if wafer_slots is None: - if type(configs) == str: - configs = yaml.safe_load(open(configs, "r")) - - ctx = core.Context(configs.get('context_file')) - optics_config_fn = configs.get('optics_config_fn') +def main_one_wafer_dummy(configs, obs_id, wafer_slot, restrict_dets_for_debug=False): + if type(configs) == str: + configs = yaml.safe_load(open(configs, "r")) + ctx = core.Context(configs.get('context_file')) + result_dir = configs.get('result_dir') + force_zero_roll = configs.get('force_zero_roll', True) + if force_zero_roll: + result_dir = result_dir + '_force_zero_roll' + + meta = ctx.get_meta(obs_id, dets={'wafer_slot': wafer_slot}) + if restrict_dets_for_debug is not False: + meta.restrict('dets', meta.dets.vals[:restrict_dets_for_debug]) + result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' + + fp_rset_dummy = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', + 'xi_err', 'eta_err', 'R2', 'redchi2']) + for det in meta.dets.vals: + fp_rset_dummy.rows.append((det, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan)) + os.makedirs(result_dir, exist_ok=True) + write_dataset(fp_rset_dummy, + filename=os.path.join(result_dir, result_filename), + address='focal_plane', + overwrite=True) + return + +def combine_pointings(pointing_result_files): + combined_dict = {} + for file in pointing_result_files: + rset = read_dataset(file, 'focal_plane') + for row in rset[:]: + if row['dets:readout_id'] not in combined_dict.keys(): + combined_dict[row['dets:readout_id']] = {} + combined_dict[row['dets:readout_id']]['xi'] = row['xi'] + combined_dict[row['dets:readout_id']]['eta'] = row['eta'] + combined_dict[row['dets:readout_id']]['gamma'] = row['gamma'] + combined_dict[row['dets:readout_id']]['xi_err'] = row['xi_err'] + combined_dict[row['dets:readout_id']]['eta_err'] = row['eta_err'] + combined_dict[row['dets:readout_id']]['R2'] = row['R2'] + combined_dict[row['dets:readout_id']]['redchi2'] = row['redchi2'] + + focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', 'xi_err', 'eta_err', 'R2', 'redchi2']) + + for det, val in combined_dict.items(): + focal_plane.rows.append((det, val['xi'], val['eta'], val['gamma'], val['xi_err'], val['eta_err'], val['R2'], val['redchi2'])) + return focal_plane + +def main_one_obs(configs, obs_id, sso_name=None, + restrict_dets_for_debug=False): + if type(configs) == str: + configs = yaml.safe_load(open(configs, "r")) + ctx = core.Context(configs.get('context_file')) + result_dir = configs.get('result_dir') + force_zero_roll = configs.get('force_zero_roll', True) + if force_zero_roll: + result_dir = result_dir + '_force_zero_roll' + optics_config_fn = configs.get('optics_config_fn') + + hit_time_threshold = configs.get('hit_time_threshold', 600) + hit_circle_r_deg = configs.get('hit_circle_r_deg', 7.0) + + if sso_name is None: + logger.info('deriving sso_name from observation tag') obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] - if sso_name is None: - known_source_names = ['moon', 'jupiter'] - for _source_name in known_source_names: - if _source_name in obs_tags: - sso_name = _source_name - if _source_name is None: - raise ValueError('sso_name is not specified') - - wafer_slots = [] - tod = ctx.get_obs(obs_id, dets=[]) - for ws in [f'ws{i}' for i in range(7)]: - hit_time = mbp.get_rough_hit_time(tod, wafer_slot=ws, sso_name=sso_name, circle_r_deg=hit_circle_r_deg, - optics_config_fn=optics_config_fn) - logger.info(f'hit_time for {ws} is {hit_time:.1f} [sec]') - if hit_time > hit_time_threshold: - wafer_slots.append(ws) - assert np.all(np.array(wafer_slots, dtype='U2') == 'ws') - - logger.info(f'wafer_slots which pointing calculated: {wafer_slots}') - for wafer_slot in wafer_slots: - main_one_wafer(configs, - obs_id, - wafer_slot, + sso_names = _get_sso_names_from_tags(ctx, obs_id) + sso_name = sso_names[0] + if len(sso_names) >= 2: + logger.info(f'sso_names of {sso_names} are found from observation tags.' + + f'Processing only {sso_name}') + + tod = ctx.get_obs(obs_id, dets=[]) + streamed_wafer_slots = ['ws{}'.format(index) for index, bit in enumerate(obs_id.split('_')[-1]) if bit == '1'] + processed_wafer_slots = [] + skipped_wafer_slots = [] + for ws in streamed_wafer_slots: + hit_time = mbp.get_rough_hit_time(tod, wafer_slot=ws, sso_name=sso_name, circle_r_deg=hit_circle_r_deg, + optics_config_fn=optics_config_fn) + logger.info(f'hit_time for {ws} is {hit_time:.1f} [sec]') + if hit_time > hit_time_threshold: + processed_wafer_slots.append(ws) + else: + skipped_wafer_slots.append(ws) + + logger.info(f'wafer_slots which pointing calculated: {processed_wafer_slots}') + for wafer_slot in processed_wafer_slots: + logger.info(f'Processing {obs_id}, {wafer_slot}') + main_one_wafer(configs=configs, + obs_id=obs_id, + wafer_slot=wafer_slot, sso_name=sso_name, - fp_hdf_file=fp_hdf_file, - fp_hdf_dir=fp_hdf_dir, - save_dir=save_dir, restrict_dets_for_debug=restrict_dets_for_debug) + + logger.info(f'create dummy hdf for non-hitting wafer: {skipped_wafer_slots}') + for wafer_slot in skipped_wafer_slots: + main_one_wafer_dummy(configs=configs, + obs_id=obs_id, + wafer_slot=wafer_slot, + restrict_dets_for_debug=restrict_dets_for_debug) + + logger.info('making combined result') + pointing_result_files = glob.glob(os.path.join(result_dir, f'focal_plane_{obs_id}_ws[0-6].hdf')) + fp_rset_full = combine_pointings(pointing_result_files) + fp_rset_full_file = os.path.join(os.path.join(result_dir, f'focal_plane_{obs_id}_all.hdf')) + write_dataset(fp_rset_full, filename=fp_rset_full_file, + address='focal_plane', overwrite=True) + +def main(configs, min_ctime=None, max_ctime=None, update_delay=None, + obs_id=None, wafer_slot=None, sso_name=None, restrict_dets_for_debug=False): + if (min_ctime is None) and (update_delay is not None): + # If min_ctime is provided it will use that.. + # Otherwise it will use update_delay to set min_ctime. + min_ctime = int(time.time()) - update_delay*86400 + + if type(configs) == str: + configs = yaml.safe_load(open(configs, "r")) + ctx = core.Context(configs.get('context_file')) + + if obs_id is None: + query_text = configs.get('query_text', None) + query_tags = configs.get('query_tags', None) + tot_query = "and " + if query_text is not None: + tot_query += f"{query_text} and " + if min_ctime is not None: + tot_query += f"timestamp>={min_ctime} and " + if max_ctime is not None: + tot_query += f"timestamp<={max_ctime} and " + tot_query = tot_query[4:-4] + if tot_query == "": + tot_query = "1" + + logger.info(f'tot_query: {tot_query}') + obs_list= ctx.obsdb.query(tot_query, query_tags) + + for obs in obs_list: + obs_id = obs['obs_id'] + logger.info(f'Processing {obs_id}') + main_one_obs(configs=configs, obs_id=obs_id, + restrict_dets_for_debug=restrict_dets_for_debug) + + elif obs_id is not None: + if wafer_slot is None: + main_one_obs(configs=configs, obs_id=obs_id, sso_name=sso_name, + restrict_dets_for_debug=restrict_dets_for_debug) + else: + main_one_wafer(configs=configs, obs_id=obs_id, wafer_slot=wafer_slot, sso_name=sso_name, + restrict_dets_for_debug=restrict_dets_for_debug) return def get_parser(): parser = argparse.ArgumentParser(description="Get updated result of pointings with tod-based results") parser.add_argument("configs", type=str, help="Path to the configuration file") - parser.add_argument("obs_id", type=str, help="Observation ID.") - parser.add_argument("--wafer_slots", nargs='*', default=None, help="Wafer slots to be processed") - parser.add_argument("--sso_name", type=str, default=None, help="Name of the Solar System Object (SSO).") - parser.add_argument("--fp_hdf_file", type=str, default=None, help="File path to the focal_plane hdf file used as a prior") - parser.add_argument("--fp_hdf_dir", type=str, default=None, - help="Directory path where focal_plane hdf file of each observation are stored. Used only fp_hdf_file is not specified.") - parser.add_argument("--save_dir", type=str, help="Directory to save the result.") - parser.add_argument("--restrict_dets_for_debug", type=int, help="Flag to restrict the number of detectors.") + parser.add_argument('--min_ctime', type=int, help="Minimum timestamp for the beginning of an observation list") + parser.add_argument('--max_ctime', type=int, help="Maximum timestamp for the beginning of an observation list") + parser.add_argument('--update-delay', type=int, help="Number of days (unit is days) in the past to start observation list.") + parser.add_argument("--obs_id", type=str, + help="Specific observation obs_id to process. If provided, overrides other filtering parameters.") + + parser.add_argument("--wafer_slot", type=str, default=None, + help="Wafer slot to be processed (e.g., 'ws0', 'ws3'). Valid only when obs_id is specified.") + + parser.add_argument("--sso_name", type=str, default=None, + help="Name of solar system object (e.g., 'moon', 'jupiter'). If not specified, get sso_name from observation tags. "\ + + "Valid only when obs_id is specified") + parser.add_argument("--restrict_dets_for_debug", type=int, default=False) return parser if __name__ == '__main__': From 2da796509a015182df4d0f33f46f9992c1bd859f Mon Sep 17 00:00:00 2001 From: tterasaki Date: Fri, 10 May 2024 05:11:48 +0000 Subject: [PATCH 11/27] removed combine_focal_plane.py which is no longer used --- .../site_pipeline/combine_focal_planes.py | 130 ------------------ 1 file changed, 130 deletions(-) delete mode 100644 sotodlib/site_pipeline/combine_focal_planes.py diff --git a/sotodlib/site_pipeline/combine_focal_planes.py b/sotodlib/site_pipeline/combine_focal_planes.py deleted file mode 100644 index b8a1184a4..000000000 --- a/sotodlib/site_pipeline/combine_focal_planes.py +++ /dev/null @@ -1,130 +0,0 @@ -import os -import re -import glob -import numpy as np - -from sotodlib.core import metadata -from sotodlib.io.metadata import write_dataset, read_dataset - -from sotodlib.site_pipeline import util -logger = util.init_logger(__name__, 'combine_focal_planes: ') - -def combine_pointings(pointing_result_files, method='highest_R2', R2_threshold=0.3, - save=False, output_dir=None, save_name=None): - combined_dict = {} - for file in pointing_result_files: - rset = read_dataset(file, 'focal_plane') - for row in rset[:]: - if row['dets:readout_id'] not in combined_dict.keys(): - combined_dict[row['dets:readout_id']] = {} - combined_dict[row['dets:readout_id']]['band'] = row['band'] - combined_dict[row['dets:readout_id']]['channel'] = row['channel'] - - combined_dict[row['dets:readout_id']]['R2'] = np.atleast_1d([]) - combined_dict[row['dets:readout_id']]['xi'] = np.atleast_1d([]) - combined_dict[row['dets:readout_id']]['eta'] = np.atleast_1d([]) - combined_dict[row['dets:readout_id']]['gamma'] = np.atleast_1d([]) - - combined_dict[row['dets:readout_id']]['R2'] = np.append(combined_dict[row['dets:readout_id']]['R2'], row['R2']) - combined_dict[row['dets:readout_id']]['xi'] = np.append(combined_dict[row['dets:readout_id']]['xi'], row['xi']) - combined_dict[row['dets:readout_id']]['eta'] = np.append(combined_dict[row['dets:readout_id']]['eta'], row['eta']) - combined_dict[row['dets:readout_id']]['gamma'] = np.append(combined_dict[row['dets:readout_id']]['gamma'], row['gamma']) - - focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'band', 'channel', 'R2', 'xi', 'eta', 'gamma']) - for det, val in combined_dict.items(): - band = int(val['band']) - channel = int(val['channel']) - - mask = val['R2'] > R2_threshold - if np.all(~mask): - xi, eta, gamma, R2 = np.nan, np.nan, np.nan, np.nan - else: - if method == 'highest_R2': - idx = np.argmax(val['R2'][mask]) - xi, eta, gamma, R2 = val['xi'][mask][idx], val['eta'][mask][idx], val['gamma'][mask][idx], val['R2'][mask][idx] - elif method == 'mean': - xi, eta, gamma = np.mean(val['xi'][mask]), np.mean(val['eta'][mask]), np.mean(val['gamma'][mask]) - R2 = np.nan - elif method == 'median': - xi, eta, gamma = np.median(val['xi'][mask]), np.median(val['eta'][mask]), np.median(val['gamma'][mask]) - R2 = np.nan - else: - raise ValueError('Not supported method. Supported methods are `highest_R2`, `mean` or `median`') - focal_plane.rows.append((det, band, channel, R2, xi, eta, gamma)) - if save: - if output_dir is None: - output_dir = os.path.join(os.getcwd(), 'combined_pointing_results') - if not os.path.exists(output_dir): - os.makedirs(output_dir) - if save_name is None: - ctimes = np.atleast_1d([]) - wafer_slots = np.atleast_1d([]) - for file in pointing_result_files: - filename = os.path.basename(file) - match = re.search('\d{10}', filename) - ctime = int(match.group(0) if match else None) - match = re.search('ws\d{1}', filename) - ws = match.group(0) - ctimes = np.append(ctimes, ctime) - wafer_slots = np.append(wafer_slots, ws) - ctimes = ctimes.astype('int') - wafer_slots = np.sort(np.unique(wafer_slots.astype('U3'))) - save_name = f'focal_plane_{ctimes.min()}_{ctimes.max()}_' + ''.join(wafer_slots) + '.hdf' - - write_dataset(focal_plane, os.path.join(output_dir, save_name), 'focal_plane', overwrite=True) - return focal_plane - -def combine_onewafer_results(pointing_dir, ws, output_dir, filename=None, - method='highest_R2', R2_threshold=0.3,): - pointing_result_files = glob.glob(os.path.join(pointing_dir, f'focal_plane*{ws}.hdf')) - if filename is None: - filename = f'focal_plane_{ws}_combined.hdf' - _ = combine_pointings(pointing_result_files, save=True, output_dir=output_dir, save_name=filename) - return - -def combine_allwafer_results(pointing_dir, output_dir, filename=None, - method='highest_R2', R2_threshold=0.3,): - pointing_result_files = glob.glob(os.path.join(pointing_dir, 'focal_plane*.hdf')) - if filename is None: - filename = f'focal_plane_combined.hdf' - _ = combine_pointings(pointing_result_files, save=True, output_dir=output_dir, save_name=filename) - return - -def make_detabase(focal_plane_file, db_file,): - scheme = metadata.ManifestScheme().add_data_field('dataset') - db = metadata.ManifestDb(scheme=scheme) - db.add_entry({'dataset': 'focal_plane'}, filename=focal_plane_file) - db.to_file(db_file) - return - -def main(pointing_dir, output_dir=None, method='highest_R2', R2_threshold=0.3,): - if output_dir is None: - output_dir = os.path.join(os.getcwd(), 'combined_results') - - logger.info('Combining each wafer resluts') - wafer_slots = [f'ws{i}' for i in range(7)] - for ws in wafer_slots: - combine_onewafer_results(pointing_dir=pointing_dir, ws=ws, - output_dir=output_dir, filename=None, - method=method, R2_threshold=R2_threshold) - - logger.info('Combining all wafer resluts') - combine_allwafer_results(pointing_dir=pointing_dir, output_dir=output_dir, filename='focal_plane_combined.hdf', - method=method, R2_threshold=R2_threshold) - - logger.info('Making a database') - focal_plane_file = os.path.join(output_dir, 'focal_plane_combined.hdf') - db_file = os.path.join(output_dir, 'focal_plane_combined.sqlite') - make_detabase(focal_plane_file, db_file,) - return - -def get_parser(): - parser = argparse.ArgumentParser(description="Combine multiple result of pointing.") - parser.add_argument('--pointing_dir', type=str, required=True, help='Directory containing pointing result files.') - parser.add_argument('--output_dir', type=str, default=None, help='Directory to save combined results. Default is "combined_results".') - parser.add_argument('--method', type=str, default='highest_R2', choices=['highest_R2', 'mean', 'median'], help='Combination method. Default is "highest_R2".') - parser.add_argument('--R2_threshold', type=float, default=0.3, help='Threshold for R2 value. Default is 0.3.') - return parser - -if __name__ == '__main__': - util.main_launcher(main, get_parser) From 8d8fae46b9c7b03368bd0bbdbb288b5ed026ad2c Mon Sep 17 00:00:00 2001 From: tterasaki Date: Thu, 30 May 2024 09:04:47 +0000 Subject: [PATCH 12/27] removed planets.close() --- sotodlib/coords/map_based_pointing.py | 1 + sotodlib/coords/planets.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/sotodlib/coords/map_based_pointing.py b/sotodlib/coords/map_based_pointing.py index 23e814713..43305b37c 100644 --- a/sotodlib/coords/map_based_pointing.py +++ b/sotodlib/coords/map_based_pointing.py @@ -32,6 +32,7 @@ def get_planet_trajectry(tod, planet, _split=20, return_model=False): If return_model is False: array: Array of quaternions representing trajectry of the planet at each timestamp. """ + print(planet) timestamps_sparse = np.linspace(tod.timestamps[0], tod.timestamps[-1], _split) planet_az_sparse = np.zeros_like(timestamps_sparse) diff --git a/sotodlib/coords/planets.py b/sotodlib/coords/planets.py index 322cd7e66..0c96c40d2 100644 --- a/sotodlib/coords/planets.py +++ b/sotodlib/coords/planets.py @@ -303,7 +303,7 @@ def _get_astrometric(source_name, timestamp, site='_default'): sf_timestamp = timescale.from_datetime( datetime.datetime.fromtimestamp(timestamp, tz=skyfield_api.utc)) astrometric = observatory.at(sf_timestamp).observe(target) - planets.close() + #planets.close() return astrometric From 59f55ce515db4315b9b852b5a7beb6bbe4ac206e Mon Sep 17 00:00:00 2001 From: tterasaki Date: Fri, 31 May 2024 04:57:22 +0000 Subject: [PATCH 13/27] tuning for master branch --- sotodlib/preprocess/processes.py | 3 +++ sotodlib/site_pipeline/update_pointing.py | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/sotodlib/preprocess/processes.py b/sotodlib/preprocess/processes.py index f10cd4d49..9d3d72ceb 100644 --- a/sotodlib/preprocess/processes.py +++ b/sotodlib/preprocess/processes.py @@ -845,6 +845,9 @@ class SourceFlags(_Preprocess): """ name = "source_flags" + def process(self, aman, proc_aman): + tod_ops.flags.get_source_flags(aman, **self.process_cfgs) + def calc_and_save(self, aman, proc_aman): center_on = self.calc_cfgs.get('center_on', 'planet') # Get source from tags diff --git a/sotodlib/site_pipeline/update_pointing.py b/sotodlib/site_pipeline/update_pointing.py index 7f1c73b31..9af4fc96f 100644 --- a/sotodlib/site_pipeline/update_pointing.py +++ b/sotodlib/site_pipeline/update_pointing.py @@ -315,7 +315,7 @@ def main_one_wafer(configs, obs_id, wafer_slot, sso_name=None, # construct pipeline from configs pipe = Pipeline(configs["process_pipe"], logger=logger) for pipe_component in pipe: - if pipe_component.name == 'compute_source_flags': + if pipe_component.name == 'source_flags': pipe_component.process_cfgs['center_on'] = sso_name # Other parameters From 5b80e503cec9c8fa68d6cddf50d7d8076885fe13 Mon Sep 17 00:00:00 2001 From: tterasaki Date: Fri, 31 May 2024 06:21:39 +0000 Subject: [PATCH 14/27] fix preprocessing stuff --- sotodlib/preprocess/processes.py | 3 --- sotodlib/site_pipeline/update_pointing.py | 2 +- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/sotodlib/preprocess/processes.py b/sotodlib/preprocess/processes.py index 7998ac517..0b236c782 100644 --- a/sotodlib/preprocess/processes.py +++ b/sotodlib/preprocess/processes.py @@ -844,9 +844,6 @@ class SourceFlags(_Preprocess): .. autofunction:: sotodlib.tod_ops.flags.get_source_flags """ name = "source_flags" - - def process(self, aman, proc_aman): - tod_ops.flags.get_source_flags(aman, **self.process_cfgs) def calc_and_save(self, aman, proc_aman): center_on = self.calc_cfgs.get('center_on', 'planet') diff --git a/sotodlib/site_pipeline/update_pointing.py b/sotodlib/site_pipeline/update_pointing.py index 9af4fc96f..2098ebbd3 100644 --- a/sotodlib/site_pipeline/update_pointing.py +++ b/sotodlib/site_pipeline/update_pointing.py @@ -316,7 +316,7 @@ def main_one_wafer(configs, obs_id, wafer_slot, sso_name=None, pipe = Pipeline(configs["process_pipe"], logger=logger) for pipe_component in pipe: if pipe_component.name == 'source_flags': - pipe_component.process_cfgs['center_on'] = sso_name + pipe_component.calc_cfgs['center_on'] = sso_name # Other parameters force_zero_roll = configs.get('force_zero_roll') From 89ec639c88d5e511501a91b6fc50da83a6ecb0f8 Mon Sep 17 00:00:00 2001 From: Elle Shaw Date: Fri, 23 May 2025 13:00:48 -0700 Subject: [PATCH 15/27] Renamed map_based_pointing and update pointing to get_brightsrc_pointing_step1 and step2. Added functionality to run as parallel processing job on a cluster. Included some documentation in site_pipeline.rst --- docs/site_pipeline.rst | 184 ++++++++++++++++++ ...ased_pointing.py => brightsrc_pointing.py} | 33 ++-- ...ing.py => get_brightsrc_pointing_step1.py} | 116 ++++++++--- ...ing.py => get_brightsrc_pointing_step2.py} | 87 +++++++-- 4 files changed, 361 insertions(+), 59 deletions(-) rename sotodlib/coords/{map_based_pointing.py => brightsrc_pointing.py} (96%) rename sotodlib/site_pipeline/{make_mapbased_pointing.py => get_brightsrc_pointing_step1.py} (75%) rename sotodlib/site_pipeline/{update_pointing.py => get_brightsrc_pointing_step2.py} (88%) diff --git a/docs/site_pipeline.rst b/docs/site_pipeline.rst index d228a5580..4f8d7481b 100644 --- a/docs/site_pipeline.rst +++ b/docs/site_pipeline.rst @@ -157,6 +157,190 @@ work, here's a more basic example that will work:: stream_id: ufm_mv9 + +get_brightsrc_pointing_part1 and get_brightsrc_pointing_part2 +---------------------------------------------------------------- + +The two-part ``get_brightsrc_pointing`` script set will will run solve for the xieta +coordinates of detectors that observe a bright source during an observation. + +It is a two part process that requires a map step and then a tod step. +To run, the scripts require config files described below. + +The code will process all wafers unless otherwise specified. +It is recommended to run with ``parallel_job: True`` in the config files if analyzing +multiple wafers at once. Otherwise, specify a wafer slot or restrict detectors in CL args. + +Recommended SLURM settings + - ``--nodes=1`` + - ``--ntasks=1`` + - ``--time=00:45:00`` + - ``--cpus-per-task=14`` + - ``--mem=150G`` + - export OMP_NUM_THREADS=1 + + +Recommended Command Line arguments: + - ``configs`` + - ``--obs_id`` + - ``--sso_name`` + +Optional Command Line arguments: + - ``--wafer_slot`` e.g. ws0 + - ``--restrict_dets_for_debug`` integer, or comma separated list of det readout_ids. + +Options to include min_ctime and max_ctime arguments, which will proces all obs +in the time frame. + +.. argparse:: + :module: sotodlib.site_pipeline.get_brightsrc_pointing_part1 or _part2 + :func: get_parser + + + +Generated results +``````````````````` +Saves results as ResultSet .hdf file in the results_dir. +ResultSet<[dets:readout_id, xi, eta, gamma, xi_err, eta_err, R2, redchi2], N rows> + +Load data with sotodlib.io.metadata.read_dataset( results.hdf, 'focal_plane') + +Configuration +````````````````` +These scripts take in a config yaml file + +Part 1 is the map-based step. Its config file should look like the following: +The parameters in these examples are used for SAT mid-freq moon observations. + +.. code-block:: yaml + + context_file: /path/to/context.yaml + query_tags: ['moon', 'jupiter', 'mars'] (alternatively specify --sso_name in kwargs + + optics_config_fn: /path/to/ufm_to_fp.yaml + single_det_maps_dir: /path/to/results/single_det_maps + results_dir: /path/to/results/map_based_results + + parallel_job: True + wafer_mask_det: 8. + res_deg: 0.3 + xieta_bs_offset: [0., 0.] + save_normal_roll: False #false for SAT, true for LAT + save_force_zero_roll: True #true for SAT, false for LAT + + hit_time_threshold: 600 #seconds + hit_circle_r_deg: 7. + + process_pipe: + - name: 'detrend' + process: + count: 2000 + method: 'linear' + - name: 'apodize' + process: + apodize_samps: 2000 + - name: 'fourier_filter' + process: + signal_name: "signal" + wrap_name: null + filt_function: "low_pass_sine2" + trim_samps: null + filter_params: + cutoff: 1.9 + width: 0.2 + - name: 'fourier_filter' + process: + signal_name: "signal" + wrap_name: null + filt_function: "high_pass_sine2" + trim_samps: 2000 + filter_params: + cutoff: 0.05 + width: 0.1 + +Part 2 is the TOD-based step. Its config file should look like the following. +The parameters in these examples are used for SAT mid-freq moon observations. + +.. code-block:: yaml + + context_file: /path/to/context.yaml + query_tags: ['moon', 'jupiter', 'mars'] (alternatively specify --sso_name in kwargs + + optics_config_fn: /path/to/ufm_to_fp.yaml + + fp_hdf_dir: /path/to/results/map_based_results from step 1 config file. + # If force_zero_roll is was True, then append _force_zero_roll to the end + result_dir: /path/to/resuls/tod_based_results + + parallel_job: True + force_zero_roll: True + + ds_factor: 40 + mask_deg: 2.5 + fit_func_name: 'gaussian2d_nonlin' + max_non_linear_order: 3 + fwhm_init_deg: 0.5 + error_estimation_method: 'force_one_redchi2' + flag_name_rms_calc: 'around_source' + flag_rms_calc_exclusive: False + + process_pipe: + - name: 'detrend' + process: + count: 2000 + method: 'linear' + - name: 'fourier_filter' + process: + signal_name: 'signal' + filt_function: 'iir_filter' + trim_samps: null + filter_params: + invert: True + - name: 'apodize' + process: + apodize_samps: 2000 + - name: 'fourier_filter' + process: + signal_name: "signal" + wrap_name: null + filt_function: "low_pass_sine2" + trim_samps: null + filter_params: + cutoff: 1.9 + width: 0.2 + - name: 'source_flags' + calc: + merge: True + max_pix: 10000000000 + source_flags_name: 'source_wide' + mask: + shape: circle + xyr: [0., 0., 5.0] + - name: 'source_flags' + calc: + merge: True + max_pix: 10000000000 + source_flags_name: 'source_narrow' + mask: + shape: circle + xyr: [0., 0., 3.0] + - name: 'reduce_flags' + process: + flags: ['source_wide', 'source_narrow'] + method: 'except' + wrap: True + new_flag: 'around_source' + - name: 'flag_turnarounds' + process: + truncate: True + - name: 'sub_polyf' + process: + method: 'legendre' + degree: 2 + mask: 'around_source' + exclusive: False + + make_read_det_match ``````````````````` This script generates the readout ID to detector ID mapping required to diff --git a/sotodlib/coords/map_based_pointing.py b/sotodlib/coords/brightsrc_pointing.py similarity index 96% rename from sotodlib/coords/map_based_pointing.py rename to sotodlib/coords/brightsrc_pointing.py index 43305b37c..5985da6fe 100644 --- a/sotodlib/coords/map_based_pointing.py +++ b/sotodlib/coords/brightsrc_pointing.py @@ -4,6 +4,7 @@ import numpy as np from scipy import interpolate from scipy.optimize import curve_fit +from joblib import Parallel, delayed from sotodlib import core from sotodlib import coords @@ -16,7 +17,8 @@ import h5py from scipy.ndimage import maximum_filter -def get_planet_trajectry(tod, planet, _split=20, return_model=False): + +def get_planet_trajectory(tod, planet, _split=20, return_model=False): """ Generate the trajectory of a given planet over a specified time range. @@ -30,7 +32,7 @@ def get_planet_trajectry(tod, planet, _split=20, return_model=False): If return_model is True: tuple: Tuple containing interpolation functions for azimuth and elevation. If return_model is False: - array: Array of quaternions representing trajectry of the planet at each timestamp. + array: Array of quaternions representing trajectory of the planet at each timestamp. """ print(planet) timestamps_sparse = np.linspace(tod.timestamps[0], tod.timestamps[-1], _split) @@ -58,19 +60,19 @@ def get_wafer_centered_sight(tod=None, planet=None, q_planet=None, q_bs=None, q_ Parameters: tod : An axis manager planet (str): The name of the planet to calculate the sightline vector. - q_planet (optional): Quaternion representing the trajectry of the planet. + q_planet (optional): Quaternion representing the trajectory of the planet. If None, it will be computed using get_planet_trajectory. Defaults to None. - q_bs (optional): Quaternion representing the trajectry of the boresight. + q_bs (optional): Quaternion representing the trajectory of the boresight. If None, it will be computed using the current boresight angles from tod. Defaults to None. q_wafer (optional): Quaternion representing the center of wafer to the center of boresight. If None, it will be computed using the median of the focal plane xi and eta from tod.focal_plane. Defaults to None. Returns: - Sightline vector for the planet trajectry centered on the center of the wafer. + Sightline vector for the planet trajectory centered on the center of the wafer. """ if q_planet is None: - q_planet = get_planet_trajectry(tod, planet) + q_planet = get_planet_trajectory(tod, planet) if q_bs is None: q_bs = quat.rotation_lonlat(tod.boresight.az, tod.boresight.el) if q_wafer is None: @@ -147,7 +149,7 @@ def get_rough_hit_time(tod, wafer_slot, sso_name, circle_r_deg=7.,optics_config_ float: Estimated rough hit time within the circular region around the wafer center. """ q_bs = quat.rotation_lonlat(tod.boresight.az, tod.boresight.el) - q_planet = get_planet_trajectry(tod, sso_name) + q_planet = get_planet_trajectory(tod, sso_name) xi_wafer, eta_wafer = get_wafer_xieta(wafer_slot, optics_config_fn=optics_config_fn, roll_bs_offset=np.median(tod.boresight.roll), wrap_to_tod=False) q_wafer = quat.rotation_xieta(xi_wafer, eta_wafer) @@ -159,13 +161,11 @@ def get_rough_hit_time(tod, wafer_slot, sso_name, circle_r_deg=7.,optics_config_ hit_time = (tod.timestamps[-1] - tod.timestamps[0]) * np.mean(np.rad2deg(r_wafer_centered) < circle_r_deg) return hit_time - def make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf, xieta_bs_offset=(0., 0.), roll_bs_offset=None, signal='signal', wafer_mask_deg=8., res_deg=0.3, cuts=None,): """ Generate boresight-centered maps from Time-Ordered Data (TOD) for each individual detector. - Parameters: tod : an axismanager object sso_name (str): Name of the planet for which the trajectory is calculated. @@ -181,9 +181,8 @@ def make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf, Returns: None """ - q_planet = get_planet_trajectry(tod, sso_name) + q_planet = get_planet_trajectory(tod, sso_name) q_bs = quat.rotation_lonlat(tod.boresight.az, tod.boresight.el) - if roll_bs_offset is None: roll_bs_offset = np.mean(tod.boresight.roll) @@ -218,16 +217,22 @@ def make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf, if cuts is None: cuts = ~tod.flags['source'] P = coords.P.for_tod(tod=tod, geom=geom, comps='T', cuts=cuts, sight=sight, threads=False) + + wT = None for di, det in enumerate(tqdm(tod.dets.vals)): det_weights = np.zeros(tod.dets.count, dtype='float32') det_weights[di] = 1. mT_weighted = P.to_map(tod=tod, signal=signal, comps='T', det_weights=det_weights) - wT = P.to_weights(tod, signal=signal, comps='T', det_weights=det_weights) + if wT is None: + wT = P.to_weights(tod, signal=signal, comps='T', det_weights=det_weights) mT = P.remove_weights(signal_map=mT_weighted, weights_map=wT, comps='T')[0] enmap.write_hdf(map_hdf, mT, address=det, - extra={'xi0': xi0, 'eta0': eta0, - 'xi_bs_offset': xi_bs_offset, 'eta_bs_offset': eta_bs_offset, 'roll_bs_offset': roll_bs_offset}) + extra={'xi0': xi0, + 'eta0': eta0, + 'xi_bs_offset': xi_bs_offset, + 'eta_bs_offset': eta_bs_offset, + 'roll_bs_offset': roll_bs_offset}) return def detect_peak_xieta(mT, filter_size=None): diff --git a/sotodlib/site_pipeline/make_mapbased_pointing.py b/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py similarity index 75% rename from sotodlib/site_pipeline/make_mapbased_pointing.py rename to sotodlib/site_pipeline/get_brightsrc_pointing_step1.py index faed09a43..2b7aea511 100644 --- a/sotodlib/site_pipeline/make_mapbased_pointing.py +++ b/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py @@ -4,12 +4,12 @@ import argparse import time import glob +from joblib import Parallel, delayed from sotodlib import core from sotodlib import coords from sotodlib import tod_ops -from sotodlib.coords import map_based_pointing as mbp -from sotodlib.site_pipeline import update_pointing as up +from sotodlib.coords import brightsrc_pointing as bsp from sotodlib.io import metadata from sotodlib.io.metadata import read_dataset, write_dataset @@ -17,7 +17,7 @@ from sotodlib.preprocess import Pipeline logger = util.init_logger(__name__, 'make_map_based_pointing: ') -def _get_sso_names_from_tags(ctx, obs_id, candidate_names=['moon', 'jupiter']): +def _get_sso_names_from_tags(ctx, obs_id, candidate_names=['moon', 'jupiter', 'mars']): obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] sso_names = [] for _name in candidate_names: @@ -58,31 +58,44 @@ def main_one_wafer(configs, obs_id, wafer_slot, sso_name=None, f'Processing only {sso_name}') # Load data - logger.info('loading data') + logger.info(f'loading meta data: {wafer_slot}') meta = ctx.get_meta(obs_id, dets={'wafer_slot': wafer_slot}) + logger.info(f'finished loading meta data: {wafer_slot}') + try: + meta.restrict('dets', meta.detcal.bg > -1) + except: + pass if restrict_dets_for_debug is not False: - meta.restrict('dets', meta.dets.vals[:restrict_dets_for_debug]) + try: + restrict_dets_for_debug = int(restrict_dets_for_debug) + meta.restrict('dets', meta.dets.vals[:restrict_dets_for_debug]) + except ValueError: + _testdets = restrict_dets_for_debug.split(',') + restrict_list = [det.split('\'')[1].strip() for det in _testdets] + meta.restrict('dets', restrict_list) + logger.info(f'loading tod data: {wafer_slot}') tod = ctx.get_obs(meta) - + logger.info(f'finished loading tod data: {wafer_slot}') # tod processing - logger.info('tod processing') + logger.info(f'tod processing {wafer_slot}') pipe = Pipeline(configs["process_pipe"], logger=logger) proc_aman, success = pipe.run(tod) - + logger.info(f'done with tod processing {wafer_slot}') # make single detecctor maps logger.info(f'Making single detector maps') os.makedirs(single_det_maps_dir, exist_ok=True) map_hdf = os.path.join(single_det_maps_dir, f'{obs_id}_{wafer_slot}.hdf') - mbp.make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf=map_hdf, + bsp.make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf=map_hdf, xieta_bs_offset=xieta_bs_offset, wafer_mask_deg=wafer_mask_deg, res_deg=res_deg) - + + #next step result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' # reconstruct pointing from single detector maps if save_normal_roll: logger.info(f'Saving map-based pointing results') - fp_rset_map_based = mbp.get_xieta_from_maps(map_hdf, save=True, + fp_rset_map_based = bsp.get_xieta_from_maps(map_hdf, save=True, output_dir=result_dir, filename=result_filename, force_zero_roll=False, @@ -91,7 +104,7 @@ def main_one_wafer(configs, obs_id, wafer_slot, sso_name=None, if save_force_zero_roll: logger.info(f'Saving map-based pointing results (force-zero-roll)') result_dir_force_zero_roll = result_dir + '_force_zero_roll' - fp_rset_map_based_force_zero_roll = mbp.get_xieta_from_maps(map_hdf, save=True, + fp_rset_map_based_force_zero_roll = bsp.get_xieta_from_maps(map_hdf, save=True, output_dir=result_dir_force_zero_roll, filename=result_filename, force_zero_roll=True, @@ -109,7 +122,13 @@ def main_one_wafer_dummy(configs, obs_id, wafer_slot, restrict_dets_for_debug=Fa meta = ctx.get_meta(obs_id, dets={'wafer_slot': wafer_slot}) if restrict_dets_for_debug is not False: - meta.restrict('dets', meta.dets.vals[:restrict_dets_for_debug]) + try: + restrict_dets_for_debug = int(restrict_dets_for_debug) + meta.restrict('dets', meta.dets.vals[:restrict_dets_for_debug]) + except ValueError: + _testdets = restrict_dets_for_debug.split(',') + restrict_list = [det.split('\'')[1].strip() for det in _testdets] + meta.restrict('dets', restrict_list) result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' fp_rset_dummy_map_based = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', 'R2']) @@ -150,6 +169,15 @@ def combine_pointings(pointing_result_files): focal_plane.rows.append((det, val['xi'], val['eta'], val['gamma'], val['R2'])) return focal_plane +def parallel_process_wafer_slot(configs, obs_id, wafer_slot, sso_name, restrict_dets_for_debug): + logger.info(f'Processing {obs_id}, {wafer_slot}') + main_one_wafer(configs=configs, + obs_id=obs_id, + wafer_slot=wafer_slot, + sso_name=sso_name, + restrict_dets_for_debug=restrict_dets_for_debug) + + def main_one_obs(configs, obs_id, sso_name=None, restrict_dets_for_debug=False): if type(configs) == str: @@ -176,26 +204,57 @@ def main_one_obs(configs, obs_id, sso_name=None, tod = ctx.get_obs(obs_id, dets=[]) streamed_wafer_slots = ['ws{}'.format(index) for index, bit in enumerate(obs_id.split('_')[-1]) if bit == '1'] processed_wafer_slots = [] + finished_wafer_slots = [] skipped_wafer_slots = [] + check_dir = result_dir + '_force_zero_roll' if save_force_zero_roll else result_dir for ws in streamed_wafer_slots: - hit_time = mbp.get_rough_hit_time(tod, wafer_slot=ws, sso_name=sso_name, circle_r_deg=hit_circle_r_deg, - optics_config_fn=optics_config_fn) + hit_time = bsp.get_rough_hit_time(tod, + wafer_slot=ws, + sso_name=sso_name, + circle_r_deg=hit_circle_r_deg, + optics_config_fn=optics_config_fn) logger.info(f'hit_time for {ws} is {hit_time:.1f} [sec]') - if hit_time > hit_time_threshold: - processed_wafer_slots.append(ws) + if hit_time >= hit_time_threshold: + if os.path.exists(os.path.join(check_dir, f'focal_plane_{obs_id}_{ws}.hdf')): + finished_wafer_slots.append(ws) + else: + processed_wafer_slots.append(ws) else: skipped_wafer_slots.append(ws) - logger.info(f'wafer_slots which pointing calculated: {processed_wafer_slots}') - for wafer_slot in processed_wafer_slots: - logger.info(f'Processing {obs_id}, {wafer_slot}') - main_one_wafer(configs=configs, - obs_id=obs_id, - wafer_slot=wafer_slot, - sso_name=sso_name, - restrict_dets_for_debug=restrict_dets_for_debug) - + logger.info(f'Found saved data for these wafer_slots: {finished_wafer_slots}') + logger.info(f'Will continue for these wafer_slots: {processed_wafer_slots}') + + if configs.get('parallel_job'): + logger.info('Continuing with parallel job') + try: + n_jobs = int(os.environ.get('SLURM_CPUS_PER_TASK', 1)) + except: + n_jobs = -1 + + logger.info('Entering wafer pool') + Parallel(n_jobs=n_jobs)( + delayed(parallel_process_wafer_slot)( + configs, + obs_id, + wafer_slot, + sso_name, + restrict_dets_for_debug, + ) + for wafer_slot in processed_wafer_slots + ) + logger.info('Exiting wafer pool') + else: + logger.info('Continuing with serial processing of wafers.') + for wafer_slot in processed_wafer_slots: + logger.info(f'Processing {obs_id}, {wafer_slot}') + main_one_wafer(configs=configs, + obs_id=obs_id, + wafer_slot=wafer_slot, + sso_name=sso_name, + restrict_dets_for_debug=restrict_dets_for_debug) + logger.info(f'create dummy hdf for non-hitting wafer: {skipped_wafer_slots}') for wafer_slot in skipped_wafer_slots: main_one_wafer_dummy(configs=configs, @@ -220,7 +279,7 @@ def main_one_obs(configs, obs_id, sso_name=None, write_dataset(fp_rset_full, filename=fp_rset_full_file, address='focal_plane', overwrite=True) - + logger.info(f'ta da! Finished with {obs_id}') return def main(configs, min_ctime=None, max_ctime=None, update_delay=None, @@ -258,6 +317,7 @@ def main(configs, min_ctime=None, max_ctime=None, update_delay=None, restrict_dets_for_debug=restrict_dets_for_debug) elif obs_id is not None: + logger.info(f'Processing {obs_id}') if wafer_slot is None: main_one_obs(configs=configs, obs_id=obs_id, sso_name=sso_name, restrict_dets_for_debug=restrict_dets_for_debug) @@ -280,7 +340,7 @@ def get_parser(): parser.add_argument("--sso_name", type=str, default=None, help="Name of solar system object (e.g., 'moon', 'jupiter'). If not specified, get sso_name from observation tags. "\ + "Valid only when obs_id is specified") - parser.add_argument("--restrict_dets_for_debug", type=int, default=False) + parser.add_argument("--restrict_dets_for_debug", type=str, default=False) return parser if __name__ == '__main__': diff --git a/sotodlib/site_pipeline/update_pointing.py b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py similarity index 88% rename from sotodlib/site_pipeline/update_pointing.py rename to sotodlib/site_pipeline/get_brightsrc_pointing_step2.py index 2098ebbd3..6081d1800 100644 --- a/sotodlib/site_pipeline/update_pointing.py +++ b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py @@ -6,11 +6,12 @@ import time import glob from tqdm import tqdm +from joblib import Parallel, delayed from scipy.optimize import curve_fit from sotodlib.core import metadata from sotodlib.io.metadata import read_dataset, write_dataset -from sotodlib.coords import map_based_pointing as mbp +from sotodlib.coords import brightsrc_pointing as bsp from sotodlib import core from sotodlib import coords from sotodlib import tod_ops @@ -334,7 +335,14 @@ def main_one_wafer(configs, obs_id, wafer_slot, sso_name=None, logger.info('loading data') meta = ctx.get_meta(obs_id, dets={'wafer_slot': wafer_slot}) if restrict_dets_for_debug is not False: - meta.restrict('dets', meta.dets.vals[:restrict_dets_for_debug]) + try: + restrict_dets_for_debug = int(restrict_dets_for_debug) + meta.restrict('dets', meta.dets.vals[:restrict_dets_for_debug]) + except ValueError: + _testdets = restrict_dets_for_debug.split(',') + restrict_list = [det.split('\'')[1].strip() for det in _testdets] + meta.restrict('dets', restrict_list) + tod = ctx.get_obs(meta) # get pointing @@ -372,7 +380,13 @@ def main_one_wafer_dummy(configs, obs_id, wafer_slot, restrict_dets_for_debug=Fa meta = ctx.get_meta(obs_id, dets={'wafer_slot': wafer_slot}) if restrict_dets_for_debug is not False: - meta.restrict('dets', meta.dets.vals[:restrict_dets_for_debug]) + try: + restrict_dets_for_debug = int(restrict_dets_for_debug) + meta.restrict('dets', meta.dets.vals[:restrict_dets_for_debug]) + except ValueError: + _testdets = restrict_dets_for_debug.split(',') + restrict_list = [det.split('\'')[1].strip() for det in _testdets] + meta.restrict('dets', restrict_list) result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' fp_rset_dummy = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', @@ -408,6 +422,14 @@ def combine_pointings(pointing_result_files): focal_plane.rows.append((det, val['xi'], val['eta'], val['gamma'], val['xi_err'], val['eta_err'], val['R2'], val['redchi2'])) return focal_plane +def parallel_process_wafer_slot(configs, obs_id, wafer_slot, sso_name, restrict_dets_for_debug): + logger.info(f'Processing {obs_id}, {wafer_slot}') + main_one_wafer(configs=configs, + obs_id=obs_id, + wafer_slot=wafer_slot, + sso_name=sso_name, + restrict_dets_for_debug=restrict_dets_for_debug) + def main_one_obs(configs, obs_id, sso_name=None, restrict_dets_for_debug=False): if type(configs) == str: @@ -434,25 +456,54 @@ def main_one_obs(configs, obs_id, sso_name=None, tod = ctx.get_obs(obs_id, dets=[]) streamed_wafer_slots = ['ws{}'.format(index) for index, bit in enumerate(obs_id.split('_')[-1]) if bit == '1'] processed_wafer_slots = [] + finished_wafer_slots = [] skipped_wafer_slots = [] + check_dir = result_dir + '_force_zero_roll' if force_zero_roll else result_dir + for ws in streamed_wafer_slots: - hit_time = mbp.get_rough_hit_time(tod, wafer_slot=ws, sso_name=sso_name, circle_r_deg=hit_circle_r_deg, - optics_config_fn=optics_config_fn) + hit_time = bsp.get_rough_hit_time(tod, + wafer_slot=ws, + sso_name=sso_name, + circle_r_deg=hit_circle_r_deg, + optics_config_fn=optics_config_fn) logger.info(f'hit_time for {ws} is {hit_time:.1f} [sec]') - if hit_time > hit_time_threshold: - processed_wafer_slots.append(ws) + if hit_time >= hit_time_threshold: + if os.path.exists(os.path.join(check_dir, f'focal_plane_{obs_id}_{ws}.hdf')): + finished_wafer_slots.append(ws) + else: + processed_wafer_slots.append(ws) else: skipped_wafer_slots.append(ws) - - logger.info(f'wafer_slots which pointing calculated: {processed_wafer_slots}') - for wafer_slot in processed_wafer_slots: - logger.info(f'Processing {obs_id}, {wafer_slot}') - main_one_wafer(configs=configs, - obs_id=obs_id, - wafer_slot=wafer_slot, - sso_name=sso_name, - restrict_dets_for_debug=restrict_dets_for_debug) - + + logger.info(f'Found saved data for these wafer_slots: {finished_wafer_slots}') + logger.info(f'Will continue for these wafer_slots: {processed_wafer_slots}') + + if configs.get('parallel_job'): + logger.info('Continuing with parallel job') + try: + n_jobs = int(os.environ.get('SLURM_CPUS_PER_TASK', 1)) + except: + n_jobs = -1 + Parallel(n_jobs=n_jobs)( + delayed(parallel_process_wafer_slot)( + configs, + obs_id, + wafer_slot, + sso_name, + restrict_dets_for_debug + ) + for wafer_slot in processed_wafer_slots + ) + else: + logger.info('Continuing with serial processing of wafers.') + for wafer_slot in processed_wafer_slots: + logger.info(f'Processing {obs_id}, {wafer_slot}') + main_one_wafer(configs=configs, + obs_id=obs_id, + wafer_slot=wafer_slot, + sso_name=sso_name, + restrict_dets_for_debug=restrict_dets_for_debug) + logger.info(f'create dummy hdf for non-hitting wafer: {skipped_wafer_slots}') for wafer_slot in skipped_wafer_slots: main_one_wafer_dummy(configs=configs, @@ -466,6 +517,7 @@ def main_one_obs(configs, obs_id, sso_name=None, fp_rset_full_file = os.path.join(os.path.join(result_dir, f'focal_plane_{obs_id}_all.hdf')) write_dataset(fp_rset_full, filename=fp_rset_full_file, address='focal_plane', overwrite=True) + logger.info(f'ta da! Finsihed with {obs_id}') def main(configs, min_ctime=None, max_ctime=None, update_delay=None, obs_id=None, wafer_slot=None, sso_name=None, restrict_dets_for_debug=False): @@ -502,6 +554,7 @@ def main(configs, min_ctime=None, max_ctime=None, update_delay=None, restrict_dets_for_debug=restrict_dets_for_debug) elif obs_id is not None: + logger.info(f'Processing {obs_id}') if wafer_slot is None: main_one_obs(configs=configs, obs_id=obs_id, sso_name=sso_name, restrict_dets_for_debug=restrict_dets_for_debug) From 954bdd5502183cb20faaa5a620ee1abebe16b8e1 Mon Sep 17 00:00:00 2001 From: Elle Shaw Date: Mon, 25 Aug 2025 14:12:29 -0700 Subject: [PATCH 16/27] Made small changes to make code be compatible with latest master branch. --- sotodlib/coords/brightsrc_pointing.py | 2 +- sotodlib/preprocess/processes.py | 3 ++- sotodlib/site_pipeline/get_brightsrc_pointing_step2.py | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/sotodlib/coords/brightsrc_pointing.py b/sotodlib/coords/brightsrc_pointing.py index 5985da6fe..8330fce1f 100644 --- a/sotodlib/coords/brightsrc_pointing.py +++ b/sotodlib/coords/brightsrc_pointing.py @@ -107,7 +107,7 @@ def get_wafer_xieta(wafer_slot, optics_config_fn, xieta_bs_offset=(0., 0.), wafer_r = np.sqrt(wafer_x**2 + wafer_y**2) wafer_theta = np.arctan2(wafer_y, wafer_x) - fp_to_sky = optics.sat_to_sky(optics.SAT_X, optics.SAT_LON) + fp_to_sky = optics.sat_to_sky(optics.SAT_R_FP, optics.SAT_R_SKY) lon = fp_to_sky(wafer_r) q1 = quat.rotation_iso(lon, 0) diff --git a/sotodlib/preprocess/processes.py b/sotodlib/preprocess/processes.py index ad274bf9d..cbfd0f7a0 100644 --- a/sotodlib/preprocess/processes.py +++ b/sotodlib/preprocess/processes.py @@ -1705,7 +1705,7 @@ def process(self, aman, proc_aman, sim=False): class ReduceFlags(_Preprocess): name = 'reduce_flags' - def process(self, aman, proc_aman): + def process(self, aman, proc_aman, sim=False): aman.flags.reduce(**self.process_cfgs) class DetcalNanCuts(_Preprocess): @@ -2506,6 +2506,7 @@ def save(self, proc_aman, flag_aman): if self.save_cfgs: proc_aman.wrap("smurfgaps", flag_aman) +_Preprocess.register(ReduceFlags) _Preprocess.register(SplitFlags) _Preprocess.register(SubtractT2P) _Preprocess.register(EstimateT2P) diff --git a/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py index 6081d1800..89354f9d3 100644 --- a/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py +++ b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py @@ -70,7 +70,7 @@ def wrapper_gaussian2d_nonlin(xieta, xi0, eta0, fwhm_xi, fwhm_eta, phi, a, *args return gaussian2d_nonlin(xieta, xi0, eta0, fwhm_xi, fwhm_eta, phi, a, nonlin_coeffs) def wrap_fp_rset(tod, fp_rset): - tod.restrict('dets', tod.dets.vals[np.in1d(tod.dets.vals, fp_rset['dets:readout_id'])]) + tod.restrict('dets', tod.dets.vals[np.isin(tod.dets.vals, fp_rset['dets:readout_id'])]) focal_plane = core.AxisManager(tod.dets) focal_plane.wrap_new('xi', shape=('dets', )) focal_plane.wrap_new('eta', shape=('dets', )) From 092b2ab8b8fc7b413b0751bdff06ce0c29de3f18 Mon Sep 17 00:00:00 2001 From: Elle Shaw Date: Fri, 29 Aug 2025 12:28:38 -0700 Subject: [PATCH 17/27] Removed offensive square brackets in a get_obs call --- sotodlib/site_pipeline/get_brightsrc_pointing_step1.py | 2 +- sotodlib/site_pipeline/get_brightsrc_pointing_step2.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py b/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py index 2b7aea511..331165ce6 100644 --- a/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py +++ b/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py @@ -201,7 +201,7 @@ def main_one_obs(configs, obs_id, sso_name=None, logger.info(f'sso_names of {sso_names} are found from observation tags.' + f'Processing only {sso_name}') - tod = ctx.get_obs(obs_id, dets=[]) + tod = ctx.get_obs(obs_id, no_signal=True) streamed_wafer_slots = ['ws{}'.format(index) for index, bit in enumerate(obs_id.split('_')[-1]) if bit == '1'] processed_wafer_slots = [] finished_wafer_slots = [] diff --git a/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py index 89354f9d3..6842cbce6 100644 --- a/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py +++ b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py @@ -453,7 +453,7 @@ def main_one_obs(configs, obs_id, sso_name=None, logger.info(f'sso_names of {sso_names} are found from observation tags.' + f'Processing only {sso_name}') - tod = ctx.get_obs(obs_id, dets=[]) + tod = ctx.get_obs(obs_id, no_signal=True) streamed_wafer_slots = ['ws{}'.format(index) for index, bit in enumerate(obs_id.split('_')[-1]) if bit == '1'] processed_wafer_slots = [] finished_wafer_slots = [] From 951de7a33dee74fcef02442599447bce891e04e2 Mon Sep 17 00:00:00 2001 From: Elle Shaw Date: Wed, 3 Sep 2025 12:05:52 -0700 Subject: [PATCH 18/27] add saturn to possible source objects --- sotodlib/site_pipeline/get_brightsrc_pointing_step1.py | 4 ++-- sotodlib/site_pipeline/get_brightsrc_pointing_step2.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py b/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py index 331165ce6..17d99f039 100644 --- a/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py +++ b/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py @@ -17,7 +17,7 @@ from sotodlib.preprocess import Pipeline logger = util.init_logger(__name__, 'make_map_based_pointing: ') -def _get_sso_names_from_tags(ctx, obs_id, candidate_names=['moon', 'jupiter', 'mars']): +def _get_sso_names_from_tags(ctx, obs_id, candidate_names=['moon', 'jupiter', 'mars', 'saturn']): obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] sso_names = [] for _name in candidate_names: @@ -225,7 +225,7 @@ def main_one_obs(configs, obs_id, sso_name=None, logger.info(f'Found saved data for these wafer_slots: {finished_wafer_slots}') logger.info(f'Will continue for these wafer_slots: {processed_wafer_slots}') - + logger.info("using filelock") if configs.get('parallel_job'): logger.info('Continuing with parallel job') try: diff --git a/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py index 6842cbce6..20ba5b754 100644 --- a/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py +++ b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py @@ -22,7 +22,7 @@ from sotodlib.preprocess import Pipeline logger = util.init_logger(__name__, 'update_pointing: ') -def _get_sso_names_from_tags(ctx, obs_id, candidate_names=['moon', 'jupiter']): +def _get_sso_names_from_tags(ctx, obs_id, candidate_names=['moon', 'jupiter', 'mars', 'saturn']): obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] sso_names = [] for _name in candidate_names: From 2a38cb3ecc46eb25785100d85dbfeafbbf7a156c Mon Sep 17 00:00:00 2001 From: Elle Shaw Date: Thu, 11 Sep 2025 07:41:41 -0700 Subject: [PATCH 19/27] removed obsolete logger entry --- sotodlib/site_pipeline/get_brightsrc_pointing_step1.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py b/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py index 17d99f039..e29c3bc76 100644 --- a/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py +++ b/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py @@ -225,7 +225,7 @@ def main_one_obs(configs, obs_id, sso_name=None, logger.info(f'Found saved data for these wafer_slots: {finished_wafer_slots}') logger.info(f'Will continue for these wafer_slots: {processed_wafer_slots}') - logger.info("using filelock") + if configs.get('parallel_job'): logger.info('Continuing with parallel job') try: From 512d4eb085bf124080157473cd669120f5c92405 Mon Sep 17 00:00:00 2001 From: Elle Shaw Date: Wed, 19 Nov 2025 14:31:42 -0800 Subject: [PATCH 20/27] fixed bug on including a list of detectors to debug with --- sotodlib/site_pipeline/get_brightsrc_pointing_step2.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py index 20ba5b754..b648010a1 100644 --- a/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py +++ b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py @@ -579,7 +579,7 @@ def get_parser(): parser.add_argument("--sso_name", type=str, default=None, help="Name of solar system object (e.g., 'moon', 'jupiter'). If not specified, get sso_name from observation tags. "\ + "Valid only when obs_id is specified") - parser.add_argument("--restrict_dets_for_debug", type=int, default=False) + parser.add_argument("--restrict_dets_for_debug", type=str, default=False) return parser if __name__ == '__main__': From e7f007fd1a7e215ef55ddbde792d0c6283ec9b25 Mon Sep 17 00:00:00 2001 From: Elle Shaw Date: Tue, 20 Jan 2026 15:19:16 -0800 Subject: [PATCH 21/27] Updated rst documentation with example NERSC submission script --- docs/site_pipeline.rst | 142 ++++++++++++++++++++++++++++------------- 1 file changed, 96 insertions(+), 46 deletions(-) diff --git a/docs/site_pipeline.rst b/docs/site_pipeline.rst index e87110890..6a103d981 100644 --- a/docs/site_pipeline.rst +++ b/docs/site_pipeline.rst @@ -197,28 +197,26 @@ if it does not exist. get_brightsrc_pointing_part1 and get_brightsrc_pointing_part2 ----------------------------------------------------------------- +--------------------------------------------------------------- The two-part ``get_brightsrc_pointing`` script set will will run solve for the xieta coordinates of detectors that observe a bright source during an observation. -It is a two part process that requires a map step and then a tod step. -To run, the scripts require config files described below. +It is a two part process that requires a map step and then a TOD step. +The scripts require the settings and preprocessing config files described below. + +For job submission and parallelization, see example NERSC slurm submission config at the end of this section. The code will process all wafers unless otherwise specified. It is recommended to run with ``parallel_job: True`` in the config files if analyzing -multiple wafers at once. Otherwise, specify a wafer slot or restrict detectors in CL args. - -Recommended SLURM settings - - ``--nodes=1`` - - ``--ntasks=1`` - - ``--time=00:45:00`` - - ``--cpus-per-task=14`` - - ``--mem=150G`` - - export OMP_NUM_THREADS=1 +multiple wafers at once. +Otherwise, specify a wafer slot or restrict detectors in command line args to debug. +.. argparse:: + :module: sotodlib.site_pipeline.get_brightsrc_pointing_part1 or _part2 + :func: get_parser -Recommended Command Line arguments: +Command Line arguments: - ``configs`` - ``--obs_id`` - ``--sso_name`` @@ -227,47 +225,55 @@ Optional Command Line arguments: - ``--wafer_slot`` e.g. ws0 - ``--restrict_dets_for_debug`` integer, or comma separated list of det readout_ids. -Options to include min_ctime and max_ctime arguments, which will proces all obs -in the time frame. - -.. argparse:: - :module: sotodlib.site_pipeline.get_brightsrc_pointing_part1 or _part2 - :func: get_parser - +There are options to include min_ctime and max_ctime arguments, which will process all observations +in the time frame. (not recommended) Generated results ``````````````````` -Saves results as ResultSet .hdf file in the results_dir. -ResultSet<[dets:readout_id, xi, eta, gamma, xi_err, eta_err, R2, redchi2], N rows> +The Step 1 map-based analysis scripts will generate the following outputs in the specified directory: -Load data with sotodlib.io.metadata.read_dataset( results.hdf, 'focal_plane') + 1. Single detector maps in ``/results/single_det_maps/_.hdf``. -Configuration -````````````````` -These scripts take in a config yaml file + * All single maps are packaged in a single hdf file, with detector readout_id as the keys in the h5py file. -Part 1 is the map-based step. Its config file should look like the following: -The parameters in these examples are used for SAT mid-freq moon observations. + 2. Fitted xi-eta focal plane position results saved as ResultSet in ``/path/to/results/map_based_results`` as specified in the Step 1 config file. Script will append 'force_zero_roll' onto the specified results_dir if True in config file. Load ResultSet with keyword 'focal_plane' + + * Contents: ``ResultSet<[dets:readout_id, xi, eta, gamma, R2], N rows>`` + + +The Step 2 TOD-based analysis scripts will use the map-based results as a starting point and then generate the finalized outputs in the specified directory: + + 1. Fitted xi-eta focal plane position results saved as ResultSet in ``/path/to/results/tod_based_results`` as specified in config file for Step-2. Script will append 'force_zero_roll' onto the specified results_dir if True in config file. Load ResultSet with keyword 'focal_plane' + + * Contents: ``ResultSet<[dets:readout_id, xi, eta, gamma, xi_err, eta_err, R2, redchi2], N rows>`` + +Configuration Files +``````````````````` +The configuration files to be input as ``configs`` in the command line arguments should have the following arguments as well as any preprocessing steps wished to be taken. Only processing steps that are agnostic of det-match can be used to do initial analyses without formalized metadata. + +The parameters in these examples could be used for SAT mid-freq moon observations. + +Step 1 Config: .. code-block:: yaml context_file: /path/to/context.yaml - query_tags: ['moon', 'jupiter', 'mars'] (alternatively specify --sso_name in kwargs + query_tags: ['moon'=1] #(alternatively specify --sso_name in kwargs) - optics_config_fn: /path/to/ufm_to_fp.yaml + optics_config_fn: '/global/cfs/cdirs/sobs/users/elleshaw/process_brightsrc/ufm_to_fp.yaml' single_det_maps_dir: /path/to/results/single_det_maps results_dir: /path/to/results/map_based_results - parallel_job: True - wafer_mask_det: 8. + parallel_job: True #For job submission. Parallel across wafers. + wafer_mask_det: 8. #mask around detector to cut TOD when source too far away. res_deg: 0.3 - xieta_bs_offset: [0., 0.] + xieta_bs_offset: [0., 0.] #Good to input xieta offset in radians. (!!! for satp2) save_normal_roll: False #false for SAT, true for LAT save_force_zero_roll: True #true for SAT, false for LAT - hit_time_threshold: 600 #seconds - hit_circle_r_deg: 7. + hit_circle_r_deg: 7. # radial mask to decide which UFMs are hit by source and should be analyzed. + hit_time_threshold: 600 #seconds, if hit_time not met then UFM does not get analyzed. process_pipe: - name: 'detrend' @@ -302,22 +308,20 @@ The parameters in these examples are used for SAT mid-freq moon observations. .. code-block:: yaml context_file: /path/to/context.yaml - query_tags: ['moon', 'jupiter', 'mars'] (alternatively specify --sso_name in kwargs - - optics_config_fn: /path/to/ufm_to_fp.yaml + query_tags: ['moon'=1] #(alternatively specify --sso_name in kwargs) + optics_config_fn: '/global/cfs/cdirs/sobs/users/elleshaw/process_brightsrc/ufm_to_fp.yaml' fp_hdf_dir: /path/to/results/map_based_results from step 1 config file. - # If force_zero_roll is was True, then append _force_zero_roll to the end - result_dir: /path/to/resuls/tod_based_results + # If force_zero_roll is was True, then append _force_zero_roll to the end. Just make sure it matches where the results from Step 1 are. + result_dir: /path/to/resuls/tod_based_results #Where you want the final Step2 results to show up. parallel_job: True - force_zero_roll: True - + force_zero_roll: True #Results will show up roatated in the xi-eta results as they are on the sky. ds_factor: 40 - mask_deg: 2.5 + mask_deg: 2.5 # size for circular mask around SSO (helps exclude focal plane reflections too) fit_func_name: 'gaussian2d_nonlin' - max_non_linear_order: 3 - fwhm_init_deg: 0.5 + max_non_linear_order: 3 #Suggested to use 1 for jupiter or sso's that do not saturate. + fwhm_init_deg: 0.5 # Lower for SATp2 error_estimation_method: 'force_one_redchi2' flag_name_rms_calc: 'around_source' flag_rms_calc_exclusive: False @@ -378,9 +382,55 @@ The parameters in these examples are used for SAT mid-freq moon observations. mask: 'around_source' exclusive: False +Example NERSC slurm job submission config file +`````````````````````````````````````````````` + +.. code-block:: yaml + #!/bin/bash -l + + #SBATCH --qos=shared + #SBATCH --constraint=cpu + #SBATCH --nodes=1 + #SBATCH --ntasks=1 + + #SBATCH --cpus-per-task=14 + #SBATCH --time=00:30:00 + #SBATCH --mem=220G`` #(may require regular queue and up to 400 Gb for extra long observations) + + export OMP_NUM_THREADS=1 + set -e + + tele=$1 + obs=$2 + map=$3 + basis=$4 + source="moon_from_moon" + + ymldir="/path/to/processing_settings_config_folder" + yfile="${ymldir}/preprocess_config_moon_${basis}_based_${tele}.yaml" + + if (($map)); then + echo submitted map job; + srun -n 1 -N 1 -c 14 python3 /path/to/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py $yfile --obs_id=${2} --sso_name="moon"; + else + echo submitted tod job; + srun -n 1 -N 1 -c 14 python3 /path/to/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py $yfile --obs_id=${2} --sso_name="moon"; + fi + + +Submit the job submission file with the following commands: + +1. For Step 1 map-based + + * ``sbatch submit_moon_job_script.sh 1 map`` + +2. For Step 2 TOD based + + * ``sbatch submit_moon_job_script.sh 0 tod`` + make_read_det_match -``````````````````` +------------------- This script generates the readout ID to detector ID mapping required to translate between the detector hardware information (ex: pixel position) and the readout IDs of the resonators used to index the SMuRF data. The script uses the From 2de0b169fb2964146746e061602a83eed10f15d9 Mon Sep 17 00:00:00 2001 From: Elle Shaw Date: Mon, 16 Feb 2026 12:38:11 -0800 Subject: [PATCH 22/27] Changes according to requested comments. Condensed ReduceFlags into CombineFlags, implemented MPI parallel processing, and other syntax updates. --- docs/site_pipeline.rst | 96 +++++++++++-------- sotodlib/coords/brightsrc_pointing.py | 17 ++-- sotodlib/preprocess/processes.py | 32 +++---- .../get_brightsrc_pointing_step1.py | 61 ++++++------ .../get_brightsrc_pointing_step2.py | 69 +++++++------ sotodlib/tod_ops/sub_polyf.py | 7 +- 6 files changed, 152 insertions(+), 130 deletions(-) diff --git a/docs/site_pipeline.rst b/docs/site_pipeline.rst index 6a103d981..ea9ea5537 100644 --- a/docs/site_pipeline.rst +++ b/docs/site_pipeline.rst @@ -196,10 +196,10 @@ The output database ``wafer_info.sqlite`` and HDF5 file if it does not exist. -get_brightsrc_pointing_part1 and get_brightsrc_pointing_part2 ---------------------------------------------------------------- +get_brightsrc_pointing_part1 +---------------------------- -The two-part ``get_brightsrc_pointing`` script set will will run solve for the xieta +The two-part ``get_brightsrc_pointing_part{}`` script set will solve for the xieta coordinates of detectors that observe a bright source during an observation. It is a two part process that requires a map step and then a TOD step. @@ -212,21 +212,14 @@ It is recommended to run with ``parallel_job: True`` in the config files if anal multiple wafers at once. Otherwise, specify a wafer slot or restrict detectors in command line args to debug. +Command Line arguments: .. argparse:: - :module: sotodlib.site_pipeline.get_brightsrc_pointing_part1 or _part2 + :module: sotodlib.site_pipeline.get_brightsrc_pointing_part1 :func: get_parser -Command Line arguments: - - ``configs`` - - ``--obs_id`` - - ``--sso_name`` - -Optional Command Line arguments: - - ``--wafer_slot`` e.g. ws0 - - ``--restrict_dets_for_debug`` integer, or comma separated list of det readout_ids. -There are options to include min_ctime and max_ctime arguments, which will process all observations -in the time frame. (not recommended) +There options to include min_ctime and max_ctime arguments, which will process all observations +in the time frame, is not recommended unless severely restricting the detectors for debugging. Generated results @@ -237,20 +230,28 @@ The Step 1 map-based analysis scripts will generate the following outputs in the * All single maps are packaged in a single hdf file, with detector readout_id as the keys in the h5py file. - 2. Fitted xi-eta focal plane position results saved as ResultSet in ``/path/to/results/map_based_results`` as specified in the Step 1 config file. Script will append 'force_zero_roll' onto the specified results_dir if True in config file. Load ResultSet with keyword 'focal_plane' + 2. Fitted xi-eta focal plane position results saved as ResultSet in ``/path/to/results/map_based_results`` + as specified in the Step 1 config file. Script will append 'force_zero_roll' onto the specified results_dir + if True in config file. Load ResultSet with keyword 'focal_plane' * Contents: ``ResultSet<[dets:readout_id, xi, eta, gamma, R2], N rows>`` -The Step 2 TOD-based analysis scripts will use the map-based results as a starting point and then generate the finalized outputs in the specified directory: +The Step 2 TOD-based analysis scripts will use the map-based results as a starting point + and then generate the finalized outputs in the specified directory: - 1. Fitted xi-eta focal plane position results saved as ResultSet in ``/path/to/results/tod_based_results`` as specified in config file for Step-2. Script will append 'force_zero_roll' onto the specified results_dir if True in config file. Load ResultSet with keyword 'focal_plane' + 1. Fitted xi-eta focal plane position results saved as ResultSet in ``/path/to/results/tod_based_results`` + as specified in config file for Step-2. Script will append 'force_zero_roll' onto the specified results_dir + if True in config file. Load ResultSet with keyword 'focal_plane' * Contents: ``ResultSet<[dets:readout_id, xi, eta, gamma, xi_err, eta_err, R2, redchi2], N rows>`` Configuration Files ``````````````````` -The configuration files to be input as ``configs`` in the command line arguments should have the following arguments as well as any preprocessing steps wished to be taken. Only processing steps that are agnostic of det-match can be used to do initial analyses without formalized metadata. +The configuration files to be input as ``configs`` in the command line should +have the following arguments as well as any preprocessing steps wished to be taken. +Only processing steps that are agnostic of det-match can be used to do +initial analyses without formalized metadata. The parameters in these examples could be used for SAT mid-freq moon observations. @@ -259,14 +260,14 @@ Step 1 Config: .. code-block:: yaml context_file: /path/to/context.yaml - query_tags: ['moon'=1] #(alternatively specify --sso_name in kwargs) + query_tags: ['moon=1'] #(alternatively specify --sso_name in kwargs) optics_config_fn: '/global/cfs/cdirs/sobs/users/elleshaw/process_brightsrc/ufm_to_fp.yaml' single_det_maps_dir: /path/to/results/single_det_maps results_dir: /path/to/results/map_based_results parallel_job: True #For job submission. Parallel across wafers. - wafer_mask_det: 8. #mask around detector to cut TOD when source too far away. + wafer_mask_det: 8. # (degrees) mask around detector to cut TOD when source too far away. res_deg: 0.3 xieta_bs_offset: [0., 0.] #Good to input xieta offset in radians. (!!! for satp2) save_normal_roll: False #false for SAT, true for LAT @@ -308,20 +309,22 @@ The parameters in these examples are used for SAT mid-freq moon observations. .. code-block:: yaml context_file: /path/to/context.yaml - query_tags: ['moon'=1] #(alternatively specify --sso_name in kwargs) + query_tags: ['moon=1'] #(alternatively specify --sso_name in kwargs) optics_config_fn: '/global/cfs/cdirs/sobs/users/elleshaw/process_brightsrc/ufm_to_fp.yaml' fp_hdf_dir: /path/to/results/map_based_results from step 1 config file. - # If force_zero_roll is was True, then append _force_zero_roll to the end. Just make sure it matches where the results from Step 1 are. + # If force_zero_roll is was True, then append _force_zero_roll to the end. + # Just make sure it matches where the results from Step 1 are. result_dir: /path/to/resuls/tod_based_results #Where you want the final Step2 results to show up. parallel_job: True force_zero_roll: True #Results will show up roatated in the xi-eta results as they are on the sky. ds_factor: 40 - mask_deg: 2.5 # size for circular mask around SSO (helps exclude focal plane reflections too) + mask_deg: 2.5 # (degrees) size for circular mask around SSO (helps exclude focal plane reflections too) fit_func_name: 'gaussian2d_nonlin' - max_non_linear_order: 3 #Suggested to use 1 for jupiter or sso's that do not saturate. - fwhm_init_deg: 0.5 # Lower for SATp2 + max_non_linear_order: 3 #Suggested to use 1 for jupiter or sso's + #that do not saturate. + fwhm_init_deg: 0.5 # (degrees) Lower for SATp2 error_estimation_method: 'force_one_redchi2' flag_name_rms_calc: 'around_source' flag_rms_calc_exclusive: False @@ -351,27 +354,28 @@ The parameters in these examples are used for SAT mid-freq moon observations. cutoff: 1.9 width: 0.2 - name: 'source_flags' + source_flags_name: 'source_wide' + save: True calc: - merge: True - max_pix: 10000000000 - source_flags_name: 'source_wide' mask: shape: circle xyr: [0., 0., 5.0] - - name: 'source_flags' - calc: merge: True max_pix: 10000000000 - source_flags_name: 'source_narrow' + - name: 'source_flags' + source_flags_name: 'source_narrow' + save: True + calc: mask: shape: circle xyr: [0., 0., 3.0] - - name: 'reduce_flags' + merge: True + max_pix: 10000000000 + - name: 'combine_flags' process: - flags: ['source_wide', 'source_narrow'] - method: 'except' - wrap: True - new_flag: 'around_source' + flag_labels: ['source_wide.moon', 'source_narrow.moon'] + method: 'except' + total_flags_label: 'around_source' - name: 'flag_turnarounds' process: truncate: True @@ -395,7 +399,7 @@ Example NERSC slurm job submission config file #SBATCH --cpus-per-task=14 #SBATCH --time=00:30:00 - #SBATCH --mem=220G`` #(may require regular queue and up to 400 Gb for extra long observations) + #SBATCH --mem=220G`` #(may need regular queue & up to 400 Gb for long obs) export OMP_NUM_THREADS=1 set -e @@ -411,10 +415,14 @@ Example NERSC slurm job submission config file if (($map)); then echo submitted map job; - srun -n 1 -N 1 -c 14 python3 /path/to/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py $yfile --obs_id=${2} --sso_name="moon"; + srun -n 1 -N 1 -c 14 python3 + /path/to/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py $yfile + --obs_id=${2} --sso_name="moon"; else echo submitted tod job; - srun -n 1 -N 1 -c 14 python3 /path/to/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py $yfile --obs_id=${2} --sso_name="moon"; + srun -n 1 -N 1 -c 14 python3 + /path/to/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py $yfile + --obs_id=${2} --sso_name="moon"; fi @@ -499,6 +507,16 @@ entries mater. det_info: true multi: true + +get_brightsrc_pointing_part2 +---------------------------- +See Part 1 for description + +.. argparse:: + :module: sotodlib.site_pipeline.get_brightsrc_pointing_part2 + :func: get_parser + + update_det_match ------------------ diff --git a/sotodlib/coords/brightsrc_pointing.py b/sotodlib/coords/brightsrc_pointing.py index 8330fce1f..2cb8605d6 100644 --- a/sotodlib/coords/brightsrc_pointing.py +++ b/sotodlib/coords/brightsrc_pointing.py @@ -1,10 +1,12 @@ +# These functions are used for fitting detector positions from bright point sources +# called by site_pipeline.get_brightsrc_pointing_step1 and site_pipeline.get_brightsrc_pointing_step2 import os import re from tqdm import tqdm import numpy as np from scipy import interpolate from scipy.optimize import curve_fit -from joblib import Parallel, delayed +#from joblib import Parallel, delayed from sotodlib import core from sotodlib import coords @@ -20,11 +22,12 @@ def get_planet_trajectory(tod, planet, _split=20, return_model=False): """ - Generate the trajectory of a given planet over a specified time range. + Generate the trajectory in horizon coordinates of a given planet over a specified time range. Parameters: - tod : An axis manager - planet (str): The name of the planet for which to generate the trajectory. + tod : An axis manager containing a timestamps field, which is used to + determine the time range and generate the trajectory. + planet (str): The name of the object for which to generate the trajectory. e.g. "moon" or "saturn" _split (int, optional): Number of points to interpolate the trajectory. Defaults to 20. return_model (bool, optional): If True, returns interpolation functions of az and el. Defaults to False. @@ -34,7 +37,6 @@ def get_planet_trajectory(tod, planet, _split=20, return_model=False): If return_model is False: array: Array of quaternions representing trajectory of the planet at each timestamp. """ - print(planet) timestamps_sparse = np.linspace(tod.timestamps[0], tod.timestamps[-1], _split) planet_az_sparse = np.zeros_like(timestamps_sparse) @@ -136,7 +138,8 @@ def get_wafer_xieta(wafer_slot, optics_config_fn, xieta_bs_offset=(0., 0.), def get_rough_hit_time(tod, wafer_slot, sso_name, circle_r_deg=7.,optics_config_fn=None): """ - Estimate the rough hit time for a axismanager, wafer_slot, and sso_name. + Estimate the rough hit time, which is the amount of time for which the source + is within some distance from the center of a wafer slot. Parameters: tod : An AxisManager object @@ -166,6 +169,8 @@ def make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf, signal='signal', wafer_mask_deg=8., res_deg=0.3, cuts=None,): """ Generate boresight-centered maps from Time-Ordered Data (TOD) for each individual detector. + This script modifies tod.focal_plane and tod.boresight + Parameters: tod : an axismanager object sso_name (str): Name of the planet for which the trajectory is calculated. diff --git a/sotodlib/preprocess/processes.py b/sotodlib/preprocess/processes.py index f6ee863ee..1b77db1da 100644 --- a/sotodlib/preprocess/processes.py +++ b/sotodlib/preprocess/processes.py @@ -1,3 +1,4 @@ +import pdb import numpy as np from operator import attrgetter import copy @@ -1784,7 +1785,7 @@ def calc_and_save(self, aman, proc_aman): source_aman.wrap(source + '_inv', RangesMatrix.ones([aman.dets.count, aman.samps.count]), [(0, 'dets'), (1, 'samps')]) - + self.save(proc_aman, source_aman) return aman, proc_aman @@ -1991,12 +1992,7 @@ def process(self, aman, proc_aman, sim=False): aman.samps.offset + aman.samps.count - trim)) proc_aman.restrict('samps', (proc_aman.samps.offset + trim, proc_aman.samps.offset + proc_aman.samps.count - trim)) - return aman, proc_aman - -class ReduceFlags(_Preprocess): - name = 'reduce_flags' - def process(self, aman, proc_aman, sim=False): - aman.flags.reduce(**self.process_cfgs) + return aman, proc_aman class DetcalNanCuts(_Preprocess): """ @@ -2524,9 +2520,9 @@ def process(self, aman, proc_aman, sim=False): aman['flags'].wrap(self.process_cfgs['total_flags_label'], total_flags) return aman, proc_aman - + class CombineFlags(_Preprocess): - """Do the conbine of relevant flags for mapping + """Do the combination of relevant flags for mapping Saves results for aman under the "flags.[total_flags_label]" field. @@ -2537,8 +2533,11 @@ class CombineFlags(_Preprocess): process: flag_labels: ['glitches.glitch_flags', 'source_flags.jupiter_inv'] total_flags_label: 'glitch_flags' - method: 'union' # You can select a method from ['union', '+', 'intersect', '*']. - #method: ['+', '*'] # Or you can pass individual method for each flags as a list. Lentgh must match the length of flag_labels. + method: 'union' # You can select a method from ['union', '+', 'intersect', '*', 'except', '-']. + #method: ['+', '*'] # Or you can pass individual method for each flags as a list. + # Length of list must match the length of flag_labels. + # If a list, the first method must be '+', as if adding the first flag set to an empty flag set. + # Operations are performed strictly from Left to Right, '*' are not performed first. """ name = "combine_flags" @@ -2548,13 +2547,12 @@ def process(self, aman, proc_aman, sim=False): if isinstance(self.process_cfgs['method'], list): if len(self.process_cfgs['flag_labels']) != len(self.process_cfgs['method']): raise ValueError("The length of method does not match to the length of flag_labels") - elif any(method not in ['+', 'union', '*', 'intersect'] for method in self.process_cfgs['method']): - raise ValueError("The method provided does not match one of '+', '*', 'union', or 'intersect'") - elif self.process_cfgs['method'] in ['+', 'union', '*', 'intersect']: + elif any(method not in ['+', 'union', '*', 'intersect', '-', 'except'] for method in self.process_cfgs['method']): + raise ValueError("One or more methods in list are not valid") + elif self.process_cfgs['method'] in ['+', 'union', '*', 'intersect', '-', 'except']: self.process_cfgs['method'] = ['+'] + (len(self.process_cfgs['flag_labels']) - 1)*[self.process_cfgs['method']] else: - raise ValueError("The method matches neither list nor the one of the ['+', 'union', '*', 'intersect']") - + raise ValueError("The method matches neither list nor the one of the valid operations") total_flags = RangesMatrix.zeros([proc_aman.dets.count, proc_aman.samps.count]) # get an empty flags with shape (Ndets,Nsamps) for i, (method, label) in enumerate(zip(self.process_cfgs['method'], self.process_cfgs['flag_labels'])): _label = attrgetter(label) @@ -2567,6 +2565,8 @@ def process(self, aman, proc_aman, sim=False): total_flags += _label(proc_aman) # The + operator is the union operator in this case elif method in ['*', 'intersect']: total_flags *= _label(proc_aman) # The * operator is the intersect operator in this case + elif method in ['-', 'except']: + total_flags *= ~ _label(proc_aman) # The - operator is the except operator in this case if 'flags' not in aman._fields: from sotodlib.core import FlagManager diff --git a/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py b/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py index e29c3bc76..c77a229d6 100644 --- a/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py +++ b/sotodlib/site_pipeline/get_brightsrc_pointing_step1.py @@ -4,7 +4,7 @@ import argparse import time import glob -from joblib import Parallel, delayed +import logging from sotodlib import core from sotodlib import coords @@ -12,10 +12,14 @@ from sotodlib.coords import brightsrc_pointing as bsp from sotodlib.io import metadata from sotodlib.io.metadata import read_dataset, write_dataset - -from sotodlib.site_pipeline import util +from sotodlib.site_pipeline.utils.pipeline import main_launcher from sotodlib.preprocess import Pipeline -logger = util.init_logger(__name__, 'make_map_based_pointing: ') +from sotodlib.utils.procs_pool import get_exec_env +from sotodlib.site_pipeline.utils.logging import init_logger as sp_init_logger + +logger = logging.getLogger("get_brightsrc_pointing_step1") +if not logger.hasHandlers(): + sp_init_logger("get_brightsrc_pointing_step1") def _get_sso_names_from_tags(ctx, obs_id, candidate_names=['moon', 'jupiter', 'mars', 'saturn']): obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] @@ -109,7 +113,7 @@ def main_one_wafer(configs, obs_id, wafer_slot, sso_name=None, filename=result_filename, force_zero_roll=True, edge_avoidance = edge_avoidance_deg*coords.DEG) - return + return f"Finished processing {obs_id}, {wafer_slot}" def main_one_wafer_dummy(configs, obs_id, wafer_slot, restrict_dets_for_debug=False): if type(configs) == str: @@ -169,15 +173,6 @@ def combine_pointings(pointing_result_files): focal_plane.rows.append((det, val['xi'], val['eta'], val['gamma'], val['R2'])) return focal_plane -def parallel_process_wafer_slot(configs, obs_id, wafer_slot, sso_name, restrict_dets_for_debug): - logger.info(f'Processing {obs_id}, {wafer_slot}') - main_one_wafer(configs=configs, - obs_id=obs_id, - wafer_slot=wafer_slot, - sso_name=sso_name, - restrict_dets_for_debug=restrict_dets_for_debug) - - def main_one_obs(configs, obs_id, sso_name=None, restrict_dets_for_debug=False): if type(configs) == str: @@ -234,17 +229,19 @@ def main_one_obs(configs, obs_id, sso_name=None, n_jobs = -1 logger.info('Entering wafer pool') - Parallel(n_jobs=n_jobs)( - delayed(parallel_process_wafer_slot)( - configs, - obs_id, - wafer_slot, - sso_name, - restrict_dets_for_debug, - ) - for wafer_slot in processed_wafer_slots - ) - logger.info('Exiting wafer pool') + rank, executor, as_completed_callable = get_exec_env(nprocs=n_jobs) + futures = [executor.submit( + main_one_wafer, + configs=configs, + obs_id=obs_id, + wafer_slot=wafer_slot, + sso_name=sso_name, + restrict_dets_for_debug=restrict_dets_for_debug, + ) + for wafer_slot in processed_wafer_slots] + for future in as_completed_callable(futures): + logger.info(future.result()) + else: logger.info('Continuing with serial processing of wafers.') for wafer_slot in processed_wafer_slots: @@ -328,20 +325,20 @@ def main(configs, min_ctime=None, max_ctime=None, update_delay=None, def get_parser(): parser = argparse.ArgumentParser(description="Process TOD data and update pointing") parser.add_argument("configs", type=str, help="Path to the configuration file") - parser.add_argument('--min_ctime', type=int, help="Minimum timestamp for the beginning of an observation list") - parser.add_argument('--max_ctime', type=int, help="Maximum timestamp for the beginning of an observation list") + parser.add_argument('--min-ctime', type=int, help="Minimum timestamp for the beginning of an observation list") + parser.add_argument('--max-ctime', type=int, help="Maximum timestamp for the beginning of an observation list") parser.add_argument('--update-delay', type=int, help="Number of days (unit is days) in the past to start observation list.") - parser.add_argument("--obs_id", type=str, + parser.add_argument("--obs-id", type=str, help="Specific observation obs_id to process. If provided, overrides other filtering parameters.") - parser.add_argument("--wafer_slot", type=str, default=None, + parser.add_argument("--wafer-slot", type=str, default=None, help="Wafer slot to be processed (e.g., 'ws0', 'ws3'). Valid only when obs_id is specified.") - parser.add_argument("--sso_name", type=str, default=None, + parser.add_argument("--sso-name", type=str, default=None, help="Name of solar system object (e.g., 'moon', 'jupiter'). If not specified, get sso_name from observation tags. "\ + "Valid only when obs_id is specified") - parser.add_argument("--restrict_dets_for_debug", type=str, default=False) + parser.add_argument("--restrict-dets-for-debug", type=str, default=False) return parser if __name__ == '__main__': - util.main_launcher(main, get_parser) + main_launcher(main, get_parser) diff --git a/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py index b648010a1..17cb2a7ff 100644 --- a/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py +++ b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py @@ -6,7 +6,7 @@ import time import glob from tqdm import tqdm -from joblib import Parallel, delayed +import logging from scipy.optimize import curve_fit from sotodlib.core import metadata @@ -18,9 +18,14 @@ import so3g from so3g.proj import quat import sotodlib.coords.planets as planets -from sotodlib.site_pipeline import util +from sotodlib.site_pipeline.utils.pipeline import main_launcher from sotodlib.preprocess import Pipeline -logger = util.init_logger(__name__, 'update_pointing: ') +from sotodlib.utils.procs_pool import get_exec_env +from sotodlib.site_pipeline.utils.logging import init_logger as sp_init_logger + +logger = logging.getLogger("get_brightsrc_pointing_step2") +if not logger.hasHandlers(): + sp_init_logger("get_brightsrc_pointing_step2") def _get_sso_names_from_tags(ctx, obs_id, candidate_names=['moon', 'jupiter', 'mars', 'saturn']): obs_tags = ctx.obsdb.get(obs_id, tags=True)['tags'] @@ -71,16 +76,14 @@ def wrapper_gaussian2d_nonlin(xieta, xi0, eta0, fwhm_xi, fwhm_eta, phi, a, *args def wrap_fp_rset(tod, fp_rset): tod.restrict('dets', tod.dets.vals[np.isin(tod.dets.vals, fp_rset['dets:readout_id'])]) + _, ind_tod, ind_rset = core.util.get_coindices(tod.dets.vals, fp_rset['dets:readout_id']) focal_plane = core.AxisManager(tod.dets) focal_plane.wrap_new('xi', shape=('dets', )) focal_plane.wrap_new('eta', shape=('dets', )) - focal_plane.wrap_new('gamma', shape=('dets', )) - - for di, det in enumerate(tod.dets.vals): - di_rset = np.where(fp_rset['dets:readout_id'] == det)[0][0] - focal_plane.xi[di] = fp_rset['xi'][di_rset] - focal_plane.eta[di] = fp_rset['eta'][di_rset] - focal_plane.gamma[di] = fp_rset['gamma'][di_rset] + focal_plane.wrap_new('gamma', shape=('dets', )) + focal_plane.xi[ind_tod] = fp_rset['xi'][ind_rset] + focal_plane.eta[ind_tod] = fp_rset['eta'][ind_rset] + focal_plane.gamma[ind_tod] = fp_rset['gamma'][ind_rset] if 'focal_plane' in tod._fields.keys(): tod.move('focal_plane', None) @@ -367,7 +370,7 @@ def main_one_wafer(configs, obs_id, wafer_slot, sso_name=None, address='focal_plane', overwrite=True) - return + return f"Finished processing {obs_id}, {wafer_slot}" def main_one_wafer_dummy(configs, obs_id, wafer_slot, restrict_dets_for_debug=False): if type(configs) == str: @@ -422,14 +425,6 @@ def combine_pointings(pointing_result_files): focal_plane.rows.append((det, val['xi'], val['eta'], val['gamma'], val['xi_err'], val['eta_err'], val['R2'], val['redchi2'])) return focal_plane -def parallel_process_wafer_slot(configs, obs_id, wafer_slot, sso_name, restrict_dets_for_debug): - logger.info(f'Processing {obs_id}, {wafer_slot}') - main_one_wafer(configs=configs, - obs_id=obs_id, - wafer_slot=wafer_slot, - sso_name=sso_name, - restrict_dets_for_debug=restrict_dets_for_debug) - def main_one_obs(configs, obs_id, sso_name=None, restrict_dets_for_debug=False): if type(configs) == str: @@ -484,16 +479,18 @@ def main_one_obs(configs, obs_id, sso_name=None, n_jobs = int(os.environ.get('SLURM_CPUS_PER_TASK', 1)) except: n_jobs = -1 - Parallel(n_jobs=n_jobs)( - delayed(parallel_process_wafer_slot)( - configs, - obs_id, - wafer_slot, - sso_name, - restrict_dets_for_debug - ) - for wafer_slot in processed_wafer_slots - ) + rank, executor, as_completed_callable = get_exec_env(nprocs=n_jobs) + futures = [executor.submit( + main_one_wafer, + configs=configs, + obs_id=obs_id, + wafer_slot=wafer_slot, + sso_name=sso_name, + restrict_dets_for_debug=restrict_dets_for_debug, + ) + for wafer_slot in processed_wafer_slots] + for future in as_completed_callable(futures): + logger.info(future.result()) else: logger.info('Continuing with serial processing of wafers.') for wafer_slot in processed_wafer_slots: @@ -567,20 +564,20 @@ def main(configs, min_ctime=None, max_ctime=None, update_delay=None, def get_parser(): parser = argparse.ArgumentParser(description="Get updated result of pointings with tod-based results") parser.add_argument("configs", type=str, help="Path to the configuration file") - parser.add_argument('--min_ctime', type=int, help="Minimum timestamp for the beginning of an observation list") - parser.add_argument('--max_ctime', type=int, help="Maximum timestamp for the beginning of an observation list") + parser.add_argument('--min-ctime', type=int, help="Minimum timestamp for the beginning of an observation list") + parser.add_argument('--max-ctime', type=int, help="Maximum timestamp for the beginning of an observation list") parser.add_argument('--update-delay', type=int, help="Number of days (unit is days) in the past to start observation list.") - parser.add_argument("--obs_id", type=str, + parser.add_argument("--obs-id", type=str, help="Specific observation obs_id to process. If provided, overrides other filtering parameters.") - parser.add_argument("--wafer_slot", type=str, default=None, + parser.add_argument("--wafer-slot", type=str, default=None, help="Wafer slot to be processed (e.g., 'ws0', 'ws3'). Valid only when obs_id is specified.") - parser.add_argument("--sso_name", type=str, default=None, + parser.add_argument("--sso-name", type=str, default=None, help="Name of solar system object (e.g., 'moon', 'jupiter'). If not specified, get sso_name from observation tags. "\ + "Valid only when obs_id is specified") - parser.add_argument("--restrict_dets_for_debug", type=str, default=False) + parser.add_argument("--restrict-dets-for-debug", type=str, default=False) return parser if __name__ == '__main__': - util.main_launcher(main, get_parser) + main_launcher(main, get_parser) diff --git a/sotodlib/tod_ops/sub_polyf.py b/sotodlib/tod_ops/sub_polyf.py index 57bf35a31..66aeec36c 100644 --- a/sotodlib/tod_ops/sub_polyf.py +++ b/sotodlib/tod_ops/sub_polyf.py @@ -44,6 +44,7 @@ def subscan_polyfilter(aman, degree, signal_name="signal", exclude_turnarounds=F """ if method not in ["polyfit", "legendre"] : raise ValueError("Only polyfit and legendre are acceptable.") + if exclude_turnarounds: if ("left_scan" not in aman.flags) or ("turnarounds" not in aman.flags): logger.warning('aman does not have left/right scan or turnarounds flag. `sotodlib.flags.get_turnaround_flags` will be ran with default parameters') @@ -113,6 +114,7 @@ def subscan_polyfilter(aman, degree, signal_name="signal", exclude_turnarounds=F elif method == "legendre": degree_corr = degree + 1 + time = np.copy(aman["timestamps"]) for start, end in subscan_indices: @@ -123,6 +125,7 @@ def subscan_polyfilter(aman, degree, signal_name="signal", exclude_turnarounds=F # Get each subscan to be filtered tod_mat = copy.deepcopy(signal[:, start:end]) + # Scale time range into [-1,1] x = np.linspace(-1, 1, tod_mat.shape[1]) dx = np.mean(np.diff(x)) @@ -163,13 +166,15 @@ def subscan_polyfilter(aman, degree, signal_name="signal", exclude_turnarounds=F tod_mat[idet,msk_indx] = interped else : pass - + + means = np.mean(tod_mat, axis=1)[:, np.newaxis] tod_mat -= means # Make model to be subtracted coeffs = np.dot(arr_legendre, tod_mat.T) model = np.dot((coeffs/norm_vector[:, np.newaxis]).T,arr_legendre)*dx + model += means signal[:,start:end] -= model From f3befaae72272685532f0bc0a1d48ef7bb8507f8 Mon Sep 17 00:00:00 2001 From: Elle Shaw Date: Tue, 24 Feb 2026 16:37:03 -0800 Subject: [PATCH 23/27] Fixed issues with lonlat quat rotations and hacky negative xi solutions --- sotodlib/coords/brightsrc_pointing.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/sotodlib/coords/brightsrc_pointing.py b/sotodlib/coords/brightsrc_pointing.py index 2cb8605d6..0ed141e21 100644 --- a/sotodlib/coords/brightsrc_pointing.py +++ b/sotodlib/coords/brightsrc_pointing.py @@ -52,7 +52,7 @@ def get_planet_trajectory(tod, planet, _split=20, return_model=False): else: planet_az = planet_az_func(tod.timestamps) planet_el = planet_el_func(tod.timestamps) - q_planet = quat.rotation_lonlat(planet_az, planet_el) + q_planet = quat.rotation_lonlat(-1 * planet_az, planet_el) return q_planet def get_wafer_centered_sight(tod=None, planet=None, q_planet=None, q_bs=None, q_wafer=None): @@ -73,18 +73,18 @@ def get_wafer_centered_sight(tod=None, planet=None, q_planet=None, q_bs=None, q_ Returns: Sightline vector for the planet trajectory centered on the center of the wafer. """ + #breakpoint() if q_planet is None: q_planet = get_planet_trajectory(tod, planet) if q_bs is None: - q_bs = quat.rotation_lonlat(tod.boresight.az, tod.boresight.el) + q_bs = quat.rotation_lonlat(-1 * tod.boresight.az, tod.boresight.el) if q_wafer is None: q_wafer = quat.rotation_xieta(np.nanmedian(tod.focal_plane.xi), np.nanmedian(tod.focal_plane.eta)) xi_wafer, eta_wafer, _ = quat.decompose_xieta(q_wafer) - q_wafer_f = quat.rotation_xieta(-xi_wafer, eta_wafer) z_to_x = quat.rotation_lonlat(0, 0) - sight = z_to_x * ~(q_bs * q_wafer_f) * q_planet + sight = z_to_x * ~(q_bs * q_wafer) * q_planet return sight def get_wafer_xieta(wafer_slot, optics_config_fn, xieta_bs_offset=(0., 0.), @@ -151,7 +151,7 @@ def get_rough_hit_time(tod, wafer_slot, sso_name, circle_r_deg=7.,optics_config_ Returns: float: Estimated rough hit time within the circular region around the wafer center. """ - q_bs = quat.rotation_lonlat(tod.boresight.az, tod.boresight.el) + q_bs = quat.rotation_lonlat(-1 * tod.boresight.az, tod.boresight.el) q_planet = get_planet_trajectory(tod, sso_name) xi_wafer, eta_wafer = get_wafer_xieta(wafer_slot, optics_config_fn=optics_config_fn, roll_bs_offset=np.median(tod.boresight.roll), wrap_to_tod=False) @@ -186,8 +186,9 @@ def make_wafer_centered_maps(tod, sso_name, optics_config_fn, map_hdf, Returns: None """ + #breakpoint() q_planet = get_planet_trajectory(tod, sso_name) - q_bs = quat.rotation_lonlat(tod.boresight.az, tod.boresight.el) + q_bs = quat.rotation_lonlat(-1 * tod.boresight.az, tod.boresight.el) if roll_bs_offset is None: roll_bs_offset = np.mean(tod.boresight.roll) @@ -419,7 +420,7 @@ def map_to_xieta(mT, edge_avoidance=1.0*coords.DEG, edge_check='nan', (np.inf, beam_sigma_init*5, np.inf),), max_nfev = 1000000) R2 = 1 - np.sum((_z - _gauss1d(_r, *popt))**2)/np.sum((_z - np.mean(_z))**2) - xi_det, eta_det, R2_det = -xi_peak, eta_peak, R2 + xi_det, eta_det, R2_det = xi_peak, eta_peak, R2 else: xi_det, eta_det, R2_det = np.nan, np.nan, np.nan return xi_det, eta_det, R2_det From 76ae8dfe5719ce2372e21ccc8d53a4c1e7d004a9 Mon Sep 17 00:00:00 2001 From: Elle Shaw Date: Wed, 4 Mar 2026 13:06:05 -0800 Subject: [PATCH 24/27] Find boresight azimuth and elevation during detector crossings. --- .../get_brightsrc_pointing_step2.py | 22 +++++++++++++------ 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py index 17cb2a7ff..34d82839c 100644 --- a/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py +++ b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py @@ -8,6 +8,10 @@ from tqdm import tqdm import logging +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt + from scipy.optimize import curve_fit from sotodlib.core import metadata from sotodlib.io.metadata import read_dataset, write_dataset @@ -154,7 +158,7 @@ def update_xieta(tod, # if focal_plane result is specified, use the information as a prior if fp_hdf_file is not None: wrap_fp_from_hdf(tod, fp_hdf_file) - + # set dets without focal_plane info to have (xi, eta, gamma) = (0, 0, 0), just to avoid error xieta_isnan = (np.isnan(tod.focal_plane.xi)) | (np.isnan(tod.focal_plane.eta)) gamma_isnan = np.isnan(tod.focal_plane.gamma) @@ -213,14 +217,17 @@ def update_xieta(tod, xieta_dict = {} for di, det in enumerate(tqdm(tod.dets.vals)): mask_di = source_flags_ds[di] + bs_az = np.nanmedian(tod.boresight.az[mask_ds][mask_di]) + bs_el = np.nanmedian(tod.boresight.el[mask_ds][mask_di]) + if np.any([xieta_isnan[di], np.all(mask_di==False), tod.rms[di]==0.]): xieta_dict[det] = {'xi': np.nan, 'eta': np.nan, 'xi_err': np.nan, 'eta_err': np.nan, - 'R2': np.nan, 'redchi2': np.nan} + 'R2': np.nan, 'redchi2': np.nan, 'az': np.nan, 'el': np.nan} else: ts = ts_ds[mask_di] d1_unix = np.median(ts) - xieta_det = np.array([tod.focal_plane.xi[di], tod.focal_plane.eta[di]]) + xieta_det = np.array([tod.focal_plane.xi[di], tod.focal_plane.eta[di]]) q_det = so3g.proj.quat.rotation_xieta(xieta_det[0], xieta_det[1]) planet = planets.SlowSource.for_named_source(sso_name, d1_unix * 1.) ra0, dec0 = planet.pos(d1_unix) @@ -232,7 +239,7 @@ def update_xieta(tod, xieta_src = xieta_src[:, mask_di] sig = sig_ds[di][mask_di] ptp_val = np.ptp(np.percentile(sig, [0.1, 99.9])) - + if fit_func_name == 'gaussian2d_nonlin': p0 = np.array([0., 0., fwhm_init_deg*coords.DEG, fwhm_init_deg*coords.DEG, 0., ptp_val]) bounds = np.array( @@ -270,16 +277,17 @@ def update_xieta(tod, xieta_det += np.array([xi_opt, eta_opt]) xieta_dict[det] = {'xi': xieta_det[0], 'eta': xieta_det[1], 'xi_err': xi_err, 'eta_err': eta_err, - 'R2': R2, 'redchi2': redchi2} + 'R2': R2, 'redchi2': redchi2, 'az' : bs_az, 'el': bs_el} except RuntimeError: xieta_dict[det] = {'xi': np.nan, 'eta': np.nan, 'xi_err': np.nan, 'eta_err': np.nan, - 'R2': np.nan, 'redchi2': np.nan} + 'R2': np.nan, 'redchi2': np.nan, 'az': np.nan, 'el': np.nan} - focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', 'xi_err', 'eta_err', 'R2', 'redchi2']) + focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', 'xi_err', 'eta_err', 'R2', 'redchi2', 'az', 'el']) for det in tod.dets.vals: focal_plane.rows.append((det, xieta_dict[det]['xi'], xieta_dict[det]['eta'], 0., xieta_dict[det]['xi_err'], xieta_dict[det]['eta_err'], xieta_dict[det]['R2'], xieta_dict[det]['redchi2'], + xieta_dict[det]['az'], xieta_dict[det]['el'], )) return focal_plane From 81ac88696a15da661398b7b0b16a444e7081cc2a Mon Sep 17 00:00:00 2001 From: Elle Shaw Date: Tue, 17 Mar 2026 20:26:47 -0700 Subject: [PATCH 25/27] propagate addition of az and el and roll to result set output in dummy wafer and final full wafer products --- .../get_brightsrc_pointing_step2.py | 24 +++++++++++-------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py index 34d82839c..ef1a77853 100644 --- a/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py +++ b/sotodlib/site_pipeline/get_brightsrc_pointing_step2.py @@ -183,7 +183,7 @@ def update_xieta(tod, max_pix=1e10, wrap='source', mask={'shape':'circle', 'xyr':[0.,0.,mask_deg]}) - + # restrict data to duration when at least one detector hit the source summed_flag = np.sum(tod.flags['source'].mask()[~xieta_isnan], axis=0).astype('bool') idx_hit = np.where(summed_flag)[0] @@ -219,10 +219,11 @@ def update_xieta(tod, mask_di = source_flags_ds[di] bs_az = np.nanmedian(tod.boresight.az[mask_ds][mask_di]) bs_el = np.nanmedian(tod.boresight.el[mask_ds][mask_di]) + bs_roll = np.nanmedian(tod.boresight.roll[mask_ds][mask_di]) if np.any([xieta_isnan[di], np.all(mask_di==False), tod.rms[di]==0.]): xieta_dict[det] = {'xi': np.nan, 'eta': np.nan, 'xi_err': np.nan, 'eta_err': np.nan, - 'R2': np.nan, 'redchi2': np.nan, 'az': np.nan, 'el': np.nan} + 'R2': np.nan, 'redchi2': np.nan, 'az': np.nan, 'el': np.nan, 'roll': np.nan} else: ts = ts_ds[mask_di] d1_unix = np.median(ts) @@ -277,17 +278,17 @@ def update_xieta(tod, xieta_det += np.array([xi_opt, eta_opt]) xieta_dict[det] = {'xi': xieta_det[0], 'eta': xieta_det[1], 'xi_err': xi_err, 'eta_err': eta_err, - 'R2': R2, 'redchi2': redchi2, 'az' : bs_az, 'el': bs_el} + 'R2': R2, 'redchi2': redchi2, 'az' : bs_az, 'el': bs_el, 'roll': bs_roll} except RuntimeError: xieta_dict[det] = {'xi': np.nan, 'eta': np.nan, 'xi_err': np.nan, 'eta_err': np.nan, - 'R2': np.nan, 'redchi2': np.nan, 'az': np.nan, 'el': np.nan} + 'R2': np.nan, 'redchi2': np.nan, 'az': np.nan, 'el': np.nan, 'roll': np.nan} - focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', 'xi_err', 'eta_err', 'R2', 'redchi2', 'az', 'el']) + focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', 'xi_err', 'eta_err', 'R2', 'redchi2', 'az', 'el', 'roll']) for det in tod.dets.vals: focal_plane.rows.append((det, xieta_dict[det]['xi'], xieta_dict[det]['eta'], 0., xieta_dict[det]['xi_err'], xieta_dict[det]['eta_err'], xieta_dict[det]['R2'], xieta_dict[det]['redchi2'], - xieta_dict[det]['az'], xieta_dict[det]['el'], + xieta_dict[det]['az'], xieta_dict[det]['el'], xieta_dict[det]['roll'], )) return focal_plane @@ -401,9 +402,9 @@ def main_one_wafer_dummy(configs, obs_id, wafer_slot, restrict_dets_for_debug=Fa result_filename = f'focal_plane_{obs_id}_{wafer_slot}.hdf' fp_rset_dummy = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', - 'xi_err', 'eta_err', 'R2', 'redchi2']) + 'xi_err', 'eta_err', 'R2', 'redchi2', 'az', 'el', 'roll']) for det in meta.dets.vals: - fp_rset_dummy.rows.append((det, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan)) + fp_rset_dummy.rows.append((det, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan)) os.makedirs(result_dir, exist_ok=True) write_dataset(fp_rset_dummy, @@ -426,11 +427,14 @@ def combine_pointings(pointing_result_files): combined_dict[row['dets:readout_id']]['eta_err'] = row['eta_err'] combined_dict[row['dets:readout_id']]['R2'] = row['R2'] combined_dict[row['dets:readout_id']]['redchi2'] = row['redchi2'] + combined_dict[row['dets:readout_id']]['az'] = row['az'] + combined_dict[row['dets:readout_id']]['el'] = row['el'] + combined_dict[row['dets:readout_id']]['roll'] = row['roll'] - focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', 'xi_err', 'eta_err', 'R2', 'redchi2']) + focal_plane = metadata.ResultSet(keys=['dets:readout_id', 'xi', 'eta', 'gamma', 'xi_err', 'eta_err', 'R2', 'redchi2', 'az', 'el', 'roll']) for det, val in combined_dict.items(): - focal_plane.rows.append((det, val['xi'], val['eta'], val['gamma'], val['xi_err'], val['eta_err'], val['R2'], val['redchi2'])) + focal_plane.rows.append((det, val['xi'], val['eta'], val['gamma'], val['xi_err'], val['eta_err'], val['R2'], val['redchi2'],val['az'], val['el'], val['roll'])) return focal_plane def main_one_obs(configs, obs_id, sso_name=None, From 5879bf3323232c6e5a3a656fb59f1749c386db7c Mon Sep 17 00:00:00 2001 From: Elle Shaw Date: Wed, 25 Mar 2026 08:46:49 -0700 Subject: [PATCH 26/27] Added ability to pass az el and roll values from get_brightsrc datasets through to the outputs of finalize focal plane. Useful for per-obs usecase, and when loaded with Receiver.load the values will be under det_boresight --- docs/site_pipeline.rst | 11 +++++-- sotodlib/coords/fp_containers.py | 33 +++++++++++++++++-- .../site_pipeline/finalize_focal_plane.py | 12 ++++--- 3 files changed, 46 insertions(+), 10 deletions(-) diff --git a/docs/site_pipeline.rst b/docs/site_pipeline.rst index ea9ea5537..89f728908 100644 --- a/docs/site_pipeline.rst +++ b/docs/site_pipeline.rst @@ -242,9 +242,10 @@ The Step 2 TOD-based analysis scripts will use the map-based results as a starti 1. Fitted xi-eta focal plane position results saved as ResultSet in ``/path/to/results/tod_based_results`` as specified in config file for Step-2. Script will append 'force_zero_roll' onto the specified results_dir - if True in config file. Load ResultSet with keyword 'focal_plane' + if True in config file. Load ResultSet with keyword 'focal_plane'. + The median boresight values from small time range the source was visible to each detector is included. - * Contents: ``ResultSet<[dets:readout_id, xi, eta, gamma, xi_err, eta_err, R2, redchi2], N rows>`` + * Contents: ``ResultSet<[dets:readout_id, xi, eta, gamma, xi_err, eta_err, R2, redchi2, az, el, roll], N rows>`` Configuration Files ``````````````````` @@ -762,6 +763,10 @@ The ``focal_plane_full`` dataset contains nine columns: - ``eta_m``: The measured eta in radians - ``gamma_m``: The measured gamma in radians. - ``weights``: The average weights of the measurements for this det. +- ``r2``: The fit weight passed in from the get_brightsrc_pointing dataset +- ``az``: The median Az value in radians from source-detector crossing +- ``el``: The median El value in radians from source-detector crossing +- ``roll``: The median Roll value in radians from source-detector crossing - ``n_point``: The number of pointing fits used for the det. - ``n_gamma``: The number of gamma fits used for this det. @@ -802,7 +807,7 @@ always be ``(1, 1, 1)`` and ``shear`` will be ``0``. ``finalize_focal_plane`` will also output a ``ManifestDb`` as a file called ``db.sqlite`` in the output directory. By default this will be indexed by ``stream_id`` and ``obs:timestamp`` and will point to the ``focal_plane`` dataset. -If you are running in ``per_obs`` mode then it wirbe indexed by ``obs_id`` and will point +If you are running in ``per_obs`` mode then it will be indexed by ``obs_id`` and will point to results associated with data observation. Be warned that in this case there will only be entries for observations with pointing fits, so design your context accordingly. diff --git a/sotodlib/coords/fp_containers.py b/sotodlib/coords/fp_containers.py index fb089c631..976e7e36d 100644 --- a/sotodlib/coords/fp_containers.py +++ b/sotodlib/coords/fp_containers.py @@ -141,6 +141,7 @@ class FocalPlane: id_strs: NDArray[np.str_] # (ndet,) avg_fp: NDArray[np.floating] # (ndim, ndet) weights: NDArray[np.floating] # (ndet,) + det_boresight: NDArray[np.floating] # (ndim, ndet) transformed: NDArray[np.floating] # (ndet, ndim) center: NDArray[np.floating] # (1, ndim) center_transformed: NDArray[np.floating] # (1, ndim) @@ -219,6 +220,7 @@ def empty(cls, template, stream_id, wafer_slot, n_aman, config=""): tot_weight = np.zeros((len(template.det_ids), 2)) avg_fp = np.full_like(template.fp, np.nan) weight = np.zeros((len(template.det_ids), 2)) + det_boresight = np.zeros((len(template.det_ids), 3)) + np.nan # az, el, roll transformed = template.fp.copy() center = template.center.copy() center_transformed = template.center.copy() @@ -232,6 +234,7 @@ def empty(cls, template, stream_id, wafer_slot, n_aman, config=""): template.id_strs, avg_fp, weight, + det_boresight, transformed, center, center_transformed, @@ -260,8 +263,17 @@ def map_by_det_id(self, aman): xi = aman.pointing.xi[msk][srt][mapping] eta = aman.pointing.eta[msk][srt][mapping] r2 = np.nan + np.zeros_like(eta) + az = np.nan + np.zeros_like(eta) + el = np.nan + np.zeros_like(eta) + roll = np.nan + np.zeros_like(eta) if "R2" in aman.pointing: r2 = aman.pointing.R2[msk][srt][mapping] + if "az" in aman.pointing: + az = aman.pointing.az[msk][srt][mapping] + if "el" in aman.pointing: + el = aman.pointing.el[msk][srt][mapping] + if "roll" in aman.pointing: + roll = aman.pointing.roll[msk][srt][mapping] if "polarization" in aman: # name of field just a placeholder for now gamma = aman.polarization.polang[msk][srt][mapping] @@ -270,14 +282,16 @@ def map_by_det_id(self, aman): else: gamma = np.full(len(xi), np.nan) fp = np.column_stack((xi, eta, gamma)) - return fp, r2, template_msk + det_boresight = np.column_stack((az, el, roll)) + return fp, r2, det_boresight, template_msk - def add_fp(self, i, fp, weights, template_msk): - if self.full_fp is None or self.tot_weight is None: + def add_fp(self, i, fp, weights, det_boresight, template_msk): + if self.full_fp is None or self.tot_weight is None or self.det_boresight is None: raise ValueError("full_fp or tot_weight not initialized") self.full_fp[template_msk, :, i] = fp * weights[:, 0][..., None] weights = np.nan_to_num(weights) self.tot_weight[template_msk] += weights + self.det_boresight[template_msk, :] = det_boresight def save(self, f, db_info, group): logger.info("Saving %s", self.stream_id) @@ -323,6 +337,9 @@ def save(self, f, db_info, group): ("gamma_m", np.float32), ("weights", np.float32), ("r2", np.float32), + ("az", np.float32), + ("el", np.float32), + ("roll", np.float32), ("n_point", np.int8), ("n_gamma", np.int8), ] @@ -333,6 +350,7 @@ def save(self, f, db_info, group): *(self.transformed.T), *(self.avg_fp.T), *(self.weights.T), + *(self.det_boresight.T), self.n_point, self.n_gamma, ), @@ -378,6 +396,14 @@ def load(cls, group, include_cm=None): np.array(fp_full["gamma_m"]), ) ) + + det_boresight = np.column_stack( + ( + np.array(fp_full["az"]), + np.array(fp_full["el"]), + np.array(fp_full["roll"]), + ) + ) # For backwards compatibility weights = np.array(fp_full["weights"]) if "r2" in fp_full.keys: @@ -424,6 +450,7 @@ def load(cls, group, include_cm=None): np.array(id_strs), avg_fp, np.array(weights), + det_boresight, transformed, center, center_transformed, diff --git a/sotodlib/site_pipeline/finalize_focal_plane.py b/sotodlib/site_pipeline/finalize_focal_plane.py index fe19a6e97..8df8a9c59 100644 --- a/sotodlib/site_pipeline/finalize_focal_plane.py +++ b/sotodlib/site_pipeline/finalize_focal_plane.py @@ -6,6 +6,10 @@ from importlib import import_module from typing import List, Optional +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt + import git import h5py import megham.transform as mt @@ -399,7 +403,7 @@ def _mk_pointing_config(telescope_flavor, tube_slot, wafer_slot, config): def _restrict_inliers(aman, focal_plane): # TODO: Use gamma as well # Map to template - fp, _, template_msk = focal_plane.map_by_det_id(aman) + fp, _, _, template_msk = focal_plane.map_by_det_id(aman) fp = fp[:, :2] inliers = np.ones(len(fp), dtype=bool) @@ -651,7 +655,7 @@ def main(): plot_dir = os.path.join(plot_dir_base, str(config["start_time"])) os.makedirs(plot_dir, exist_ok=True) logger.info("Working on batch containing: %s", str(obs_ids)) - + # Setup db and Receiver db, base, group = _create_db( dbpath, @@ -771,7 +775,7 @@ def main(): _restrict_inliers(aman, focal_plane) # Mapping to template - fp, r2, template_msk = focal_plane.map_by_det_id(aman) + fp, r2, det_boresight, template_msk = focal_plane.map_by_det_id(aman) focal_plane.template.add_wafer_info(aman, template_msk) # Try an initial alignment and get weights @@ -808,7 +812,7 @@ def main(): # Store weighted values weights = np.column_stack((weights, r2)) - focal_plane.add_fp(i, fp, weights, template_msk) + focal_plane.add_fp(i, fp, weights, det_boresight, template_msk) n_obs += 1 From 3dbb7edaf4cb6100898f230994ab95476d0e0096 Mon Sep 17 00:00:00 2001 From: Elle Shaw Date: Thu, 26 Mar 2026 15:46:02 -0700 Subject: [PATCH 27/27] Pre-fill detector pointing with central values from obs_info if per detector pointing not included in inputs --- sotodlib/coords/fp_containers.py | 6 +++--- sotodlib/site_pipeline/finalize_focal_plane.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/sotodlib/coords/fp_containers.py b/sotodlib/coords/fp_containers.py index 976e7e36d..3e8a4ab92 100644 --- a/sotodlib/coords/fp_containers.py +++ b/sotodlib/coords/fp_containers.py @@ -263,9 +263,9 @@ def map_by_det_id(self, aman): xi = aman.pointing.xi[msk][srt][mapping] eta = aman.pointing.eta[msk][srt][mapping] r2 = np.nan + np.zeros_like(eta) - az = np.nan + np.zeros_like(eta) - el = np.nan + np.zeros_like(eta) - roll = np.nan + np.zeros_like(eta) + az = np.deg2rad(aman.obs_info.az_center) * np.ones_like(eta) + el = np.deg2rad(aman.obs_info.el_center) * np.ones_like(eta) + roll = np.deg2rad(aman.obs_info.roll_center) * np.ones_like(eta) if "R2" in aman.pointing: r2 = aman.pointing.R2[msk][srt][mapping] if "az" in aman.pointing: diff --git a/sotodlib/site_pipeline/finalize_focal_plane.py b/sotodlib/site_pipeline/finalize_focal_plane.py index 8b8f0d826..0a0209cfb 100644 --- a/sotodlib/site_pipeline/finalize_focal_plane.py +++ b/sotodlib/site_pipeline/finalize_focal_plane.py @@ -663,7 +663,7 @@ def main(): plot_dir = os.path.join(plot_dir_base, str(config["start_time"])) os.makedirs(plot_dir, exist_ok=True) logger.info("Working on batch containing: %s", str(obs_ids)) - + # Setup db and Receiver db, base, group = _create_db( dbpath,