From 05f55ea7aa7390600f6254e5091faca5a3fe03b2 Mon Sep 17 00:00:00 2001 From: ykyohei <38639108+ykyohei@users.noreply.github.com> Date: Mon, 13 Apr 2026 20:54:37 -0400 Subject: [PATCH 1/4] Fix az mask of azss legendre fit --- sotodlib/tod_ops/azss.py | 97 +++++++++++++++++++++++++++++----------- tests/test_azss.py | 38 +++++++++++++--- 2 files changed, 105 insertions(+), 30 deletions(-) diff --git a/sotodlib/tod_ops/azss.py b/sotodlib/tod_ops/azss.py index 5cb94267f..05b01bfb7 100644 --- a/sotodlib/tod_ops/azss.py +++ b/sotodlib/tod_ops/azss.py @@ -148,28 +148,75 @@ def fit_azss(az, azss_stats, max_mode, fit_range=None, overwrite=False): Model fit for each detector size ndets x n_samps """ bin_width = azss_stats.binned_az[1] - azss_stats.binned_az[0] - m = ~np.isnan(azss_stats.binned_signal[0]) # masks bins without counts - if np.count_nonzero(m) < max_mode + 1: - raise ValueError('Number of valid bins is smaller than mode of Legendre function') - if fit_range is None: - az_min = np.min(azss_stats.binned_az[m]) - bin_width / 2 - az_max = np.max(azss_stats.binned_az[m]) + bin_width / 2 - else: - az_min, az_max = fit_range[0], fit_range[1] + model = np.zeros((azss_stats.dets.count, len(az))) + + m = ~np.isnan(azss_stats.binned_signal) + # Number of valid bins cannot be smaller than mode of Legendre function + valid_dets = np.sum(m, axis=1) >= max_mode + 1 + if 'bad_dets' in azss_stats: + valid_dets = np.logical_and(valid_dets, ~azss_stats['bad_dets']) - x_legendre = (2 * az - (az_min+az_max)) / (az_max - az_min) - x_legendre_bin_centers = (2 * azss_stats.binned_az - (az_min+az_max)) / (az_max - az_min) - x_legendre_bin_centers = np.where(~m, np.nan, x_legendre_bin_centers) - if ('coeffs' in azss_stats) and not overwrite: - return L.legval(x_legendre, azss_stats.coeffs.T) + if sum(valid_dets) == 0: + logger.info('All the detectors have low az coverage and cannot make model') + return model - coeffs = L.legfit(x_legendre_bin_centers[m], azss_stats.binned_signal[:, m].T, max_mode) - coeffs = coeffs.T - binned_model = L.legval(x_legendre_bin_centers, coeffs.T) - binned_model = np.where(~m, np.nan, binned_model) - sum_of_squares = np.sum(((azss_stats.binned_signal[:, m] - binned_model[:, m])**2), axis=-1) - redchi2s = sum_of_squares/azss_stats.uniform_binned_signal_sigma**2 / (len(x_legendre_bin_centers[m]) - max_mode - 1) + mask = ~np.isnan(azss_stats.binned_signal[valid_dets, :]) + is_uniform = np.all(mask == mask[0, :]) + + coeffs = np.zeros((azss_stats.dets.count, max_mode + 1)) + binned_model = np.full((azss_stats.dets.count, azss_stats.bin_az_samps.count), np.nan) + sum_of_squares = np.full(azss_stats.dets.count, np.nan) + redchi2s = np.full(azss_stats.dets.count, np.nan) + if fit_range is not None: + az_min, az_max = fit_range[0], fit_range[1] + x_legendre = (2 * az - (az_min+az_max)) / (az_max - az_min) + x_legendre_bin_centers = (2 * azss_stats.binned_az - (az_min+az_max)) / (az_max - az_min) + x_legendre_bin_centers = np.where(~m, np.nan, x_legendre_bin_centers) + + if is_uniform: + m = mask[0, :] + if fit_range is None: + az_min = np.min(azss_stats.binned_az[m]) - bin_width / 2 + az_max = np.max(azss_stats.binned_az[m]) + bin_width / 2 + x_legendre = (2 * az - (az_min+az_max)) / (az_max - az_min) + x_legendre_bin_centers = (2 * azss_stats.binned_az - (az_min+az_max)) / (az_max - az_min) + x_legendre_bin_centers = np.where(~m, np.nan, x_legendre_bin_centers) + + if ('coeffs' in azss_stats) and not overwrite: + return L.legval(x_legendre, azss_stats.coeffs.T) + + coeffs[valid_dets, :] = L.legfit(x_legendre_bin_centers[m], azss_stats.binned_signal[:, m][valid_dets, :].T, max_mode).T + assert ~np.any(np.isnan(coeffs)) + model = L.legval(x_legendre, coeffs.T) + binned_model = L.legval(x_legendre_bin_centers, coeffs.T) + binned_model = np.where(~m, np.nan, binned_model) + sum_of_squares = np.sum(((azss_stats.binned_signal[:, m] - binned_model[:, m])**2), axis=-1) + redchi2s = sum_of_squares/azss_stats.uniform_binned_signal_sigma**2 / (len(x_legendre_bin_centers[m]) - max_mode - 1) + + else: + for i, m, binned_signal in zip(np.where(valid_dets)[0], mask, azss_stats.binned_signal[valid_dets, :]): + if fit_range is None: + az_min = np.min(azss_stats.binned_az[m]) - bin_width / 2 + az_max = np.max(azss_stats.binned_az[m]) + bin_width / 2 + x_legendre = (2 * az - (az_min+az_max)) / (az_max - az_min) + x_legendre_bin_centers = (2 * azss_stats.binned_az - (az_min+az_max)) / (az_max - az_min) + x_legendre_bin_centers = np.where(~m, np.nan, x_legendre_bin_centers) + if ('coeffs' in azss_stats) and not overwrite: + model[i, :] = L.legval(x_legendre, azss_stats.coeffs[:, i]) + else: + _coeffs = L.legfit(x_legendre_bin_centers[m], azss_stats.binned_signal[:, m][i, :].T, max_mode) + _binned_model = L.legval(x_legendre_bin_centers, _coeffs) + _binned_model = np.where(~m, np.nan, _binned_model) + _sum_of_squares = np.sum(((azss_stats.binned_signal[:, m][i, :] - _binned_model[m])**2), axis=-1) + _redchi2s = _sum_of_squares/azss_stats.uniform_binned_signal_sigma[i]**2 / (len(x_legendre_bin_centers[m]) - max_mode - 1) + coeffs[i, :] = _coeffs + binned_model[i, :] = _binned_model + sum_of_squares[i] = _sum_of_squares + redchi2s[i] = _redchi2s + model[i, :] = L.legval(x_legendre, _coeffs) + if ('coeffs' in azss_stats) and not overwrite: + return model mode_names = np.array([f'legendre{mode}' for mode in range(max_mode + 1)], dtype=' 0 if 'bad_dets' in azss_stats: - valid_dets = ~azss_stats['bad_dets'] - else: - valid_dets = np.ones(aman.dets.count, dtype=bool) + valid_dets = np.logical_and(valid_dets, ~azss_stats['bad_dets']) if sum(valid_dets) == 0: logger.info('All the detectors have low az coverage and cannot make model') return model @@ -380,9 +427,9 @@ def get_azss_model(aman, azss_stats, az=None, method='interpolate', f_template = interp1d(azss_stats.binned_az[m], azss_stats.binned_signal[:, m][valid_dets, :], fill_value='extrapolate') model[valid_dets, :] = f_template(az) else: - for i, (m, binned_signal) in enumerate(zip(mask, azss_stats.binned_signal[valid_dets, :])): + for i, m, binned_signal in zip(np.where(valid_dets)[0], mask, azss_stats.binned_signal[valid_dets, :]): f_template = interp1d(azss_stats.binned_az[m], binned_signal[m], fill_value='extrapolate') - model[valid_dets, :][i] = f_template(az) + model[i, :] = f_template(az) if np.any(~np.isfinite(model)): logger.warning('azss model has nan. set zero to nan but this may make glitch') diff --git a/tests/test_azss.py b/tests/test_azss.py index 36b9cf364..33b67fdb1 100644 --- a/tests/test_azss.py +++ b/tests/test_azss.py @@ -52,7 +52,7 @@ def get_scan(n_scans=33, scan_accel=0.025, scanrate=0.025, def make_fake_azss_tod(max_mode=20, noise_amp=1, n_scans=10, - ndets=2, input_coeffs=None): + ndets=3, input_coeffs=None): """ Makes an axis manager with azimuth synchronous signal in it, populated via legendre polynomials plus gaussian noise. @@ -94,16 +94,18 @@ def get_coeff_metric(tod): """ Evaluates fit is working by comparing coefficients in to out. """ - print(tod.input_coeffs[0]) - print(tod.azss_stats.coeffs[0]) + print(tod.input_coeffs) + print(tod.azss_stats.coeffs) outmetric_num = (tod.azss_stats.coeffs - tod.input_coeffs)**2 outmetric_denom = (tod.input_coeffs)**2 return np.median(100*(outmetric_num/outmetric_denom)) class AzssTest(unittest.TestCase): - "Test the Azimuth Synchronous Signal fitting functions" def test_fit(self): + """ + Test the Azimuth Synchronous Signal fitting functions + """ max_mode = 10 tod = make_fake_azss_tod(noise_amp=0, n_scans=50, max_mode=max_mode) azss_stats, model_sig_tod = azss.get_azss( @@ -119,6 +121,30 @@ def test_fit(self): self.assertTrue(~np.any(np.isnan(tod.signal))) self.assertTrue(np.std(tod.signal) > np.std(tod.signal - model_sig_tod)) + def test_fit_with_flags(self): + """ + Test the Azimuth Synchronous Signal fitting functions with flags + """ + max_mode = 10 + tod = make_fake_azss_tod(noise_amp=0, n_scans=50, max_mode=max_mode) + + mask = np.zeros((tod.dets.count, tod.samps.count), dtype=bool) + # one detector has partial az coverage, one detector has zero az coverage + mask[0, tod.boresight.az > np.percentile(tod.boresight.az, 95)] = True + mask[-1, :] = True + flags = RangesMatrix.from_mask(mask) + + azss_stats, model_sig_tod = azss.get_azss( + tod, + method='fit', + max_mode=max_mode, + azrange=None, + bins=100, + flags=flags, + ) + self.assertTrue(~np.any(np.isnan(tod.signal))) + self.assertTrue(np.std(tod.signal[~mask]) > np.std(tod.signal[~mask] - model_sig_tod[~mask])) + def test_interpolate(self): """ Test the interpolation method of Azimuth Synchronous Signal subtraction. @@ -144,7 +170,9 @@ def test_interpolate_with_flags(self): tod = make_fake_azss_tod(noise_amp=0, n_scans=50, max_mode=max_mode) mask = np.zeros((tod.dets.count, tod.samps.count), dtype=bool) - mask[-1, :] = True # one detector has low az coverage + # one detector has partial az coverage, one detector has zero az coverage + mask[0, tod.boresight.az > np.percentile(tod.boresight.az, 95)] = True + mask[-1, :] = True flags = RangesMatrix.from_mask(mask) azss_stats, model_sig_tod = azss.get_azss( From f07ba349cfbb81e50d6e5cd5d4b6a328bbd0a5b8 Mon Sep 17 00:00:00 2001 From: ykyohei <38639108+ykyohei@users.noreply.github.com> Date: Mon, 13 Apr 2026 21:29:54 -0400 Subject: [PATCH 2/4] Fix and add test --- sotodlib/tod_ops/azss.py | 3 +-- tests/test_azss.py | 6 ++++++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/sotodlib/tod_ops/azss.py b/sotodlib/tod_ops/azss.py index 05b01bfb7..ef44ca28e 100644 --- a/sotodlib/tod_ops/azss.py +++ b/sotodlib/tod_ops/azss.py @@ -187,7 +187,6 @@ def fit_azss(az, azss_stats, max_mode, fit_range=None, overwrite=False): return L.legval(x_legendre, azss_stats.coeffs.T) coeffs[valid_dets, :] = L.legfit(x_legendre_bin_centers[m], azss_stats.binned_signal[:, m][valid_dets, :].T, max_mode).T - assert ~np.any(np.isnan(coeffs)) model = L.legval(x_legendre, coeffs.T) binned_model = L.legval(x_legendre_bin_centers, coeffs.T) binned_model = np.where(~m, np.nan, binned_model) @@ -203,7 +202,7 @@ def fit_azss(az, azss_stats, max_mode, fit_range=None, overwrite=False): x_legendre_bin_centers = (2 * azss_stats.binned_az - (az_min+az_max)) / (az_max - az_min) x_legendre_bin_centers = np.where(~m, np.nan, x_legendre_bin_centers) if ('coeffs' in azss_stats) and not overwrite: - model[i, :] = L.legval(x_legendre, azss_stats.coeffs[:, i]) + model[i, :] = L.legval(x_legendre, azss_stats.coeffs[i, :]) else: _coeffs = L.legfit(x_legendre_bin_centers[m], azss_stats.binned_signal[:, m][i, :].T, max_mode) _binned_model = L.legval(x_legendre_bin_centers, _coeffs) diff --git a/tests/test_azss.py b/tests/test_azss.py index 33b67fdb1..a6ef32138 100644 --- a/tests/test_azss.py +++ b/tests/test_azss.py @@ -121,6 +121,9 @@ def test_fit(self): self.assertTrue(~np.any(np.isnan(tod.signal))) self.assertTrue(np.std(tod.signal) > np.std(tod.signal - model_sig_tod)) + model = azss.fit_azss(tod.boresight.az, azss_stats, max_mode) + self.assertTrue(np.all(model_sig_tod == model)) + def test_fit_with_flags(self): """ Test the Azimuth Synchronous Signal fitting functions with flags @@ -145,6 +148,9 @@ def test_fit_with_flags(self): self.assertTrue(~np.any(np.isnan(tod.signal))) self.assertTrue(np.std(tod.signal[~mask]) > np.std(tod.signal[~mask] - model_sig_tod[~mask])) + model = azss.fit_azss(tod.boresight.az, azss_stats, max_mode) + self.assertTrue(np.all(model_sig_tod == model)) + def test_interpolate(self): """ Test the interpolation method of Azimuth Synchronous Signal subtraction. From e6aebfaa36e2a6b1e0664f15826b5ac44d55414e Mon Sep 17 00:00:00 2001 From: ykyohei <38639108+ykyohei@users.noreply.github.com> Date: Tue, 28 Apr 2026 14:23:35 -0400 Subject: [PATCH 3/4] refactor to reduce redundancy --- sotodlib/tod_ops/azss.py | 174 +++++++++++++++++++++++++-------------- 1 file changed, 111 insertions(+), 63 deletions(-) diff --git a/sotodlib/tod_ops/azss.py b/sotodlib/tod_ops/azss.py index ef44ca28e..ffa4a6d84 100644 --- a/sotodlib/tod_ops/azss.py +++ b/sotodlib/tod_ops/azss.py @@ -43,6 +43,19 @@ def _check_azcoverage(aman, flags, az=None, coverage_threshold=0.95, return coverages < coverage_threshold, coverages +def _valid_dets_mask(azss_stats, min_valid_bins=1): + """ + Boolean mask of detectors usable for fitting / interpolation. + + A detector is considered valid if it has at least ``min_valid_bins`` + non-NaN bins and is not flagged in ``azss_stats['bad_dets']``, if present. + """ + valid = (~np.isnan(azss_stats.binned_signal)).sum(axis=1) >= min_valid_bins + if 'bad_dets' in azss_stats: + valid &= ~azss_stats['bad_dets'] + return valid + + def bin_by_az(aman, signal=None, az=None, azrange=None, bins=100, flags=None, apodize_edges=True, apodize_edges_samps=1600, apodize_flags=True, apodize_flags_samps=200): @@ -120,6 +133,45 @@ def bin_by_az(aman, signal=None, az=None, azrange=None, bins=100, flags=None, range=azrange, bins=bins, flags=flags, weight_for_signal=weight_for_signal) return binning_dict + +def _legendre_x(az, binned_az, m, fit_range, bin_width): + """ + Build Legendre x-coordinates (rescaled to [-1, 1]) for both the sample + azimuth array and the bin centers, given a 1D mask m of valid bins. + + If fit_range is None the range is inferred from the bin centers + selected by m (extended by half a bin on each side). + """ + if fit_range is None: + az_min = binned_az[m].min() - bin_width / 2 + az_max = binned_az[m].max() + bin_width / 2 + else: + az_min, az_max = fit_range[0], fit_range[1] + span = az_max - az_min + mid = az_min + az_max + x_samp = (2 * az - mid) / span + x_bin = np.where(m, (2 * binned_az - mid) / span, np.nan) + return x_samp, x_bin + + +def _fit_one_det(binned_signal_i, x_bin, m, sigma_i, max_mode): + """ + Fit a Legendre polynomial to one detector's binned signal. + + Returns + ------- + coeffs : ndarray, shape (max_mode + 1,) + binned_model : ndarray, shape (nbins,) (NaN at masked bins) + sum_of_squares : float + redchi2 : float + """ + coeffs = L.legfit(x_bin[m], binned_signal_i[m], max_mode) + binned_model = np.where(m, L.legval(x_bin, coeffs), np.nan) + ssq = np.sum((binned_signal_i[m] - binned_model[m])**2) + redchi2 = ssq / sigma_i**2 / (m.sum() - max_mode - 1) + return coeffs, binned_model, ssq, redchi2 + + def fit_azss(az, azss_stats, max_mode, fit_range=None, overwrite=False): """ Function for fitting Legendre polynomials to signal binned in azimuth. @@ -148,73 +200,71 @@ def fit_azss(az, azss_stats, max_mode, fit_range=None, overwrite=False): Model fit for each detector size ndets x n_samps """ bin_width = azss_stats.binned_az[1] - azss_stats.binned_az[0] + ndets = azss_stats.dets.count + nbins = azss_stats.bin_az_samps.count - model = np.zeros((azss_stats.dets.count, len(az))) + m_2d = ~np.isnan(azss_stats.binned_signal) - m = ~np.isnan(azss_stats.binned_signal) - # Number of valid bins cannot be smaller than mode of Legendre function - valid_dets = np.sum(m, axis=1) >= max_mode + 1 - if 'bad_dets' in azss_stats: - valid_dets = np.logical_and(valid_dets, ~azss_stats['bad_dets']) + # Number of valid bins cannot be smaller than mode of Legendre function. + valid_dets = _valid_dets_mask(azss_stats, min_valid_bins=max_mode + 1) + + model = np.zeros((ndets, len(az))) + coeffs = np.zeros((ndets, max_mode + 1)) + binned_model = np.full((ndets, nbins), np.nan) + sum_of_squares = np.full(ndets, np.nan) + redchi2s = np.full(ndets, np.nan) - if sum(valid_dets) == 0: - logger.info('All the detectors have low az coverage and cannot make model') + if not valid_dets.any(): + logger.warning('All the detectors have low az coverage and cannot make model') return model - mask = ~np.isnan(azss_stats.binned_signal[valid_dets, :]) - is_uniform = np.all(mask == mask[0, :]) + use_cached = ('coeffs' in azss_stats) and not overwrite - coeffs = np.zeros((azss_stats.dets.count, max_mode + 1)) - binned_model = np.full((azss_stats.dets.count, azss_stats.bin_az_samps.count), np.nan) - sum_of_squares = np.full(azss_stats.dets.count, np.nan) - redchi2s = np.full(azss_stats.dets.count, np.nan) - if fit_range is not None: - az_min, az_max = fit_range[0], fit_range[1] - x_legendre = (2 * az - (az_min+az_max)) / (az_max - az_min) - x_legendre_bin_centers = (2 * azss_stats.binned_az - (az_min+az_max)) / (az_max - az_min) - x_legendre_bin_centers = np.where(~m, np.nan, x_legendre_bin_centers) + is_uniform = np.array_equal(m_2d[valid_dets].any(axis=0), + m_2d[valid_dets].all(axis=0)) if is_uniform: - m = mask[0, :] - if fit_range is None: - az_min = np.min(azss_stats.binned_az[m]) - bin_width / 2 - az_max = np.max(azss_stats.binned_az[m]) + bin_width / 2 - x_legendre = (2 * az - (az_min+az_max)) / (az_max - az_min) - x_legendre_bin_centers = (2 * azss_stats.binned_az - (az_min+az_max)) / (az_max - az_min) - x_legendre_bin_centers = np.where(~m, np.nan, x_legendre_bin_centers) - - if ('coeffs' in azss_stats) and not overwrite: - return L.legval(x_legendre, azss_stats.coeffs.T) - - coeffs[valid_dets, :] = L.legfit(x_legendre_bin_centers[m], azss_stats.binned_signal[:, m][valid_dets, :].T, max_mode).T - model = L.legval(x_legendre, coeffs.T) - binned_model = L.legval(x_legendre_bin_centers, coeffs.T) - binned_model = np.where(~m, np.nan, binned_model) - sum_of_squares = np.sum(((azss_stats.binned_signal[:, m] - binned_model[:, m])**2), axis=-1) - redchi2s = sum_of_squares/azss_stats.uniform_binned_signal_sigma**2 / (len(x_legendre_bin_centers[m]) - max_mode - 1) - + m = m_2d[valid_dets][0] + x_samp, x_legendre_bin_centers = _legendre_x( + az, azss_stats.binned_az, m, fit_range, bin_width) + + if use_cached: + return L.legval(x_samp, azss_stats.coeffs.T) + + # Vectorized Legendre fit over all valid detectors at once. + c = L.legfit(x_legendre_bin_centers[m], + azss_stats.binned_signal[valid_dets][:, m].T, + max_mode).T # shape: (n_valid, max_mode + 1) + coeffs[valid_dets] = c + bm = L.legval(x_legendre_bin_centers, c.T) + binned_model[valid_dets] = np.where(m, bm, np.nan) + model[valid_dets] = L.legval(x_samp, c.T) + + diff = azss_stats.binned_signal[valid_dets][:, m] - bm[:, m] + sum_of_squares[valid_dets] = np.sum(diff**2, axis=-1) + redchi2s[valid_dets] = ( + sum_of_squares[valid_dets] + / azss_stats.uniform_binned_signal_sigma[valid_dets]**2 + / (m.sum() - max_mode - 1) + ) else: - for i, m, binned_signal in zip(np.where(valid_dets)[0], mask, azss_stats.binned_signal[valid_dets, :]): - if fit_range is None: - az_min = np.min(azss_stats.binned_az[m]) - bin_width / 2 - az_max = np.max(azss_stats.binned_az[m]) + bin_width / 2 - x_legendre = (2 * az - (az_min+az_max)) / (az_max - az_min) - x_legendre_bin_centers = (2 * azss_stats.binned_az - (az_min+az_max)) / (az_max - az_min) - x_legendre_bin_centers = np.where(~m, np.nan, x_legendre_bin_centers) - if ('coeffs' in azss_stats) and not overwrite: - model[i, :] = L.legval(x_legendre, azss_stats.coeffs[i, :]) - else: - _coeffs = L.legfit(x_legendre_bin_centers[m], azss_stats.binned_signal[:, m][i, :].T, max_mode) - _binned_model = L.legval(x_legendre_bin_centers, _coeffs) - _binned_model = np.where(~m, np.nan, _binned_model) - _sum_of_squares = np.sum(((azss_stats.binned_signal[:, m][i, :] - _binned_model[m])**2), axis=-1) - _redchi2s = _sum_of_squares/azss_stats.uniform_binned_signal_sigma[i]**2 / (len(x_legendre_bin_centers[m]) - max_mode - 1) - coeffs[i, :] = _coeffs - binned_model[i, :] = _binned_model - sum_of_squares[i] = _sum_of_squares - redchi2s[i] = _redchi2s - model[i, :] = L.legval(x_legendre, _coeffs) - if ('coeffs' in azss_stats) and not overwrite: + x_legendre_bin_centers = None # set inside the loop + for i in np.where(valid_dets)[0]: + m_i = m_2d[i] + x_samp, x_legendre_bin_centers = _legendre_x( + az, azss_stats.binned_az, m_i, fit_range, bin_width) + if use_cached: + model[i] = L.legval(x_samp, azss_stats.coeffs[i]) + continue + (coeffs[i], + binned_model[i], + sum_of_squares[i], + redchi2s[i]) = _fit_one_det( + azss_stats.binned_signal[i], x_legendre_bin_centers, m_i, + azss_stats.uniform_binned_signal_sigma[i], max_mode) + model[i] = L.legval(x_samp, coeffs[i]) + + if use_cached: return model mode_names = np.array([f'legendre{mode}' for mode in range(max_mode + 1)], dtype=' 0 - if 'bad_dets' in azss_stats: - valid_dets = np.logical_and(valid_dets, ~azss_stats['bad_dets']) - if sum(valid_dets) == 0: - logger.info('All the detectors have low az coverage and cannot make model') + valid_dets = _valid_dets_mask(azss_stats, min_valid_bins=1) + if not valid_dets.any(): + logger.warning('All the detectors have low az coverage and cannot make model') return model mask = ~np.isnan(azss_stats.binned_signal[valid_dets, :]) From 3c9b29ff1f97d7195c0fa2a37ee8ebe0b138702b Mon Sep 17 00:00:00 2001 From: ykyohei <38639108+ykyohei@users.noreply.github.com> Date: Tue, 28 Apr 2026 14:55:33 -0400 Subject: [PATCH 4/4] fix --- sotodlib/tod_ops/azss.py | 13 +++++++------ tests/test_azss.py | 1 + 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/sotodlib/tod_ops/azss.py b/sotodlib/tod_ops/azss.py index ffa4a6d84..18481d4b4 100644 --- a/sotodlib/tod_ops/azss.py +++ b/sotodlib/tod_ops/azss.py @@ -248,7 +248,6 @@ def fit_azss(az, azss_stats, max_mode, fit_range=None, overwrite=False): / (m.sum() - max_mode - 1) ) else: - x_legendre_bin_centers = None # set inside the loop for i in np.where(valid_dets)[0]: m_i = m_2d[i] x_samp, x_legendre_bin_centers = _legendre_x( @@ -467,15 +466,17 @@ def get_azss_model(aman, azss_stats, az=None, method='interpolate', logger.warning('All the detectors have low az coverage and cannot make model') return model - mask = ~np.isnan(azss_stats.binned_signal[valid_dets, :]) - is_uniform = np.all(mask == mask[0, :]) + m_2d = ~np.isnan(azss_stats.binned_signal) + is_uniform = np.array_equal(m_2d[valid_dets].any(axis=0), + m_2d[valid_dets].all(axis=0)) if is_uniform: - m = mask[0, :] + m = m_2d[valid_dets][0] f_template = interp1d(azss_stats.binned_az[m], azss_stats.binned_signal[:, m][valid_dets, :], fill_value='extrapolate') model[valid_dets, :] = f_template(az) else: - for i, m, binned_signal in zip(np.where(valid_dets)[0], mask, azss_stats.binned_signal[valid_dets, :]): - f_template = interp1d(azss_stats.binned_az[m], binned_signal[m], fill_value='extrapolate') + for i in np.where(valid_dets)[0]: + m = m_2d[i] + f_template = interp1d(azss_stats.binned_az[m], azss_stats.binned_signal[i][m], fill_value='extrapolate') model[i, :] = f_template(az) if np.any(~np.isfinite(model)): diff --git a/tests/test_azss.py b/tests/test_azss.py index a6ef32138..256956771 100644 --- a/tests/test_azss.py +++ b/tests/test_azss.py @@ -148,6 +148,7 @@ def test_fit_with_flags(self): self.assertTrue(~np.any(np.isnan(tod.signal))) self.assertTrue(np.std(tod.signal[~mask]) > np.std(tod.signal[~mask] - model_sig_tod[~mask])) + # check consistency of model made with cached legendre coeffs model = azss.fit_azss(tod.boresight.az, azss_stats, max_mode) self.assertTrue(np.all(model_sig_tod == model))