From 8d2c132922940a23cd8e3847d01325adb34bf27f Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Sat, 30 May 2026 21:47:59 +1200 Subject: [PATCH 01/22] Close TPF file handle before cache removal to fix Windows deletion --- tessreduce/tessreduce.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 2246b25..0ca0abf 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -478,6 +478,7 @@ def get_TESS(self,ra=None,dec=None,name=None,size=None,sector=None, tpf = tess.download(quality_bitmask=quality_bitmask,cutout_size=size,download_dir=cache_dir) if not cache: try: + tpf.hdu.close() os.remove(tpf.path) if self.verbose > 0: print('Cache removed') From 8260f429cf4a68a214e30a76d6e9bb0125c8e681 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Sun, 31 May 2026 17:13:12 +1200 Subject: [PATCH 02/22] Delete TPF cache file after reduction completes rather than immediately after download --- tessreduce/tessreduce.py | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 0ca0abf..d28654c 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -223,6 +223,7 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect self._quality_bitmask = quality_bitmask self._smooth_motion = smooth_motion self._timing = timing + self._cache_path = None # Offline Paths if catalogue_path is None: @@ -477,13 +478,9 @@ def get_TESS(self,ra=None,dec=None,name=None,size=None,sector=None, # Download tpf = tess.download(quality_bitmask=quality_bitmask,cutout_size=size,download_dir=cache_dir) if not cache: - try: - tpf.hdu.close() - os.remove(tpf.path) - if self.verbose > 0: - print('Cache removed') - except OSError: - print(f'Failed to remove: {tpf.path}') + self._cache_path = tpf.path + else: + self._cache_path = None # Check to ensure it succeeded if tpf is None: @@ -2818,6 +2815,15 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, except Exception: print(traceback.format_exc()) + if self._cache_path is not None: + try: + os.remove(self._cache_path) + if self.verbose > 0: + print('Cache removed') + except OSError: + print(f'Failed to remove cache: {self._cache_path}') + self._cache_path = None + def external_photometry(self,size=50,phot=None): """ From 1c331e82be8bc9768d2fa2ba7d692d43cd386396 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 4 Jun 2026 09:39:53 +1200 Subject: [PATCH 03/22] Use SLURM-aware CPU detection to fix parallelism on HPC multiprocessing.cpu_count() returns the total node CPU count, not the cores allocated by SLURM. With n_jobs=-1 joblib would spawn workers for every node CPU, causing severe oversubscription on jobs with a smaller allocation and making parallel steps appear serial. _available_cores() now checks SLURM_CPUS_PER_TASK first, then os.sched_getaffinity(0) (cgroup-aware on Linux), before falling back to cpu_count(). Also removes a duplicate self.num_cores assignment that was overwriting the resolved value. --- tessreduce/tessreduce.py | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index d28654c..0c5bc96 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -22,6 +22,26 @@ import multiprocessing from joblib import Parallel, delayed +def _available_cores(): + """Return the number of CPU cores available to this process. + + Checks in priority order: + 1. SLURM_CPUS_PER_TASK - cores allocated by SLURM scheduler + 2. os.sched_getaffinity - respects cgroup/affinity limits (Linux) + 3. multiprocessing.cpu_count - total node CPUs (fallback) + """ + slurm = os.environ.get('SLURM_CPUS_PER_TASK') + if slurm is not None: + try: + return int(slurm) + except ValueError: + pass + try: + return len(os.sched_getaffinity(0)) + except AttributeError: + pass + return multiprocessing.cpu_count() + from .catalog_tools import * from .calibration_tools import * from .ground_tools import ground @@ -207,8 +227,8 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect self.imaging = imaging self.parallel = parallel self._col_offset = col_offset - if isinstance(num_cores, str): - self.num_cores = multiprocessing.cpu_count() + if num_cores == -1 or isinstance(num_cores, str): + self.num_cores = _available_cores() else: self.num_cores = num_cores self._assign_phot_method(phot_method) @@ -231,7 +251,6 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect elif catalogue_path is False: catalogue_path = None self._catalogue_path = catalogue_path - self.num_cores = num_cores self.imaging = imaging self._prf_path = prf_path self._vector_path = vector_path From 0587f2969f709c56c649139f94a7781d8856b606 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 4 Jun 2026 10:22:17 +1200 Subject: [PATCH 04/22] Vectorise component loops and parallelise serial frame operations helpers.py: - fix_background_anomalies: replace labeled==i component loops with np.bincount/np.isin vectorisation in sep_validate, bad_bkg_mask, and _fit_residual; precompute distance array per SEP object to avoid redundant sqrt calls in annulus search - grad_clip_fill_bkg: replace five Python component loops (size accumulation, size filtering, overlap ratio calculation, ratio filtering, small-region removal) with np.bincount + np.isin - blend_dynamic_background: extract _blend_frame helper and dispatch all T frames via Parallel(n_jobs=n_jobs); add n_jobs parameter - parallel_strap_fit: extract _strap_fit_col helper and dispatch strap columns via Parallel(n_jobs=n_jobs); add n_jobs parameter - regional_stats_mask: extract _clip_region helper and dispatch regions via Parallel(n_jobs=n_jobs); add n_jobs parameter - Add _shift_one and _shift_ref_one module-level helpers for picklable per-frame shift operations tessreduce.py: - shift_images: parallelise both median and normal branches via Parallel(n_jobs=self.num_cores) using the new shift helpers - _bkg_median: replace list comprehension with np.nanmedian(axis=(1,2)) - blend_dynamic_background call: pass n_jobs=self.num_cores background_separator.py: - Replace nested list comprehension for yx_all coordinate array with np.mgrid + np.column_stack (two locations) --- tessreduce/background_separator.py | 6 +- tessreduce/helpers.py | 269 ++++++++++++++++------------- tessreduce/tessreduce.py | 29 +++- 3 files changed, 171 insertions(+), 133 deletions(-) diff --git a/tessreduce/background_separator.py b/tessreduce/background_separator.py index 60c35ed..5de049f 100644 --- a/tessreduce/background_separator.py +++ b/tessreduce/background_separator.py @@ -336,7 +336,8 @@ def fit(self, method: str = 'local_pca', **kwargs) -> np.ndarray: ) # (T,) ceiling = self.flux + (k * frame_noise)[:, np.newaxis, np.newaxis] - yx_all = np.array([[y, x] for y in range(self._X) for x in range(self._Y)]) + _yy, _xx = np.mgrid[:self._X, :self._Y] + yx_all = np.column_stack([_yy.ravel(), _xx.ravel()]) for t in range(self._T): bad = bkg[t] > ceiling[t] if not bad.any(): @@ -705,7 +706,8 @@ def floor_negative_flux(self, noise_floor=0.5, max_iter=3): noise = 1.4826 * mad # (X, Y) threshold = -noise_floor * noise # (X, Y) — always <= 0 - yx_all = np.array([[y, x] for y in range(self._X) for x in range(self._Y)]) + _yy, _xx = np.mgrid[:self._X, :self._Y] + yx_all = np.column_stack([_yy.ravel(), _xx.ravel()]) for t in range(self._T): frame = self.flux[t].copy() diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index 4e93ba5..864848c 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -166,6 +166,16 @@ def unknown_mask(image): return mask +def _shift_one(frame, s): + if np.nansum(abs(frame)) > 0: + return shift(frame, [s[0], s[1]], mode='nearest', order=5) + return frame + +def _shift_ref_one(ref, frame_for_check, s): + if np.nansum(abs(frame_for_check)) > 0: + return shift(ref, [-s[1], -s[0]], mode='nearest', order=5) + return frame_for_check + def parallel_bkg3(data,mask): data = deepcopy(data) data[mask] = np.nan @@ -1104,20 +1114,27 @@ def Extract_fits(pixelfile): print('OSError ',pixelfile) return -def regional_stats_mask(image,size=90,sigma=3,iters=10): +def _clip_region(image, rx, ry, sigma, iters): + m, me, s = sigma_clipped_stats(image[ry, rx], maxiters=iters) + cut = (image[rx, ry] >= me + sigma * s) | (image[rx, ry] <= me - sigma * s) + return rx[cut], ry[cut] + +def regional_stats_mask(image, size=90, sigma=3, iters=10, n_jobs=1): if size < 30: print('!!! Region size is small !!!') sx, sy = image.shape X, Y = np.ogrid[0:sx, 0:sy] - regions = sy//size * (X//size) + Y//size - max_reg = np.max(regions) + regions = sy // size * (X // size) + Y // size + max_reg = int(np.max(regions)) + + region_pixels = [np.where(regions == i) for i in range(max_reg + 1)] clip = np.zeros_like(image) - for i in range(max_reg+1): - rx,ry = np.where(regions == i) - m,me, s = sigma_clipped_stats(image[ry,rx],maxiters=iters) - cut_ind = np.where((image[rx,ry] >= me+sigma*s) | (image[rx,ry] <= me-sigma*s)) - clip[rx[cut_ind],ry[cut_ind]] = 1 + results = Parallel(n_jobs=n_jobs)( + delayed(_clip_region)(image, rx, ry, sigma, iters) + for rx, ry in region_pixels) + for rx_cut, ry_cut in results: + clip[rx_cut, ry_cut] = 1 return clip @@ -1201,63 +1218,59 @@ def grad_clip_fill_bkg(bkg,sigma=3,max_size=1000): ap = (abs(a) - a_med) > 3*a_std ap = fftconvolve(ap,np.ones((3,3)),mode='same') > 0.8 - b_labeled, b_objects = label(bp) - a_labeled, a_objects = label(ap) + b_labeled, b_objects = label(bp) + a_labeled, a_objects = label(ap) - b_obj_size = [] - for i in range(b_objects): - b_obj_size += [np.sum(b_labeled==i)] - b_obj_size = np.array(b_obj_size) + # Component sizes via bincount (preserves original loop range 0..n_objects-1) + b_obj_size = np.bincount(b_labeled.ravel(), minlength=b_objects + 1)[:b_objects] + a_obj_size = np.bincount(a_labeled.ravel(), minlength=a_objects + 1)[:a_objects] - a_obj_size = [] - for i in range(a_objects): - a_obj_size += [np.sum(a_labeled==i)] - a_obj_size = np.array(a_obj_size) + # Zero out components outside the valid size range + bad_a = np.where((a_obj_size >= max_size) | (a_obj_size <= 9))[0] + if len(bad_a) > 0: + a_labeled[np.isin(a_labeled, bad_a)] = 0 - for i in range(a_objects): - if (a_obj_size[i] >= max_size) | (a_obj_size[i] <= 9): - a_labeled[a_labeled==i] = 0 + bad_b = np.where((b_obj_size >= max_size) | (b_obj_size <= 9))[0] + if len(bad_b) > 0: + b_labeled[np.isin(b_labeled, bad_b)] = 0 - for i in range(b_objects): - if (b_obj_size[i] >= max_size) | (b_obj_size[i] <= 9): - b_labeled[b_labeled==i] = 0 - - - overlap = (a_labeled>0) & (b_labeled>0) - y,x = np.where(overlap) + overlap = (a_labeled > 0) & (b_labeled > 0) + y, x = np.where(overlap) + good_a = np.unique(a_labeled[y, x]) + good_b = np.unique(b_labeled[y, x]) - good_a = np.unique(a_labeled[y,x]) - good_b = np.unique(b_labeled[y,x]) + # Overlap ratios via bincount — counts per label, then fraction overlapping + _a_max = int(a_labeled.max()) + 1 if a_labeled.max() > 0 else 1 + _b_max = int(b_labeled.max()) + 1 if b_labeled.max() > 0 else 1 + a_lab_flat = a_labeled.ravel() + b_lab_flat = b_labeled.ravel() + ov_flat = overlap.ravel().astype(float) + + a_sizes = np.bincount(a_lab_flat, minlength=_a_max) + a_overlap_counts = np.bincount(a_lab_flat, weights=ov_flat, minlength=_a_max) + a_ratio = a_overlap_counts[good_a] / np.maximum(a_sizes[good_a], 1) + + b_sizes = np.bincount(b_lab_flat, minlength=_b_max) + b_overlap_counts = np.bincount(b_lab_flat, weights=ov_flat, minlength=_b_max) + b_ratio = b_overlap_counts[good_b] / np.maximum(b_sizes[good_b], 1) + + suppress_a = good_a[a_ratio < 0.2] + if len(suppress_a) > 0: + a_labeled[np.isin(a_labeled, suppress_a)] = 0 + + suppress_b = good_b[b_ratio < 0.2] + if len(suppress_b) > 0: + b_labeled[np.isin(b_labeled, suppress_b)] = 0 - a_ratio = [] - for ind in good_a: - eh = a_labeled == ind - eh2 = eh * overlap - ratio = np.sum(eh2) / np.sum(eh) - a_ratio += [ratio] - a_ratio = np.array(a_ratio) - - - b_ratio = [] - for ind in good_b: - eh = b_labeled == ind - eh2 = eh * overlap - ratio = np.sum(eh2) / np.sum(eh) - b_ratio += [ratio] - b_ratio = np.array(b_ratio) - - - for i in good_a[a_ratio<0.2]: - a_labeled[a_labeled==i] = 0 - for i in good_b[b_ratio<0.2]: - b_labeled[b_labeled==i] = 0 c = (a_labeled + b_labeled) > 0 - c_labeled, c_objects = label(c==0) - for i in range(c_objects): - if np.sum(c_labeled==i) < 10: - c[c_labeled==i] = 1 + c_labeled, c_objects = label(c == 0) + if c_objects > 0: + c_sizes = np.bincount(c_labeled.ravel(), minlength=c_objects + 1)[:c_objects] + small = np.where(c_sizes < 10)[0] + if len(small) > 0: + c[np.isin(c_labeled, small)] = 1 #points = fftconvolve(c,np.ones((5,5)),mode='same') points = c>0#oints > 0.8 @@ -1323,29 +1336,30 @@ def fit_strap(data,mask): p[np.isnan(p)] = p2[np.isnan(p)] return p -def parallel_strap_fit(frame,frame_bkg,frame_err,mask,repeats=3,tol=3): +def _strap_fit_col(col_data, norm_col): + d = abs(np.gradient(col_data)) + m, med, std = sigma_clipped_stats(d, maxiters=10) + nm = (d > med + std) * 1 + nm = np.convolve(nm, np.ones(3), mode='same') + nm = (nm == 0) + q = fit_strap(norm_col, nm) + if len(col_data) > 110: + q = savgol_filter(q, 101, 1) + else: + q[:] = np.nanmedian(q) + return q + +def parallel_strap_fit(frame, frame_bkg, frame_err, mask, repeats=3, tol=3, n_jobs=1): norm = frame / frame_bkg - sind = np.where(np.nansum(mask,axis=0)>0)[0] + sind = np.where(np.nansum(mask, axis=0) > 0)[0] qe = np.ones_like(frame) - for i in sind: - y = frame[:,i] - d = abs(np.gradient(y)) - m, med, std = sigma_clipped_stats(d,maxiters=10) - nm = (d > med + std) * 1 - nm = np.convolve(nm,np.ones(3),mode='same') - nm = (nm == 0) - #for r in range(repeats): - q = fit_strap(norm[:,i],nm) - #q /= frame_bkg[:,i] - if len(y) > 110: - q = savgol_filter(q,101,1) - else: - mq = np.nanmedian(q) - q[:] = mq - - qe[:,i] = q - - #qe[:,np.nanmean((qe-1),axis=) < 5e-3] = 1 + if len(sind) == 0: + return qe + results = Parallel(n_jobs=n_jobs)( + delayed(_strap_fit_col)(frame[:, i], norm[:, i]) + for i in sind) + for col, q in zip(sind, results): + qe[:, col] = q return qe @@ -1520,10 +1534,10 @@ def _block_sigma(r, box): if ap_mask.sum() == 0 or snr_map[ap_mask].mean() <= sep_snr_thresh: continue cx, cy = obj['x'], obj['y'] + dist = np.sqrt((xx - cx)**2 + (yy - cy)**2) true_r = None for r in range(2, 20): - ann = (np.sqrt((xx - cx)**2 + (yy - cy)**2) >= r - 0.5) & \ - (np.sqrt((xx - cx)**2 + (yy - cy)**2) < r + 0.5) + ann = (dist >= r - 0.5) & (dist < r + 0.5) if ann.sum() == 0: break if lap_sub[ann].mean() < noise: @@ -1533,8 +1547,7 @@ def _block_sigma(r, box): true_r = 19 if true_r is None or true_r < 2 or true_r > 5: continue - circ = np.sqrt((xx - cx)**2 + (yy - cy)**2) <= true_r - sep_mask |= circ + sep_mask |= dist <= true_r lap_med = np.nanmedian(lap_abs) lap_mad = np.nanmedian(np.abs(lap_abs - lap_med)) @@ -1547,15 +1560,18 @@ def _block_sigma(r, box): edge_border[0, :] = True; edge_border[-1, :] = True edge_border[:, 0] = True; edge_border[:, -1] = True labeled, n_comp = label(flagged_coarse) - smooth_mask = np.zeros((NY, NX), dtype=bool) - for comp_id in range(1, n_comp + 1): - comp = labeled == comp_id - if (comp & edge_border).any() and not (comp & sep_mask).any(): - smooth_mask |= comp - elif (comp & is_sharp).any() and (comp & sep_mask).any(): - sharp_mask |= comp - else: - smooth_mask |= comp + if n_comp > 0: + lab_flat = labeled.ravel() + touches_edge = np.zeros(n_comp + 1, dtype=bool) + touches_sep_lbl = np.zeros(n_comp + 1, dtype=bool) + touches_sharp_lbl = np.zeros(n_comp + 1, dtype=bool) + np.bitwise_or.at(touches_edge, lab_flat, edge_border.ravel()) + np.bitwise_or.at(touches_sep_lbl, lab_flat, sep_mask.ravel()) + np.bitwise_or.at(touches_sharp_lbl, lab_flat, is_sharp.ravel()) + # sharp: first condition False AND touches_sharp AND touches_sep + is_sharp_lbl = ~(touches_edge & ~touches_sep_lbl) & touches_sharp_lbl & touches_sep_lbl + is_sharp_lbl[0] = False + sharp_mask |= is_sharp_lbl[labeled] else: sharp_mask |= flagged_coarse & is_sharp @@ -1618,11 +1634,12 @@ def _block_sigma(r, box): lap_mad = np.nanmedian(np.abs(lap_abs - lap_med)) high_lap = lap_abs > lap_med + bad_bkg_sigma * 1.4826 * lap_mad labeled_lap, n_lap = label(high_lap) - bad_bkg_mask = np.zeros((NY, NX), dtype=bool) - for cid in range(1, n_lap + 1): - comp = labeled_lap == cid - if comp.sum() >= bad_bkg_min_area: - bad_bkg_mask |= comp + if n_lap > 0: + areas = np.bincount(labeled_lap.ravel(), minlength=n_lap + 1) + large = np.flatnonzero(areas[1:] >= bad_bkg_min_area) + 1 + bad_bkg_mask = np.isin(labeled_lap, large) if len(large) > 0 else np.zeros((NY, NX), dtype=bool) + else: + bad_bkg_mask = np.zeros((NY, NX), dtype=bool) return fixed, excess, sharp_mask, bad_bkg_mask @@ -1663,11 +1680,14 @@ def _fit_residual(residual, exclude_mask): lap_mad = np.nanmedian(np.abs(lap_abs - lap_med)) is_sharp = lap_abs > lap_med + 3 * 1.4826 * lap_mad labeled_c, n_c = label(flagged) - for cid in range(1, n_c + 1): - comp = labeled_c == cid - n_sp = (comp & is_sharp).sum() - if n_sp / comp.sum() >= 0.3: - corr[comp] = 0.0 + if n_c > 0: + lab_flat = labeled_c.ravel() + comp_sizes = np.bincount(lab_flat, minlength=n_c + 1) + sharp_counts = np.bincount(lab_flat, weights=is_sharp.ravel(), minlength=n_c + 1) + sharp_frac = sharp_counts / np.maximum(comp_sizes, 1) + suppress = np.flatnonzero(sharp_frac[1:] >= 0.3) + 1 + if len(suppress) > 0: + corr[np.isin(labeled_c, suppress)] = 0.0 return corr residuals = flux - bkg_fixed @@ -1776,7 +1796,22 @@ def orbit_ref_subtract(flux, times_mjd, sector=None, camera=None, return result, segments, orbit_refs -def blend_dynamic_background(bkg_new, bkg_prev, flux, sigma=2.0, sharp_masks=None): +def _blend_frame(bkg_new_i, bkg_prev_i, delta_i, resid_prev_i, sigma, sharp_mask_i, gauss_kernel): + _, _, scale = sigma_clipped_stats(resid_prev_i) + w = np.clip(delta_i / (sigma * scale + 1e-10), 0, 1) + w_zero = (w == 0) + if w_zero.any(): + w_med = median_filter(w_zero.astype(float), size=7) + diff_mask = w_zero & ~(w_med > 0.5) + w[w_med > 0.5] = 0 + if diff_mask.any(): + w[diff_mask] = np.nan + w = interpolate_replace_nans(w, gauss_kernel) + if sharp_mask_i is not None: + w[sharp_mask_i] = 0.0 + return (1 - w) * bkg_new_i + w * bkg_prev_i + +def blend_dynamic_background(bkg_new, bkg_prev, flux, sigma=2.0, sharp_masks=None, n_jobs=1): """Per-pixel blend bkg_new toward bkg_prev based on residual quality. For each pixel, compares |flux - bkg_new| vs |flux - bkg_prev|. Where @@ -1803,25 +1838,15 @@ def blend_dynamic_background(bkg_new, bkg_prev, flux, sigma=2.0, sharp_masks=Non flux = np.array(flux) resid_new = np.abs(flux - bkg_new) resid_prev = np.abs(flux - bkg_prev) - delta = resid_new - resid_prev # positive = bkg_new is worse - _gauss_kernel = Gaussian2DKernel(1.0) - - out = bkg_new.copy() - for i in range(bkg_new.shape[0]): - _, _, scale = sigma_clipped_stats(resid_prev[i]) - w = np.clip(delta[i] / (sigma * scale + 1e-10), 0, 1) - - w_zero = (w == 0) - if w_zero.any(): - w_med = median_filter(w_zero.astype(float), size=7) - diff_mask = w_zero & ~(w_med > 0.5) - w[w_med > 0.5] = 0 - if diff_mask.any(): - w[diff_mask] = np.nan - w = interpolate_replace_nans(w, _gauss_kernel) - - if sharp_masks is not None: - w[sharp_masks[i]] = 0.0 - out[i] = (1 - w) * bkg_new[i] + w * bkg_prev[i] - return out + delta = resid_new - resid_prev + gauss_kernel = Gaussian2DKernel(1.0) + T = bkg_new.shape[0] + + sharp_list = [sharp_masks[i] if sharp_masks is not None else None for i in range(T)] + + results = Parallel(n_jobs=n_jobs)( + delayed(_blend_frame)(bkg_new[i], bkg_prev[i], delta[i], resid_prev[i], + sigma, sharp_list[i], gauss_kernel) + for i in range(T)) + return np.array(results) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 0c5bc96..e137043 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -854,7 +854,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa bkg_s1 = np.array(bkg_smth) bkg_smth = Parallel(n_jobs=self.num_cores)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*new_mask) if blend_dynamic: - bkg_smth = blend_dynamic_background(bkg_smth, bkg_s1, flux) + bkg_smth = blend_dynamic_background(bkg_smth, bkg_s1, flux, n_jobs=self.num_cores) _times['residual surface rerun'] = time.perf_counter() - _t else: @@ -878,7 +878,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa bkg_pre_fix = np.array(self.bkg) from .adaptive_background import get_tessvectors, _interpolate_angles _df = get_tessvectors(self.sector, self.tpf.camera, data_path=self._vector_path) - _bkg_median = np.array([np.nanmedian(self.bkg[i]) for i in range(len(self.bkg))]) + _bkg_median = np.nanmedian(self.bkg, axis=(1, 2)) if _df is not None: _earth_angle, _moon_angle = _interpolate_angles(self.mjd, _df) _high_bkg_frames = (_earth_angle < 30.0) | (_moon_angle < 30.0) | (_bkg_median > 300.0) @@ -1541,20 +1541,31 @@ def shift_images(self,median=False): """ + from .helpers import _shift_one, _shift_ref_one shifted = self.flux.copy() nans = ~np.isfinite(shifted) shifted[nans] = 0. if median: - for i in range(len(shifted)): - if np.nansum(abs(shifted[i])) > 0: - shifted[i] = shift(self.ref,[-self.shift[i,1],-self.shift[i,0]], mode='nearest',order=5) + if self.parallel: + result = Parallel(n_jobs=self.num_cores)( + delayed(_shift_ref_one)(self.ref, shifted[i], self.shift[i]) + for i in range(len(shifted))) + shifted = np.array(result) + else: + for i in range(len(shifted)): + shifted[i] = _shift_ref_one(self.ref, shifted[i], self.shift[i]) self.flux -= shifted else: - for i in range(len(shifted)): - if np.nansum(abs(shifted[i])) > 0: - shifted[i] = shift(shifted[i],[self.shift[i,0],self.shift[i,1]],mode='nearest',order=5)#mode='constant',cval=np.nan) - self.flux = shifted + if self.parallel: + result = Parallel(n_jobs=self.num_cores)( + delayed(_shift_one)(shifted[i], self.shift[i]) + for i in range(len(shifted))) + self.flux = np.array(result) + else: + for i in range(len(shifted)): + shifted[i] = _shift_one(shifted[i], self.shift[i]) + self.flux = shifted def bin_data(self,lc=None,time_bin=6/24,frames = None): """ From 4fd2f1bfcc2e6679f9c20821ff68256ffef6dbc1 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 4 Jun 2026 10:34:05 +1200 Subject: [PATCH 05/22] Fix index swap, component off-by-one, and QE infinity in pipeline MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit helpers.py: - _clip_region: fix image[rx,ry] → image[ry,rx] in the threshold check so it is consistent with the stats computation; previously stats were computed on image[ry,rx] (y,x) but the mask was applied to image[rx,ry] (x,y), clipping the wrong pixels - grad_clip_fill_bkg: fix component size filter to cover all labeled regions 1..n_objects; previous [:n_objects] slice excluded the last real component and included the background (label 0), allowing the final artifact to bypass size filtering tessreduce.py: - _calc_qe: replace zero-only guard with np.isfinite check so infinities produced by division by zero background are also converted to NaN before sigma-clipping; suppress the runtime warning with np.errstate --- tessreduce/helpers.py | 14 +++++++------- tessreduce/tessreduce.py | 5 +++-- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index 864848c..e905649 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -1116,7 +1116,7 @@ def Extract_fits(pixelfile): def _clip_region(image, rx, ry, sigma, iters): m, me, s = sigma_clipped_stats(image[ry, rx], maxiters=iters) - cut = (image[rx, ry] >= me + sigma * s) | (image[rx, ry] <= me - sigma * s) + cut = (image[ry, rx] >= me + sigma * s) | (image[ry, rx] <= me - sigma * s) return rx[cut], ry[cut] def regional_stats_mask(image, size=90, sigma=3, iters=10, n_jobs=1): @@ -1221,16 +1221,16 @@ def grad_clip_fill_bkg(bkg,sigma=3,max_size=1000): b_labeled, b_objects = label(bp) a_labeled, a_objects = label(ap) - # Component sizes via bincount (preserves original loop range 0..n_objects-1) - b_obj_size = np.bincount(b_labeled.ravel(), minlength=b_objects + 1)[:b_objects] - a_obj_size = np.bincount(a_labeled.ravel(), minlength=a_objects + 1)[:a_objects] + # Component sizes via bincount — labels 1..n_objects, index 0 is background + b_obj_size = np.bincount(b_labeled.ravel(), minlength=b_objects + 1) + a_obj_size = np.bincount(a_labeled.ravel(), minlength=a_objects + 1) - # Zero out components outside the valid size range - bad_a = np.where((a_obj_size >= max_size) | (a_obj_size <= 9))[0] + # Zero out components outside the valid size range (skip label 0 = background) + bad_a = np.where((a_obj_size[1:] >= max_size) | (a_obj_size[1:] <= 9))[0] + 1 if len(bad_a) > 0: a_labeled[np.isin(a_labeled, bad_a)] = 0 - bad_b = np.where((b_obj_size >= max_size) | (b_obj_size <= 9))[0] + bad_b = np.where((b_obj_size[1:] >= max_size) | (b_obj_size[1:] <= 9))[0] + 1 if len(bad_b) > 0: b_labeled[np.isin(b_labeled, bad_b)] = 0 diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index e137043..99b71bf 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -671,8 +671,9 @@ def _calc_qe(self):#,flux_e): ''' time = deepcopy(self.mjd) strap_data = (self.flux) * ((self.mask&4) > 0)*(~self.mask&1) - qe = strap_data/self.bkg - qe[qe == 0] = np.nan + with np.errstate(divide='ignore', invalid='ignore'): + qe = strap_data / self.bkg + qe[~np.isfinite(qe)] = np.nan m,med,std = sigma_clipped_stats(qe,axis=1,sigma_upper=2) qes = np.ones_like(qe) qes[:,:,:] = med[:,np.newaxis,:] From 31e1e6afcc6768d76eec9db0874fa79e8da9425e Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 4 Jun 2026 10:40:55 +1200 Subject: [PATCH 06/22] Fix PSF position shift inversion, kernel accumulation, and ref baseline MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit psf_photom.py: - psf_position: invert finite check so invalid (NaN/inf) ext_shift is zeroed rather than valid shifts being discarded - minimize_psf_flux: convolve a local copy of self.psf with the kernel instead of mutating self.psf in-place; previously each optimizer iteration re-convolved the already-convolved PSF, progressively broadening it through the fit tessreduce.py: - reduce: replace nanmin with nanpercentile(1) for reference baseline removal in difference imaging mode; nanmin is dominated by a single outlier negative pixel (cosmic ray, bad pixel), biasing every difference frame by that amount helpers.py: - _clip_region: fix image[rx,ry] → image[ry,rx] in threshold check so mask is applied to the same pixels as the sigma-clipped stats - grad_clip_fill_bkg: fix component size arrays to cover labels 1..n_objects (exclude background label 0, include last component) --- tessreduce/psf_photom.py | 7 +++---- tessreduce/tessreduce.py | 2 +- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/tessreduce/psf_photom.py b/tessreduce/psf_photom.py index d2e8073..f27baf7 100644 --- a/tessreduce/psf_photom.py +++ b/tessreduce/psf_photom.py @@ -162,7 +162,7 @@ def psf_position(self,image,error,limx=0.8,limy=0.8,ext_shift=[0,0]):#,surface=F error = np.ones_like(image) #brightloc = if np.nansum(image) > 0: - if np.isfinite(ext_shift).all(): + if not np.isfinite(ext_shift).all(): ext_shift[0] = 0; ext_shift[1] = 0 normimage = image / np.nansum(image) # normalise the image @@ -215,10 +215,9 @@ def minimize_psf_flux(self,coeff,image,error=None,surface=True,order=2,kernel=No s = polynomial_surface(xx,yy,plane_coeff,order) else: s = 0 - if kernel is not None: - self.psf = fftconvolve(self.psf, kernel, mode='same') + psf = fftconvolve(self.psf, kernel, mode='same') if kernel is not None else self.psf - res = np.nansum((image - self.psf*coeff[0] - s)**2/error) + res = np.nansum((image - psf*coeff[0] - s)**2/error) return res def psf_flux(self,image,error=None,ext_shift=None,surface=True,poly_order=3,kernel=None): diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 99b71bf..d88893b 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -2763,7 +2763,7 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, self.ref = deepcopy(self.flux[self.ref_ind]) elif self._ref_type.lower() == 'stack': self.stack_ref() - self.ref -= np.nanmin(self.ref) + self.ref -= np.nanpercentile(self.ref, 1) self.flux -= self.ref # self.ref -= self.bkg[self.ref_ind] From 1cfb0872a00a0bd57dd070c95148c94fab49018a Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 4 Jun 2026 11:43:02 +1200 Subject: [PATCH 07/22] Switch all Parallel calls to prefer='threads' to fix serial fallback on HPC On OzSTAR and similar HPC systems, joblib's default loky backend (process-based) can fail silently due to process-spawning restrictions, security policies, or memory limits, causing every Parallel() call to fall back to sequential execution with no error or warning. All heavy functions (Smooth_bkg, shift, inpaint_biharmonic, griddata, Background2D, gaussian_filter, sigma_clipped_stats) are implemented in C extensions that release the GIL, so thread-based parallelism gives genuine speedup without requiring process spawning. Changed 33 Parallel() calls across tessreduce.py, helpers.py, sep_aligner.py, adaptive_background.py, lastpercent.py, background.py, and rescale_straps.py. Also fixes the standalone _fit_residual_surface call which was still using n_jobs=-1 without SLURM awareness. --- tessreduce/adaptive_background.py | 6 ++-- tessreduce/background.py | 2 +- tessreduce/helpers.py | 8 +++--- tessreduce/lastpercent.py | 4 +-- tessreduce/rescale_straps.py | 2 +- tessreduce/sep_aligner.py | 2 +- tessreduce/tessreduce.py | 46 +++++++++++++++---------------- 7 files changed, 35 insertions(+), 35 deletions(-) diff --git a/tessreduce/adaptive_background.py b/tessreduce/adaptive_background.py index b636cf2..76aaf1c 100644 --- a/tessreduce/adaptive_background.py +++ b/tessreduce/adaptive_background.py @@ -236,7 +236,7 @@ def _compute_dev(cw, data_metric=data_metric, segments=segments, gw=gw): dev[s:e] = median_filter(diff, size=(dw, 1, 1), mode='reflect') return dev - all_devs = Parallel(n_jobs=n_jobs)(delayed(_compute_dev)(cw) for cw in scales) + all_devs = Parallel(n_jobs=n_jobs, prefer="threads")(delayed(_compute_dev)(cw) for cw in scales) _frame_mean = data_metric.mean(axis=(1, 2)) _clipped = _frame_mean[np.isfinite(_frame_mean)] @@ -280,7 +280,7 @@ def _compute_norm(i_dev, bright_mask=bright_mask, segments=segments, norm_scale = norm_scale * bright_mask return norm_scale - scale_norms = Parallel(n_jobs=n_jobs)( + scale_norms = Parallel(n_jobs=n_jobs, prefer="threads")( delayed(_compute_norm)((i, dev)) for i, dev in enumerate(all_devs) ) @@ -377,7 +377,7 @@ def _norm01(x): def _smooth_seg(w, seg=seg_data): return w, median_filter(seg, size=(w, 1, 1), mode='reflect') - for w, smoothed_w in Parallel(n_jobs=n_jobs)(delayed(_smooth_seg)(w) for w in seg_levels): + for w, smoothed_w in Parallel(n_jobs=n_jobs, prefer="threads")(delayed(_smooth_seg)(w) for w in seg_levels): result[s:e][seg_wins == w] = smoothed_w[seg_wins == w] if sigma_clip is not None: diff --git a/tessreduce/background.py b/tessreduce/background.py index 751c806..987a3a2 100755 --- a/tessreduce/background.py +++ b/tessreduce/background.py @@ -80,7 +80,7 @@ def _smooth_wrapper(self): if self.parallel: num_cores = multiprocessing.cpu_count() - bkg_smth = Parallel(n_jobs=num_cores)( + bkg_smth = Parallel(n_jobs=num_cores, prefer="threads")( delayed(Smooth_bkg)(frame) for frame in flux * m) else: bkg_smth = np.zeros_like(flux) * np.nan diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index e905649..104f57d 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -197,7 +197,7 @@ def parallel_background2d(cube, box_size=5, filter_size=3, sigma=3, maxiters=5, from astropy.stats import SigmaClip sc = SigmaClip(sigma=sigma, maxiters=maxiters) estimator = MedianBackground() - return np.array(Parallel(n_jobs=n_jobs)( + return np.array(Parallel(n_jobs=n_jobs, prefer="threads")( delayed(_background2d_frame)(frame, box_size, filter_size, sc, estimator, mask) for frame in cube)) @@ -1130,7 +1130,7 @@ def regional_stats_mask(image, size=90, sigma=3, iters=10, n_jobs=1): region_pixels = [np.where(regions == i) for i in range(max_reg + 1)] clip = np.zeros_like(image) - results = Parallel(n_jobs=n_jobs)( + results = Parallel(n_jobs=n_jobs, prefer="threads")( delayed(_clip_region)(image, rx, ry, sigma, iters) for rx, ry in region_pixels) for rx_cut, ry_cut in results: @@ -1355,7 +1355,7 @@ def parallel_strap_fit(frame, frame_bkg, frame_err, mask, repeats=3, tol=3, n_jo qe = np.ones_like(frame) if len(sind) == 0: return qe - results = Parallel(n_jobs=n_jobs)( + results = Parallel(n_jobs=n_jobs, prefer="threads")( delayed(_strap_fit_col)(frame[:, i], norm[:, i]) for i in sind) for col, q in zip(sind, results): @@ -1844,7 +1844,7 @@ def blend_dynamic_background(bkg_new, bkg_prev, flux, sigma=2.0, sharp_masks=Non sharp_list = [sharp_masks[i] if sharp_masks is not None else None for i in range(T)] - results = Parallel(n_jobs=n_jobs)( + results = Parallel(n_jobs=n_jobs, prefer="threads")( delayed(_blend_frame)(bkg_new[i], bkg_prev[i], delta[i], resid_prev[i], sigma, sharp_list[i], gauss_kernel) for i in range(T)) diff --git a/tessreduce/lastpercent.py b/tessreduce/lastpercent.py index 2a9486c..38dacb4 100644 --- a/tessreduce/lastpercent.py +++ b/tessreduce/lastpercent.py @@ -55,7 +55,7 @@ def _find_bkg_cor(tess,cores): coord = np.c_[y,x] cors = np.zeros_like(tess.ref) - cor = Parallel(n_jobs=cores)(delayed(_parallel_correlation) + cor = Parallel(n_jobs=cores, prefer="threads")(delayed(_parallel_correlation) (tess.flux[:,coord[i,0],coord[i,1]], tess.bkg[:,coord[i,0],coord[i,1]], cors,coord[i],30) for i in range(len(coord))) @@ -156,7 +156,7 @@ def multi_correlation_cor(tess, limit=0.8, cores=7): if len(y) == 0: return flux, bkg - results = Parallel(n_jobs=cores)( + results = Parallel(n_jobs=cores, prefer="threads")( delayed(_correct_pixel_correlation)( tess.flux[:, y[i], x[i]], tess.bkg[:, y[i], x[i]], diff --git a/tessreduce/rescale_straps.py b/tessreduce/rescale_straps.py index 85dfc08..6d91832 100755 --- a/tessreduce/rescale_straps.py +++ b/tessreduce/rescale_straps.py @@ -123,7 +123,7 @@ def correct_straps(Image,mask,av_size=5,parallel=True): if parallel: num_cores = multiprocessing.cpu_count() x = np.arange(0,len(breaks),dtype=int) - qe = np.array(Parallel(n_jobs=num_cores)(delayed(calc_strap_factor)(i,breaks,size,av_size,normals,data) for i in x)) + qe = np.array(Parallel(n_jobs=num_cores, prefer="threads")(delayed(calc_strap_factor)(i,breaks,size,av_size,normals,data) for i in x)) qe = np.nanmedian(qe,axis=0) qe[np.isnan(qe)] = 1 else: diff --git a/tessreduce/sep_aligner.py b/tessreduce/sep_aligner.py index 4b58475..c465c62 100644 --- a/tessreduce/sep_aligner.py +++ b/tessreduce/sep_aligner.py @@ -703,7 +703,7 @@ def run(self, time: Optional[np.ndarray] = None, self._sub_ref, cores, positions, np.asarray(weights), p['core_half']) results = Parallel(n_jobs=self.n_jobs, verbose=verbose, - backend='loky')( + prefer="threads")( delayed(_align_one_frame)( t, self.flux[t], ref_comp, w_cols, diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index d88893b..4a3bb13 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -125,9 +125,9 @@ def _fit_frame(residual): except Exception: return np.full_like(residual, np.nanmedian(residual[~transient_mask])) - n_jobs = -1 + n_jobs = _available_cores() residuals = flux - bkg - corrections = Parallel(n_jobs=n_jobs)(delayed(_fit_frame)(residuals[i]) for i in range(flux.shape[0])) + corrections = Parallel(n_jobs=n_jobs, prefer="threads")(delayed(_fit_frame)(residuals[i]) for i in range(flux.shape[0])) bkg += np.array(corrections) return bkg @@ -651,7 +651,7 @@ def psf_source_mask(self,sigma=5): data = (self._flux_aligned - self.ref) #* mask if self.parallel: try: - m = Parallel(n_jobs=self.num_cores)(delayed(par_psf_source_mask)(frame,self.prf,sigma) for frame in data) + m = Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_source_mask)(frame,self.prf,sigma) for frame in data) m = np.array(m) except: m = np.ones_like(data) @@ -778,7 +778,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa bkg_smth = np.zeros_like(flux) * np.nan if self.parallel: _t = time.perf_counter() - bkg_smth = Parallel(n_jobs=self.num_cores)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*m) + bkg_smth = Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*m) _times['initial smooth background'] = time.perf_counter() - _t if rerun_negative: _t = time.perf_counter() @@ -800,7 +800,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa else: m[over_sub] = 1 self._bkgmask = m - bkg_smth = Parallel(n_jobs=self.num_cores)(delayed(Smooth_bkg)(frame,gauss_smooth,interpolate) for frame in flux*m) + bkg_smth = Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(Smooth_bkg)(frame,gauss_smooth,interpolate) for frame in flux*m) _times['negative over-subtraction rerun'] = time.perf_counter() - _t @@ -853,7 +853,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa new_mask = abs(new_mask - 1) self._bkgmask = new_mask bkg_s1 = np.array(bkg_smth) - bkg_smth = Parallel(n_jobs=self.num_cores)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*new_mask) + bkg_smth = Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*new_mask) if blend_dynamic: bkg_smth = blend_dynamic_background(bkg_smth, bkg_s1, flux, n_jobs=self.num_cores) _times['residual surface rerun'] = time.perf_counter() - _t @@ -996,7 +996,7 @@ def _bkg_round_3(self,iters=5): kern = np.ones((1,3,3)) dist_mask = convolve(dist_mask,kern) > 0 if self.parallel: - bkg_3 = Parallel(n_jobs=self.num_cores)(delayed(parallel_bkg3)(self.bkg[i],dist_mask[i]) + bkg_3 = Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(parallel_bkg3)(self.bkg[i],dist_mask[i]) for i in np.arange(len(dist_mask))) else: bkg_3 = np.zeros_like(self.bkg) @@ -1020,7 +1020,7 @@ def _clip_background(self,sigma=5,ideal_size=90): """ if self.parallel: - bkg_clip = Parallel(n_jobs=self.num_cores)(delayed(clip_background)(self.bkg[i],self.mask,sigma,ideal_size) + bkg_clip = Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(clip_background)(self.bkg[i],self.mask,sigma,ideal_size) for i in np.arange(len(self.bkg))) else: bkg_clip = np.zeros_like(self.bkg) @@ -1046,7 +1046,7 @@ def _grad_bkg_clip(self,sigma=3,max_size=1000): """ if self.parallel: - bkg_clip = Parallel(n_jobs=self.num_cores)(delayed(grad_clip_fill_bkg)(self.bkg[i],sigma,max_size) + bkg_clip = Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(grad_clip_fill_bkg)(self.bkg[i],sigma,max_size) for i in np.arange(len(self.bkg))) else: bkg_clip = np.zeros_like(self.bkg) @@ -1402,7 +1402,7 @@ def centroids_shifts_starfind(self,plot=None,savename=None): self._dat_sources = s.to_pandas() if self.parallel: - shifts = Parallel(n_jobs=self.num_cores)( + shifts = Parallel(n_jobs=self.num_cores, prefer="threads")( delayed(Calculate_shifts)(frame,mx,my,finder) for frame in f) shifts = np.array(shifts) else: @@ -1474,7 +1474,7 @@ def fit_shift(self,smooth=True,plot=None,savename=None): if self.parallel: ind = np.arange(len(f)) - shifts = Parallel(n_jobs=self.num_cores)( + shifts = Parallel(n_jobs=self.num_cores, prefer="threads")( #delayed(difference_shifts)(f[i],m,self.eflux[i],eref) for i in ind) delayed(difference_shifts)(f[i],m) for i in ind) shifts = np.array(shifts) @@ -1548,7 +1548,7 @@ def shift_images(self,median=False): shifted[nans] = 0. if median: if self.parallel: - result = Parallel(n_jobs=self.num_cores)( + result = Parallel(n_jobs=self.num_cores, prefer="threads")( delayed(_shift_ref_one)(self.ref, shifted[i], self.shift[i]) for i in range(len(shifted))) shifted = np.array(result) @@ -1559,7 +1559,7 @@ def shift_images(self,median=False): else: if self.parallel: - result = Parallel(n_jobs=self.num_cores)( + result = Parallel(n_jobs=self.num_cores, prefer="threads")( delayed(_shift_one)(shifted[i], self.shift[i]) for i in range(len(shifted))) self.flux = np.array(result) @@ -2267,7 +2267,7 @@ def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2): raise ValueError(m) inds = np.arange(0,len(xpos)) if self.parallel: - prfs, cutouts, ecutouts = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_initialise)(self.flux,self.tpf.camera,self.tpf.ccd, + prfs, cutouts, ecutouts = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_initialise)(self.flux,self.tpf.camera,self.tpf.ccd, self.tpf.sector,self.tpf.column,self.tpf.row, size,[xpos[i],ypos[i]],time_ind) for i in inds)) else: @@ -2280,7 +2280,7 @@ def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2): cutouts = np.array(cutouts) print('made cutouts') if self.parallel: - flux, pos = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_full)(cutouts[i],prfs[i],self.shift[i],xlim,ylim) for i in inds)) + flux, pos = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_full)(cutouts[i],prfs[i],self.shift[i],xlim,ylim) for i in inds)) else: flux = [] pos = [] @@ -2366,14 +2366,14 @@ def psf_photutils(self,xPix=None,yPix=None,size=5,local_bkg=False,epsf=None, eflux = np.zeros(len(self.flux)) * np.nan psfphot2 = PSFPhotometry(epsf, fit_shape, finder=None,aperture_radius=1.5, xy_bounds=(0.05),localbkg_estimator=localbkg_estimator) - f,ef = zip(*Parallel(n_jobs=self.num_cores)(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot2,init) for i in np.arange(len(cutouts)))) + f,ef = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot2,init) for i in np.arange(len(cutouts)))) f = np.array(f).flatten() ef = np.array(ef).flatten() phot = phot.to_pandas() pos = phot[['x_fit','y_fit']].values + np.array([xPix,yPix]) - size//2 epos = phot[['x_err','y_err']].values else: - f,ef,pos,epos = zip(*Parallel(n_jobs=self.num_cores)(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot,init,True) for i in np.arange(len(cutouts)))) + f,ef,pos,epos = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot,init,True) for i in np.arange(len(cutouts)))) pos['x_fit'] += xPix - size//2 pos['y_fit'] += xPix - size//2 @@ -2437,7 +2437,7 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa prf, cutouts, ecutouts = self._psf_initialise(size,(xPix,yPix),ref=(not diff)) # gather base PRF and the array of cutouts data inds = np.arange(len(cutouts)) base = create_psf(prf,size) - flux, eflux, pos = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_full)(cutouts[i],base,self.shift[i]) for i in inds)) + flux, eflux, pos = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_full)(cutouts[i],base,self.shift[i]) for i in inds)) #prf, cutouts = self._psf_initialise(size,(xPix,yPix)) # gather base PRF and the array of cutouts data #xShifts = [] @@ -2483,9 +2483,9 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa if self.parallel: inds = np.arange(len(cutouts)) if self.delta_kernel is not None: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) else: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) else: for i in range(len(cutouts)): flux += [par_psf_flux(cutouts[i],ecutouts[i],base,self.shift[i])] @@ -2501,9 +2501,9 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa if self.parallel: inds = np.arange(len(cutouts)) if self.delta_kernel is not None: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) else: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) else: for i in range(len(cutouts)): flux += [par_psf_flux(cutouts[i],ecutouts[i],base,self.shift[i])] @@ -2548,7 +2548,7 @@ def kernel_matching(self,size=7,diff=True): mask = self.mask == 1 if self.parallel: - d, kernel = zip(*Parallel(n_jobs=self.num_cores)(delayed(parallel_delta_diff)(frame,self.ref,mask,size) for frame in flux)) + d, kernel = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(parallel_delta_diff)(frame,self.ref,mask,size) for frame in flux)) else: d = [] kernel = [] From abbe4db95ee68a455ea496906947fb94c8c93f62 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 4 Jun 2026 11:56:45 +1200 Subject: [PATCH 08/22] Switch from prefer='threads' to backend='multiprocessing' (fork-based) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit prefer='threads' caused high system CPU time on OzSTAR because all threads share one address space: simultaneous malloc/free calls from numpy temporary arrays contend on the glibc allocator lock, and cross-socket memory access on NUMA hardware requires OS arbitration — both appearing as system rather than user CPU time. backend='multiprocessing' uses Linux fork() directly. Each worker gets its own address space via copy-on-write, eliminating allocator contention and NUMA cross-talk. Unlike the loky backend (which uses a forkserver that was failing silently on OzSTAR), the multiprocessing backend uses a direct fork() call that works reliably in SLURM environments. --- tessreduce/adaptive_background.py | 6 ++--- tessreduce/background.py | 2 +- tessreduce/helpers.py | 8 +++--- tessreduce/lastpercent.py | 4 +-- tessreduce/rescale_straps.py | 2 +- tessreduce/sep_aligner.py | 2 +- tessreduce/tessreduce.py | 44 +++++++++++++++---------------- 7 files changed, 34 insertions(+), 34 deletions(-) diff --git a/tessreduce/adaptive_background.py b/tessreduce/adaptive_background.py index 76aaf1c..64a9dd8 100644 --- a/tessreduce/adaptive_background.py +++ b/tessreduce/adaptive_background.py @@ -236,7 +236,7 @@ def _compute_dev(cw, data_metric=data_metric, segments=segments, gw=gw): dev[s:e] = median_filter(diff, size=(dw, 1, 1), mode='reflect') return dev - all_devs = Parallel(n_jobs=n_jobs, prefer="threads")(delayed(_compute_dev)(cw) for cw in scales) + all_devs = Parallel(n_jobs=n_jobs, backend="multiprocessing")(delayed(_compute_dev)(cw) for cw in scales) _frame_mean = data_metric.mean(axis=(1, 2)) _clipped = _frame_mean[np.isfinite(_frame_mean)] @@ -280,7 +280,7 @@ def _compute_norm(i_dev, bright_mask=bright_mask, segments=segments, norm_scale = norm_scale * bright_mask return norm_scale - scale_norms = Parallel(n_jobs=n_jobs, prefer="threads")( + scale_norms = Parallel(n_jobs=n_jobs, backend="multiprocessing")( delayed(_compute_norm)((i, dev)) for i, dev in enumerate(all_devs) ) @@ -377,7 +377,7 @@ def _norm01(x): def _smooth_seg(w, seg=seg_data): return w, median_filter(seg, size=(w, 1, 1), mode='reflect') - for w, smoothed_w in Parallel(n_jobs=n_jobs, prefer="threads")(delayed(_smooth_seg)(w) for w in seg_levels): + for w, smoothed_w in Parallel(n_jobs=n_jobs, backend="multiprocessing")(delayed(_smooth_seg)(w) for w in seg_levels): result[s:e][seg_wins == w] = smoothed_w[seg_wins == w] if sigma_clip is not None: diff --git a/tessreduce/background.py b/tessreduce/background.py index 987a3a2..706b608 100755 --- a/tessreduce/background.py +++ b/tessreduce/background.py @@ -80,7 +80,7 @@ def _smooth_wrapper(self): if self.parallel: num_cores = multiprocessing.cpu_count() - bkg_smth = Parallel(n_jobs=num_cores, prefer="threads")( + bkg_smth = Parallel(n_jobs=num_cores, backend="multiprocessing")( delayed(Smooth_bkg)(frame) for frame in flux * m) else: bkg_smth = np.zeros_like(flux) * np.nan diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index 104f57d..fca9ac7 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -197,7 +197,7 @@ def parallel_background2d(cube, box_size=5, filter_size=3, sigma=3, maxiters=5, from astropy.stats import SigmaClip sc = SigmaClip(sigma=sigma, maxiters=maxiters) estimator = MedianBackground() - return np.array(Parallel(n_jobs=n_jobs, prefer="threads")( + return np.array(Parallel(n_jobs=n_jobs, backend="multiprocessing")( delayed(_background2d_frame)(frame, box_size, filter_size, sc, estimator, mask) for frame in cube)) @@ -1130,7 +1130,7 @@ def regional_stats_mask(image, size=90, sigma=3, iters=10, n_jobs=1): region_pixels = [np.where(regions == i) for i in range(max_reg + 1)] clip = np.zeros_like(image) - results = Parallel(n_jobs=n_jobs, prefer="threads")( + results = Parallel(n_jobs=n_jobs, backend="multiprocessing")( delayed(_clip_region)(image, rx, ry, sigma, iters) for rx, ry in region_pixels) for rx_cut, ry_cut in results: @@ -1355,7 +1355,7 @@ def parallel_strap_fit(frame, frame_bkg, frame_err, mask, repeats=3, tol=3, n_jo qe = np.ones_like(frame) if len(sind) == 0: return qe - results = Parallel(n_jobs=n_jobs, prefer="threads")( + results = Parallel(n_jobs=n_jobs, backend="multiprocessing")( delayed(_strap_fit_col)(frame[:, i], norm[:, i]) for i in sind) for col, q in zip(sind, results): @@ -1844,7 +1844,7 @@ def blend_dynamic_background(bkg_new, bkg_prev, flux, sigma=2.0, sharp_masks=Non sharp_list = [sharp_masks[i] if sharp_masks is not None else None for i in range(T)] - results = Parallel(n_jobs=n_jobs, prefer="threads")( + results = Parallel(n_jobs=n_jobs, backend="multiprocessing")( delayed(_blend_frame)(bkg_new[i], bkg_prev[i], delta[i], resid_prev[i], sigma, sharp_list[i], gauss_kernel) for i in range(T)) diff --git a/tessreduce/lastpercent.py b/tessreduce/lastpercent.py index 38dacb4..134c0ee 100644 --- a/tessreduce/lastpercent.py +++ b/tessreduce/lastpercent.py @@ -55,7 +55,7 @@ def _find_bkg_cor(tess,cores): coord = np.c_[y,x] cors = np.zeros_like(tess.ref) - cor = Parallel(n_jobs=cores, prefer="threads")(delayed(_parallel_correlation) + cor = Parallel(n_jobs=cores, backend="multiprocessing")(delayed(_parallel_correlation) (tess.flux[:,coord[i,0],coord[i,1]], tess.bkg[:,coord[i,0],coord[i,1]], cors,coord[i],30) for i in range(len(coord))) @@ -156,7 +156,7 @@ def multi_correlation_cor(tess, limit=0.8, cores=7): if len(y) == 0: return flux, bkg - results = Parallel(n_jobs=cores, prefer="threads")( + results = Parallel(n_jobs=cores, backend="multiprocessing")( delayed(_correct_pixel_correlation)( tess.flux[:, y[i], x[i]], tess.bkg[:, y[i], x[i]], diff --git a/tessreduce/rescale_straps.py b/tessreduce/rescale_straps.py index 6d91832..31d3957 100755 --- a/tessreduce/rescale_straps.py +++ b/tessreduce/rescale_straps.py @@ -123,7 +123,7 @@ def correct_straps(Image,mask,av_size=5,parallel=True): if parallel: num_cores = multiprocessing.cpu_count() x = np.arange(0,len(breaks),dtype=int) - qe = np.array(Parallel(n_jobs=num_cores, prefer="threads")(delayed(calc_strap_factor)(i,breaks,size,av_size,normals,data) for i in x)) + qe = np.array(Parallel(n_jobs=num_cores, backend="multiprocessing")(delayed(calc_strap_factor)(i,breaks,size,av_size,normals,data) for i in x)) qe = np.nanmedian(qe,axis=0) qe[np.isnan(qe)] = 1 else: diff --git a/tessreduce/sep_aligner.py b/tessreduce/sep_aligner.py index c465c62..0f14abe 100644 --- a/tessreduce/sep_aligner.py +++ b/tessreduce/sep_aligner.py @@ -703,7 +703,7 @@ def run(self, time: Optional[np.ndarray] = None, self._sub_ref, cores, positions, np.asarray(weights), p['core_half']) results = Parallel(n_jobs=self.n_jobs, verbose=verbose, - prefer="threads")( + backend="multiprocessing")( delayed(_align_one_frame)( t, self.flux[t], ref_comp, w_cols, diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 4a3bb13..01b9fb8 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -127,7 +127,7 @@ def _fit_frame(residual): n_jobs = _available_cores() residuals = flux - bkg - corrections = Parallel(n_jobs=n_jobs, prefer="threads")(delayed(_fit_frame)(residuals[i]) for i in range(flux.shape[0])) + corrections = Parallel(n_jobs=n_jobs, backend="multiprocessing")(delayed(_fit_frame)(residuals[i]) for i in range(flux.shape[0])) bkg += np.array(corrections) return bkg @@ -651,7 +651,7 @@ def psf_source_mask(self,sigma=5): data = (self._flux_aligned - self.ref) #* mask if self.parallel: try: - m = Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_source_mask)(frame,self.prf,sigma) for frame in data) + m = Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_source_mask)(frame,self.prf,sigma) for frame in data) m = np.array(m) except: m = np.ones_like(data) @@ -778,7 +778,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa bkg_smth = np.zeros_like(flux) * np.nan if self.parallel: _t = time.perf_counter() - bkg_smth = Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*m) + bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*m) _times['initial smooth background'] = time.perf_counter() - _t if rerun_negative: _t = time.perf_counter() @@ -800,7 +800,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa else: m[over_sub] = 1 self._bkgmask = m - bkg_smth = Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(Smooth_bkg)(frame,gauss_smooth,interpolate) for frame in flux*m) + bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(Smooth_bkg)(frame,gauss_smooth,interpolate) for frame in flux*m) _times['negative over-subtraction rerun'] = time.perf_counter() - _t @@ -853,7 +853,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa new_mask = abs(new_mask - 1) self._bkgmask = new_mask bkg_s1 = np.array(bkg_smth) - bkg_smth = Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*new_mask) + bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*new_mask) if blend_dynamic: bkg_smth = blend_dynamic_background(bkg_smth, bkg_s1, flux, n_jobs=self.num_cores) _times['residual surface rerun'] = time.perf_counter() - _t @@ -996,7 +996,7 @@ def _bkg_round_3(self,iters=5): kern = np.ones((1,3,3)) dist_mask = convolve(dist_mask,kern) > 0 if self.parallel: - bkg_3 = Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(parallel_bkg3)(self.bkg[i],dist_mask[i]) + bkg_3 = Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(parallel_bkg3)(self.bkg[i],dist_mask[i]) for i in np.arange(len(dist_mask))) else: bkg_3 = np.zeros_like(self.bkg) @@ -1020,7 +1020,7 @@ def _clip_background(self,sigma=5,ideal_size=90): """ if self.parallel: - bkg_clip = Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(clip_background)(self.bkg[i],self.mask,sigma,ideal_size) + bkg_clip = Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(clip_background)(self.bkg[i],self.mask,sigma,ideal_size) for i in np.arange(len(self.bkg))) else: bkg_clip = np.zeros_like(self.bkg) @@ -1046,7 +1046,7 @@ def _grad_bkg_clip(self,sigma=3,max_size=1000): """ if self.parallel: - bkg_clip = Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(grad_clip_fill_bkg)(self.bkg[i],sigma,max_size) + bkg_clip = Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(grad_clip_fill_bkg)(self.bkg[i],sigma,max_size) for i in np.arange(len(self.bkg))) else: bkg_clip = np.zeros_like(self.bkg) @@ -1402,7 +1402,7 @@ def centroids_shifts_starfind(self,plot=None,savename=None): self._dat_sources = s.to_pandas() if self.parallel: - shifts = Parallel(n_jobs=self.num_cores, prefer="threads")( + shifts = Parallel(n_jobs=self.num_cores, backend="multiprocessing")( delayed(Calculate_shifts)(frame,mx,my,finder) for frame in f) shifts = np.array(shifts) else: @@ -1474,7 +1474,7 @@ def fit_shift(self,smooth=True,plot=None,savename=None): if self.parallel: ind = np.arange(len(f)) - shifts = Parallel(n_jobs=self.num_cores, prefer="threads")( + shifts = Parallel(n_jobs=self.num_cores, backend="multiprocessing")( #delayed(difference_shifts)(f[i],m,self.eflux[i],eref) for i in ind) delayed(difference_shifts)(f[i],m) for i in ind) shifts = np.array(shifts) @@ -1548,7 +1548,7 @@ def shift_images(self,median=False): shifted[nans] = 0. if median: if self.parallel: - result = Parallel(n_jobs=self.num_cores, prefer="threads")( + result = Parallel(n_jobs=self.num_cores, backend="multiprocessing")( delayed(_shift_ref_one)(self.ref, shifted[i], self.shift[i]) for i in range(len(shifted))) shifted = np.array(result) @@ -1559,7 +1559,7 @@ def shift_images(self,median=False): else: if self.parallel: - result = Parallel(n_jobs=self.num_cores, prefer="threads")( + result = Parallel(n_jobs=self.num_cores, backend="multiprocessing")( delayed(_shift_one)(shifted[i], self.shift[i]) for i in range(len(shifted))) self.flux = np.array(result) @@ -2267,7 +2267,7 @@ def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2): raise ValueError(m) inds = np.arange(0,len(xpos)) if self.parallel: - prfs, cutouts, ecutouts = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_initialise)(self.flux,self.tpf.camera,self.tpf.ccd, + prfs, cutouts, ecutouts = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_initialise)(self.flux,self.tpf.camera,self.tpf.ccd, self.tpf.sector,self.tpf.column,self.tpf.row, size,[xpos[i],ypos[i]],time_ind) for i in inds)) else: @@ -2280,7 +2280,7 @@ def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2): cutouts = np.array(cutouts) print('made cutouts') if self.parallel: - flux, pos = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_full)(cutouts[i],prfs[i],self.shift[i],xlim,ylim) for i in inds)) + flux, pos = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_full)(cutouts[i],prfs[i],self.shift[i],xlim,ylim) for i in inds)) else: flux = [] pos = [] @@ -2366,14 +2366,14 @@ def psf_photutils(self,xPix=None,yPix=None,size=5,local_bkg=False,epsf=None, eflux = np.zeros(len(self.flux)) * np.nan psfphot2 = PSFPhotometry(epsf, fit_shape, finder=None,aperture_radius=1.5, xy_bounds=(0.05),localbkg_estimator=localbkg_estimator) - f,ef = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot2,init) for i in np.arange(len(cutouts)))) + f,ef = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot2,init) for i in np.arange(len(cutouts)))) f = np.array(f).flatten() ef = np.array(ef).flatten() phot = phot.to_pandas() pos = phot[['x_fit','y_fit']].values + np.array([xPix,yPix]) - size//2 epos = phot[['x_err','y_err']].values else: - f,ef,pos,epos = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot,init,True) for i in np.arange(len(cutouts)))) + f,ef,pos,epos = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot,init,True) for i in np.arange(len(cutouts)))) pos['x_fit'] += xPix - size//2 pos['y_fit'] += xPix - size//2 @@ -2437,7 +2437,7 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa prf, cutouts, ecutouts = self._psf_initialise(size,(xPix,yPix),ref=(not diff)) # gather base PRF and the array of cutouts data inds = np.arange(len(cutouts)) base = create_psf(prf,size) - flux, eflux, pos = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_full)(cutouts[i],base,self.shift[i]) for i in inds)) + flux, eflux, pos = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_full)(cutouts[i],base,self.shift[i]) for i in inds)) #prf, cutouts = self._psf_initialise(size,(xPix,yPix)) # gather base PRF and the array of cutouts data #xShifts = [] @@ -2483,9 +2483,9 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa if self.parallel: inds = np.arange(len(cutouts)) if self.delta_kernel is not None: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) else: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) else: for i in range(len(cutouts)): flux += [par_psf_flux(cutouts[i],ecutouts[i],base,self.shift[i])] @@ -2501,9 +2501,9 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa if self.parallel: inds = np.arange(len(cutouts)) if self.delta_kernel is not None: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) else: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) else: for i in range(len(cutouts)): flux += [par_psf_flux(cutouts[i],ecutouts[i],base,self.shift[i])] @@ -2548,7 +2548,7 @@ def kernel_matching(self,size=7,diff=True): mask = self.mask == 1 if self.parallel: - d, kernel = zip(*Parallel(n_jobs=self.num_cores, prefer="threads")(delayed(parallel_delta_diff)(frame,self.ref,mask,size) for frame in flux)) + d, kernel = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(parallel_delta_diff)(frame,self.ref,mask,size) for frame in flux)) else: d = [] kernel = [] From f65e29d3be5ac7f0a7e15083e78a6df7d3273763 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 4 Jun 2026 12:21:47 +1200 Subject: [PATCH 09/22] Add verbose=1 to all Parallel calls for HPC diagnostics --- tessreduce/adaptive_background.py | 6 ++--- tessreduce/background.py | 2 +- tessreduce/helpers.py | 8 +++--- tessreduce/lastpercent.py | 4 +-- tessreduce/rescale_straps.py | 2 +- tessreduce/sep_aligner.py | 2 +- tessreduce/tessreduce.py | 44 +++++++++++++++---------------- 7 files changed, 34 insertions(+), 34 deletions(-) diff --git a/tessreduce/adaptive_background.py b/tessreduce/adaptive_background.py index 64a9dd8..2c09ab5 100644 --- a/tessreduce/adaptive_background.py +++ b/tessreduce/adaptive_background.py @@ -236,7 +236,7 @@ def _compute_dev(cw, data_metric=data_metric, segments=segments, gw=gw): dev[s:e] = median_filter(diff, size=(dw, 1, 1), mode='reflect') return dev - all_devs = Parallel(n_jobs=n_jobs, backend="multiprocessing")(delayed(_compute_dev)(cw) for cw in scales) + all_devs = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)(delayed(_compute_dev)(cw) for cw in scales) _frame_mean = data_metric.mean(axis=(1, 2)) _clipped = _frame_mean[np.isfinite(_frame_mean)] @@ -280,7 +280,7 @@ def _compute_norm(i_dev, bright_mask=bright_mask, segments=segments, norm_scale = norm_scale * bright_mask return norm_scale - scale_norms = Parallel(n_jobs=n_jobs, backend="multiprocessing")( + scale_norms = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( delayed(_compute_norm)((i, dev)) for i, dev in enumerate(all_devs) ) @@ -377,7 +377,7 @@ def _norm01(x): def _smooth_seg(w, seg=seg_data): return w, median_filter(seg, size=(w, 1, 1), mode='reflect') - for w, smoothed_w in Parallel(n_jobs=n_jobs, backend="multiprocessing")(delayed(_smooth_seg)(w) for w in seg_levels): + for w, smoothed_w in Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)(delayed(_smooth_seg)(w) for w in seg_levels): result[s:e][seg_wins == w] = smoothed_w[seg_wins == w] if sigma_clip is not None: diff --git a/tessreduce/background.py b/tessreduce/background.py index 706b608..c398831 100755 --- a/tessreduce/background.py +++ b/tessreduce/background.py @@ -80,7 +80,7 @@ def _smooth_wrapper(self): if self.parallel: num_cores = multiprocessing.cpu_count() - bkg_smth = Parallel(n_jobs=num_cores, backend="multiprocessing")( + bkg_smth = Parallel(n_jobs=num_cores, backend="multiprocessing", verbose=1)( delayed(Smooth_bkg)(frame) for frame in flux * m) else: bkg_smth = np.zeros_like(flux) * np.nan diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index fca9ac7..60cc21d 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -197,7 +197,7 @@ def parallel_background2d(cube, box_size=5, filter_size=3, sigma=3, maxiters=5, from astropy.stats import SigmaClip sc = SigmaClip(sigma=sigma, maxiters=maxiters) estimator = MedianBackground() - return np.array(Parallel(n_jobs=n_jobs, backend="multiprocessing")( + return np.array(Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( delayed(_background2d_frame)(frame, box_size, filter_size, sc, estimator, mask) for frame in cube)) @@ -1130,7 +1130,7 @@ def regional_stats_mask(image, size=90, sigma=3, iters=10, n_jobs=1): region_pixels = [np.where(regions == i) for i in range(max_reg + 1)] clip = np.zeros_like(image) - results = Parallel(n_jobs=n_jobs, backend="multiprocessing")( + results = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( delayed(_clip_region)(image, rx, ry, sigma, iters) for rx, ry in region_pixels) for rx_cut, ry_cut in results: @@ -1355,7 +1355,7 @@ def parallel_strap_fit(frame, frame_bkg, frame_err, mask, repeats=3, tol=3, n_jo qe = np.ones_like(frame) if len(sind) == 0: return qe - results = Parallel(n_jobs=n_jobs, backend="multiprocessing")( + results = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( delayed(_strap_fit_col)(frame[:, i], norm[:, i]) for i in sind) for col, q in zip(sind, results): @@ -1844,7 +1844,7 @@ def blend_dynamic_background(bkg_new, bkg_prev, flux, sigma=2.0, sharp_masks=Non sharp_list = [sharp_masks[i] if sharp_masks is not None else None for i in range(T)] - results = Parallel(n_jobs=n_jobs, backend="multiprocessing")( + results = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( delayed(_blend_frame)(bkg_new[i], bkg_prev[i], delta[i], resid_prev[i], sigma, sharp_list[i], gauss_kernel) for i in range(T)) diff --git a/tessreduce/lastpercent.py b/tessreduce/lastpercent.py index 134c0ee..624b9d8 100644 --- a/tessreduce/lastpercent.py +++ b/tessreduce/lastpercent.py @@ -55,7 +55,7 @@ def _find_bkg_cor(tess,cores): coord = np.c_[y,x] cors = np.zeros_like(tess.ref) - cor = Parallel(n_jobs=cores, backend="multiprocessing")(delayed(_parallel_correlation) + cor = Parallel(n_jobs=cores, backend="multiprocessing", verbose=1)(delayed(_parallel_correlation) (tess.flux[:,coord[i,0],coord[i,1]], tess.bkg[:,coord[i,0],coord[i,1]], cors,coord[i],30) for i in range(len(coord))) @@ -156,7 +156,7 @@ def multi_correlation_cor(tess, limit=0.8, cores=7): if len(y) == 0: return flux, bkg - results = Parallel(n_jobs=cores, backend="multiprocessing")( + results = Parallel(n_jobs=cores, backend="multiprocessing", verbose=1)( delayed(_correct_pixel_correlation)( tess.flux[:, y[i], x[i]], tess.bkg[:, y[i], x[i]], diff --git a/tessreduce/rescale_straps.py b/tessreduce/rescale_straps.py index 31d3957..c8efdb7 100755 --- a/tessreduce/rescale_straps.py +++ b/tessreduce/rescale_straps.py @@ -123,7 +123,7 @@ def correct_straps(Image,mask,av_size=5,parallel=True): if parallel: num_cores = multiprocessing.cpu_count() x = np.arange(0,len(breaks),dtype=int) - qe = np.array(Parallel(n_jobs=num_cores, backend="multiprocessing")(delayed(calc_strap_factor)(i,breaks,size,av_size,normals,data) for i in x)) + qe = np.array(Parallel(n_jobs=num_cores, backend="multiprocessing", verbose=1)(delayed(calc_strap_factor)(i,breaks,size,av_size,normals,data) for i in x)) qe = np.nanmedian(qe,axis=0) qe[np.isnan(qe)] = 1 else: diff --git a/tessreduce/sep_aligner.py b/tessreduce/sep_aligner.py index 0f14abe..4f27ac3 100644 --- a/tessreduce/sep_aligner.py +++ b/tessreduce/sep_aligner.py @@ -703,7 +703,7 @@ def run(self, time: Optional[np.ndarray] = None, self._sub_ref, cores, positions, np.asarray(weights), p['core_half']) results = Parallel(n_jobs=self.n_jobs, verbose=verbose, - backend="multiprocessing")( + backend="multiprocessing", verbose=1)( delayed(_align_one_frame)( t, self.flux[t], ref_comp, w_cols, diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 01b9fb8..9db216f 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -127,7 +127,7 @@ def _fit_frame(residual): n_jobs = _available_cores() residuals = flux - bkg - corrections = Parallel(n_jobs=n_jobs, backend="multiprocessing")(delayed(_fit_frame)(residuals[i]) for i in range(flux.shape[0])) + corrections = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)(delayed(_fit_frame)(residuals[i]) for i in range(flux.shape[0])) bkg += np.array(corrections) return bkg @@ -651,7 +651,7 @@ def psf_source_mask(self,sigma=5): data = (self._flux_aligned - self.ref) #* mask if self.parallel: try: - m = Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_source_mask)(frame,self.prf,sigma) for frame in data) + m = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_source_mask)(frame,self.prf,sigma) for frame in data) m = np.array(m) except: m = np.ones_like(data) @@ -778,7 +778,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa bkg_smth = np.zeros_like(flux) * np.nan if self.parallel: _t = time.perf_counter() - bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*m) + bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*m) _times['initial smooth background'] = time.perf_counter() - _t if rerun_negative: _t = time.perf_counter() @@ -800,7 +800,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa else: m[over_sub] = 1 self._bkgmask = m - bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(Smooth_bkg)(frame,gauss_smooth,interpolate) for frame in flux*m) + bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(Smooth_bkg)(frame,gauss_smooth,interpolate) for frame in flux*m) _times['negative over-subtraction rerun'] = time.perf_counter() - _t @@ -853,7 +853,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa new_mask = abs(new_mask - 1) self._bkgmask = new_mask bkg_s1 = np.array(bkg_smth) - bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*new_mask) + bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*new_mask) if blend_dynamic: bkg_smth = blend_dynamic_background(bkg_smth, bkg_s1, flux, n_jobs=self.num_cores) _times['residual surface rerun'] = time.perf_counter() - _t @@ -996,7 +996,7 @@ def _bkg_round_3(self,iters=5): kern = np.ones((1,3,3)) dist_mask = convolve(dist_mask,kern) > 0 if self.parallel: - bkg_3 = Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(parallel_bkg3)(self.bkg[i],dist_mask[i]) + bkg_3 = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(parallel_bkg3)(self.bkg[i],dist_mask[i]) for i in np.arange(len(dist_mask))) else: bkg_3 = np.zeros_like(self.bkg) @@ -1020,7 +1020,7 @@ def _clip_background(self,sigma=5,ideal_size=90): """ if self.parallel: - bkg_clip = Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(clip_background)(self.bkg[i],self.mask,sigma,ideal_size) + bkg_clip = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(clip_background)(self.bkg[i],self.mask,sigma,ideal_size) for i in np.arange(len(self.bkg))) else: bkg_clip = np.zeros_like(self.bkg) @@ -1046,7 +1046,7 @@ def _grad_bkg_clip(self,sigma=3,max_size=1000): """ if self.parallel: - bkg_clip = Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(grad_clip_fill_bkg)(self.bkg[i],sigma,max_size) + bkg_clip = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(grad_clip_fill_bkg)(self.bkg[i],sigma,max_size) for i in np.arange(len(self.bkg))) else: bkg_clip = np.zeros_like(self.bkg) @@ -1402,7 +1402,7 @@ def centroids_shifts_starfind(self,plot=None,savename=None): self._dat_sources = s.to_pandas() if self.parallel: - shifts = Parallel(n_jobs=self.num_cores, backend="multiprocessing")( + shifts = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)( delayed(Calculate_shifts)(frame,mx,my,finder) for frame in f) shifts = np.array(shifts) else: @@ -1474,7 +1474,7 @@ def fit_shift(self,smooth=True,plot=None,savename=None): if self.parallel: ind = np.arange(len(f)) - shifts = Parallel(n_jobs=self.num_cores, backend="multiprocessing")( + shifts = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)( #delayed(difference_shifts)(f[i],m,self.eflux[i],eref) for i in ind) delayed(difference_shifts)(f[i],m) for i in ind) shifts = np.array(shifts) @@ -1548,7 +1548,7 @@ def shift_images(self,median=False): shifted[nans] = 0. if median: if self.parallel: - result = Parallel(n_jobs=self.num_cores, backend="multiprocessing")( + result = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)( delayed(_shift_ref_one)(self.ref, shifted[i], self.shift[i]) for i in range(len(shifted))) shifted = np.array(result) @@ -1559,7 +1559,7 @@ def shift_images(self,median=False): else: if self.parallel: - result = Parallel(n_jobs=self.num_cores, backend="multiprocessing")( + result = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)( delayed(_shift_one)(shifted[i], self.shift[i]) for i in range(len(shifted))) self.flux = np.array(result) @@ -2267,7 +2267,7 @@ def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2): raise ValueError(m) inds = np.arange(0,len(xpos)) if self.parallel: - prfs, cutouts, ecutouts = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_initialise)(self.flux,self.tpf.camera,self.tpf.ccd, + prfs, cutouts, ecutouts = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_initialise)(self.flux,self.tpf.camera,self.tpf.ccd, self.tpf.sector,self.tpf.column,self.tpf.row, size,[xpos[i],ypos[i]],time_ind) for i in inds)) else: @@ -2280,7 +2280,7 @@ def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2): cutouts = np.array(cutouts) print('made cutouts') if self.parallel: - flux, pos = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_full)(cutouts[i],prfs[i],self.shift[i],xlim,ylim) for i in inds)) + flux, pos = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_full)(cutouts[i],prfs[i],self.shift[i],xlim,ylim) for i in inds)) else: flux = [] pos = [] @@ -2366,14 +2366,14 @@ def psf_photutils(self,xPix=None,yPix=None,size=5,local_bkg=False,epsf=None, eflux = np.zeros(len(self.flux)) * np.nan psfphot2 = PSFPhotometry(epsf, fit_shape, finder=None,aperture_radius=1.5, xy_bounds=(0.05),localbkg_estimator=localbkg_estimator) - f,ef = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot2,init) for i in np.arange(len(cutouts)))) + f,ef = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot2,init) for i in np.arange(len(cutouts)))) f = np.array(f).flatten() ef = np.array(ef).flatten() phot = phot.to_pandas() pos = phot[['x_fit','y_fit']].values + np.array([xPix,yPix]) - size//2 epos = phot[['x_err','y_err']].values else: - f,ef,pos,epos = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot,init,True) for i in np.arange(len(cutouts)))) + f,ef,pos,epos = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot,init,True) for i in np.arange(len(cutouts)))) pos['x_fit'] += xPix - size//2 pos['y_fit'] += xPix - size//2 @@ -2437,7 +2437,7 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa prf, cutouts, ecutouts = self._psf_initialise(size,(xPix,yPix),ref=(not diff)) # gather base PRF and the array of cutouts data inds = np.arange(len(cutouts)) base = create_psf(prf,size) - flux, eflux, pos = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_full)(cutouts[i],base,self.shift[i]) for i in inds)) + flux, eflux, pos = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_full)(cutouts[i],base,self.shift[i]) for i in inds)) #prf, cutouts = self._psf_initialise(size,(xPix,yPix)) # gather base PRF and the array of cutouts data #xShifts = [] @@ -2483,9 +2483,9 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa if self.parallel: inds = np.arange(len(cutouts)) if self.delta_kernel is not None: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) else: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) else: for i in range(len(cutouts)): flux += [par_psf_flux(cutouts[i],ecutouts[i],base,self.shift[i])] @@ -2501,9 +2501,9 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa if self.parallel: inds = np.arange(len(cutouts)) if self.delta_kernel is not None: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) else: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) else: for i in range(len(cutouts)): flux += [par_psf_flux(cutouts[i],ecutouts[i],base,self.shift[i])] @@ -2548,7 +2548,7 @@ def kernel_matching(self,size=7,diff=True): mask = self.mask == 1 if self.parallel: - d, kernel = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing")(delayed(parallel_delta_diff)(frame,self.ref,mask,size) for frame in flux)) + d, kernel = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(parallel_delta_diff)(frame,self.ref,mask,size) for frame in flux)) else: d = [] kernel = [] From 0c18012c3e5d8ea61cf39a22c11e85d1ae1204fb Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 4 Jun 2026 12:39:24 +1200 Subject: [PATCH 10/22] Add core-detection and SLURM env diagnostics to _available_cores and __init__ --- tessreduce/tessreduce.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 9db216f..cce63e7 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -33,14 +33,20 @@ def _available_cores(): slurm = os.environ.get('SLURM_CPUS_PER_TASK') if slurm is not None: try: - return int(slurm) + n = int(slurm) + print(f'[tessreduce] _available_cores: SLURM_CPUS_PER_TASK={slurm} → {n} cores') + return n except ValueError: pass try: - return len(os.sched_getaffinity(0)) + n = len(os.sched_getaffinity(0)) + print(f'[tessreduce] _available_cores: sched_getaffinity → {n} cores') + return n except AttributeError: pass - return multiprocessing.cpu_count() + n = multiprocessing.cpu_count() + print(f'[tessreduce] _available_cores: cpu_count fallback → {n} cores') + return n from .catalog_tools import * from .calibration_tools import * @@ -245,6 +251,13 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect self._timing = timing self._cache_path = None + # SLURM environment diagnostics + _slurm_vars = ['SLURM_CPUS_PER_TASK', 'SLURM_NTASKS', 'SLURM_NTASKS_PER_NODE', + 'SLURM_JOB_CPUS_PER_NODE', 'SLURM_CPUS_ON_NODE'] + _slurm_env = {k: os.environ.get(k, 'not set') for k in _slurm_vars} + print(f'[tessreduce] SLURM env: {_slurm_env}') + print(f'[tessreduce] num_cores resolved to: {self.num_cores}') + # Offline Paths if catalogue_path is None: catalogue_path = os.getcwd() From 56cf828ed6674deff3346354bdd2ee96936e1ad8 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 4 Jun 2026 12:55:11 +1200 Subject: [PATCH 11/22] Convert closures to module-level functions to fix multiprocessing pickling backend='multiprocessing' uses standard pickle to dispatch jobs even with fork, which cannot serialize closures. The affected functions were silently falling back to sequential execution. helpers.py: - Extract _process closure from fix_background_anomalies to module-level _fix_bkg_frame, passing per-frame 2D slices instead of the full 3D cube to avoid large-array pickling - Extract _fit_residual closure to module-level _fit_residual_bkg with explicit res_box and n_sigma parameters - Both Parallel calls now use backend='multiprocessing' explicitly tessreduce.py: - Extract _fit_frame closure from _fit_residual_surface to module-level _fit_bkg_surface_frame with explicit parameters --- tessreduce/helpers.py | 423 +++++++++++++++++++-------------------- tessreduce/tessreduce.py | 40 ++-- 2 files changed, 228 insertions(+), 235 deletions(-) diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index 60cc21d..c66397b 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -1405,6 +1405,193 @@ def parallel_photutils(cutout,e_cutout,psf_phot,init_params=None,return_pos=Fals return np.array([np.nan]), np.array([np.nan]) +def _fix_bkg_frame(bkg_i, flux_i, bkgmask_i, prev_i, high_bkg_i, + strap, strap_cols, has_straps, src_mask_2d, + eff_box, disk, yy, xx, NY, NX, + n_sigma, anom_box, anom_box_fine, sep_thresh, sep_snr_thresh, + sep_validate, gauss_smooth, bad_bkg_sigma, bad_bkg_min_area): + from photutils.background import Background2D, MedianBackground + data_src = np.isnan(bkgmask_i) if bkgmask_i is not None else src_mask_2d + phot_mask = strap | data_src + frame = bkg_i.copy() + excess = np.zeros(len(strap_cols)) + if has_straps and flux_i is not None: + _, excess, _ = sigma_clipped_stats((flux_i - frame)[:, strap_cols], axis=0) + frame[:, strap_cols] -= excess + try: + bkg2d = Background2D(frame, box_size=eff_box, filter_size=3, + mask=phot_mask, bkg_estimator=MedianBackground(), + exclude_percentile=50) + trend = bkg2d.background + except Exception: + trend = np.full_like(frame, np.nanmedian(frame)) + resid = frame - trend + valid_mask = ~phot_mask + + def _block_sigma(r, box): + sigma = np.full((NY, NX), np.inf, dtype=float) + row_starts = [min(r0, NY - box) for r0 in range(0, NY, box)] + col_starts = [min(c0, NX - box) for c0 in range(0, NX, box)] + for r0 in row_starts: + for c0 in col_starts: + r1, c1 = r0 + box, c0 + box + vals = r[r0:r1, c0:c1][valid_mask[r0:r1, c0:c1]] + if vals.size >= 4: + m = np.nanmedian(vals) + sigma[r0:r1, c0:c1] = 1.4826 * np.nanmedian(np.abs(vals - m)) + return sigma + + is_high_bkg = bool(high_bkg_i) + sharp_mask = np.zeros((NY, NX), dtype=bool) + if not is_high_bkg: + sigma_coarse = _block_sigma(resid, anom_box) + flagged_coarse = (np.abs(resid) > n_sigma * sigma_coarse) & ~phot_mask + import sep as _sep + lap_abs = np.abs(laplace(frame)).astype(np.float64) + _bkg = _sep.Background(lap_abs) + lap_sub = lap_abs - _bkg.back() + lap_err = _bkg.rms() + snr_map = lap_sub / (lap_err + 1e-10) + try: + objects = _sep.extract(lap_sub, thresh=sep_thresh, err=lap_err) + except Exception: + objects = [] + sep_mask = np.zeros((NY, NX), dtype=bool) + noise = np.nanmedian(lap_err) + for obj in objects: + ap_mask = np.zeros((NY, NX), dtype=bool) + _sep.mask_ellipse(ap_mask, obj['x'], obj['y'], obj['a'], obj['b'], obj['theta'], r=3.0) + if ap_mask.sum() == 0 or snr_map[ap_mask].mean() <= sep_snr_thresh: + continue + cx, cy = obj['x'], obj['y'] + dist = np.sqrt((xx - cx)**2 + (yy - cy)**2) + true_r = None + for r in range(2, 20): + ann = (dist >= r - 0.5) & (dist < r + 0.5) + if ann.sum() == 0: + break + if lap_sub[ann].mean() < noise: + true_r = r - 1 + break + else: + true_r = 19 + if true_r is None or true_r < 2 or true_r > 5: + continue + sep_mask |= dist <= true_r + lap_med = np.nanmedian(lap_abs) + lap_mad = np.nanmedian(np.abs(lap_abs - lap_med)) + is_sharp = lap_abs > lap_med + 3 * 1.4826 * lap_mad + sharp_mask |= sep_mask + if sep_validate: + edge_border = np.zeros((NY, NX), dtype=bool) + edge_border[0, :] = True; edge_border[-1, :] = True + edge_border[:, 0] = True; edge_border[:, -1] = True + labeled, n_comp = label(flagged_coarse) + if n_comp > 0: + lab_flat = labeled.ravel() + touches_edge = np.zeros(n_comp + 1, dtype=bool) + touches_sep_lbl = np.zeros(n_comp + 1, dtype=bool) + touches_sharp_lbl = np.zeros(n_comp + 1, dtype=bool) + np.bitwise_or.at(touches_edge, lab_flat, edge_border.ravel()) + np.bitwise_or.at(touches_sep_lbl, lab_flat, sep_mask.ravel()) + np.bitwise_or.at(touches_sharp_lbl, lab_flat, is_sharp.ravel()) + is_sharp_lbl = ~(touches_edge & ~touches_sep_lbl) & touches_sharp_lbl & touches_sep_lbl + is_sharp_lbl[0] = False + sharp_mask |= is_sharp_lbl[labeled] + else: + sharp_mask |= flagged_coarse & is_sharp + smooth_mask = flagged_coarse & ~sharp_mask + any_anom = sharp_mask.any() or smooth_mask.any() + if any_anom: + eff_fine = max(min(anom_box_fine, min(NY, NX) // 2), 4) + try: + bkg2d_fine = Background2D(frame, box_size=eff_fine, filter_size=3, + mask=phot_mask | sharp_mask, bkg_estimator=MedianBackground(), + exclude_percentile=50) + trend_fine = bkg2d_fine.background + except Exception: + trend_fine = trend + if sharp_mask.any(): + frame[sharp_mask & ~phot_mask] = trend_fine[sharp_mask & ~phot_mask] + resid_fine = frame - trend_fine + sigma_fine = _block_sigma(resid_fine, anom_box_fine) + flagged_fine = (np.abs(resid_fine) > n_sigma * sigma_fine) & ~phot_mask + confirmed = smooth_mask & flagged_fine + if confirmed.any(): + dilated = binary_dilation(confirmed, structure=disk) & ~phot_mask + frame[dilated] = trend_fine[dilated] + if prev_i is not None and flux_i is not None: + resid_new = np.abs(flux_i - frame) + resid_prev = np.abs(flux_i - prev_i) + delta = resid_new - resid_prev + _, _, scale = sigma_clipped_stats(resid_prev) + w = np.clip(delta / (scale + 1e-10), 0, 1) + w_zero = (w == 0) + if w_zero.any(): + w_med = median_filter(w_zero.astype(float), size=7) + diff_mask = w_zero & ~(w_med > 0.5) + w[w_med > 0.5] = 0 + if diff_mask.any(): + w[diff_mask] = np.nan + _, idx = distance_transform_edt(np.isnan(w), return_indices=True) + w[diff_mask] = w[tuple(idx[:, diff_mask])] + w[sharp_mask] = 0.0 + w[phot_mask] = 0.0 + frame = (1 - w) * frame + w * prev_i + smooth_sigma = 2.0 if is_high_bkg else gauss_smooth + fixed = gaussian_filter(frame, sigma=smooth_sigma) + if has_straps and flux_i is not None: + _, excess, _ = sigma_clipped_stats((flux_i - fixed)[:, strap_cols], axis=0) + fixed[:, strap_cols] += excess + lap_abs = np.abs(laplace(fixed)) + lap_med = np.nanmedian(lap_abs) + lap_mad = np.nanmedian(np.abs(lap_abs - lap_med)) + high_lap = lap_abs > lap_med + bad_bkg_sigma * 1.4826 * lap_mad + labeled_lap, n_lap = label(high_lap) + if n_lap > 0: + areas = np.bincount(labeled_lap.ravel(), minlength=n_lap + 1) + large = np.flatnonzero(areas[1:] >= bad_bkg_min_area) + 1 + bad_bkg_mask = np.isin(labeled_lap, large) if len(large) > 0 else np.zeros((NY, NX), dtype=bool) + else: + bad_bkg_mask = np.zeros((NY, NX), dtype=bool) + return fixed, excess, sharp_mask, bad_bkg_mask + + +def _fit_residual_bkg(residual, exclude_mask, res_box, n_sigma): + from photutils.background import Background2D, MedianBackground + from astropy.stats import SigmaClip + sc = SigmaClip(sigma=3.0, maxiters=5) + finite_vals = residual[~exclude_mask & np.isfinite(residual)] + med = np.nanmedian(finite_vals) + std = np.nanstd(finite_vals) + transient_mask = exclude_mask | (np.abs(residual - med) > 5 * std) + try: + b = Background2D(residual, box_size=res_box, filter_size=3, + sigma_clip=sc, bkg_estimator=MedianBackground(), + mask=transient_mask, fill_value=0.0) + corr = b.background + except Exception: + corr = np.full_like(residual, np.nanmedian(residual[~transient_mask])) + corr_resid = residual - corr + _, corr_med, corr_std = sigma_clipped_stats(corr_resid[~exclude_mask]) + flagged = np.abs(corr_resid) > n_sigma * corr_std + if flagged.any(): + lap_abs = np.abs(laplace(corr_resid)) + lap_med = np.nanmedian(lap_abs) + lap_mad = np.nanmedian(np.abs(lap_abs - lap_med)) + is_sharp = lap_abs > lap_med + 3 * 1.4826 * lap_mad + labeled_c, n_c = label(flagged) + if n_c > 0: + lab_flat = labeled_c.ravel() + comp_sizes = np.bincount(lab_flat, minlength=n_c + 1) + sharp_counts = np.bincount(lab_flat, weights=is_sharp.ravel(), minlength=n_c + 1) + sharp_frac = sharp_counts / np.maximum(comp_sizes, 1) + suppress = np.flatnonzero(sharp_frac[1:] >= 0.3) + 1 + if len(suppress) > 0: + corr[np.isin(labeled_c, suppress)] = 0.0 + return corr + + def fix_background_anomalies(bkg, mask, flux=None, bkg_prev=None, bkgmask=None, n_sigma=5.0, box_size=16, anom_box=30, anom_box_fine=4, dilate_r=2, gauss_smooth=2, sep_thresh=3.0, sep_snr_thresh=2.0, sep_validate=True, @@ -1468,233 +1655,41 @@ def fix_background_anomalies(bkg, mask, flux=None, bkg_prev=None, bkgmask=None, yr, xr = np.ogrid[-dilate_r:dilate_r+1, -dilate_r:dilate_r+1] disk = xr**2 + yr**2 <= dilate_r**2 yy, xx = np.mgrid[:NY, :NX] - def _process(i): - if bkgmask_arr is not None: - data_src = np.isnan(bkgmask_arr[i]) if bkgmask_arr.ndim == 3 else np.isnan(bkgmask_arr) - else: - data_src = src_mask_2d - phot_mask = strap | data_src - - frame = bkg[i].copy() - excess = np.zeros(len(strap_cols)) - if has_straps and flux is not None: - _, excess, _ = sigma_clipped_stats((flux[i] - frame)[:, strap_cols], axis=0) - frame[:, strap_cols] -= excess + def _bkgmask_i(i): + if bkgmask_arr is None: + return None + return bkgmask_arr[i] if bkgmask_arr.ndim == 3 else bkgmask_arr - try: - bkg2d = Background2D(frame, box_size=eff_box, filter_size=3, - mask=phot_mask, bkg_estimator=MedianBackground(), - exclude_percentile=50) - trend = bkg2d.background - except Exception: - trend = np.full_like(frame, np.nanmedian(frame)) - - resid = frame - trend - - valid_mask = ~phot_mask - - def _block_sigma(r, box): - sigma = np.full((NY, NX), np.inf, dtype=float) - row_starts = [min(r0, NY - box) for r0 in range(0, NY, box)] - col_starts = [min(c0, NX - box) for c0 in range(0, NX, box)] - for r0 in row_starts: - for c0 in col_starts: - r1, c1 = r0 + box, c0 + box - vals = r[r0:r1, c0:c1][valid_mask[r0:r1, c0:c1]] - if vals.size >= 4: - m = np.nanmedian(vals) - sigma[r0:r1, c0:c1] = 1.4826 * np.nanmedian(np.abs(vals - m)) - return sigma - - is_high_bkg = (high_bkg_frames is not None) and bool(high_bkg_frames[i]) - - sharp_mask = np.zeros((NY, NX), dtype=bool) - - if not is_high_bkg: - sigma_coarse = _block_sigma(resid, anom_box) - flagged_coarse = (np.abs(resid) > n_sigma * sigma_coarse) & ~phot_mask - - # Build SEP source mask on |Laplacian(frame)|. - import sep as _sep - lap_abs = np.abs(laplace(frame)).astype(np.float64) - _bkg = _sep.Background(lap_abs) - lap_sub = lap_abs - _bkg.back() - lap_err = _bkg.rms() - snr_map = lap_sub / (lap_err + 1e-10) - try: - objects = _sep.extract(lap_sub, thresh=sep_thresh, err=lap_err) - except Exception: - objects = [] - sep_mask = np.zeros((NY, NX), dtype=bool) - noise = np.nanmedian(lap_err) - for obj in objects: - ap_mask = np.zeros((NY, NX), dtype=bool) - _sep.mask_ellipse(ap_mask, obj['x'], obj['y'], obj['a'], obj['b'], obj['theta'], r=3.0) - if ap_mask.sum() == 0 or snr_map[ap_mask].mean() <= sep_snr_thresh: - continue - cx, cy = obj['x'], obj['y'] - dist = np.sqrt((xx - cx)**2 + (yy - cy)**2) - true_r = None - for r in range(2, 20): - ann = (dist >= r - 0.5) & (dist < r + 0.5) - if ann.sum() == 0: - break - if lap_sub[ann].mean() < noise: - true_r = r - 1 - break - else: - true_r = 19 - if true_r is None or true_r < 2 or true_r > 5: - continue - sep_mask |= dist <= true_r - - lap_med = np.nanmedian(lap_abs) - lap_mad = np.nanmedian(np.abs(lap_abs - lap_med)) - is_sharp = lap_abs > lap_med + 3 * 1.4826 * lap_mad - - sharp_mask |= sep_mask - - if sep_validate: - edge_border = np.zeros((NY, NX), dtype=bool) - edge_border[0, :] = True; edge_border[-1, :] = True - edge_border[:, 0] = True; edge_border[:, -1] = True - labeled, n_comp = label(flagged_coarse) - if n_comp > 0: - lab_flat = labeled.ravel() - touches_edge = np.zeros(n_comp + 1, dtype=bool) - touches_sep_lbl = np.zeros(n_comp + 1, dtype=bool) - touches_sharp_lbl = np.zeros(n_comp + 1, dtype=bool) - np.bitwise_or.at(touches_edge, lab_flat, edge_border.ravel()) - np.bitwise_or.at(touches_sep_lbl, lab_flat, sep_mask.ravel()) - np.bitwise_or.at(touches_sharp_lbl, lab_flat, is_sharp.ravel()) - # sharp: first condition False AND touches_sharp AND touches_sep - is_sharp_lbl = ~(touches_edge & ~touches_sep_lbl) & touches_sharp_lbl & touches_sep_lbl - is_sharp_lbl[0] = False - sharp_mask |= is_sharp_lbl[labeled] - else: - sharp_mask |= flagged_coarse & is_sharp - - smooth_mask = flagged_coarse & ~sharp_mask - - # Fit fine Background2D excluding sharp anomaly pixels. - any_anom = sharp_mask.any() or smooth_mask.any() - if any_anom: - eff_fine = max(min(anom_box_fine, min(NY, NX) // 2), 4) - try: - bkg2d_fine = Background2D(frame, box_size=eff_fine, filter_size=3, - mask=phot_mask | sharp_mask, bkg_estimator=MedianBackground(), - exclude_percentile=50) - trend_fine = bkg2d_fine.background - except Exception: - trend_fine = trend - - if sharp_mask.any(): - frame[sharp_mask & ~phot_mask] = trend_fine[sharp_mask & ~phot_mask] - - resid_fine = frame - trend_fine - sigma_fine = _block_sigma(resid_fine, anom_box_fine) - flagged_fine = (np.abs(resid_fine) > n_sigma * sigma_fine) & ~phot_mask - confirmed = smooth_mask & flagged_fine - if confirmed.any(): - dilated = binary_dilation(confirmed, structure=disk) & ~phot_mask - frame[dilated] = trend_fine[dilated] - - # Blend toward bkg_prev before smoothing so smoothing is not overridden. - if bkg_prev is not None and flux is not None: - flux_frame = flux[i] - prev_frame = np.asarray(bkg_prev)[i] - resid_new = np.abs(flux_frame - frame) - resid_prev = np.abs(flux_frame - prev_frame) - delta = resid_new - resid_prev - _, _, scale = sigma_clipped_stats(resid_prev) - w = np.clip(delta / (scale + 1e-10), 0, 1) - w_zero = (w == 0) - if w_zero.any(): - w_med = median_filter(w_zero.astype(float), size=7) - diff_mask = w_zero & ~(w_med > 0.5) - w[w_med > 0.5] = 0 - if diff_mask.any(): - w[diff_mask] = np.nan - _, idx = distance_transform_edt(np.isnan(w), return_indices=True) - w[diff_mask] = w[tuple(idx[:, diff_mask])] - w[sharp_mask] = 0.0 - w[phot_mask] = 0.0 - frame = (1 - w) * frame + w * prev_frame - - smooth_sigma = 2.0 if is_high_bkg else gauss_smooth - fixed = gaussian_filter(frame, sigma=smooth_sigma) - - if has_straps and flux is not None: - _, excess, _ = sigma_clipped_stats((flux[i] - fixed)[:, strap_cols], axis=0) - fixed[:, strap_cols] += excess - - lap_abs = np.abs(laplace(fixed)) - lap_med = np.nanmedian(lap_abs) - lap_mad = np.nanmedian(np.abs(lap_abs - lap_med)) - high_lap = lap_abs > lap_med + bad_bkg_sigma * 1.4826 * lap_mad - labeled_lap, n_lap = label(high_lap) - if n_lap > 0: - areas = np.bincount(labeled_lap.ravel(), minlength=n_lap + 1) - large = np.flatnonzero(areas[1:] >= bad_bkg_min_area) + 1 - bad_bkg_mask = np.isin(labeled_lap, large) if len(large) > 0 else np.zeros((NY, NX), dtype=bool) - else: - bad_bkg_mask = np.zeros((NY, NX), dtype=bool) - - return fixed, excess, sharp_mask, bad_bkg_mask - - results = Parallel(n_jobs=n_jobs, prefer='threads')(delayed(_process)(i) for i in range(T)) + results = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( + delayed(_fix_bkg_frame)( + bkg[i], + flux[i] if flux is not None else None, + _bkgmask_i(i), + np.asarray(bkg_prev)[i] if bkg_prev is not None else None, + bool(high_bkg_frames[i]) if high_bkg_frames is not None else False, + strap, strap_cols, has_straps, src_mask_2d, + eff_box, disk, yy, xx, NY, NX, + n_sigma, anom_box, anom_box_fine, sep_thresh, sep_snr_thresh, + sep_validate, gauss_smooth, bad_bkg_sigma, bad_bkg_min_area, + ) + for i in range(T) + ) bkg_fixed = np.array([r[0] for r in results]) excesses = [r[1] for r in results] sharp_masks = np.array([r[2] for r in results]) bad_bkg_masks = np.array([r[3] for r in results]) # (T, NY, NX) bool if flux is not None and bkgmask is not None: - from astropy.stats import SigmaClip - sc = SigmaClip(sigma=3.0, maxiters=5) res_box = min(20, min(NY, NX) // 2) res_box = max(res_box, 4) - def _fit_residual(residual, exclude_mask): - finite_vals = residual[~exclude_mask & np.isfinite(residual)] - med = np.nanmedian(finite_vals) - std = np.nanstd(finite_vals) - transient_mask = exclude_mask | (np.abs(residual - med) > 5 * std) - try: - b = Background2D(residual, box_size=res_box, filter_size=3, - sigma_clip=sc, bkg_estimator=MedianBackground(), - mask=transient_mask, fill_value=0.0) - corr = b.background - except Exception: - corr = np.full_like(residual, np.nanmedian(residual[~transient_mask])) - - # Apply Laplacian sharp/smooth classification to the correction. - # Only suppress corrections in sharp (point-like) regions — smooth - # background gradients should be applied. - corr_resid = residual - corr - _, corr_med, corr_std = sigma_clipped_stats(corr_resid[~exclude_mask]) - flagged = np.abs(corr_resid) > n_sigma * corr_std - if flagged.any(): - lap_abs = np.abs(laplace(corr_resid)) - lap_med = np.nanmedian(lap_abs) - lap_mad = np.nanmedian(np.abs(lap_abs - lap_med)) - is_sharp = lap_abs > lap_med + 3 * 1.4826 * lap_mad - labeled_c, n_c = label(flagged) - if n_c > 0: - lab_flat = labeled_c.ravel() - comp_sizes = np.bincount(lab_flat, minlength=n_c + 1) - sharp_counts = np.bincount(lab_flat, weights=is_sharp.ravel(), minlength=n_c + 1) - sharp_frac = sharp_counts / np.maximum(comp_sizes, 1) - suppress = np.flatnonzero(sharp_frac[1:] >= 0.3) + 1 - if len(suppress) > 0: - corr[np.isin(labeled_c, suppress)] = 0.0 - return corr - residuals = flux - bkg_fixed - corrections = Parallel(n_jobs=n_jobs, prefer='threads')( - delayed(_fit_residual)( + corrections = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( + delayed(_fit_residual_bkg)( residuals[i], - np.isnan(bkgmask_arr[i]) if bkgmask_arr.ndim == 3 else np.isnan(bkgmask_arr) + np.isnan(bkgmask_arr[i]) if bkgmask_arr.ndim == 3 else np.isnan(bkgmask_arr), + res_box, n_sigma, ) for i in range(T) ) bkg_fixed += np.array(corrections) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index cce63e7..1c66b5a 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -68,6 +68,22 @@ def _available_cores(): # set the package directory so we can load in a file later package_directory = os.path.dirname(os.path.abspath(__file__)) + '/' +def _fit_bkg_surface_frame(residual, exclude_mask, box_size, filter_size, sigma): + from photutils.background import Background2D, MedianBackground + from astropy.stats import SigmaClip + sc = SigmaClip(sigma=sigma, maxiters=5) + finite_vals = residual[~exclude_mask & np.isfinite(residual)] + med = np.nanmedian(finite_vals) + std = np.nanstd(finite_vals) + transient_mask = exclude_mask | (residual > med + 5 * std) + try: + b = Background2D(residual, box_size=box_size, filter_size=filter_size, + sigma_clip=sc, bkg_estimator=MedianBackground(), + mask=transient_mask, fill_value=0.0) + return b.background + except Exception: + return np.full_like(residual, np.nanmedian(residual[~transient_mask])) + fig_width_pt = 240.0 # Get this from LaTeX using \showthe\columnwidth inches_per_pt = 1.0/72.27 # Convert pt to inches golden_mean = (np.sqrt(5)-1.0)/2.0 # Aesthetic ratio @@ -105,35 +121,17 @@ def _subtract_residual_surface(bkg, flux, bkgmask, box_size=20, filter_size=5, s from astropy.stats import SigmaClip from joblib import Parallel, delayed - sc = SigmaClip(sigma=sigma, maxiters=5) - estimator = MedianBackground() - bkgmask = np.asarray(bkgmask) if bkgmask.ndim == 3: exclude_mask = np.any(np.isnan(bkgmask), axis=0) else: exclude_mask = np.isnan(bkgmask) - def _fit_frame(residual): - finite_vals = residual[~exclude_mask & np.isfinite(residual)] - med = np.nanmedian(finite_vals) - std = np.nanstd(finite_vals) - transient_mask = exclude_mask | (residual > med + 5 * std) - try: - b = Background2D(residual, - box_size=box_size, - filter_size=filter_size, - sigma_clip=sc, - bkg_estimator=estimator, - mask=transient_mask, - fill_value=0.0) - return b.background - except Exception: - return np.full_like(residual, np.nanmedian(residual[~transient_mask])) - n_jobs = _available_cores() residuals = flux - bkg - corrections = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)(delayed(_fit_frame)(residuals[i]) for i in range(flux.shape[0])) + corrections = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( + delayed(_fit_bkg_surface_frame)(residuals[i], exclude_mask, box_size, filter_size, sigma) + for i in range(flux.shape[0])) bkg += np.array(corrections) return bkg From 8ac1904c685998277622d3e2dfcb59b4e1450bc5 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 4 Jun 2026 13:02:51 +1200 Subject: [PATCH 12/22] Add progress prints at each key pipeline stage --- tessreduce/tessreduce.py | 29 ++++++++++++++++++----------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 1c66b5a..b994639 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -789,8 +789,10 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa bkg_smth = np.zeros_like(flux) * np.nan if self.parallel: _t = time.perf_counter() + print('smooth background...') bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*m) _times['initial smooth background'] = time.perf_counter() - _t + print('smooth background done') if rerun_negative: _t = time.perf_counter() if self._use_error_image: @@ -811,8 +813,10 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa else: m[over_sub] = 1 self._bkgmask = m + print('smooth background rerun (negative correction)...') bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(Smooth_bkg)(frame,gauss_smooth,interpolate) for frame in flux*m) _times['negative over-subtraction rerun'] = time.perf_counter() - _t + print('smooth background rerun done') if rerun_diff: @@ -864,6 +868,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa new_mask = abs(new_mask - 1) self._bkgmask = new_mask bkg_s1 = np.array(bkg_smth) + print('smooth background rerun (residual surface)...') bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*new_mask) if blend_dynamic: bkg_smth = blend_dynamic_background(bkg_smth, bkg_s1, flux, n_jobs=self.num_cores) @@ -896,6 +901,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa _high_bkg_frames = (_earth_angle < 30.0) | (_moon_angle < 30.0) | (_bkg_median > 300.0) else: _high_bkg_frames = _bkg_median > 300.0 + print('fixing background anomalies...') self.bkg, _sharp_masks, self.bad_bkg = fix_background_anomalies(self.bkg, self.mask, flux=deepcopy(self.flux), bkg_prev=bkg_pre_fix if blend_dynamic else None, @@ -904,8 +910,10 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa high_bkg_frames=_high_bkg_frames, n_jobs=self.num_cores) _times['anomaly fixing'] = time.perf_counter() - _t + print('fixing background anomalies done') _t = time.perf_counter() + print('adaptive temporal smoothing...') from .adaptive_background import AdaptiveBackground smoother = AdaptiveBackground(self.bkg, self.mjd, sector=self.sector, camera=self.tpf.camera, data_path=self._vector_path,n_jobs=self.num_cores) @@ -913,6 +921,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa smoothed = smoother.smooth(method='savgol').smoothed self.bkg = smoothed _times['adaptive temporal smoothing'] = time.perf_counter() - _t + print('adaptive temporal smoothing done') # Store data-driven sources from _bkgmask as bit 8, preserving the catalogue mask (bit 1) if rerun_diff: @@ -2680,14 +2689,14 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, # calculate the background #self.flux -= self.ref + print('background pass 1...') _t = time.perf_counter() self.background(rerun_negative=True) - #self.flux += self.ref self.flux -= self.bkg _times['background (pass 1)'] = time.perf_counter() - _t + print('background pass 1 done') if np.isnan(self.bkg).all(): - # check to see if the background worked raise ValueError('bkg all nans') # flux = strip_units(self.flux) @@ -2740,7 +2749,8 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, else: self.shift = np.zeros((len(self.flux),2)) _times['alignment'] = time.perf_counter() - _t - + print('alignment done') + if not self.diff: if self.align: _t = time.perf_counter() @@ -2800,27 +2810,24 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, temp[:,:,:] = self.mask self.mask = temp | moving_mask _times['difference imaging setup'] = time.perf_counter() - _t + print('diff setup done') - if self.verbose > 0: - print('remade mask') - # background - if self.verbose > 0: - print('background') self.bkg_orig = deepcopy(self.bkg) + print('background pass 2...') _t = time.perf_counter() self.background(calc_qe = False,strap_iso = False,source_hunt=self._sourcehunt, gauss_smooth=self._bkg_gauss_sigma,interpolate=False, rerun_negative=False,rerun_diff=True,blend_dynamic=True) - # self._grad_bkg_clip() self.flux -= self.bkg _times['background (pass 2)'] = time.perf_counter() - _t + print('background pass 2 done') if self.corr_correction: _t = time.perf_counter() - if self.verbose > 0: - print('background correlation correction') + print('correlation correction...') self.correlation_corrector() _times['correlation correction'] = time.perf_counter() - _t + print('correlation correction done') if self.kernel_match: self.kernel_matching(diff=self.diff) if self.verbose > 0: From 790be067fd9fda5e004d7af54f957ba58fb62160 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 4 Jun 2026 13:14:35 +1200 Subject: [PATCH 13/22] Fix duplicate verbose keyword in sep_aligner Parallel call --- tessreduce/sep_aligner.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tessreduce/sep_aligner.py b/tessreduce/sep_aligner.py index 4f27ac3..0f14abe 100644 --- a/tessreduce/sep_aligner.py +++ b/tessreduce/sep_aligner.py @@ -703,7 +703,7 @@ def run(self, time: Optional[np.ndarray] = None, self._sub_ref, cores, positions, np.asarray(weights), p['core_half']) results = Parallel(n_jobs=self.n_jobs, verbose=verbose, - backend="multiprocessing", verbose=1)( + backend="multiprocessing")( delayed(_align_one_frame)( t, self.flux[t], ref_comp, w_cols, From ae6185869764f79fb1235c1a1f297fc0d74da497 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 4 Jun 2026 13:30:30 +1200 Subject: [PATCH 14/22] Fix parallel_background2d: create SigmaClip/MedianBackground inside worker to avoid pickle failure --- tessreduce/helpers.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index c66397b..4c81b03 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -182,23 +182,21 @@ def parallel_bkg3(data,mask): estimate = inpaint.inpaint_biharmonic(data,mask) return estimate -def _background2d_frame(frame, box_size, filter_size, sc, estimator, mask=None): - from photutils.background import Background2D +def _background2d_frame(frame, box_size, filter_size, sigma, maxiters, mask=None): + from photutils.background import Background2D, MedianBackground + from astropy.stats import SigmaClip + sc = SigmaClip(sigma=sigma, maxiters=maxiters) try: return Background2D(frame, box_size=box_size, filter_size=filter_size, - sigma_clip=sc, bkg_estimator=estimator, + sigma_clip=sc, bkg_estimator=MedianBackground(), mask=mask, fill_value=0.0, exclude_percentile=50).background except Exception: return np.full_like(frame, np.nanmedian(frame)) def parallel_background2d(cube, box_size=5, filter_size=3, sigma=3, maxiters=5, n_jobs=-1, mask=None): - from photutils.background import MedianBackground - from astropy.stats import SigmaClip - sc = SigmaClip(sigma=sigma, maxiters=maxiters) - estimator = MedianBackground() return np.array(Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( - delayed(_background2d_frame)(frame, box_size, filter_size, sc, estimator, mask) + delayed(_background2d_frame)(frame, box_size, filter_size, sigma, maxiters, mask) for frame in cube)) def Smooth_bkg(data, gauss_smooth=0, interpolate=False, extrapolate=True): From 615f0a02f392b2c4b22538ddb5b2f9407e804261 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 4 Jun 2026 20:45:46 +1200 Subject: [PATCH 15/22] Fix blend_dynamic_background: create Gaussian2DKernel inside worker to fix pickle failure Gaussian2DKernel was passed as an argument to each joblib worker. Like SigmaClip and MedianBackground, it contains numpy dtype dispatch wrappers that cannot be pickled by standard pickle, causing silent fallback to serial execution for all 3372 frames (~50s overhead in residual surface rerun). Creating it inside _blend_frame avoids any pickling of the object. --- tessreduce/helpers.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index 4c81b03..a5e5d09 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -1789,7 +1789,7 @@ def orbit_ref_subtract(flux, times_mjd, sector=None, camera=None, return result, segments, orbit_refs -def _blend_frame(bkg_new_i, bkg_prev_i, delta_i, resid_prev_i, sigma, sharp_mask_i, gauss_kernel): +def _blend_frame(bkg_new_i, bkg_prev_i, delta_i, resid_prev_i, sigma, sharp_mask_i): _, _, scale = sigma_clipped_stats(resid_prev_i) w = np.clip(delta_i / (sigma * scale + 1e-10), 0, 1) w_zero = (w == 0) @@ -1799,7 +1799,7 @@ def _blend_frame(bkg_new_i, bkg_prev_i, delta_i, resid_prev_i, sigma, sharp_mask w[w_med > 0.5] = 0 if diff_mask.any(): w[diff_mask] = np.nan - w = interpolate_replace_nans(w, gauss_kernel) + w = interpolate_replace_nans(w, Gaussian2DKernel(1.0)) if sharp_mask_i is not None: w[sharp_mask_i] = 0.0 return (1 - w) * bkg_new_i + w * bkg_prev_i @@ -1832,14 +1832,12 @@ def blend_dynamic_background(bkg_new, bkg_prev, flux, sigma=2.0, sharp_masks=Non resid_new = np.abs(flux - bkg_new) resid_prev = np.abs(flux - bkg_prev) delta = resid_new - resid_prev - gauss_kernel = Gaussian2DKernel(1.0) T = bkg_new.shape[0] - sharp_list = [sharp_masks[i] if sharp_masks is not None else None for i in range(T)] results = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( delayed(_blend_frame)(bkg_new[i], bkg_prev[i], delta[i], resid_prev[i], - sigma, sharp_list[i], gauss_kernel) + sigma, sharp_list[i]) for i in range(T)) return np.array(results) From a06ae75c0dadcfaad2798d9e3915f69fc9470a8c Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Tue, 9 Jun 2026 21:10:09 +1200 Subject: [PATCH 16/22] Fix Background2D crash on heavily masked fields with exclude_percentile=50 fallback --- tessreduce/tessreduce.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index b994639..d337dae 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -836,13 +836,17 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa estimator = MedianBackground() sc = SigmaClip(sigma=5, maxiters=5) - b = Background2D(std, - box_size=5, - filter_size=3, - sigma_clip=sc, - bkg_estimator=estimator, - fill_value=0.0) - std_sub = std - b.background + try: + b = Background2D(std, + box_size=5, + filter_size=3, + sigma_clip=sc, + bkg_estimator=estimator, + exclude_percentile=50, + fill_value=0.0) + std_sub = std - b.background + except ValueError: + std_sub = std _,smed,sstd = sigma_clipped_stats(std_sub) resid_mask = (std_sub > smed + 3*sstd) * 1.0 ny, nx = resid_mask.shape From e335285e95778d3d3cba4e8b2c274502e093c161 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Fri, 10 Jul 2026 10:23:45 +1000 Subject: [PATCH 17/22] Make joblib backend and verbosity configurable instead of hardcoded backend='loky' is now a parameter (default) threaded through every Parallel call, since raw 'multiprocessing' breaks in Jupyter/IPython (spawned workers try to re-run the ipykernel launcher). Callers on HPC/SLURM can still pass backend='multiprocessing'. verbose is now a tiered scheme: 0 silent, 1 (default) stage announcements, 2 adds joblib's per-task Parallel output. Also fixes a TypeError in diff_lc where wcs.all_world2pix returned a length-1 array that newer NumPy refuses to cast directly to int/float. --- tessreduce/adaptive_background.py | 20 ++++- tessreduce/background.py | 6 +- tessreduce/helpers.py | 40 +++++---- tessreduce/lastpercent.py | 4 +- tessreduce/rescale_straps.py | 4 +- tessreduce/sep_aligner.py | 5 +- tessreduce/tessreduce.py | 138 ++++++++++++++++++------------ 7 files changed, 132 insertions(+), 85 deletions(-) diff --git a/tessreduce/adaptive_background.py b/tessreduce/adaptive_background.py index 2c09ab5..4f985ea 100644 --- a/tessreduce/adaptive_background.py +++ b/tessreduce/adaptive_background.py @@ -128,6 +128,8 @@ def adaptive_medfilt_3d( per_pixel_norm=True, n_levels=7, n_jobs=1, + backend='loky', + verbose=0, metric='deviation', coarse_windows=(11, 21, 51, 101), combined_weight=0.5, @@ -236,7 +238,7 @@ def _compute_dev(cw, data_metric=data_metric, segments=segments, gw=gw): dev[s:e] = median_filter(diff, size=(dw, 1, 1), mode='reflect') return dev - all_devs = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)(delayed(_compute_dev)(cw) for cw in scales) + all_devs = Parallel(n_jobs=n_jobs, backend=backend, verbose=verbose)(delayed(_compute_dev)(cw) for cw in scales) _frame_mean = data_metric.mean(axis=(1, 2)) _clipped = _frame_mean[np.isfinite(_frame_mean)] @@ -280,7 +282,7 @@ def _compute_norm(i_dev, bright_mask=bright_mask, segments=segments, norm_scale = norm_scale * bright_mask return norm_scale - scale_norms = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( + scale_norms = Parallel(n_jobs=n_jobs, backend=backend, verbose=verbose)( delayed(_compute_norm)((i, dev)) for i, dev in enumerate(all_devs) ) @@ -377,7 +379,7 @@ def _norm01(x): def _smooth_seg(w, seg=seg_data): return w, median_filter(seg, size=(w, 1, 1), mode='reflect') - for w, smoothed_w in Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)(delayed(_smooth_seg)(w) for w in seg_levels): + for w, smoothed_w in Parallel(n_jobs=n_jobs, backend=backend, verbose=verbose)(delayed(_smooth_seg)(w) for w in seg_levels): result[s:e][seg_wins == w] = smoothed_w[seg_wins == w] if sigma_clip is not None: @@ -544,13 +546,15 @@ class AdaptiveBackground: >>> smoothed = ab.smoothed """ - def __init__(self, data, time, sector, camera, data_path=None, n_jobs=-1, block_size=5): + def __init__(self, data, time, sector, camera, data_path=None, n_jobs=-1, block_size=5, backend='loky', verbose=0): self.data = np.asarray(data, dtype=np.float32) self.time = np.asarray(time, dtype=float) self.sector = int(sector) self.camera = int(camera) self.data_path = data_path self.n_jobs = n_jobs + self.backend = backend + self.verbose = verbose self.block_size = int(block_size) self.smoothed = None @@ -578,6 +582,8 @@ def smooth( per_pixel_norm=True, n_levels=7, n_jobs=None, + backend=None, + verbose=None, metric='deviation', coarse_windows=(11, 21, 51, 101), combined_weight=0.5, @@ -618,6 +624,10 @@ def smooth( else: if n_jobs is None: n_jobs = self.n_jobs + if backend is None: + backend = self.backend + if verbose is None: + verbose = self.verbose if block_size is None: block_size = self.block_size self.smoothed, self.windows, self.variability, self._windows_pre_smooth = ( @@ -633,6 +643,8 @@ def smooth( per_pixel_norm=per_pixel_norm, n_levels=n_levels, n_jobs=n_jobs, + backend=backend, + verbose=verbose, metric=metric, coarse_windows=coarse_windows, combined_weight=combined_weight, diff --git a/tessreduce/background.py b/tessreduce/background.py index c398831..aacacce 100755 --- a/tessreduce/background.py +++ b/tessreduce/background.py @@ -10,13 +10,15 @@ class Background(): - def __init__(self, flux, mask, buffer=3, extrapolate=True, parallel=True): + def __init__(self, flux, mask, buffer=3, extrapolate=True, parallel=True, backend='loky', verbose=0): self.flux = flux self.mask = mask self.buffer = buffer self.extrapolate = extrapolate self.parallel = parallel + self.backend = backend + self.verbose = verbose self.size = self._check_size() self.smooth_bkg = np.zeros_like(flux) @@ -80,7 +82,7 @@ def _smooth_wrapper(self): if self.parallel: num_cores = multiprocessing.cpu_count() - bkg_smth = Parallel(n_jobs=num_cores, backend="multiprocessing", verbose=1)( + bkg_smth = Parallel(n_jobs=num_cores, backend=self.backend, verbose=self.verbose)( delayed(Smooth_bkg)(frame) for frame in flux * m) else: bkg_smth = np.zeros_like(flux) * np.nan diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index a5e5d09..b814f52 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -194,8 +194,8 @@ def _background2d_frame(frame, box_size, filter_size, sigma, maxiters, mask=None except Exception: return np.full_like(frame, np.nanmedian(frame)) -def parallel_background2d(cube, box_size=5, filter_size=3, sigma=3, maxiters=5, n_jobs=-1, mask=None): - return np.array(Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( +def parallel_background2d(cube, box_size=5, filter_size=3, sigma=3, maxiters=5, n_jobs=-1, mask=None, backend='loky', verbose=0): + return np.array(Parallel(n_jobs=n_jobs, backend=backend, verbose=verbose)( delayed(_background2d_frame)(frame, box_size, filter_size, sigma, maxiters, mask) for frame in cube)) @@ -255,9 +255,9 @@ def Smooth_bkg(data, gauss_smooth=0, interpolate=False, extrapolate=True): gauss_smooth = gauss_smooth * 4 estimate = gaussian_filter(estimate,gauss_smooth) else: - estimate = np.zeros_like(data) * np.nan + estimate = np.zeros_like(data) else: - estimate = np.zeros_like(data) #* np.nan + estimate = np.zeros_like(data) return estimate @@ -1117,7 +1117,7 @@ def _clip_region(image, rx, ry, sigma, iters): cut = (image[ry, rx] >= me + sigma * s) | (image[ry, rx] <= me - sigma * s) return rx[cut], ry[cut] -def regional_stats_mask(image, size=90, sigma=3, iters=10, n_jobs=1): +def regional_stats_mask(image, size=90, sigma=3, iters=10, n_jobs=1, backend='loky', verbose=0): if size < 30: print('!!! Region size is small !!!') sx, sy = image.shape @@ -1128,7 +1128,7 @@ def regional_stats_mask(image, size=90, sigma=3, iters=10, n_jobs=1): region_pixels = [np.where(regions == i) for i in range(max_reg + 1)] clip = np.zeros_like(image) - results = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( + results = Parallel(n_jobs=n_jobs, backend=backend, verbose=verbose)( delayed(_clip_region)(image, rx, ry, sigma, iters) for rx, ry in region_pixels) for rx_cut, ry_cut in results: @@ -1347,13 +1347,13 @@ def _strap_fit_col(col_data, norm_col): q[:] = np.nanmedian(q) return q -def parallel_strap_fit(frame, frame_bkg, frame_err, mask, repeats=3, tol=3, n_jobs=1): +def parallel_strap_fit(frame, frame_bkg, frame_err, mask, repeats=3, tol=3, n_jobs=1, backend='loky', verbose=0): norm = frame / frame_bkg sind = np.where(np.nansum(mask, axis=0) > 0)[0] qe = np.ones_like(frame) if len(sind) == 0: return qe - results = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( + results = Parallel(n_jobs=n_jobs, backend=backend, verbose=verbose)( delayed(_strap_fit_col)(frame[:, i], norm[:, i]) for i in sind) for col, q in zip(sind, results): @@ -1558,10 +1558,12 @@ def _block_sigma(r, box): def _fit_residual_bkg(residual, exclude_mask, res_box, n_sigma): from photutils.background import Background2D, MedianBackground from astropy.stats import SigmaClip + if (~exclude_mask).sum() < 4: + return np.zeros_like(residual) sc = SigmaClip(sigma=3.0, maxiters=5) finite_vals = residual[~exclude_mask & np.isfinite(residual)] - med = np.nanmedian(finite_vals) - std = np.nanstd(finite_vals) + med = np.nanmedian(finite_vals) if finite_vals.size > 0 else 0.0 + std = np.nanstd(finite_vals) if finite_vals.size > 0 else 0.0 transient_mask = exclude_mask | (np.abs(residual - med) > 5 * std) try: b = Background2D(residual, box_size=res_box, filter_size=3, @@ -1569,7 +1571,8 @@ def _fit_residual_bkg(residual, exclude_mask, res_box, n_sigma): mask=transient_mask, fill_value=0.0) corr = b.background except Exception: - corr = np.full_like(residual, np.nanmedian(residual[~transient_mask])) + valid = residual[~transient_mask] + corr = np.full_like(residual, np.nanmedian(valid) if valid.size > 0 else 0.0) corr_resid = residual - corr _, corr_med, corr_std = sigma_clipped_stats(corr_resid[~exclude_mask]) flagged = np.abs(corr_resid) > n_sigma * corr_std @@ -1594,7 +1597,7 @@ def fix_background_anomalies(bkg, mask, flux=None, bkg_prev=None, bkgmask=None, box_size=16, anom_box=30, anom_box_fine=4, dilate_r=2, gauss_smooth=2, sep_thresh=3.0, sep_snr_thresh=2.0, sep_validate=True, high_bkg_frames=None, high_bkg_thresh=200.0, - bad_bkg_sigma=10.0, bad_bkg_min_area=100, n_jobs=-1): + bad_bkg_sigma=10.0, bad_bkg_min_area=100, n_jobs=-1, backend='loky', verbose=0): """ Fix anomalies (asteroids, cosmic rays) in a background cube. @@ -1659,7 +1662,7 @@ def _bkgmask_i(i): return None return bkgmask_arr[i] if bkgmask_arr.ndim == 3 else bkgmask_arr - results = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( + results = Parallel(n_jobs=n_jobs, backend=backend, verbose=verbose)( delayed(_fix_bkg_frame)( bkg[i], flux[i] if flux is not None else None, @@ -1683,7 +1686,7 @@ def _bkgmask_i(i): res_box = max(res_box, 4) residuals = flux - bkg_fixed - corrections = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( + corrections = Parallel(n_jobs=n_jobs, backend=backend, verbose=verbose)( delayed(_fit_residual_bkg)( residuals[i], np.isnan(bkgmask_arr[i]) if bkgmask_arr.ndim == 3 else np.isnan(bkgmask_arr), @@ -1800,11 +1803,14 @@ def _blend_frame(bkg_new_i, bkg_prev_i, delta_i, resid_prev_i, sigma, sharp_mask if diff_mask.any(): w[diff_mask] = np.nan w = interpolate_replace_nans(w, Gaussian2DKernel(1.0)) + w = np.where(np.isfinite(w), w, 0.0) if sharp_mask_i is not None: w[sharp_mask_i] = 0.0 - return (1 - w) * bkg_new_i + w * bkg_prev_i + result = (1 - w) * bkg_new_i + w * bkg_prev_i + result = np.where(np.isfinite(result), result, bkg_prev_i) + return result -def blend_dynamic_background(bkg_new, bkg_prev, flux, sigma=2.0, sharp_masks=None, n_jobs=1): +def blend_dynamic_background(bkg_new, bkg_prev, flux, sigma=2.0, sharp_masks=None, n_jobs=1, backend='loky', verbose=0): """Per-pixel blend bkg_new toward bkg_prev based on residual quality. For each pixel, compares |flux - bkg_new| vs |flux - bkg_prev|. Where @@ -1835,7 +1841,7 @@ def blend_dynamic_background(bkg_new, bkg_prev, flux, sigma=2.0, sharp_masks=Non T = bkg_new.shape[0] sharp_list = [sharp_masks[i] if sharp_masks is not None else None for i in range(T)] - results = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( + results = Parallel(n_jobs=n_jobs, backend=backend, verbose=verbose)( delayed(_blend_frame)(bkg_new[i], bkg_prev[i], delta[i], resid_prev[i], sigma, sharp_list[i]) for i in range(T)) diff --git a/tessreduce/lastpercent.py b/tessreduce/lastpercent.py index 624b9d8..9721ba0 100644 --- a/tessreduce/lastpercent.py +++ b/tessreduce/lastpercent.py @@ -55,7 +55,7 @@ def _find_bkg_cor(tess,cores): coord = np.c_[y,x] cors = np.zeros_like(tess.ref) - cor = Parallel(n_jobs=cores, backend="multiprocessing", verbose=1)(delayed(_parallel_correlation) + cor = Parallel(n_jobs=cores, backend=getattr(tess, 'backend', 'loky'), verbose=(1 if getattr(tess, 'verbose', 1) >= 2 else 0))(delayed(_parallel_correlation) (tess.flux[:,coord[i,0],coord[i,1]], tess.bkg[:,coord[i,0],coord[i,1]], cors,coord[i],30) for i in range(len(coord))) @@ -156,7 +156,7 @@ def multi_correlation_cor(tess, limit=0.8, cores=7): if len(y) == 0: return flux, bkg - results = Parallel(n_jobs=cores, backend="multiprocessing", verbose=1)( + results = Parallel(n_jobs=cores, backend=getattr(tess, 'backend', 'loky'), verbose=(1 if getattr(tess, 'verbose', 1) >= 2 else 0))( delayed(_correct_pixel_correlation)( tess.flux[:, y[i], x[i]], tess.bkg[:, y[i], x[i]], diff --git a/tessreduce/rescale_straps.py b/tessreduce/rescale_straps.py index c8efdb7..13f4bc0 100755 --- a/tessreduce/rescale_straps.py +++ b/tessreduce/rescale_straps.py @@ -108,7 +108,7 @@ def calc_strap_factor(i,breaks,size,av_size,normals,data): qe[:,normals[b]+1+j] = factor return qe -def correct_straps(Image,mask,av_size=5,parallel=True): +def correct_straps(Image,mask,av_size=5,parallel=True,backend='loky',verbose=0): data = deepcopy(Image) mask = deepcopy(mask) av_size = int(av_size) @@ -123,7 +123,7 @@ def correct_straps(Image,mask,av_size=5,parallel=True): if parallel: num_cores = multiprocessing.cpu_count() x = np.arange(0,len(breaks),dtype=int) - qe = np.array(Parallel(n_jobs=num_cores, backend="multiprocessing", verbose=1)(delayed(calc_strap_factor)(i,breaks,size,av_size,normals,data) for i in x)) + qe = np.array(Parallel(n_jobs=num_cores, backend=backend, verbose=verbose)(delayed(calc_strap_factor)(i,breaks,size,av_size,normals,data) for i in x)) qe = np.nanmedian(qe,axis=0) qe[np.isnan(qe)] = 1 else: diff --git a/tessreduce/sep_aligner.py b/tessreduce/sep_aligner.py index 0f14abe..5c856b3 100644 --- a/tessreduce/sep_aligner.py +++ b/tessreduce/sep_aligner.py @@ -559,6 +559,7 @@ def __init__(self, ref: np.ndarray, flux: np.ndarray, n_jobs: int = -1, + backend: str = 'loky', thresh: float = 3.0, sat_frac: float = 0.7, ell_max: float = 0.5, @@ -581,6 +582,7 @@ def __init__(self, self.ref = ref self.flux = flux self.n_jobs = n_jobs + self.backend = backend self.thresh = thresh self.sat_frac = sat_frac self.ell_max = ell_max @@ -645,6 +647,7 @@ def from_tessreduce(cls, tr, **kwargs) -> 'SepAligner': inst = cls(ref=np.asarray(tr.ref), flux=np.asarray(tr.flux), n_jobs=kwargs.pop('n_jobs', n_jobs), + backend=kwargs.pop('backend', getattr(tr, 'backend', 'loky')), pixel_mask=kwargs.pop('pixel_mask', tr_mask), source_mask=kwargs.pop('source_mask', tr_source_mask), **kwargs) @@ -703,7 +706,7 @@ def run(self, time: Optional[np.ndarray] = None, self._sub_ref, cores, positions, np.asarray(weights), p['core_half']) results = Parallel(n_jobs=self.n_jobs, verbose=verbose, - backend="multiprocessing")( + backend=self.backend)( delayed(_align_one_frame)( t, self.flux[t], ref_comp, w_cols, diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index d337dae..3015ac9 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -22,7 +22,7 @@ import multiprocessing from joblib import Parallel, delayed -def _available_cores(): +def _available_cores(verbose=1): """Return the number of CPU cores available to this process. Checks in priority order: @@ -34,18 +34,21 @@ def _available_cores(): if slurm is not None: try: n = int(slurm) - print(f'[tessreduce] _available_cores: SLURM_CPUS_PER_TASK={slurm} → {n} cores') + if verbose > 0: + print(f'[tessreduce] _available_cores: SLURM_CPUS_PER_TASK={slurm} → {n} cores') return n except ValueError: pass try: n = len(os.sched_getaffinity(0)) - print(f'[tessreduce] _available_cores: sched_getaffinity → {n} cores') + if verbose > 0: + print(f'[tessreduce] _available_cores: sched_getaffinity → {n} cores') return n except AttributeError: pass n = multiprocessing.cpu_count() - print(f'[tessreduce] _available_cores: cpu_count fallback → {n} cores') + if verbose > 0: + print(f'[tessreduce] _available_cores: cpu_count fallback → {n} cores') return n from .catalog_tools import * @@ -89,7 +92,7 @@ def _fit_bkg_surface_frame(residual, exclude_mask, box_size, filter_size, sigma) golden_mean = (np.sqrt(5)-1.0)/2.0 # Aesthetic ratio fig_width = fig_width_pt*inches_per_pt # width in inches -def _subtract_residual_surface(bkg, flux, bkgmask, box_size=20, filter_size=5, sigma=3.0): +def _subtract_residual_surface(bkg, flux, bkgmask, box_size=20, filter_size=5, sigma=3.0, backend='loky', verbose=0): """ Fit and subtract a smooth 2D residual surface per frame using photutils Background2D. Operates on (flux - bkg) to capture large-scale structure @@ -129,7 +132,7 @@ def _subtract_residual_surface(bkg, flux, bkgmask, box_size=20, filter_size=5, s n_jobs = _available_cores() residuals = flux - bkg - corrections = Parallel(n_jobs=n_jobs, backend="multiprocessing", verbose=1)( + corrections = Parallel(n_jobs=n_jobs, backend=backend, verbose=verbose)( delayed(_fit_bkg_surface_frame)(residuals[i], exclude_mask, box_size, filter_size, sigma) for i in range(flux.shape[0])) bkg += np.array(corrections) @@ -141,7 +144,7 @@ class tessreduce(): def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sector=None, reduce=True,align=True,diff=True,corr_correction=False,kernel_match=False,calibrate=True,sourcehunt=True, - phot_method='aperture',imaging=False,parallel=True,num_cores=-1,diagnostic_plot=False,plot=True, + phot_method='aperture',imaging=False,parallel=True,num_cores=-1,backend='loky',diagnostic_plot=False,plot=True, savename=None,quality_bitmask='hard',cache_dir=None,cache=True,catalogue_path=False, shift_method='sep_core',use_error_image=False,prf_path=None,verbose=1,col_offset=0, bkg_temporal_window=501,ref_ind=None,ref_type='stack',ref_time_window=2,vector_path=None, @@ -189,6 +192,10 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect Perform computation with parallel processing using 'num_cores'. The default is True. num_cores : int, optional Number of cores to run parallel process on. The default is -1 which uses max system cores. + backend : str, optional + joblib backend used for parallel processing. 'loky' (default) is robust in interactive + sessions such as Jupyter/IPython. Use 'multiprocessing' on servers/clusters (e.g. SLURM/OzStar) + where it is required. diagnostic_plot : bool, optional During reduction, plot figures which outline various calculation steps, such as the image shifts over time or the zeropoint calculation. The default is False. plot : bool, optional @@ -204,7 +211,8 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect prf_path : str, optional Path to local TESS PRF files. The default is currently a specific location on the OzStar supercomputer. verbose : int, optional - Controls the level of verbosity, 0 is none, 1 is verbose. The default is 1. + Controls the level of verbosity. 0 is silent, 1 (default) prints reduction stage + announcements, 2 additionally prints joblib's per-task Parallel output. timing : bool, optional Print execution time reports for major pipeline blocks in background() and reduce(). The default is False. @@ -230,9 +238,10 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect self._center_mask = center_mask self.imaging = imaging self.parallel = parallel + self.backend = backend self._col_offset = col_offset if num_cores == -1 or isinstance(num_cores, str): - self.num_cores = _available_cores() + self.num_cores = _available_cores(verbose=verbose) else: self.num_cores = num_cores self._assign_phot_method(phot_method) @@ -250,11 +259,12 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect self._cache_path = None # SLURM environment diagnostics - _slurm_vars = ['SLURM_CPUS_PER_TASK', 'SLURM_NTASKS', 'SLURM_NTASKS_PER_NODE', - 'SLURM_JOB_CPUS_PER_NODE', 'SLURM_CPUS_ON_NODE'] - _slurm_env = {k: os.environ.get(k, 'not set') for k in _slurm_vars} - print(f'[tessreduce] SLURM env: {_slurm_env}') - print(f'[tessreduce] num_cores resolved to: {self.num_cores}') + if verbose > 0: + _slurm_vars = ['SLURM_CPUS_PER_TASK', 'SLURM_NTASKS', 'SLURM_NTASKS_PER_NODE', + 'SLURM_JOB_CPUS_PER_NODE', 'SLURM_CPUS_ON_NODE'] + _slurm_env = {k: os.environ.get(k, 'not set') for k in _slurm_vars} + print(f'[tessreduce] SLURM env: {_slurm_env}') + print(f'[tessreduce] num_cores resolved to: {self.num_cores}') # Offline Paths if catalogue_path is None: @@ -359,6 +369,11 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect if reduce: self.reduce() + @property + def _joblib_verbose(self): + """joblib Parallel verbosity: 0 unless self.verbose requests joblib output (>=2).""" + return 1 if self.verbose >= 2 else 0 + def check_coord(self): """ Checks if target coordinate / name input is valid. @@ -662,7 +677,7 @@ def psf_source_mask(self,sigma=5): data = (self._flux_aligned - self.ref) #* mask if self.parallel: try: - m = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_source_mask)(frame,self.prf,sigma) for frame in data) + m = Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(par_psf_source_mask)(frame,self.prf,sigma) for frame in data) m = np.array(m) except: m = np.ones_like(data) @@ -789,10 +804,12 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa bkg_smth = np.zeros_like(flux) * np.nan if self.parallel: _t = time.perf_counter() - print('smooth background...') - bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*m) + if self.verbose > 0: + print('smooth background...') + bkg_smth = Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*m) _times['initial smooth background'] = time.perf_counter() - _t - print('smooth background done') + if self.verbose > 0: + print('smooth background done') if rerun_negative: _t = time.perf_counter() if self._use_error_image: @@ -813,11 +830,12 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa else: m[over_sub] = 1 self._bkgmask = m - print('smooth background rerun (negative correction)...') - bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(Smooth_bkg)(frame,gauss_smooth,interpolate) for frame in flux*m) + if self.verbose > 0: + print('smooth background rerun (negative correction)...') + bkg_smth = Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(Smooth_bkg)(frame,gauss_smooth,interpolate) for frame in flux*m) _times['negative over-subtraction rerun'] = time.perf_counter() - _t - print('smooth background rerun done') - + if self.verbose > 0: + print('smooth background rerun done') if rerun_diff: _t = time.perf_counter() @@ -872,10 +890,11 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa new_mask = abs(new_mask - 1) self._bkgmask = new_mask bkg_s1 = np.array(bkg_smth) - print('smooth background rerun (residual surface)...') - bkg_smth = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*new_mask) + if self.verbose > 0: + print('smooth background rerun (residual surface)...') + bkg_smth = Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*new_mask) if blend_dynamic: - bkg_smth = blend_dynamic_background(bkg_smth, bkg_s1, flux, n_jobs=self.num_cores) + bkg_smth = blend_dynamic_background(bkg_smth, bkg_s1, flux, n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose) _times['residual surface rerun'] = time.perf_counter() - _t else: @@ -884,7 +903,8 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa bkg_smth[i] = Smooth_bkg((flux*m)[i],0,interpolate) _times['initial smooth background'] = time.perf_counter() - _t else: - print('Small tpf, using percentile cut background') + if self.verbose > 0: + print('Small tpf, using percentile cut background') self.small_background() bkg_smth = self.bkg @@ -905,27 +925,31 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa _high_bkg_frames = (_earth_angle < 30.0) | (_moon_angle < 30.0) | (_bkg_median > 300.0) else: _high_bkg_frames = _bkg_median > 300.0 - print('fixing background anomalies...') + if self.verbose > 0: + print('fixing background anomalies...') self.bkg, _sharp_masks, self.bad_bkg = fix_background_anomalies(self.bkg, self.mask, flux=deepcopy(self.flux), bkg_prev=bkg_pre_fix if blend_dynamic else None, bkgmask=self._bkgmask, gauss_smooth=gauss_smooth, high_bkg_frames=_high_bkg_frames, - n_jobs=self.num_cores) + n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose) _times['anomaly fixing'] = time.perf_counter() - _t - print('fixing background anomalies done') + if self.verbose > 0: + print('fixing background anomalies done') _t = time.perf_counter() - print('adaptive temporal smoothing...') + if self.verbose > 0: + print('adaptive temporal smoothing...') from .adaptive_background import AdaptiveBackground smoother = AdaptiveBackground(self.bkg, self.mjd, sector=self.sector, camera=self.tpf.camera, - data_path=self._vector_path,n_jobs=self.num_cores) + data_path=self._vector_path,n_jobs=self.num_cores,backend=self.backend,verbose=self._joblib_verbose) if smoother._df is not None: smoothed = smoother.smooth(method='savgol').smoothed self.bkg = smoothed _times['adaptive temporal smoothing'] = time.perf_counter() - _t - print('adaptive temporal smoothing done') + if self.verbose > 0: + print('adaptive temporal smoothing done') # Store data-driven sources from _bkgmask as bit 8, preserving the catalogue mask (bit 1) if rerun_diff: @@ -939,7 +963,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa bkg2d_mask = np.any(bkg2d_mask, axis=0) bkg_corr = parallel_background2d(f, box_size=9, filter_size=3, sigma=3, maxiters=5, n_jobs=self.num_cores, - mask=bkg2d_mask) + mask=bkg2d_mask, backend=self.backend, verbose=self._joblib_verbose) self.bkg += bkg_corr bkgmask_arr = np.asarray(self._bkgmask) @@ -1020,7 +1044,7 @@ def _bkg_round_3(self,iters=5): kern = np.ones((1,3,3)) dist_mask = convolve(dist_mask,kern) > 0 if self.parallel: - bkg_3 = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(parallel_bkg3)(self.bkg[i],dist_mask[i]) + bkg_3 = Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(parallel_bkg3)(self.bkg[i],dist_mask[i]) for i in np.arange(len(dist_mask))) else: bkg_3 = np.zeros_like(self.bkg) @@ -1044,7 +1068,7 @@ def _clip_background(self,sigma=5,ideal_size=90): """ if self.parallel: - bkg_clip = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(clip_background)(self.bkg[i],self.mask,sigma,ideal_size) + bkg_clip = Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(clip_background)(self.bkg[i],self.mask,sigma,ideal_size) for i in np.arange(len(self.bkg))) else: bkg_clip = np.zeros_like(self.bkg) @@ -1070,7 +1094,7 @@ def _grad_bkg_clip(self,sigma=3,max_size=1000): """ if self.parallel: - bkg_clip = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(grad_clip_fill_bkg)(self.bkg[i],sigma,max_size) + bkg_clip = Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(grad_clip_fill_bkg)(self.bkg[i],sigma,max_size) for i in np.arange(len(self.bkg))) else: bkg_clip = np.zeros_like(self.bkg) @@ -1426,7 +1450,7 @@ def centroids_shifts_starfind(self,plot=None,savename=None): self._dat_sources = s.to_pandas() if self.parallel: - shifts = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)( + shifts = Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)( delayed(Calculate_shifts)(frame,mx,my,finder) for frame in f) shifts = np.array(shifts) else: @@ -1498,7 +1522,7 @@ def fit_shift(self,smooth=True,plot=None,savename=None): if self.parallel: ind = np.arange(len(f)) - shifts = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)( + shifts = Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)( #delayed(difference_shifts)(f[i],m,self.eflux[i],eref) for i in ind) delayed(difference_shifts)(f[i],m) for i in ind) shifts = np.array(shifts) @@ -1572,7 +1596,7 @@ def shift_images(self,median=False): shifted[nans] = 0. if median: if self.parallel: - result = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)( + result = Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)( delayed(_shift_ref_one)(self.ref, shifted[i], self.shift[i]) for i in range(len(shifted))) shifted = np.array(result) @@ -1583,7 +1607,7 @@ def shift_images(self,median=False): else: if self.parallel: - result = Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)( + result = Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)( delayed(_shift_one)(shifted[i], self.shift[i]) for i in range(len(shifted))) self.flux = np.array(result) @@ -1792,12 +1816,12 @@ def diff_lc(self,time=None,x=None,y=None,ra=None,dec=None,tar_ap=3, if (ra is not None) & (dec is not None) & (self.tpf is not None): x,y = self.wcs.all_world2pix(ra,dec,0) - x = int(np.round(x,0)) - y = int(np.round(y,0)) + x = int(np.round(np.ravel(x)[0],0)) + y = int(np.round(np.ravel(y)[0],0)) elif (x is None) & (y is None): x,y = self.wcs.all_world2pix(self.ra,self.dec,0) - x = int(np.round(x,0)) - y = int(np.round(y,0)) + x = int(np.round(np.ravel(x)[0],0)) + y = int(np.round(np.ravel(y)[0],0)) ap_tar = np.zeros_like(data[0]) ap_sky = np.zeros_like(data[0]) @@ -1933,8 +1957,8 @@ def dif_diag_plot(self,ap_tar,ap_sky,lc=None,sky=None,data=None): maxind = np.where((np.nanmax(lc[1]) == lc[1]))[0] try: maxind = maxind[0] - except: - pass + except (IndexError, ValueError): + return d = data[maxind] nonan1 = np.isfinite(d) nonan2 = np.isfinite(d*ap) @@ -2291,7 +2315,7 @@ def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2): raise ValueError(m) inds = np.arange(0,len(xpos)) if self.parallel: - prfs, cutouts, ecutouts = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_initialise)(self.flux,self.tpf.camera,self.tpf.ccd, + prfs, cutouts, ecutouts = zip(*Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(par_psf_initialise)(self.flux,self.tpf.camera,self.tpf.ccd, self.tpf.sector,self.tpf.column,self.tpf.row, size,[xpos[i],ypos[i]],time_ind) for i in inds)) else: @@ -2304,7 +2328,7 @@ def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2): cutouts = np.array(cutouts) print('made cutouts') if self.parallel: - flux, pos = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_full)(cutouts[i],prfs[i],self.shift[i],xlim,ylim) for i in inds)) + flux, pos = zip(*Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(par_psf_full)(cutouts[i],prfs[i],self.shift[i],xlim,ylim) for i in inds)) else: flux = [] pos = [] @@ -2390,14 +2414,14 @@ def psf_photutils(self,xPix=None,yPix=None,size=5,local_bkg=False,epsf=None, eflux = np.zeros(len(self.flux)) * np.nan psfphot2 = PSFPhotometry(epsf, fit_shape, finder=None,aperture_radius=1.5, xy_bounds=(0.05),localbkg_estimator=localbkg_estimator) - f,ef = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot2,init) for i in np.arange(len(cutouts)))) + f,ef = zip(*Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot2,init) for i in np.arange(len(cutouts)))) f = np.array(f).flatten() ef = np.array(ef).flatten() phot = phot.to_pandas() pos = phot[['x_fit','y_fit']].values + np.array([xPix,yPix]) - size//2 epos = phot[['x_err','y_err']].values else: - f,ef,pos,epos = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot,init,True) for i in np.arange(len(cutouts)))) + f,ef,pos,epos = zip(*Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(parallel_photutils)(cutouts[i],ecutouts[i],psfphot,init,True) for i in np.arange(len(cutouts)))) pos['x_fit'] += xPix - size//2 pos['y_fit'] += xPix - size//2 @@ -2461,7 +2485,7 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa prf, cutouts, ecutouts = self._psf_initialise(size,(xPix,yPix),ref=(not diff)) # gather base PRF and the array of cutouts data inds = np.arange(len(cutouts)) base = create_psf(prf,size) - flux, eflux, pos = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_full)(cutouts[i],base,self.shift[i]) for i in inds)) + flux, eflux, pos = zip(*Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(par_psf_full)(cutouts[i],base,self.shift[i]) for i in inds)) #prf, cutouts = self._psf_initialise(size,(xPix,yPix)) # gather base PRF and the array of cutouts data #xShifts = [] @@ -2507,9 +2531,9 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa if self.parallel: inds = np.arange(len(cutouts)) if self.delta_kernel is not None: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) else: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) else: for i in range(len(cutouts)): flux += [par_psf_flux(cutouts[i],ecutouts[i],base,self.shift[i])] @@ -2525,9 +2549,9 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa if self.parallel: inds = np.arange(len(cutouts)) if self.delta_kernel is not None: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order,self.delta_kernel[i]) for i in inds)) else: - flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) + flux, eflux = zip(*Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(par_psf_flux)(cutouts[i],ecutouts[i],base,self.shift[i],bkg_poly_order) for i in inds)) else: for i in range(len(cutouts)): flux += [par_psf_flux(cutouts[i],ecutouts[i],base,self.shift[i])] @@ -2572,7 +2596,7 @@ def kernel_matching(self,size=7,diff=True): mask = self.mask == 1 if self.parallel: - d, kernel = zip(*Parallel(n_jobs=self.num_cores, backend="multiprocessing", verbose=1)(delayed(parallel_delta_diff)(frame,self.ref,mask,size) for frame in flux)) + d, kernel = zip(*Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(parallel_delta_diff)(frame,self.ref,mask,size) for frame in flux)) else: d = [] kernel = [] @@ -2728,7 +2752,7 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, elif self._shift_method == 'sep_core': from .sep_aligner import SepAligner aligner = SepAligner.from_tessreduce(self) - aligner.run() + aligner.run(verbose=self._joblib_verbose) if self._smooth_motion: aligner.smooth_shift(time=self.mjd, gap_thresh=0.5, # days From 681a504c8ba4abe4452634de477de3e862641209 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Fri, 10 Jul 2026 15:12:18 +1000 Subject: [PATCH 18/22] Update PS1/SkyMapper-to-TESS synthetic magnitude coefficients via calibrimbore The hardcoded color-term coefficients in PS1_to_TESS_mag/SM_to_TESS_mag were stale relative to what calibrimbore currently produces for the same tess.dat bandpass and g-r color cuts (i/z terms off by ~20-28%). Re-derived via sauron(band='tess.dat', system='ps1'|'skymapper', gr_lims=[-.5,.8]). --- tessreduce/catalog_tools.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/tessreduce/catalog_tools.py b/tessreduce/catalog_tools.py index f0d3828..857855a 100755 --- a/tessreduce/catalog_tools.py +++ b/tessreduce/catalog_tools.py @@ -270,8 +270,9 @@ def PS1_to_TESS_mag(PS1,ebv = 0): z = mag2flux(PS1.zmag.values - ez,zp) y = mag2flux(PS1.ymag.values - ey,zp) - cr = 0.25582823; ci = 0.27609407; cz = 0.35809516 - cy = 0.11244277; cp = 0.00049096 + # re-derived via calibrimbore (sauron(band='tess.dat',system='ps1',gr_lims=[-.5,.8])) + cr = 0.23693349; ci = 0.34742086; cz = 0.28054228 + cy = 0.13630491; cp = 0.00036188 t = (cr*r + ci*i + cz*z + cy*y)*(g/i)**cp t = -2.5*np.log10(t) + zp + et @@ -295,8 +296,9 @@ def SM_to_TESS_mag(SM,ebv = 0): i = mag2flux(SM.imag.values - ei,zp) z = mag2flux(SM.zmag.values - ez,zp) - cr = 0.25825435; ci = 0.35298213 - cz = 0.39388206; cp = -0.00170817 + # re-derived via calibrimbore (sauron(band='tess.dat',system='skymapper',gr_lims=[-.5,.8])) + cr = 0.25323566; ci = 0.42107837 + cz = 0.32697125; cp = -0.01408441 t = (cr*r + ci*i + cz*z)*(g/i)**cp t = -2.5*np.log10(t) + zp + et From 13a0fab7ea04c198f122c1b0b1983c0a20b23574 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Fri, 10 Jul 2026 15:12:34 +1000 Subject: [PATCH 19/22] Add scene-modelling PSF photometry, two new calibration pathways, tiered verbose logging, and fix a NumPy 2.x crash in diff_lc scene_photom.py: shared linear scene-fit engine (bucketed PRF cache, 2D polynomial background surface instead of a flat constant, PSF-derivative columns to absorb poor difference-image subtraction residuals, catalog- informed ridge-prior crowding for neighbours, vectorized one-shot solve across an entire time series instead of a per-frame nonlinear fit). Wired in as a new tessreduce.scene_photometry() method, alongside (not replacing) psf_photometry(). field_calibration.py: two additive calibration pathways built on the same scene-fit engine -- field_calibrate_scene(catalog='ps1'|'gaia'). 'ps1' reuses the existing Tonry-locus extinction correction and PS1/SkyMapper synthetic- TESS-mag reconstruction; 'gaia' queries Gaia DR3 via astroquery (matching tessreduce's existing catalog convention) and calibrates directly against Rp with the correct Vega->AB offset. Both use magnitude-binned robust zeropoint combining and formal per-star error gating instead of field_calibrate()'s coarse sanity check; failures warn and fall back to zp_ab=20.6 rather than raising. Neither pathway touches the existing field_calibrate(). verbose is now tiered: 1 (default) prints top-level reduce() stage announcements with timing, 2 adds background()'s internal sub-steps (also timed), 3 adds joblib's per-task Parallel output and core-selection diagnostics. diff_lc: wcs.all_world2pix can return 0-d/length-1 arrays depending on input type, which newer NumPy refuses to cast directly to int/float; flatten via np.ravel(...)[0] first. --- tessreduce/field_calibration.py | 367 +++++++++++++++++++++++++++++++ tessreduce/scene_photom.py | 245 +++++++++++++++++++++ tessreduce/tessreduce.py | 304 ++++++++++++++++++++++---- tests/test_scene_photom.py | 369 ++++++++++++++++++++++++++++++++ 4 files changed, 1247 insertions(+), 38 deletions(-) create mode 100644 tessreduce/field_calibration.py create mode 100644 tessreduce/scene_photom.py create mode 100644 tests/test_scene_photom.py diff --git a/tessreduce/field_calibration.py b/tessreduce/field_calibration.py new file mode 100644 index 0000000..c504aa6 --- /dev/null +++ b/tessreduce/field_calibration.py @@ -0,0 +1,367 @@ +""" +Scene-modelling photometric calibration. + +Two additive calibration pathways, both built on the shared scene-fit engine +in `scene_photom.py` (real pixel-space deblending, formal error propagation), +sitting alongside (not replacing) `tessreduce.py`'s existing field_calibrate(): + + calibrate_ps1_skymapper(tess, ...) -- PS1/SkyMapper, with the existing + Tonry-locus extinction correction and multi-band synthetic-TESS-mag + reconstruction (both already AB-consistent). + calibrate_gaia(tess, ...) -- Gaia DR3 (queried via astroquery, + matching tessreduce's existing catalog-fetch convention), calibrated + directly against Rp with the standard Vega->AB offset applied. No + extinction step, matching the simpler single-band reference. + +On any failure (too few good stars, fits don't converge) a warning is raised +and a fixed zp_ab=20.6 fallback is used, rather than raising an exception. +""" +import warnings +import numpy as np +import pandas as pd +from astropy.coordinates import SkyCoord, Angle + +from . import scene_photom as sp +from .catalog_tools import Get_Catalogue, PS1_to_TESS_mag, SM_to_TESS_mag +from .calibration_tools import Tonry_reduce + +# Gaia DR3 Rp Vega -> AB offset, from synthetic photometry (pysynphot Vega +# spectrum through the actual Gaia3 Rp passband, via calibrimbore's +# get_pb_zpt): zp_AB - zp_Vega = 0.379. The previously-used literature value +# of 0.152 (Casagrande & VandenBerg 2018) was ~0.23 mag too small and was the +# dominant cause of the zp offset between this pathway and PS1/SkyMapper. +GAIA_RP_AB_OFFSET = 0.379 + +_ZP_FALLBACK = 20.6 +_ZP_FALLBACK_ERR = 0.1 + + +def _zp_to_scale(zp): + """AB zeropoint (mag) -> linear flux-scale factor.""" + return 10 ** (0.4 * np.asarray(zp, dtype=float)) + + +def _scale_to_zp(f): + """Linear flux-scale factor -> AB zeropoint (mag).""" + return 2.5 * np.log10(f) + + +def _zp_fallback(reason): + warnings.warn(f'field_calibration: {reason} -- falling back to zp_ab={_ZP_FALLBACK}') + return _ZP_FALLBACK, _ZP_FALLBACK_ERR + + +def _binned_zeropoint(mag, zp, bin_width=0.5, min_bin_n=5): + """Magnitude-binned, log-space-safe robust zeropoint combine. + + Bins calibration stars by magnitude, takes a robust (median/1.4826*MAD) + zeropoint per bin in LINEAR flux-scale space (avoiding the bias of + averaging a log quantity directly), then inverse-variance-weights the + per-bin values back into a single zeropoint. Falls back to a global + robust combine if there are too few usable bins. + """ + mag = np.asarray(mag, dtype=float) + zp = np.asarray(zp, dtype=float) + finite = np.isfinite(mag) & np.isfinite(zp) + mag, zp = mag[finite], zp[finite] + + def _global(): + scale = _zp_to_scale(zp) + med = np.median(scale) + mad = 1.4826 * np.median(np.abs(scale - med)) + return _scale_to_zp(med), (mad / med) * 1.0857, None + + if len(zp) < min_bin_n: + return _global() + + lo, hi = np.floor(mag.min() / bin_width) * bin_width, np.ceil(mag.max() / bin_width) * bin_width + edges = np.arange(lo, hi + bin_width, bin_width) + if len(edges) < 3: + return _global() + + bin_scale, bin_mad, bin_n, bin_centers = [], [], [], [] + for i in range(len(edges) - 1): + sel = (mag >= edges[i]) & (mag < edges[i + 1]) + n = sel.sum() + if n < 2: + continue + scale = _zp_to_scale(zp[sel]) + med = np.median(scale) + mad = 1.4826 * np.median(np.abs(scale - med)) + bin_scale.append(med) + bin_mad.append(max(mad, 1e-6 * med)) + bin_n.append(n) + bin_centers.append(0.5 * (edges[i] + edges[i + 1])) + + if len(bin_scale) < 2: + return _global() + + bin_scale = np.array(bin_scale) + bin_mad = np.array(bin_mad) + bin_n = np.array(bin_n) + se = bin_mad / np.sqrt(bin_n) + weights = 1.0 / se ** 2 + combined_scale = np.sum(weights * bin_scale) / np.sum(weights) + zp_ab = _scale_to_zp(combined_scale) + zp_scatter = float(np.average(bin_mad / bin_scale, weights=bin_n)) * 1.0857 + bins = pd.DataFrame({'mag_center': bin_centers, 'scale': bin_scale, + 'mad': bin_mad, 'n': bin_n}) + return zp_ab, zp_scatter, bins + + +def _select_isolated(ra, dec, mag, ra_all, dec_all, mag_all, + iso_radius_pix=4.0, pix_scale=21.0, delta_mag=2.0): + """Boolean mask of candidates with no brighter Gaia/catalogue neighbour + within `iso_radius_pix` (converted to arcsec via `pix_scale`). + + Uses a flat-sky approximation (valid at these small separations), same as + TESSELLATE's `_select_isolated`. + """ + iso_arcsec = iso_radius_pix * pix_scale + ra = np.asarray(ra); dec = np.asarray(dec); mag = np.asarray(mag) + ra_all = np.asarray(ra_all); dec_all = np.asarray(dec_all); mag_all = np.asarray(mag_all) + cosdec = np.cos(np.deg2rad(dec_all)) + + keep = np.ones(len(ra), dtype=bool) + for i in range(len(ra)): + dra = (ra_all - ra[i]) * cosdec + ddec = dec_all - dec[i] + sep = np.hypot(dra, ddec) * 3600.0 + neighbour = (sep > 0.5) & (sep < iso_arcsec) & (mag_all < mag[i] + delta_mag) + if neighbour.any(): + keep[i] = False + return keep + + +def _err_keep_mask(e_flux, flux, max_err_factor=3.0): + """Quality gate on formal per-star flux error, replacing the coarse + mag<=tmag+1 sanity check used by tessreduce's legacy field_calibrate().""" + frac_err = np.abs(e_flux / np.where(flux > 0, flux, np.nan)) + med = np.nanmedian(frac_err) + return np.isfinite(frac_err) & (frac_err < max_err_factor * med) & (flux > 0) + + +def _star_stamp(image, xpix, ypix, stamp_size): + half = stamp_size // 2 + xi, yi = int(round(xpix)), int(round(ypix)) + ny, nx = image.shape + if xi - half < 0 or yi - half < 0 or xi + half + 1 > nx or yi + half + 1 > ny: + return None, None, None + stamp = image[yi - half:yi + half + 1, xi - half:xi + half + 1] + return stamp, xpix - xi, ypix - yi + + +def _fit_star(image, prf, xpix, ypix, stamp_size, poly_order, + neighbour_xy=None, flux_bounds=None): + """Single-star scene fit on a single image (the calibration reference + frame), optionally deblending catalogue neighbours. Returns + (flux, e_flux) or (None, None) if the stamp is unusable.""" + stamp, x_sub, y_sub = _star_stamp(image, xpix, ypix, stamp_size) + if stamp is None or not np.all(np.isfinite(stamp)): + return None, None + cent = (stamp_size - 1) / 2.0 + + neighbour_dxdy = [] + if neighbour_xy: + xi, yi = int(round(xpix)), int(round(ypix)) + for nx, ny in neighbour_xy: + neighbour_dxdy.append((nx - xi, ny - yi)) + + A, info = sp.build_design_matrix(prf, cent, (x_sub, y_sub), neighbour_dxdy, stamp_size, + poly_order=poly_order, include_psf_derivatives=False) + bounds = None + if flux_bounds is not None: + lo = np.full(A.shape[1], -np.inf) + hi = np.full(A.shape[1], np.inf) + lo[0], hi[0] = flux_bounds + bounds = (lo, hi) + coeffs, cov, dof = sp.fit_scene_frame(A, stamp.ravel(), flux_bounds=bounds) + flux = coeffs[0] + e_flux = np.sqrt(max(cov[0, 0], 0.0)) + return flux, e_flux + + +def _calibrate_common(tess, ra, dec, mag, mag_col_label, + mag_lo, mag_hi, iso_radius_pix=4.0, stamp_size=9, + poly_order=2, refine_iter=3, refine_tol=1e-3, + var_tol_mag=0.5, max_err_factor=3.0, max_zp_err=0.1, + zp_bin_width=0.5, zp_bin_min_n=5, edge_margin=5): + """Shared calibration engine: star selection + scene-fit photometry + + robust combine + iterative refinement, fed either PS1/SkyMapper `tmag` + (already extinction-corrected & AB) or Gaia `rp_ab`. + + `ra`, `dec`, `mag` are the FULL catalogue (for isolation checks and + neighbour deblending); the calibration sample is the `mag_lo mag_lo) & (mag < mag_hi) + edge_ok = ((x_all > edge_margin) & (x_all < nx - edge_margin) + & (y_all > edge_margin) & (y_all < ny - edge_margin)) + cal_idx = np.where(sel_range & edge_ok & np.isfinite(x_all) & np.isfinite(y_all))[0] + if len(cal_idx) < 10: + zp_ab, zp_err = _zp_fallback(f'too few {mag_col_label} calibration candidates ({len(cal_idx)})') + return zp_ab, zp_err, None + + iso_keep = _select_isolated(ra[cal_idx], dec[cal_idx], mag[cal_idx], ra, dec, mag, + iso_radius_pix=iso_radius_pix) + iso_idx = cal_idx[iso_keep] + if len(iso_idx) < 10: + zp_ab, zp_err = _zp_fallback(f'too few isolated {mag_col_label} stars ({len(iso_idx)})') + return zp_ab, zp_err, None + + prf_cam, prf_ccd, prf_sector = tess.tpf.camera, tess.tpf.ccd, tess.tpf.sector + + def _prf_at(xpix, ypix): + col = int(np.clip(tess.tpf.column - int(tess.size//2) + xpix + 45, 45, 2090)) + row = int(np.clip(tess.tpf.row - int(tess.size//2) + ypix + 1, 1, 2040)) + return sp._prf_cache(prf_cam, prf_ccd, prf_sector, col, row, tess._prf_path) + + # ---- stage 1: isolated-star fits, no neighbour deblending needed ---- + fluxes, e_fluxes, mags, xs, ys = [], [], [], [], [] + for i in iso_idx: + prf = _prf_at(x_all[i], y_all[i]) + flux, e_flux = _fit_star(ref, prf, x_all[i], y_all[i], stamp_size, poly_order) + if flux is None: + continue + fluxes.append(flux); e_fluxes.append(e_flux); mags.append(mag[i]) + xs.append(x_all[i]); ys.append(y_all[i]) + + fluxes = np.array(fluxes); e_fluxes = np.array(e_fluxes); mags = np.array(mags) + if len(fluxes) < 10: + zp_ab, zp_err = _zp_fallback(f'too few successful {mag_col_label} stage-1 fits ({len(fluxes)})') + return zp_ab, zp_err, None + + keep = _err_keep_mask(e_fluxes, fluxes, max_err_factor=max_err_factor) + if keep.sum() < 10: + zp_ab, zp_err = _zp_fallback(f'too few good-quality {mag_col_label} stage-1 fits ({keep.sum()})') + return zp_ab, zp_err, None + + zp_star = mags[keep] + 2.5 * np.log10(fluxes[keep]) + zp_ab, zp_err, bins = _binned_zeropoint(mags[keep], zp_star, + bin_width=zp_bin_width, min_bin_n=zp_bin_min_n) + if not np.isfinite(zp_ab): + zp_ab, zp_err = _zp_fallback(f'{mag_col_label} stage-1 zeropoint combine failed') + return zp_ab, zp_err, None + + if zp_err is not None and zp_err > max_zp_err: + # stage-2 iterative refinement: re-fit ALL in-range stars (not just + # isolated ones), deblending catalogue neighbours, with flux bounded + # around the current zeropoint-predicted value. + for _ in range(refine_iter): + fluxes2, e_fluxes2, mags2 = [], [], [] + tol = max(3 * zp_err, var_tol_mag) + for i in cal_idx: + neighbour_sel = (np.hypot(x_all - x_all[i], y_all - y_all[i]) < stamp_size) & (np.arange(len(x_all)) != i) + neighbour_xy = list(zip(x_all[neighbour_sel], y_all[neighbour_sel])) if neighbour_sel.any() else None + expected = 10 ** (-0.4 * (mag[i] - zp_ab)) + bounds = (expected * 10 ** (-0.4 * tol), expected * 10 ** (0.4 * tol)) + prf = _prf_at(x_all[i], y_all[i]) + flux, e_flux = _fit_star(ref, prf, x_all[i], y_all[i], stamp_size, poly_order, + neighbour_xy=neighbour_xy, flux_bounds=bounds) + if flux is None: + continue + fluxes2.append(flux); e_fluxes2.append(e_flux); mags2.append(mag[i]) + fluxes2 = np.array(fluxes2); e_fluxes2 = np.array(e_fluxes2); mags2 = np.array(mags2) + if len(fluxes2) < 10: + break + keep2 = _err_keep_mask(e_fluxes2, fluxes2, max_err_factor=max_err_factor) + if keep2.sum() < 10: + break + zp_star2 = mags2[keep2] + 2.5 * np.log10(fluxes2[keep2]) + zp_new, zp_err_new, bins = _binned_zeropoint(mags2[keep2], zp_star2, + bin_width=zp_bin_width, min_bin_n=zp_bin_min_n) + if not np.isfinite(zp_new): + break + converged = abs(zp_new - zp_ab) < refine_tol + zp_ab, zp_err = zp_new, zp_err_new + if converged: + break + + return zp_ab, zp_err, bins + + +def calibrate_ps1_skymapper(tess, mag_lo=None, mag_hi=None, iso_radius_pix=4.0, + stamp_size=9, poly_order=2, plot=False, **kwargs): + """Pathway 1: PS1 (dec>-30) or SkyMapper (dec<=-30), with the Tonry-locus + extinction correction and multi-band synthetic-TESS-mag reconstruction + reused exactly from tessreduce's existing calibration machinery -- both + already produce an AB-consistent `tmag`. + """ + if tess.dec < -30: + table = Get_Catalogue(tess.tpf, Catalog='skymapper') + system = 'skymapper' + if table is None: + zp_ab, zp_err = _zp_fallback('SkyMapper catalogue unavailable') + tess.cat = None + return zp_ab, zp_err + else: + table = Get_Catalogue(tess.tpf, Catalog='ps1') + system = 'ps1' + + try: + ebv, dat = Tonry_reduce(table, plot=plot, system=system) + except ValueError: + zp_ab, zp_err = _zp_fallback(f'Tonry extinction fit failed for {system}') + tess.cat = table + tess.ebv = 0.0 + return zp_ab, zp_err + tess.ebv = float(np.atleast_1d(ebv)[0]) + + table = PS1_to_TESS_mag(table, ebv=tess.ebv) if system == 'ps1' else SM_to_TESS_mag(table, ebv=tess.ebv) + x, y = tess.wcs.all_world2pix(table.RAJ2000.values, table.DEJ2000.values, 0) + table['col'] = x + table['row'] = y + tess.cat = table + + mag_lo = 8.5 if mag_lo is None else mag_lo + mag_hi = 16.0 if mag_hi is None else mag_hi + + zp_ab, zp_err, bins = _calibrate_common( + tess, table.RAJ2000.values, table.DEJ2000.values, table['tmag'].values, 'tmag', + mag_lo=mag_lo, mag_hi=mag_hi, iso_radius_pix=iso_radius_pix, + stamp_size=stamp_size, poly_order=poly_order, **kwargs) + return zp_ab, zp_err + + +def calibrate_gaia(tess, mag_lo=11.0, mag_hi=15.5, iso_radius_pix=4.0, + stamp_size=9, poly_order=2, plot=False, **kwargs): + """Pathway 2: Gaia DR3, queried via astroquery (matching tessreduce's + existing catalogue-fetch convention, not TESSELLATE's duckdb+local-CSV + approach). Calibrated directly against Gaia Rp with the standard + Vega->AB offset; no extinction correction (matching TESSELLATE, since + this pathway doesn't reconstruct a synthetic multi-band TESS magnitude). + """ + from astroquery.vizier import Vizier + Vizier.ROW_LIMIT = -1 + c1 = SkyCoord(tess.tpf.ra, tess.tpf.dec, frame='icrs', unit='deg') + pix_scale = 21.0 + rad = Angle(np.max(tess.tpf.shape[1:]) * pix_scale + 60, 'arcsec') + result = Vizier.query_region(c1, catalog=['I/355/gaiadr3'], radius=rad, + column_filters={'Gmag': '<19'}) + if result is None or len(result) == 0: + zp_ab, zp_err = _zp_fallback('Gaia DR3 query returned no sources') + tess.cat = None + return zp_ab, zp_err + + table = result['I/355/gaiadr3'].to_pandas() + if 'RPmag' not in table.columns: + zp_ab, zp_err = _zp_fallback('Gaia DR3 query missing RPmag column') + tess.cat = None + return zp_ab, zp_err + table = table[np.isfinite(table['RPmag'])].reset_index(drop=True) + table['rp_ab'] = table['RPmag'] + GAIA_RP_AB_OFFSET + + x, y = tess.wcs.all_world2pix(table['RA_ICRS'].values, table['DE_ICRS'].values, 0) + table['col'] = x + table['row'] = y + tess.cat = table + + zp_ab, zp_err, bins = _calibrate_common( + tess, table['RA_ICRS'].values, table['DE_ICRS'].values, table['rp_ab'].values, 'rp_ab', + mag_lo=mag_lo, mag_hi=mag_hi, iso_radius_pix=iso_radius_pix, + stamp_size=stamp_size, poly_order=poly_order, **kwargs) + return zp_ab, zp_err diff --git a/tessreduce/scene_photom.py b/tessreduce/scene_photom.py new file mode 100644 index 0000000..c235311 --- /dev/null +++ b/tessreduce/scene_photom.py @@ -0,0 +1,245 @@ +""" +Shared scene-modelling PSF photometry engine. + +Catalog-agnostic, pure numpy/scipy: builds a linear design matrix (target PSF +column, optional PSF-derivative columns, optional neighbour PSF columns, and a +2D polynomial background surface) at a FIXED sub-pixel position, then solves +for flux either per-frame (calibration star fitting) or for an entire time +series in a single vectorized linear solve (forced/scene photometry). + +Because position is fixed, every basis column here is frame-independent, so +the one-shot vectorized solve in `fit_scene_lightcurve` only has to invert one +small (n_columns x n_columns) matrix regardless of how many frames there are. +""" +import numpy as np +from scipy.optimize import minimize, lsq_linear + +_PRF_CACHE = {} + + +def _prf_cache(camera, ccd, sector, col, row, prf_path, bucket=100): + """Bucketed cache for TESS_PRF objects, keyed on a coarse pixel grid. + + Mirrors the sector-dependent `localdatadir` convention already used + elsewhere in tessreduce (Sectors1_2_3 vs Sectors4+), and TESSELLATE's + bucketed-cache approach: one PRF object is built per ~100x100 pixel + region per (camera, ccd, sector, prf_path), not per star. + """ + from PRF import TESS_PRF + + col = int(np.clip(col, 45, 2090)) + row = int(np.clip(row, 1, 2040)) + cb = int(np.clip((col // bucket) * bucket + bucket // 2, 45, 2090)) + rb = int(np.clip((row // bucket) * bucket + bucket // 2, 1, 2040)) + localdatadir = None + if prf_path is not None: + subdir = 'Sectors4+' if sector >= 4 else 'Sectors1_2_3' + localdatadir = f'{prf_path}/{subdir}' + key = (camera, ccd, sector, cb, rb, prf_path) + prf = _PRF_CACHE.get(key) + if prf is None: + prf = TESS_PRF(camera, ccd, sector, cb, rb, localdatadir=localdatadir) + _PRF_CACHE[key] = prf + return prf + + +def polynomial_columns(stamp_size, order): + """Linear 2D polynomial background basis, one flattened column per term. + + Uses the same triangular index scheme (i+j <= order) as the existing + nonlinear `polynomial_surface` in psf_photom.py, so `poly_order` means the + same thing in both APIs. Evaluated on a fixed, stamp-centered pixel grid + so the columns never depend on frame. order=0 returns a single all-ones + column (equivalent to a flat background). + """ + cent = (stamp_size - 1) / 2.0 + yy, xx = np.mgrid[0:stamp_size, 0:stamp_size] + xx = (xx - cent).astype(float) + yy = (yy - cent).astype(float) + + cols = [] + for i in range(order + 1): + for j in range(order + 1 - i): + cols.append(((xx ** i) * (yy ** j)).ravel()) + return cols + + +def prf_column(prf, cent, dx, dy, stamp_size): + """Evaluate + normalize + flatten the PRF at a fixed sub-pixel offset.""" + npix = stamp_size * stamp_size + p = prf.locate(cent + dx, cent + dy, (stamp_size, stamp_size)) + s = np.nansum(p) + return (p / s).ravel() if (np.isfinite(s) and s > 0) else np.zeros(npix) + + +def prf_derivative_columns(prf, cent, dx, dy, stamp_size, eps=0.01): + """Finite-difference spatial derivatives of the PRF at a fixed position. + + These absorb the sharp, PSF-shaped residual left by a subpixel-misaligned + or imperfectly kernel-matched difference-image subtraction near a source + (often a dipole pattern) -- something a smooth polynomial background + cannot reproduce. Still frame-independent (position is fixed), so this + doesn't break the vectorized one-shot solve. + """ + dpdx = (prf_column(prf, cent, dx + eps, dy, stamp_size) + - prf_column(prf, cent, dx - eps, dy, stamp_size)) / (2 * eps) + dpdy = (prf_column(prf, cent, dx, dy + eps, stamp_size) + - prf_column(prf, cent, dx, dy - eps, stamp_size)) / (2 * eps) + return dpdx, dpdy + + +def build_design_matrix(prf, cent, target_dxdy, neighbour_dxdys, stamp_size, + poly_order=2, include_psf_derivatives=True, + derivatives_for='target'): + """Stack [target PSF, target PSF-derivatives, neighbour PSFs, + neighbour PSF-derivatives, polynomial background] into a design matrix. + + Returns (A, column_info) where column_info is a dict describing which + column index corresponds to what, so callers can locate the target flux + column, neighbour columns, etc. + """ + cols = [] + info = {'target': 0, 'neighbours': [], 'derivatives': {}, 'background': []} + + cols.append(prf_column(prf, cent, target_dxdy[0], target_dxdy[1], stamp_size)) + + if include_psf_derivatives: + dpdx, dpdy = prf_derivative_columns(prf, cent, target_dxdy[0], target_dxdy[1], stamp_size) + info['derivatives']['target'] = (len(cols), len(cols) + 1) + cols.append(dpdx) + cols.append(dpdy) + + neighbour_dxdys = neighbour_dxdys or [] + for dx, dy in neighbour_dxdys: + info['neighbours'].append(len(cols)) + cols.append(prf_column(prf, cent, dx, dy, stamp_size)) + if include_psf_derivatives and derivatives_for == 'all': + dpdx, dpdy = prf_derivative_columns(prf, cent, dx, dy, stamp_size) + info['derivatives'][len(info['neighbours']) - 1] = (len(cols), len(cols) + 1) + cols.append(dpdx) + cols.append(dpdy) + + bg_cols = polynomial_columns(stamp_size, poly_order) + info['background'] = list(range(len(cols), len(cols) + len(bg_cols))) + cols.extend(bg_cols) + + A = np.column_stack(cols) + return A, info + + +def fit_scene_position(build_A, data_pix, tol_x, tol_y, x0=(0.0, 0.0)): + """2-D nonlinear position refinement. + + `build_A(dx, dy)` must return a design matrix for that trial position. + The objective is the unconstrained linear least-squares residual at each + trial position -- i.e. we search over position, but flux/background are + always solved for in closed form at each trial (ported from TESSELLATE's + `_scene_fit_worker` position-search step, generalized to the polynomial + background design matrix used here). + """ + good = np.isfinite(data_pix) + + def _chi2(p): + dx, dy = p + A = build_A(dx, dy) + Ag = A[good] + if Ag.shape[0] < Ag.shape[1] + 1: + return np.inf + sol, *_ = np.linalg.lstsq(Ag, data_pix[good], rcond=None) + resid = data_pix[good] - Ag @ sol + return float(np.nansum(resid ** 2)) + + if tol_x <= 0 and tol_y <= 0: + return 0.0, 0.0 + + opt = minimize(_chi2, list(x0), method='L-BFGS-B', + bounds=[(-tol_x, tol_x), (-tol_y, tol_y)]) + return float(opt.x[0]), float(opt.x[1]) + + +def fit_scene_frame(A, data_pix, flux_bounds=None): + """Single-frame scene solve. + + Bounded (`scipy.optimize.lsq_linear`) when `flux_bounds` is given + (calibration use -- constrain each source's flux to a catalog-predicted + range), otherwise a plain unconstrained `np.linalg.lstsq`. Also returns + the formal parameter covariance `s2 * inv(A.T@A)` for per-star error + propagation (ported from TESSELLATE's `_scene_fit_worker`). + + Returns (coeffs, cov, dof). + """ + good = np.isfinite(data_pix) & np.all(np.isfinite(A), axis=1) + Ag = A[good] + dg = data_pix[good] + + if flux_bounds is not None: + res = lsq_linear(Ag, dg, bounds=flux_bounds, method='trf', max_iter=200) + coeffs = res.x + else: + coeffs, *_ = np.linalg.lstsq(Ag, dg, rcond=None) + + resid = dg - Ag @ coeffs + dof = max(Ag.shape[0] - Ag.shape[1], 1) + s2 = float(np.nansum(resid ** 2) / dof) + cov = s2 * np.linalg.inv(Ag.T @ Ag) + return coeffs, cov, dof + + +def fit_scene_lightcurve(stamp_cube, A, prior_flux=None, prior_strength=None): + """Vectorized one-shot scene solve across an entire time series. + + Because `A` is frame-independent (position fixed), the whole time series + is solved in a single matrix product rather than one nonlinear fit per + frame. Optionally accepts a Tikhonov/ridge prior (`prior_flux`, + `prior_strength`, both length-M arrays with zeros on unregularized + columns) which pulls specific columns (e.g. crowded neighbours) toward a + catalog-predicted value without breaking the closed-form vectorized + solve -- unlike a hard per-frame bound, which would require iterating. + + Returns (flux, e_flux, background, coeffs) where `flux`/`e_flux` are the + target (column 0) flux and its formal per-frame error, and `background` + is the fitted per-frame amplitude of the first background column (the + constant term). + """ + stamp_cube = np.asarray(stamp_cube, dtype=float) + nfr = stamp_cube.shape[0] + npix = stamp_cube.shape[1] * stamp_cube.shape[2] + D = stamp_cube.reshape(nfr, npix) + + good = np.all(np.isfinite(D), axis=0) & np.all(np.isfinite(A), axis=1) + if good.sum() < A.shape[1] + 1: + raise RuntimeError('Too few finite pixels for scene photometry.') + Ag = A[good] + Dg = D[:, good] + + ATA = Ag.T @ Ag + if prior_strength is not None: + prior_strength = np.asarray(prior_strength, dtype=float) + prior_flux = np.asarray(prior_flux, dtype=float) + # a single non-finite entry here would silently NaN-poison every + # column of the solve (one bad row spreads through the matrix + # inverse) -- fail loudly instead. Callers should zero out + # prior_strength for any column with an unknown/non-finite prior. + if not (np.all(np.isfinite(prior_strength)) and + np.all(np.isfinite(prior_flux[prior_strength != 0]))): + raise ValueError('fit_scene_lightcurve: prior_flux/prior_strength must be finite ' + 'wherever prior_strength is nonzero.') + ATA = ATA + np.diag(prior_strength) + # prior_flux may be non-finite wherever prior_strength is exactly 0 + # (no prior for that column) -- np.where avoids 0*nan=nan poisoning + # the RHS despite the zero weight. + pull = np.where(prior_strength != 0, prior_strength * prior_flux, 0.0) + rhs = Ag.T @ Dg.T + pull[:, None] + else: + rhs = Ag.T @ Dg.T + + ATA_inv = np.linalg.inv(ATA) + coeffs = ATA_inv @ rhs # (M, n_frames) + + flux = coeffs[0] + resid = Dg.T - Ag @ coeffs + med = np.median(resid, axis=0) + sigma = 1.4826 * np.median(np.abs(resid - med), axis=0) + e_flux = sigma * np.sqrt(ATA_inv[0, 0]) + + return flux, e_flux, coeffs diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 3015ac9..b544c9e 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -34,20 +34,20 @@ def _available_cores(verbose=1): if slurm is not None: try: n = int(slurm) - if verbose > 0: + if verbose >= 3: print(f'[tessreduce] _available_cores: SLURM_CPUS_PER_TASK={slurm} → {n} cores') return n except ValueError: pass try: n = len(os.sched_getaffinity(0)) - if verbose > 0: + if verbose >= 3: print(f'[tessreduce] _available_cores: sched_getaffinity → {n} cores') return n except AttributeError: pass n = multiprocessing.cpu_count() - if verbose > 0: + if verbose >= 3: print(f'[tessreduce] _available_cores: cpu_count fallback → {n} cores') return n @@ -211,8 +211,9 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect prf_path : str, optional Path to local TESS PRF files. The default is currently a specific location on the OzStar supercomputer. verbose : int, optional - Controls the level of verbosity. 0 is silent, 1 (default) prints reduction stage - announcements, 2 additionally prints joblib's per-task Parallel output. + Controls the level of verbosity. 0 is silent, 1 (default) prints top-level reduction + stage announcements, 2 additionally prints background() sub-step announcements (with + per-step timing), 3 additionally prints joblib's per-task Parallel output. timing : bool, optional Print execution time reports for major pipeline blocks in background() and reduce(). The default is False. @@ -259,7 +260,7 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect self._cache_path = None # SLURM environment diagnostics - if verbose > 0: + if verbose >= 3: _slurm_vars = ['SLURM_CPUS_PER_TASK', 'SLURM_NTASKS', 'SLURM_NTASKS_PER_NODE', 'SLURM_JOB_CPUS_PER_NODE', 'SLURM_CPUS_ON_NODE'] _slurm_env = {k: os.environ.get(k, 'not set') for k in _slurm_vars} @@ -371,8 +372,8 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect @property def _joblib_verbose(self): - """joblib Parallel verbosity: 0 unless self.verbose requests joblib output (>=2).""" - return 1 if self.verbose >= 2 else 0 + """joblib Parallel verbosity: 0 unless self.verbose requests joblib output (>=3).""" + return 1 if self.verbose >= 3 else 0 def check_coord(self): """ @@ -804,12 +805,12 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa bkg_smth = np.zeros_like(flux) * np.nan if self.parallel: _t = time.perf_counter() - if self.verbose > 0: + if self.verbose >= 2: print('smooth background...') bkg_smth = Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*m) _times['initial smooth background'] = time.perf_counter() - _t - if self.verbose > 0: - print('smooth background done') + if self.verbose >= 2: + print(f'smooth background done ({_times["initial smooth background"]:.2f}s)') if rerun_negative: _t = time.perf_counter() if self._use_error_image: @@ -830,12 +831,12 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa else: m[over_sub] = 1 self._bkgmask = m - if self.verbose > 0: + if self.verbose >= 2: print('smooth background rerun (negative correction)...') bkg_smth = Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(Smooth_bkg)(frame,gauss_smooth,interpolate) for frame in flux*m) _times['negative over-subtraction rerun'] = time.perf_counter() - _t - if self.verbose > 0: - print('smooth background rerun done') + if self.verbose >= 2: + print(f'smooth background rerun done ({_times["negative over-subtraction rerun"]:.2f}s)') if rerun_diff: _t = time.perf_counter() @@ -890,12 +891,14 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa new_mask = abs(new_mask - 1) self._bkgmask = new_mask bkg_s1 = np.array(bkg_smth) - if self.verbose > 0: + if self.verbose >= 2: print('smooth background rerun (residual surface)...') bkg_smth = Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*new_mask) if blend_dynamic: bkg_smth = blend_dynamic_background(bkg_smth, bkg_s1, flux, n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose) _times['residual surface rerun'] = time.perf_counter() - _t + if self.verbose >= 2: + print(f'smooth background rerun (residual surface) done ({_times["residual surface rerun"]:.2f}s)') else: _t = time.perf_counter() @@ -903,8 +906,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa bkg_smth[i] = Smooth_bkg((flux*m)[i],0,interpolate) _times['initial smooth background'] = time.perf_counter() - _t else: - if self.verbose > 0: - print('Small tpf, using percentile cut background') + print('Small tpf, using percentile cut background') self.small_background() bkg_smth = self.bkg @@ -925,7 +927,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa _high_bkg_frames = (_earth_angle < 30.0) | (_moon_angle < 30.0) | (_bkg_median > 300.0) else: _high_bkg_frames = _bkg_median > 300.0 - if self.verbose > 0: + if self.verbose >= 2: print('fixing background anomalies...') self.bkg, _sharp_masks, self.bad_bkg = fix_background_anomalies(self.bkg, self.mask, flux=deepcopy(self.flux), @@ -935,11 +937,11 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa high_bkg_frames=_high_bkg_frames, n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose) _times['anomaly fixing'] = time.perf_counter() - _t - if self.verbose > 0: - print('fixing background anomalies done') + if self.verbose >= 2: + print(f'fixing background anomalies done ({_times["anomaly fixing"]:.2f}s)') _t = time.perf_counter() - if self.verbose > 0: + if self.verbose >= 2: print('adaptive temporal smoothing...') from .adaptive_background import AdaptiveBackground smoother = AdaptiveBackground(self.bkg, self.mjd, sector=self.sector, camera=self.tpf.camera, @@ -948,8 +950,8 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa smoothed = smoother.smooth(method='savgol').smoothed self.bkg = smoothed _times['adaptive temporal smoothing'] = time.perf_counter() - _t - if self.verbose > 0: - print('adaptive temporal smoothing done') + if self.verbose >= 2: + print(f'adaptive temporal smoothing done ({_times["adaptive temporal smoothing"]:.2f}s)') # Store data-driven sources from _bkgmask as bit 8, preserving the catalogue mask (bit 1) if rerun_diff: @@ -2564,6 +2566,158 @@ def psf_photometry(self,xPix,yPix,size=7,snap='brightest',ext_shift=True,plot=Fa return flux, eflux + def scene_photometry(self,xPix,yPix,size=9,bkg_poly_order=2,snap='ref', + neighbours=True,neighbour_prior_mag_tol=0.5, + include_psf_derivatives=True,diff=None): + """ + Scene-modelling PSF photometry. + + Fits a fixed (or once-refined) position with a single vectorized linear + solve across the *entire* time series, using a 2D polynomial background + surface (instead of a flat constant) and optional PSF-derivative columns + that absorb poor difference-image subtraction residuals right at the + target position. This is an additional, much faster alternative to + psf_photometry() for the fixed/near-fixed position case -- it does not + replace psf_photometry(), whose per-frame position refit (snap='all') + still covers genuine per-frame motion that a fixed-basis linear model + cannot represent. + + Parameters + ---------- + xPix, yPix : float + Pixel location of the target, in the same pixel frame used + elsewhere (e.g. psf_photometry, field_calibrate's catalogue col/row). + size : int, optional + Cutout side length in pixels (should be odd). The default is 9. + bkg_poly_order : int, optional + Order of the 2D polynomial background surface. 0 reproduces a flat + background. The default is 2. + snap : {'ref','brightest','fixed'}, optional + How the fixed position is chosen. 'ref' (default) refines position + once on the reference frame; 'brightest' refines on the highest + flux-weighted frame in the cutout; 'fixed' skips position + refinement entirely (uses (xPix,yPix) exactly). + neighbours : bool, optional + If True, nearby sources from the calibration catalogue (self.cat, + populated by field_calibrate()/field_calibrate_scene()) within the + cutout are modelled as extra PSF columns for deblending, with a + ridge prior pulling their fitted flux toward the catalogue- + predicted value (see neighbour_prior_mag_tol). If no calibration + has been run yet, neighbours are modelled unconstrained and a + warning is raised. The default is True. + neighbour_prior_mag_tol : float, optional + Tolerance (mag) controlling how strongly a neighbour's fitted flux + is pulled toward its catalogue-predicted value: smaller values pull + harder. The default is 0.5. + include_psf_derivatives : bool, optional + Add PSF spatial-derivative columns at the target position to absorb + sharp subtraction-residual structure from imperfect difference-image + alignment/kernel-matching. The default is True. + diff : bool, optional + If True use the already-differenced flux cube; if False add self.ref + back before extracting the cutout. Defaults to self.diff. + + Returns + ------- + flux, eflux : array_like + Flux and formal per-frame flux error, length = n frames. + """ + from . import scene_photom as sp + + if diff is None: + diff = self.diff + + xPix = int(np.round(xPix)) + yPix = int(np.round(yPix)) + + col = self.tpf.column - int(self.size//2) + xPix + 45 + row = self.tpf.row - int(self.size//2) + yPix + 1 + col = int(np.clip(col,45,2090)) + row = int(np.clip(row,10,2040)) + prf = sp._prf_cache(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row,self._prf_path) + + half = size // 2 + if diff: + cutout = self.flux[:,yPix-half:yPix+half+1,xPix-half:xPix+half+1] + else: + cutout = (self.flux+self.ref)[:,yPix-half:yPix+half+1,xPix-half:xPix+half+1] + + cent = (size-1) / 2.0 + + # -- position: fixed, or refined once on a single high-SNR frame -- # + if snap == 'ref': + ref_ind = self.ref_ind + elif snap == 'brightest': + weight = np.abs(np.nansum(cutout[:,half-1:half+2,half-1:half+2],axis=(1,2))) + weight[~np.isfinite(weight)] = 0 + ref_ind = int(np.argmax(weight)) + elif snap == 'fixed': + ref_ind = None + else: + raise ValueError("snap must be one of 'ref','brightest','fixed' for scene_photometry") + + if ref_ind is not None: + def build_A(dx,dy): + A, _ = sp.build_design_matrix(prf,cent,(dx,dy),[],size, + poly_order=bkg_poly_order, + include_psf_derivatives=include_psf_derivatives) + return A + target_dxdy = sp.fit_scene_position(build_A,cutout[ref_ind].ravel(),tol_x=1.0,tol_y=1.0) + else: + target_dxdy = (0.0,0.0) + + # -- neighbours, informed by the calibration catalogue -- # + neighbour_dxdys = [] + expected_flux = np.array([]) + if neighbours: + cat = getattr(self,'cat',None) + zp = getattr(self,'zp',None) + if (cat is not None) and (zp is not None): + mag_col = 'tmag' if 'tmag' in cat.columns else ('rp_ab' if 'rp_ab' in cat.columns else None) + if mag_col is not None: + dcol = cat['col'].values - xPix + drow = cat['row'].values - yPix + in_stamp = (np.abs(dcol) <= half) & (np.abs(drow) <= half) + is_target = np.hypot(dcol,drow) <= 0.5 + sel = in_stamp & ~is_target + neighbour_dxdys = list(zip(dcol[sel],drow[sel])) + if len(neighbour_dxdys) > 0: + if diff: + # on a differenced cube a static neighbour's true + # amplitude is ~0 (it subtracts out against the + # reference), not its catalogue absolute flux -- + # pulling toward the absolute flux here would bias + # the fit. Prior toward zero residual instead. + expected_flux = np.zeros(sel.sum()) + else: + expected_flux = 10**(-0.4*(cat[mag_col].values[sel]-zp)) + else: + warnings.warn('scene_photometry: neighbours=True but no calibration catalogue/zeropoint ' + 'is set yet (run field_calibrate() or field_calibrate_scene() first); ' + 'proceeding with neighbours unconstrained.') + + A, info = sp.build_design_matrix(prf,cent,target_dxdy,neighbour_dxdys,size, + poly_order=bkg_poly_order, + include_psf_derivatives=include_psf_derivatives) + + prior_flux = None + prior_strength = None + if len(neighbour_dxdys) > 0 and len(expected_flux) == len(neighbour_dxdys): + prior_flux = np.zeros(A.shape[1]) + prior_strength = np.zeros(A.shape[1]) + for k,idx in enumerate(info['neighbours']): + # a non-finite catalogue magnitude must not poison the shared + # closed-form solve (one NaN row contaminates every column via + # the matrix inverse) -- leave that neighbour unconstrained + # (still deblended via its own column, just with no prior). + if np.isfinite(expected_flux[k]): + prior_flux[idx] = expected_flux[k] + prior_strength[idx] = 1.0 / (neighbour_prior_mag_tol**2) + + flux, eflux, coeffs = sp.fit_scene_lightcurve(cutout,A,prior_flux=prior_flux,prior_strength=prior_strength) + return flux, eflux + + def orbit_ref_subtract(self): """ Subtract a per-orbit reference from each orbit's frames. @@ -2689,7 +2843,7 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, self.get_ref(ref_start,ref_stop) _times['reference frame'] = time.perf_counter() - _t if self.verbose > 0: - print('made reference') + print(f'made reference ({_times["reference frame"]:.2f}s)') # make source mask _t = time.perf_counter() if mask is None: @@ -2699,30 +2853,32 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, if frac < 0.05: print('!!!WARNING!!! mask is too dense, lowering mask_scale to 0.5, and raising maglim to 15. Background quality will be reduced.') self.make_mask(catalogue_path=self._catalogue_path,maglim=15,strapsize=7,scale=0.2) - if self.verbose > 0: - print('made source mask') + _mask_msg = 'made source mask' else: self.mask = mask - if self.verbose > 0: - print('assigned source mask') + _mask_msg = 'assigned source mask' if moving_mask is not None: moving_mask = moving_mask > 0 temp = np.zeros_like(self.flux,dtype=int) temp[:,:,:] = self.mask self.mask = temp | moving_mask _times['source mask'] = time.perf_counter() - _t + if self.verbose > 0: + print(f'{_mask_msg} ({_times["source mask"]:.2f}s)') # calculate background for each frame if self.verbose > 0: print('calculating background') # calculate the background #self.flux -= self.ref - print('background pass 1...') + if self.verbose > 0: + print('background pass 1...') _t = time.perf_counter() self.background(rerun_negative=True) self.flux -= self.bkg _times['background (pass 1)'] = time.perf_counter() - _t - print('background pass 1 done') + if self.verbose > 0: + print(f'background pass 1 done ({_times["background (pass 1)"]:.2f}s)') if np.isnan(self.bkg).all(): raise ValueError('bkg all nans') @@ -2777,7 +2933,8 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, else: self.shift = np.zeros((len(self.flux),2)) _times['alignment'] = time.perf_counter() - _t - print('alignment done') + if self.verbose > 0: + print(f'alignment done ({_times["alignment"]:.2f}s)') if not self.diff: if self.align: @@ -2786,7 +2943,7 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, self.flux[np.nansum(self.tpf.flux.value,axis=(1,2))==0] = np.nan _times['image shifting'] = time.perf_counter() - _t if self.verbose > 0: - print('images shifted') + print(f'images shifted ({_times["image shifting"]:.2f}s)') #if self.kernel_match: # self.kernel_matching(diff=False) @@ -2838,24 +2995,29 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, temp[:,:,:] = self.mask self.mask = temp | moving_mask _times['difference imaging setup'] = time.perf_counter() - _t - print('diff setup done') + if self.verbose > 0: + print(f'diff setup done ({_times["difference imaging setup"]:.2f}s)') self.bkg_orig = deepcopy(self.bkg) - print('background pass 2...') + if self.verbose > 0: + print('background pass 2...') _t = time.perf_counter() self.background(calc_qe = False,strap_iso = False,source_hunt=self._sourcehunt, gauss_smooth=self._bkg_gauss_sigma,interpolate=False, rerun_negative=False,rerun_diff=True,blend_dynamic=True) self.flux -= self.bkg _times['background (pass 2)'] = time.perf_counter() - _t - print('background pass 2 done') + if self.verbose > 0: + print(f'background pass 2 done ({_times["background (pass 2)"]:.2f}s)') if self.corr_correction: _t = time.perf_counter() - print('correlation correction...') + if self.verbose > 0: + print('correlation correction...') self.correlation_corrector() _times['correlation correction'] = time.perf_counter() - _t - print('correlation correction done') + if self.verbose > 0: + print(f'correlation correction done ({_times["correlation correction"]:.2f}s)') if self.kernel_match: self.kernel_matching(diff=self.diff) if self.verbose > 0: @@ -2868,14 +3030,19 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, if self.calibrate: _t = time.perf_counter() - print('field calibration') + if self.verbose > 0: + print('field calibration...') self.field_calibrate() _times['field calibration'] = time.perf_counter() - _t + if self.verbose > 0: + print(f'field calibration done ({_times["field calibration"]:.2f}s)') if self._create_lc: _t = time.perf_counter() self.lc, self.sky = self.diff_lc(plot=self.plot,diff=self.diff,tar_ap=tar_ap,sky_in=sky_in,sky_out=sky_out) _times['light curve'] = time.perf_counter() - _t + if self.verbose > 0: + print(f'light curve extracted ({_times["light curve"]:.2f}s)') if self.imaging: # if self.verbose > 0: @@ -3778,6 +3945,67 @@ def field_calibrate(self,zp_single=True,plot=None,savename=None): return + def field_calibrate_scene(self,catalog='ps1',mag_lo=None,mag_hi=None,plot=False,**kwargs): + """ + Scene-modelling photometric calibration. + + Additive alternative to field_calibrate(): both pathways use the shared + scene-fit engine in scene_photom.py (real pixel-space deblending, formal + per-star error propagation, magnitude-binned robust zeropoint combining) + instead of field_calibrate()'s catalog-magnitude-space crowding correction + and coarse sanity-check gating. Sets self.zp/self.zp_e/self.tzp/self.tzp_e + exactly as field_calibrate() does, so diff_lc/light-curve unit conversion + works unchanged with either pathway. Does not touch field_calibrate(). + + Parameters + ---------- + catalog : {'ps1','gaia'}, optional + 'ps1' (default) uses PS1 (dec>-30) or SkyMapper (dec<=-30), with the + existing Tonry-locus extinction correction and multi-band synthetic- + TESS-mag reconstruction -- both already AB. 'gaia' uses Gaia DR3 + (queried via astroquery, matching tessreduce's existing catalogue + convention), calibrated directly against Rp with the standard + Vega->AB offset applied; no extinction step for this pathway. + mag_lo, mag_hi : float, optional + Magnitude range for the calibration sample. Defaults are pathway- + specific (8.5-16.0 for PS1/SkyMapper tmag, 11.0-15.5 for Gaia rp_ab). + plot : bool, optional + Passed through to the extinction fit (PS1/SkyMapper pathway only). + **kwargs + Forwarded to field_calibration._calibrate_common (e.g. + iso_radius_pix, stamp_size, poly_order, refine_iter, max_zp_err). + + On failure (too few good stars, fits don't converge) a warning is + raised and self.zp/self.tzp fall back to 20.6 (self.zp_e/self.tzp_e=0.1), + rather than raising an exception. + + Assigns + ------- + self.zp/tzp : float + AB photometric zeropoint. + self.zp_e/tzp_e : float + Error in the photometric zeropoint. + self.cat : pandas.DataFrame + Calibration catalogue with pixel position (col,row) and either + tmag (PS1/SkyMapper pathway) or rp_ab (Gaia pathway). + """ + from . import field_calibration as fc + + if catalog == 'ps1': + zp_ab, zp_err = fc.calibrate_ps1_skymapper(self,mag_lo=mag_lo,mag_hi=mag_hi,plot=plot,**kwargs) + elif catalog == 'gaia': + mag_lo = 11.0 if mag_lo is None else mag_lo + mag_hi = 15.5 if mag_hi is None else mag_hi + zp_ab, zp_err = fc.calibrate_gaia(self,mag_lo=mag_lo,mag_hi=mag_hi,plot=plot,**kwargs) + else: + raise ValueError(f"catalog must be 'ps1' or 'gaia', got {catalog!r}") + + self.zp = zp_ab + self.zp_e = zp_err + self.tzp = zp_ab + self.tzp_e = zp_err + return + def to_mag(self,zp=None,zp_e=0): """ Convert the TESS lc into magnitude space. diff --git a/tests/test_scene_photom.py b/tests/test_scene_photom.py new file mode 100644 index 0000000..e1c9be4 --- /dev/null +++ b/tests/test_scene_photom.py @@ -0,0 +1,369 @@ +import unittest +import numpy as np +from numpy.testing import assert_allclose + +import matplotlib +matplotlib.use('Agg') + +from tessreduce.scene_photom import ( + polynomial_columns, + prf_column, + prf_derivative_columns, + build_design_matrix, + fit_scene_position, + fit_scene_frame, + fit_scene_lightcurve, +) + +STAMP_SIZE = 9 +CENT = (STAMP_SIZE - 1) / 2.0 + + +class _GaussianPRF: + """Minimal stand-in for PRF.TESS_PRF with a .locate(x, y, shape) method.""" + + def __init__(self, sigma=1.2, skew=0.0): + self.sigma = sigma + self.skew = skew + + def locate(self, x, y, shape): + ny, nx = shape + yy, xx = np.mgrid[0:ny, 0:nx] + r2 = (xx - x) ** 2 + (yy - y) ** 2 + g = np.exp(-(r2 / (2 * self.sigma ** 2))) + if self.skew: + g = g * (1 + self.skew * (xx - x)) + return np.clip(g, 0, None) + + +class TestPolynomialColumns(unittest.TestCase): + + def test_order_zero_is_flat(self): + cols = polynomial_columns(STAMP_SIZE, order=0) + self.assertEqual(len(cols), 1) + assert_allclose(cols[0], np.ones(STAMP_SIZE * STAMP_SIZE)) + + def test_term_count_matches_triangular_scheme(self): + for order in range(4): + cols = polynomial_columns(STAMP_SIZE, order=order) + expected = (order + 1) * (order + 2) // 2 + self.assertEqual(len(cols), expected) + + def test_columns_are_frame_independent_shape(self): + cols = polynomial_columns(STAMP_SIZE, order=2) + for c in cols: + self.assertEqual(c.shape, (STAMP_SIZE * STAMP_SIZE,)) + + +class TestDesignMatrix(unittest.TestCase): + + def test_shape_without_derivatives_or_neighbours(self): + prf = _GaussianPRF() + A, info = build_design_matrix(prf, CENT, (0.0, 0.0), [], STAMP_SIZE, + poly_order=0, include_psf_derivatives=False) + self.assertEqual(A.shape, (STAMP_SIZE * STAMP_SIZE, 2)) # target + flat bg + self.assertEqual(info['target'], 0) + self.assertEqual(info['background'], [1]) + self.assertEqual(info['neighbours'], []) + + def test_shape_with_derivatives(self): + prf = _GaussianPRF() + A, info = build_design_matrix(prf, CENT, (0.0, 0.0), [], STAMP_SIZE, + poly_order=0, include_psf_derivatives=True) + self.assertEqual(A.shape, (STAMP_SIZE * STAMP_SIZE, 4)) # target + 2 deriv + bg + self.assertIn('target', info['derivatives']) + + def test_shape_with_neighbours(self): + prf = _GaussianPRF() + A, info = build_design_matrix(prf, CENT, (0.0, 0.0), [(2.0, 0.0), (-1.5, 1.0)], + STAMP_SIZE, poly_order=1, include_psf_derivatives=False) + # target(1) + 2 neighbours + poly_order=1 -> 3 bg terms + self.assertEqual(A.shape, (STAMP_SIZE * STAMP_SIZE, 1 + 2 + 3)) + self.assertEqual(len(info['neighbours']), 2) + + def test_symmetric_psf_orthogonal_to_its_own_derivative(self): + """A symmetric PSF is even; its spatial derivative is odd. Their inner + product over a symmetric stamp should vanish -- this is why a pure + dipole subtraction-residual can't bias flux regardless of whether the + derivative columns are modelled (see TestDifferenceResidualHandling).""" + prf = _GaussianPRF(skew=0.0) + target = prf_column(prf, CENT, 0.0, 0.0, STAMP_SIZE) + dpdx, dpdy = prf_derivative_columns(prf, CENT, 0.0, 0.0, STAMP_SIZE) + self.assertAlmostEqual(np.dot(target, dpdx), 0.0, places=8) + self.assertAlmostEqual(np.dot(target, dpdy), 0.0, places=8) + + +class TestFitSceneLightcurve(unittest.TestCase): + + def setUp(self): + self.rng = np.random.default_rng(0) + self.prf = _GaussianPRF() + self.nfr = 60 + + def _make_cube(self, A, col_fluxes, noise_sigma=2.0): + """col_fluxes: dict {col_index: (n_frames,) array} of per-frame amplitudes.""" + nfr = self.nfr + data = np.zeros((nfr, A.shape[0])) + for idx, series in col_fluxes.items(): + data += np.outer(series, A[:, idx]) + data += self.rng.normal(0, noise_sigma, size=data.shape) + return data.reshape(nfr, STAMP_SIZE, STAMP_SIZE) + + def test_recovers_flux_and_background(self): + A, info = build_design_matrix(self.prf, CENT, (0.0, 0.0), [], STAMP_SIZE, + poly_order=0, include_psf_derivatives=False) + true_flux = 1000 + 200 * np.sin(np.linspace(0, 6, self.nfr)) + true_bg = 50 + 5 * np.cos(np.linspace(0, 3, self.nfr)) + cube = self._make_cube(A, {0: true_flux, info['background'][0]: true_bg}) + + flux, e_flux, coeffs = fit_scene_lightcurve(cube, A) + assert_allclose(flux, true_flux, atol=5 * np.median(e_flux)) + assert_allclose(coeffs[info['background'][0]], true_bg, atol=5 * np.median(e_flux)) + + def test_poly_order_zero_matches_flat_background(self): + """order=0 should be numerically identical to a single flat column.""" + A0, _ = build_design_matrix(self.prf, CENT, (0.0, 0.0), [], STAMP_SIZE, + poly_order=0, include_psf_derivatives=False) + cols = polynomial_columns(STAMP_SIZE, order=0) + assert_allclose(A0[:, 1], cols[0]) + + def test_error_scales_with_injected_noise(self): + A, info = build_design_matrix(self.prf, CENT, (0.0, 0.0), [], STAMP_SIZE, + poly_order=0, include_psf_derivatives=False) + true_flux = np.full(self.nfr, 1000.0) + + cube_lo = self._make_cube(A, {0: true_flux}, noise_sigma=1.0) + cube_hi = self._make_cube(A, {0: true_flux}, noise_sigma=8.0) + + _, e_lo, _ = fit_scene_lightcurve(cube_lo, A) + _, e_hi, _ = fit_scene_lightcurve(cube_hi, A) + self.assertLess(np.median(e_lo), np.median(e_hi)) + + +class TestDifferenceResidualHandling(unittest.TestCase): + """Verify the PSF-derivative columns absorb poor-subtraction residuals + near the target, per the plan's difference-imaging requirement.""" + + def setUp(self): + self.rng = np.random.default_rng(1) + self.nfr = 60 + + def test_derivative_columns_reduce_noise_inflation(self): + # slightly asymmetric PRF, closer to a real TESS PRF than a pure Gaussian + prf = _GaussianPRF(skew=0.08) + true_flux = 1000 + 100 * np.sin(np.linspace(0, 6, self.nfr)) + true_bg = 50.0 + shift_amp = 40.0 + + A_no, info_no = build_design_matrix(prf, CENT, (0.0, 0.0), [], STAMP_SIZE, + poly_order=0, include_psf_derivatives=False) + A_wd, info_wd = build_design_matrix(prf, CENT, (0.0, 0.0), [], STAMP_SIZE, + poly_order=0, include_psf_derivatives=True) + dpdx, _ = prf_derivative_columns(prf, CENT, 0.0, 0.0, STAMP_SIZE) + + target_col = A_no[:, 0] + bg_col = A_no[:, info_no['background'][0]] + + data = (np.outer(true_flux, target_col) + true_bg * bg_col + + shift_amp * dpdx) + data += self.rng.normal(0, 2.0, size=data.shape) + cube = data.reshape(self.nfr, STAMP_SIZE, STAMP_SIZE) + + flux_no, e_no, _ = fit_scene_lightcurve(cube, A_no) + flux_wd, e_wd, coeffs_wd = fit_scene_lightcurve(cube, A_wd) + + # modelling the residual should not inflate (and should typically + # reduce) the robust per-frame noise estimate + self.assertLessEqual(np.median(e_wd), np.median(e_no) * 1.05) + + # the fitted derivative coefficient should recover the injected + # residual amplitude -- this is the diagnostic value, not just noise + # absorption + deriv_idx = info_wd['derivatives']['target'][0] + recovered = np.median(coeffs_wd[deriv_idx]) + self.assertAlmostEqual(recovered, shift_amp, delta=0.3 * shift_amp) + + def test_condition_number_stable_with_derivatives(self): + prf = _GaussianPRF() + A, _ = build_design_matrix(prf, CENT, (0.0, 0.0), [], STAMP_SIZE, + poly_order=2, include_psf_derivatives=True) + cond = np.linalg.cond(A.T @ A) + self.assertLess(cond, 1e8) + + +class TestCatalogInformedCrowding(unittest.TestCase): + """Ridge-prior handling of a tightly-blended neighbour, informed by a + catalog-predicted flux -- keeps the vectorized closed-form solve while + stabilizing near-degenerate PRF columns.""" + + def setUp(self): + self.rng = np.random.default_rng(2) + self.prf = _GaussianPRF() + self.nfr = 30 + + def _blended_setup(self): + neighbour_dxdy = [(0.15, 0.1)] # extremely close -> near-degenerate + A, info = build_design_matrix(self.prf, CENT, (0.0, 0.0), neighbour_dxdy, + STAMP_SIZE, poly_order=0, include_psf_derivatives=False) + return A, info + + def test_unconstrained_blend_is_ill_conditioned(self): + A, _ = self._blended_setup() + cond = np.linalg.cond(A.T @ A) + self.assertGreater(cond, 1e4) + + def test_nonfinite_prior_raises_instead_of_poisoning_solve(self): + """A non-finite catalog magnitude must not silently NaN every column + of the solve via the matrix inverse -- it should raise instead. + Regression test for a real bug found on live TESS data: one NaN + catalog mag among the selected neighbours produced all-NaN flux for + every frame and every snap mode, with no error raised.""" + A, info = self._blended_setup() + cube = self.rng.normal(1000, 5, size=(self.nfr, STAMP_SIZE, STAMP_SIZE)) + + prior_flux = np.zeros(A.shape[1]) + prior_strength = np.zeros(A.shape[1]) + prior_flux[info['neighbours'][0]] = np.nan # e.g. missing catalog mag + prior_strength[info['neighbours'][0]] = 50.0 + + with self.assertRaises(ValueError): + fit_scene_lightcurve(cube, A, prior_flux=prior_flux, prior_strength=prior_strength) + + def test_zero_strength_ignores_nonfinite_prior_flux(self): + """If a column's prior_strength is exactly 0, a non-finite prior_flux + for that column must be tolerated (it contributes nothing).""" + A, info = self._blended_setup() + cube = self.rng.normal(1000, 5, size=(self.nfr, STAMP_SIZE, STAMP_SIZE)) + + prior_flux = np.zeros(A.shape[1]) + prior_strength = np.zeros(A.shape[1]) + prior_flux[info['neighbours'][0]] = np.nan + prior_strength[info['neighbours'][0]] = 0.0 + + flux, e_flux, coeffs = fit_scene_lightcurve(cube, A, prior_flux=prior_flux, + prior_strength=prior_strength) + self.assertTrue(np.all(np.isfinite(flux))) + + def test_ridge_prior_stabilizes_target_flux(self): + A, info = self._blended_setup() + true_target_flux = np.full(self.nfr, 1000.0) + true_neighbour_flux = 800.0 + true_bg = 50.0 + + target_col = A[:, 0] + neighbour_col = A[:, info['neighbours'][0]] + bg_col = A[:, info['background'][0]] + + data = (np.outer(true_target_flux, target_col) + true_neighbour_flux * neighbour_col + + true_bg * bg_col) + data += self.rng.normal(0, 3.0, size=data.shape) + cube = data.reshape(self.nfr, STAMP_SIZE, STAMP_SIZE) + + flux_unc, e_unc, _ = fit_scene_lightcurve(cube, A) + + prior_flux = np.zeros(A.shape[1]) + prior_strength = np.zeros(A.shape[1]) + prior_flux[info['neighbours'][0]] = true_neighbour_flux + prior_strength[info['neighbours'][0]] = 50.0 + + flux_prior, e_prior, _ = fit_scene_lightcurve(cube, A, prior_flux=prior_flux, + prior_strength=prior_strength) + + # ridge prior should not make the fit less stable than unconstrained + self.assertLessEqual(np.std(flux_prior), np.std(flux_unc) * 1.2) + # and should keep the target flux closer to truth on average + self.assertLess(abs(np.median(flux_prior) - 1000.0), + abs(np.median(flux_unc) - 1000.0) + 50.0) + + def test_prior_negligible_when_weak(self): + """A very weak prior_strength shouldn't meaningfully drag a + well-separated neighbour's flux away from its unconstrained fit.""" + prf = self.prf + A, info = build_design_matrix(prf, CENT, (0.0, 0.0), [(4.0, 0.0)], STAMP_SIZE, + poly_order=0, include_psf_derivatives=False) + target_col = A[:, 0] + neighbour_col = A[:, info['neighbours'][0]] + bg_col = A[:, info['background'][0]] + + true_neighbour_flux = 800.0 + data = (np.outer(np.full(self.nfr, 1000.0), target_col) + + true_neighbour_flux * neighbour_col + 50.0 * bg_col) + data += self.rng.normal(0, 1.0, size=data.shape) + cube = data.reshape(self.nfr, STAMP_SIZE, STAMP_SIZE) + + flux_unc, _, coeffs_unc = fit_scene_lightcurve(cube, A) + + prior_flux = np.zeros(A.shape[1]) + prior_strength = np.zeros(A.shape[1]) + # deliberately "wrong" prior, but with ~zero strength + prior_flux[info['neighbours'][0]] = 200.0 + prior_strength[info['neighbours'][0]] = 1e-6 + + flux_weak, _, coeffs_weak = fit_scene_lightcurve(cube, A, prior_flux=prior_flux, + prior_strength=prior_strength) + assert_allclose(coeffs_unc[info['neighbours'][0]], coeffs_weak[info['neighbours'][0]], + rtol=1e-2) + + +class TestFitScenePosition(unittest.TestCase): + + def test_recovers_known_subpixel_shift(self): + prf = _GaussianPRF() + true_dx, true_dy = 0.3, -0.2 + target_col = prf_column(prf, CENT, true_dx, true_dy, STAMP_SIZE) + bg_col = np.ones(STAMP_SIZE * STAMP_SIZE) + + rng = np.random.default_rng(3) + data = target_col * 1000 + bg_col * 50 + rng.normal(0, 1.0, target_col.shape) + + def build_A(dx, dy): + A, _ = build_design_matrix(prf, CENT, (dx, dy), [], STAMP_SIZE, + poly_order=0, include_psf_derivatives=False) + return A + + dx_fit, dy_fit = fit_scene_position(build_A, data, tol_x=1.0, tol_y=1.0) + self.assertAlmostEqual(dx_fit, true_dx, delta=0.15) + self.assertAlmostEqual(dy_fit, true_dy, delta=0.15) + + def test_zero_tolerance_returns_fixed_position(self): + prf = _GaussianPRF() + + def build_A(dx, dy): + A, _ = build_design_matrix(prf, CENT, (dx, dy), [], STAMP_SIZE, + poly_order=0, include_psf_derivatives=False) + return A + + dx_fit, dy_fit = fit_scene_position(build_A, np.ones(STAMP_SIZE * STAMP_SIZE), + tol_x=0.0, tol_y=0.0) + self.assertEqual((dx_fit, dy_fit), (0.0, 0.0)) + + +class TestFitSceneFrame(unittest.TestCase): + + def test_bounded_solve_recovers_flux(self): + prf = _GaussianPRF() + A, _ = build_design_matrix(prf, CENT, (0.0, 0.0), [], STAMP_SIZE, + poly_order=0, include_psf_derivatives=False) + rng = np.random.default_rng(4) + data = A[:, 0] * 500 + A[:, 1] * 20 + rng.normal(0, 1.0, A.shape[0]) + + bounds = ([0, -np.inf], [np.inf, np.inf]) + coeffs, cov, dof = fit_scene_frame(A, data, flux_bounds=bounds) + self.assertAlmostEqual(coeffs[0], 500, delta=20) + self.assertGreater(dof, 0) + self.assertEqual(cov.shape, (A.shape[1], A.shape[1])) + + def test_unbounded_solve_matches_lstsq(self): + prf = _GaussianPRF() + A, _ = build_design_matrix(prf, CENT, (0.0, 0.0), [], STAMP_SIZE, + poly_order=0, include_psf_derivatives=False) + rng = np.random.default_rng(5) + data = A[:, 0] * 300 + A[:, 1] * 10 + rng.normal(0, 0.5, A.shape[0]) + + coeffs, cov, dof = fit_scene_frame(A, data, flux_bounds=None) + expected, *_ = np.linalg.lstsq(A, data, rcond=None) + assert_allclose(coeffs, expected, atol=1e-8) + + +if __name__ == '__main__': + unittest.main() From a779e1a335494c8d628d9e0eab9c55f7de90073f Mon Sep 17 00:00:00 2001 From: Hugh Roxburgh Date: Mon, 3 Aug 2026 14:00:25 +0800 Subject: [PATCH 20/22] remove tpf dependency --- tessreduce/cat_mask.py | 10 +-- tessreduce/catalog_tools.py | 42 +++++------ tessreduce/tessreduce.py | 142 ++++++++++++++++++++++++------------ 3 files changed, 123 insertions(+), 71 deletions(-) diff --git a/tessreduce/cat_mask.py b/tessreduce/cat_mask.py index a4050f9..b908ce8 100755 --- a/tessreduce/cat_mask.py +++ b/tessreduce/cat_mask.py @@ -296,7 +296,7 @@ def Strap_mask(Image, col, size=4, flux_cube=None): big_strap = fftconvolve(strap_mask, np.ones((size, size)), mode='same') > .5 return big_strap -def Cat_mask(tpf,catalogue_path=None,maglim=19,scale=1,strapsize=3,ref=None,sigma=3,col_offset=0,flux_cube=None): +def Cat_mask(ra,dec,shape,wcs,flux,column,catalogue_path=None,maglim=19,scale=1,strapsize=3,ref=None,sigma=3,col_offset=0,flux_cube=None): """ Make a source mask from the PS1 and Gaia catalogs. @@ -330,14 +330,14 @@ def Cat_mask(tpf,catalogue_path=None,maglim=19,scale=1,strapsize=3,ref=None,sigm if catalogue_path is not None: gaia = external_load_cat(catalogue_path,maglim) - coords = tpf.wcs.all_world2pix(gaia['ra'],gaia['dec'], 0) + coords = wcs.all_world2pix(gaia['ra'],gaia['dec'], 0) gaia['x'] = coords[0] gaia['y'] = coords[1] else: - gp,gm = Get_Gaia(tpf,magnitude_limit=maglim) + gp,gm = Get_Gaia(ra,dec,wcs,shape,magnitude_limit=maglim) gaia = pd.DataFrame(np.array([gp[:,0],gp[:,1],gm]).T,columns=['x','y','mag']) - image = tpf.flux[10] + image = flux[10] image = strip_units(image) NY, NX = image.shape @@ -360,7 +360,7 @@ def Cat_mask(tpf,catalogue_path=None,maglim=19,scale=1,strapsize=3,ref=None,sigm sat = (np.nansum(sat,axis=0) > 0).astype(int) * 2 # assign 2 bit if strapsize > 0: - strap = Strap_mask(image,tpf.column+col_offset,strapsize,flux_cube=flux_cube).astype(int) * 4 # assign 4 bit + strap = Strap_mask(image,column+col_offset,strapsize,flux_cube=flux_cube).astype(int) * 4 # assign 4 bit else: strap = np.zeros_like(image,dtype=int) diff --git a/tessreduce/catalog_tools.py b/tessreduce/catalog_tools.py index f0d3828..2ab2b24 100755 --- a/tessreduce/catalog_tools.py +++ b/tessreduce/catalog_tools.py @@ -13,7 +13,7 @@ from .helpers import * -def Get_Catalogue(tpf, Catalog = 'gaia'): +def Get_Catalogue(ra,dec,shape, Catalog = 'gaia'): """ Get the coordinates and mag of all sources in the field of view from a specified catalogue. @@ -32,10 +32,10 @@ def Get_Catalogue(tpf, Catalog = 'gaia'): coords array coordinates of sources Gmag array Gmags of sources """ - c1 = SkyCoord(tpf.ra, tpf.dec, frame='icrs', unit='deg') + c1 = SkyCoord(ra, dec, frame='icrs', unit='deg') # Use pixel scale for query size pix_scale = 21.0 - rad = Angle(np.max(tpf.shape[1:]) * pix_scale + 60, "arcsec") + rad = Angle(np.max(shape[1:]) * pix_scale + 60, "arcsec") # We are querying with a diameter as the radius, overfilling by 2x. from astroquery.vizier import Vizier Vizier.ROW_LIMIT = -1 @@ -52,10 +52,10 @@ def Get_Catalogue(tpf, Catalog = 'gaia'): raise ValueError(f"{catalog} not recognised as a catalog. Available options: 'gaia', 'dist','ps1'") if Catalog == 'gaia': result = Vizier.query_region(c1, catalog=[catalog], - radius=Angle(np.max(tpf.shape[1:]) * pix_scale + 60, "arcsec"),column_filters={'Gmag':'<19'}) + radius=Angle(np.max(shape[1:]) * pix_scale + 60, "arcsec"),column_filters={'Gmag':'<19'}) elif Catalog == 'ps1': result = Vizier.query_region(c1, catalog=[catalog], - radius=Angle(np.max(tpf.shape[1:]) * pix_scale + 60, "arcsec")) + radius=Angle(np.max(shape[1:]) * pix_scale + 60, "arcsec")) no_targets_found_message = ValueError('Either no sources were found in the query region ' 'or Vizier is unavailable') @@ -208,7 +208,7 @@ def Get_Gaia_External(ra,dec,size,wcsObj,magnitude_limit = 18, Offset = 10): #Jmag = Jmag[ind] return radecs, Tmag, source -def Get_Gaia(tpf, magnitude_limit = 18, Offset = 10): +def Get_Gaia(ra,dec,wcs,shape, magnitude_limit = 18, Offset = 10): """ Get the coordinates and mag of all gaia sources in the field of view. @@ -230,18 +230,18 @@ def Get_Gaia(tpf, magnitude_limit = 18, Offset = 10): 'ymag','e_ymag','yKmag','e_yKmag','tmag','gaiaid','gaiamag','gaiadist','gaiadist_u','gaiadist_l', 'row','col'] - result = Get_Catalogue(tpf, Catalog = 'gaia') + result = Get_Catalogue(ra,dec,shape, Catalog = 'gaia') result = result[result.Gmag < magnitude_limit] if len(result) == 0: raise no_targets_found_message radecs = np.vstack([result['RA_ICRS'], result['DE_ICRS']]).T - coords = tpf.wcs.all_world2pix(radecs, 0) ## TODO, is origin supposed to be zero or one? + coords = wcs.all_world2pix(radecs, 0) ## TODO, is origin supposed to be zero or one? Gmag = result['Gmag'].values RPmag = result['RPmag'].values #Jmag = result['Jmag'] ind = (((coords[:,0] >= -10) & (coords[:,1] >= -10)) & - ((coords[:,0] < (tpf.shape[2] + 10)) & (coords[:,1] < (tpf.shape[1] + 10)))) + ((coords[:,0] < (shape[2] + 10)) & (coords[:,1] < (shape[1] + 10)))) coords = coords[ind] Gmag = Gmag[ind] RPmag = RPmag[ind] @@ -305,7 +305,7 @@ def SM_to_TESS_mag(SM,ebv = 0): -def Get_PS1(tpf, magnitude_limit = 20, Offset = 10): +def Get_PS1(ra,dec,wcs,shape, magnitude_limit = 20, Offset = 10): """ Get the coordinates and mag of all PS1 sources in the field of view. @@ -322,7 +322,7 @@ def Get_PS1(tpf, magnitude_limit = 20, Offset = 10): coords array coordinates of sources Gmag array Gmags of sources """ - result = Get_Catalogue(tpf, Catalog = 'ps1') + result = Get_Catalogue(ra,dec,shape, Catalog = 'ps1') result = result[np.isfinite(result.rmag) & np.isfinite(result.imag)]# & np.isfinite(result.zmag)& np.isfinite(result.ymag)] result = PS1_to_TESS_mag(result) @@ -331,11 +331,11 @@ def Get_PS1(tpf, magnitude_limit = 20, Offset = 10): if len(result) == 0: raise no_targets_found_message radecs = np.vstack([result['RAJ2000'], result['DEJ2000']]).T - coords = tpf.wcs.all_world2pix(radecs, 0) ## TODO, is origin supposed to be zero or one? + coords = wcs.all_world2pix(radecs, 0) ## TODO, is origin supposed to be zero or one? Tessmag = result['tmag'].values #Jmag = result['Jmag'] ind = (((coords[:,0] >= -Offset) & (coords[:,1] >= -Offset)) & - ((coords[:,0] < (tpf.shape[1] + Offset)) & (coords[:,1] < (tpf.shape[2] + Offset)))) + ((coords[:,0] < (shape[1] + Offset)) & (coords[:,1] < (shape[2] + Offset)))) coords = coords[ind] Tessmag = Tessmag[ind] #Jmag = Jmag[ind] @@ -378,7 +378,7 @@ def Skymapper_df(sm): return df -def Unified_catalog(tpf,magnitude_limit=18,offset=10): +def Unified_catalog(ra,dec,wcs,shape,magnitude_limit=18,offset=10): """ Find all sources present in the TESS field from PS!, and Gaia. Catalogs are cross matched through distance, and Gaia distances are assigned from Gaia ID. @@ -398,11 +398,11 @@ def Unified_catalog(tpf,magnitude_limit=18,offset=10): pd.options.mode.chained_assignment = None # need to look at how the icrs coords are offset from J2000 # Get gaia catalogs - gaia = Get_Catalogue(tpf, Catalog = 'gaia') - gaiadist = Get_Catalogue(tpf, Catalog = 'dist') + gaia = Get_Catalogue(ra,dec,shape, Catalog = 'gaia') + gaiadist = Get_Catalogue(ra,dec,shape, Catalog = 'dist') # Get PS1 and structure it - if tpf.dec > -30: - ps1 = Get_Catalogue(tpf, Catalog = 'ps1') + if dec > -30: + ps1 = Get_Catalogue(ra,dec,shape, Catalog = 'ps1') ps1 = ps1[np.isfinite(ps1.rmag) & np.isfinite(ps1.imag)]# & np.isfinite(result.zmag)& np.isfinite(result.ymag)] ps1 = PS1_to_TESS_mag(ps1) keep = ['objID','RAJ2000', 'DEJ2000','e_RAJ2000','e_DEJ2000','gmag', 'e_gmag', 'gKmag', @@ -412,7 +412,7 @@ def Unified_catalog(tpf,magnitude_limit=18,offset=10): 'tmag'] result = ps1[keep] else: - sm = Get_Catalogue(tpf, Catalog = 'skymapper') + sm = Get_Catalogue(ra,dec,shape, Catalog = 'skymapper') sm = Skymapper_df(sm) sm = sm[np.isfinite(sm.rmag) & np.isfinite(sm.imag)]# & np.isfinite(result.zmag)& np.isfinite(result.ymag)] sm = PS1_to_TESS_mag(sm) @@ -478,12 +478,12 @@ def Unified_catalog(tpf,magnitude_limit=18,offset=10): print(no_targets_found_message) radecs = np.vstack([result['RAJ2000'], result['DEJ2000']]).T - coords = tpf.wcs.all_world2pix(radecs, 0) + coords = wcs.all_world2pix(radecs, 0) result['row'] = coords[:,1] result['col'] = coords[:,0] #Jmag = result['Jmag'] ind = (((coords[:,0] >= -offset) & (coords[:,1] >= -offset)) & - ((coords[:,0] < (tpf.shape[1] + offset)) & (coords[:,1] < (tpf.shape[2] + offset)))) + ((coords[:,0] < (shape[1] + offset)) & (coords[:,1] < (shape[2] + offset)))) result = result.iloc[ind] return result diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 3015ac9..3daae84 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -143,6 +143,7 @@ def _subtract_residual_surface(bkg, flux, bkgmask, box_size=20, filter_size=5, s class tessreduce(): def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sector=None, + flux=None,mjd=None,wcs=None,camera=None,ccd=None,shifts=None, reduce=True,align=True,diff=True,corr_correction=False,kernel_match=False,calibrate=True,sourcehunt=True, phot_method='aperture',imaging=False,parallel=True,num_cores=-1,backend='loky',diagnostic_plot=False,plot=True, savename=None,quality_bitmask='hard',cache_dir=None,cache=True,catalogue_path=False, @@ -281,16 +282,22 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect self.diagnostic_plot = diagnostic_plot self.savename = savename + # Optionally given + self.flux = flux + self.mjd = mjd + self.wcs = wcs + self.camera = camera + self.ccd = ccd + self.shift = shifts + # Calculated self.mask = None #self._over_sub = None - self.shift = None self.bkg = None - self.flux = None + # self.flux = None self.delta_kernel = None self.ref = None self.ref_ind = None - self.wcs = None self.qe = None self.lc = None self.sky = None @@ -357,6 +364,40 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect if not self.sector: self.sector = 999 + self.column = tpf.column + self.row = tpf.row + self.camera = self.tpf.camera + self.ccd = self.tpf.ccd + + # -- Allow for cube to be given directly, but require MJD, WCS, Sector to be given as well -- # + elif self.flux is not None: + if self.mjd is None: + m = 'If flux is given, the MJD must also be given.' + raise ValueError(m) + if self.wcs is None: + m = 'If flux is given, the WCS must also be given.' + raise ValueError(m) + if self.sector is None: + m = 'If flux is given, the sector must also be given.' + raise ValueError(m) + if self._force_ref_ind is None: + m = 'If flux is given, the reference frame index must also be given.' + raise ValueError(m) + if self.camera is None: + m = 'If flux is given, the camera must also be given.' + raise ValueError(m) + if self.ccd is None: + m = 'If flux is given, the CCD must also be given.' + raise ValueError(m) + + self.size = self.flux.shape[1] + self.ra,self.dec = self.wcs.all_pix2world(self.size//2,self.size//2,0) + self.eflux = None + self.column = 0 + self.row = 0 + self.rawflux = deepcopy(self.flux) + + # Retrieve TPF elif self.check_coord(): if self.verbose>0: @@ -406,7 +447,7 @@ def _get_gaia(self,maglim=21): """ # Get dataframe from Gaia around cutout - result = Get_Catalogue(self.tpf, Catalog = 'gaia') + result = Get_Catalogue(self.ra,self.dec,self.flux.shape, Catalog = 'gaia') result = result[result.Gmag < maglim] result = result.rename(columns={'RA_ICRS': 'ra', 'DE_ICRS': 'dec', @@ -541,6 +582,11 @@ def get_TESS(self,ra=None,dec=None,name=None,size=None,sector=None, self.eflux = None self.wcs = tpf.wcs self.mjd = tpf.time.mjd + self.column = tpf.column + self.row = tpf.row + self.camera = self.tpf.camera + self.ccd = self.tpf.ccd + def make_mask(self,catalogue_path=None,maglim=19,scale=1,strapsize=6,useref=False): """ @@ -580,9 +626,9 @@ def make_mask(self,catalogue_path=None,maglim=19,scale=1,strapsize=6,useref=Fals # Generate mask from source catalogue if useref: - mask, cat = Cat_mask(self.tpf,catalogue_path,maglim,scale,strapsize,ref=self.ref,col_offset=self._col_offset) + mask, cat = Cat_mask(self.ra,self.dec,self.flux.shape,self.wcs,self.flux,self.column,catalogue_path,maglim,scale,strapsize,ref=self.ref,col_offset=self._col_offset) else: - mask, cat = Cat_mask(self.tpf,catalogue_path,maglim,scale,strapsize,col_offset=self._col_offset) + mask, cat = Cat_mask(self.ra,self.dec,self.flux.shape,self.wcs,self.flux,self.column,catalogue_path,maglim,scale,strapsize,ref=self.ref,col_offset=self._col_offset) # Generate sky background as the inverse of mask sky = ((mask & 1)+1 == 1) * 1. @@ -642,8 +688,8 @@ def psf_source_mask(self,sigma=5): from PRF import TESS_PRF - col = self.tpf.column + int(self.size//2) # find column and row, when specifying location on a *say* 90x90 px cutout - row = self.tpf.row + int(self.size//2) + col = self.column + int(self.size//2) # find column and row, when specifying location on a *say* 90x90 px cutout + row = self.row + int(self.size//2) col += 45 # add on the non-science columns row += 1 # add on the non-science row @@ -656,16 +702,16 @@ def psf_source_mask(self,sigma=5): if self._catalogue_path is not None: if self.sector < 4: - prf = TESS_PRF(self.tpf.camera,self.tpf.ccd,self.sector, + prf = TESS_PRF(self.camera,self.ccd,self.sector, col,row, localdatadir=f'{self._prf_path}/Sectors1_2_3') else: - prf = TESS_PRF(self.tpf.camera,self.tpf.ccd,self.sector, + prf = TESS_PRF(self.camera,self.ccd,self.sector, col,row, localdatadir=f'{self._prf_path}/Sectors4+') else: try: - prf = TESS_PRF(self.tpf.camera,self.tpf.ccd,self.sector, + prf = TESS_PRF(self.camera,self.ccd,self.sector, col,row) except Exception as e: print(f'Warning: could not load PRF (network error?): {e}') @@ -918,7 +964,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa # if calc_qe: bkg_pre_fix = np.array(self.bkg) from .adaptive_background import get_tessvectors, _interpolate_angles - _df = get_tessvectors(self.sector, self.tpf.camera, data_path=self._vector_path) + _df = get_tessvectors(self.sector, self.camera, data_path=self._vector_path) _bkg_median = np.nanmedian(self.bkg, axis=(1, 2)) if _df is not None: _earth_angle, _moon_angle = _interpolate_angles(self.mjd, _df) @@ -942,7 +988,7 @@ def background(self,gauss_smooth=None,calc_qe=True,strap_iso=True,source_hunt=Fa if self.verbose > 0: print('adaptive temporal smoothing...') from .adaptive_background import AdaptiveBackground - smoother = AdaptiveBackground(self.bkg, self.mjd, sector=self.sector, camera=self.tpf.camera, + smoother = AdaptiveBackground(self.bkg, self.mjd, sector=self.sector, camera=self.camera, data_path=self._vector_path,n_jobs=self.num_cores,backend=self.backend,verbose=self._joblib_verbose) if smoother._df is not None: smoothed = smoother.smooth(method='savgol').smoothed @@ -1435,8 +1481,8 @@ def centroids_shifts_starfind(self,plot=None,savename=None): mean, med, std = sigma_clipped_stats(m, sigma=3.0) - prf = TESS_PRF(self.tpf.camera,self.tpf.ccd,self.tpf.sector, - self.tpf.column+self.flux.shape[2]/2,self.tpf.row+self.flux.shape[1]/2) + prf = TESS_PRF(self.camera,self.ccd,self.sector, + self.column+self.flux.shape[2]/2,self.row+self.flux.shape[1]/2) self.prf = prf.locate(5,5,(11,11)) finder = StarFinder(2*std,kernel=self.prf,exclude_border=True) @@ -1468,7 +1514,7 @@ def centroids_shifts_starfind(self,plot=None,savename=None): self.shift = meds #smooth if plot: - t = self.tpf.time.mjd + t = self.mjd ind = np.where(np.diff(t) > .5)[0] smooth[ind,:] = np.nan plt.figure(figsize=(1.5*fig_width,1*fig_width)) @@ -1541,7 +1587,7 @@ def fit_shift(self,smooth=True,plot=None,savename=None): self.shift = shifts if plot: - t = self.tpf.time.mjd + t = self.mjd ind = np.where(np.diff(t) > .5)[0] shifts[ind,:] = np.nan plt.figure(figsize=(1.5*fig_width,1*fig_width)) @@ -1559,7 +1605,7 @@ def fit_shift(self,smooth=True,plot=None,savename=None): def plot_shifts(self,savename=None): import matplotlib.pyplot as plt - t = self.tpf.time.mjd + t = self.mjd shifts = self.shift ind = np.where(np.diff(t) > .5)[0] shifts[ind,:] = np.nan @@ -1699,7 +1745,7 @@ def bin_flux(self,flux=None,time_bin=6/24,frames = None): if flux is None: flux = self.flux - t = self.tpf.time.mjd + t = self.mjd if time_bin is None: bin_size = int(frames) @@ -1814,7 +1860,7 @@ def diff_lc(self,time=None,x=None,y=None,ra=None,dec=None,tar_ap=3, print(Warning('sky_out must be odd, adding 1')) sky_in += 1 - if (ra is not None) & (dec is not None) & (self.tpf is not None): + if (ra is not None) & (dec is not None): x,y = self.wcs.all_world2pix(ra,dec,0) x = int(np.round(np.ravel(x)[0],0)) y = int(np.round(np.ravel(y)[0],0)) @@ -1894,8 +1940,8 @@ def diff_lc(self,time=None,x=None,y=None,ra=None,dec=None,tar_ap=3, mask = self.orbit_segments == seg tar[mask] += delta - if self.tpf is not None: - time = self.tpf.time.mjd + time = self.mjd + lc = np.array([time, tar, tar_err]) sky = np.array([time, sky_med, sky_std]) @@ -2248,8 +2294,8 @@ def _psf_initialise(self,cutoutSize,loc,time_ind=None,ref=False): time_ind = np.arange(0,len(self.flux)) - col = self.tpf.column - int(self.size//2) + loc[0] # find column and row, when specifying location on a *say* 90x90 px cutout - row = self.tpf.row - int(self.size//2) + loc[1] + col = self.column - int(self.size//2) + loc[0] # find column and row, when specifying location on a *say* 90x90 px cutout + row = self.row - int(self.size//2) + loc[1] if isinstance(loc[0], (float, np.floating, np.float32, np.float64)): loc[0] = int(np.round(loc[0],0)) @@ -2265,9 +2311,9 @@ def _psf_initialise(self,cutoutSize,loc,time_ind=None,ref=False): row = np.max([row,10]) col = np.max([col,45]) try: - prf = TESS_PRF(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row) # initialise psf kernel + prf = TESS_PRF(self.camera,self.ccd,self.sector,col,row) # initialise psf kernel except: - print(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row) + print(self.camera,self.ccd,self.sector,col,row) raise ValueError if ref: cutout = (self.flux+self.ref)[time_ind,loc[1]-cutoutSize//2:loc[1]+1+cutoutSize//2,loc[0]-cutoutSize//2:loc[0]+1+cutoutSize//2] # gather cutouts @@ -2315,8 +2361,8 @@ def moving_psf_photometry(self,xpos,ypos,size=5,time_ind=None,xlim=2,ylim=2): raise ValueError(m) inds = np.arange(0,len(xpos)) if self.parallel: - prfs, cutouts, ecutouts = zip(*Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(par_psf_initialise)(self.flux,self.tpf.camera,self.tpf.ccd, - self.tpf.sector,self.tpf.column,self.tpf.row, + prfs, cutouts, ecutouts = zip(*Parallel(n_jobs=self.num_cores, backend=self.backend, verbose=self._joblib_verbose)(delayed(par_psf_initialise)(self.flux,self.camera,self.ccd, + self.sector,self.column,self.row, size,[xpos[i],ypos[i]],time_ind) for i in inds)) else: prfs = [] @@ -2365,8 +2411,8 @@ def psf_photutils(self,xPix=None,yPix=None,size=5,local_bkg=False,epsf=None, if epsf is None: if self.epsf is None: - col = self.tpf.column - int(self.size//2) + yPix # find column and row, when specifying location on a *say* 90x90 px cutout - row = self.tpf.row - int(self.size//2) + xPix + col = self.column - int(self.size//2) + yPix # find column and row, when specifying location on a *say* 90x90 px cutout + row = self.row - int(self.size//2) + xPix col += 45 # add on the non-science columns row += 1 # add on the non-science row if col > 2090: @@ -2374,7 +2420,7 @@ def psf_photutils(self,xPix=None,yPix=None,size=5,local_bkg=False,epsf=None, if row > 2040: row = 2040 - self.epsf = simulate_epsf(self.tpf.camera,self.tpf.ccd,self.tpf.sector,col,row) + self.epsf = simulate_epsf(self.camera,self.ccd,self.sector,col,row) epsf = self.epsf if local_bkg: @@ -2575,8 +2621,8 @@ def orbit_ref_subtract(self): Updates self.flux in place and stores self.orbit_segments. """ - sector = (self.tpf.sector if self.tpf is not None else None) or self.sector - camera = (self.tpf.camera if self.tpf is not None else None) or self.camera + sector = self.sector + camera = self.camera flux, segments, orbit_refs = orbit_ref_subtract(deepcopy(self.flux), self.mjd, sector=sector, camera=camera, vector_path=self._vector_path) @@ -2740,7 +2786,7 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, raise ValueError('flux all nans') _t = time.perf_counter() - if self.align: + if self.align and self.shift is None: if self.verbose > 0: print('aligning images') @@ -2774,16 +2820,20 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, # self.centroids_shifts_starfind() # elif self._shift_method == 'minimize': # self.fit_shift() - else: + elif self.shift is None: self.shift = np.zeros((len(self.flux),2)) + _times['alignment'] = time.perf_counter() - _t print('alignment done') + rawcube = self.rawflux if hasattr(self,'rawflux') else self.tpf.flux.value + _bad_frames = np.nansum(rawcube,axis=(1,2))==0 + if not self.diff: if self.align: _t = time.perf_counter() self.shift_images() - self.flux[np.nansum(self.tpf.flux.value,axis=(1,2))==0] = np.nan + self.flux[_bad_frames] = np.nan _times['image shifting'] = time.perf_counter() - _t if self.verbose > 0: print('images shifted') @@ -2795,7 +2845,8 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, print('!!Re-running for difference image!!') # reseting to do diffim _t = time.perf_counter() - self.flux = strip_units(self.tpf.flux) + + self.flux = rawcube self.flux = self.flux / self.qe if self.align: @@ -2806,7 +2857,8 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, self._flux_aligned = deepcopy(self.flux) if test_seed is not None: self.flux += test_seed - self.flux[np.nansum(self.tpf.flux.value,axis=(1,2))==0] = np.nan + rawcube = self.rawflux if hasattr(self,'rawflux') else self.tpf.flux.value + self.flux[_bad_frames] = np.nan # subtract reference if self._ref_type.lower() == 'single': self.ref = deepcopy(self.flux[self.ref_ind]) @@ -2943,7 +2995,7 @@ def make_lc(self,aperture = None,bin_size=0,zeropoint=None,scale='counts',clip = # hack solution for new lightkurve flux = strip_units(self.flux) - t = self.tpf.time.mjd + t = self.mjd if type(aperture) == type(None): aper = np.zeros_like(flux[0]) @@ -3401,13 +3453,13 @@ def isolated_star_lcs(self): if self.dec < -30: if self.verbose > 0: print('target is below -30 dec, calibrating to SkyMapper photometry.') - table = Get_Catalogue(self.tpf,Catalog='skymapper') + table = Get_Catalogue(self.ra,self.dec,self.flux.shape,Catalog='skymapper') table = Skymapper_df(table) system = 'skymapper' else: if self.verbose > 0: print('target is above -30 dec, calibrating to PS1 photometry.') - table = Get_Catalogue(self.tpf,Catalog='ps1') + table = Get_Catalogue(self.ra,self.dec,self.flux.shape,Catalog='ps1') system = 'ps1' if self.diff: @@ -3539,7 +3591,7 @@ def field_calibrate(self,zp_single=True,plot=None,savename=None): if self.dec < -30: if self.verbose > 0: print('target is below -30 dec, calibrating to SkyMapper photometry.') - table = Get_Catalogue(self.tpf,Catalog='skymapper') + table = Get_Catalogue(self.ra,self.dec,self.flux.shape,Catalog='skymapper') system = 'skymapper' if table is None: print('WARNING: SkyMapper unavailable, skipping field calibration.') @@ -3547,7 +3599,7 @@ def field_calibrate(self,zp_single=True,plot=None,savename=None): else: if self.verbose > 0: print('target is above -30 dec, calibrating to PS1 photometry.') - table = Get_Catalogue(self.tpf,Catalog='ps1') + table = Get_Catalogue(self.ra,self.dec,self.flux.shape,Catalog='ps1') system = 'ps1' x,y = self.wcs.all_world2pix(table.RAJ2000.values,table.DEJ2000.values,0) table['col'] = x @@ -3736,7 +3788,7 @@ def field_calibrate(self,zp_single=True,plot=None,savename=None): plt.gca().xaxis.set_major_locator(plt.MaxNLocator(6)) plt.subplot(122) - plt.plot(self.tpf.time.mjd[mask],mzp[mask],'.',alpha=0.5) + plt.plot(self.mjd[mask],mzp[mask],'.',alpha=0.5) #plt.axhspan(averager.mean-averager.stdev,averager.mean+averager.stdev,alpha=0.3,color='C1') #plt.axhline(averager.mean,color='C1') #plt.axhspan(med-std,med+std,alpha=0.3,color='C1') @@ -3759,7 +3811,7 @@ def field_calibrate(self,zp_single=True,plot=None,savename=None): else: zp = np.nanmedian(zp,axis=0) - mzp,stdzp = smooth_zp(zp, self.tpf.time.mjd) + mzp,stdzp = smooth_zp(zp, self.mjd) compare = (abs(mzp-20.44) > 2).any() if compare: From 0114c56e48186c348fa3f1a33abb564276657ff2 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 6 Aug 2026 10:19:24 +1200 Subject: [PATCH 21/22] Store tpf-derived values as class attributes and drop the tpf object Extract flux_raw, quality, column, row, camera, ccd, mjd, wcs at load time (both the direct-tpf path and get_TESS) and close/discard the tpf object afterward so its flux cube isn't kept duplicated in memory alongside self.flux. Smooth_motion and the reference-frame quality check no longer depend on a live tpf object, which also fixes a crash in the flux-array-only init path (self.tpf is None there) since reduce()'s default alignment methods called Smooth_motion(..., self.tpf). Also fixes make_mask(useref=False) incorrectly passing ref=self.ref to Cat_mask, which made it behave like useref=True whenever self.ref was already set. --- tessreduce/helpers.py | 22 +++++++++++--------- tessreduce/tessreduce.py | 44 ++++++++++++++++++++++++++++++---------- 2 files changed, 45 insertions(+), 21 deletions(-) diff --git a/tessreduce/helpers.py b/tessreduce/helpers.py index b814f52..c911f85 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -384,18 +384,20 @@ def difference_shifts(image,ref):#,eimage,eref): s = np.zeros((2)) * np.nan return s -def Smooth_motion(Centroids,tpf,skernel=25): +def Smooth_motion(Centroids,mjd,flux,skernel=25): """ - Calculate the smoothed centroid shift + Calculate the smoothed centroid shift Parameters ---------- Centroids : array centroid shifts from all frames + mjd : array + time of each frame, used to find orbit gaps - TPF : lightkurve targetpixelfile - tpf + flux : array + flux cube, used to identify padded/empty frames Returns ------- @@ -411,11 +413,11 @@ def Smooth_motion(Centroids,tpf,skernel=25): # skernel = 25 try: try: - split = np.where(np.diff(tpf.time.mjd) > 0.5)[0][0] + 1 + split = np.where(np.diff(mjd) > 0.5)[0][0] + 1 # ugly, but who cares - ind1 = np.nansum(tpf.flux[:split],axis=(1,2)) + ind1 = np.nansum(flux[:split],axis=(1,2)) ind1 = np.where(ind1 != 0)[0] - ind2 = np.nansum(tpf.flux[split:],axis=(1,2)) + ind2 = np.nansum(flux[split:],axis=(1,2)) ind2 = np.where(ind2 != 0)[0] + split smoothed[ind1,0] = savgol_filter(Centroids[ind1,0],skernel,3) smoothed[ind2,0] = savgol_filter(Centroids[ind2,0],skernel,3) @@ -423,11 +425,11 @@ def Smooth_motion(Centroids,tpf,skernel=25): smoothed[ind1,1] = savgol_filter(Centroids[ind1,1],skernel,3) smoothed[ind2,1] = savgol_filter(Centroids[ind2,1],skernel,3) except: - split = np.where(np.diff(tpf.time.mjd) > 0.5)[0][0] + 1 + split = np.where(np.diff(mjd) > 0.5)[0][0] + 1 # ugly, but who cares - ind1 = np.nansum(tpf.flux[:split],axis=(1,2)) + ind1 = np.nansum(flux[:split],axis=(1,2)) ind1 = np.where(ind1 != 0)[0] - ind2 = np.nansum(tpf.flux[split:],axis=(1,2)) + ind2 = np.nansum(flux[split:],axis=(1,2)) ind2 = np.where(ind2 != 0)[0] + split smoothed[ind1,0] = savgol_filter(Centroids[ind1,0],skernel//2+1,3) smoothed[ind2,0] = savgol_filter(Centroids[ind2,0],skernel//2+1,3) diff --git a/tessreduce/tessreduce.py b/tessreduce/tessreduce.py index 3daae84..584be6b 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -290,11 +290,15 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect self.ccd = ccd self.shift = shifts - # Calculated + # Calculated self.mask = None #self._over_sub = None self.bkg = None # self.flux = None + self.flux_raw = None + self.quality = None + self.column = None + self.row = None self.delta_kernel = None self.ref = None self.ref_ind = None @@ -368,6 +372,16 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect self.row = tpf.row self.camera = self.tpf.camera self.ccd = self.tpf.ccd + self.quality = self.tpf.quality + self.flux_raw = deepcopy(self.flux) + + # All values needed downstream are now class attributes, so drop the + # tpf object rather than keep its flux cube duplicated in memory. + try: + self.tpf.hdu.close() + except Exception: + pass + self.tpf = None # -- Allow for cube to be given directly, but require MJD, WCS, Sector to be given as well -- # elif self.flux is not None: @@ -395,7 +409,8 @@ def __init__(self,ra=None,dec=None,name=None,obs_list=None,tpf=None,size=90,sect self.eflux = None self.column = 0 self.row = 0 - self.rawflux = deepcopy(self.flux) + self.quality = np.zeros(len(self.flux),dtype=int) + self.flux_raw = deepcopy(self.flux) # Retrieve TPF @@ -586,6 +601,16 @@ def get_TESS(self,ra=None,dec=None,name=None,size=None,sector=None, self.row = tpf.row self.camera = self.tpf.camera self.ccd = self.tpf.ccd + self.quality = self.tpf.quality + self.flux_raw = deepcopy(self.flux) + + # All values needed downstream are now class attributes, so drop the + # tpf object rather than keep its flux cube duplicated in memory. + try: + self.tpf.hdu.close() + except Exception: + pass + self.tpf = None def make_mask(self,catalogue_path=None,maglim=19,scale=1,strapsize=6,useref=False): @@ -628,7 +653,7 @@ def make_mask(self,catalogue_path=None,maglim=19,scale=1,strapsize=6,useref=Fals if useref: mask, cat = Cat_mask(self.ra,self.dec,self.flux.shape,self.wcs,self.flux,self.column,catalogue_path,maglim,scale,strapsize,ref=self.ref,col_offset=self._col_offset) else: - mask, cat = Cat_mask(self.ra,self.dec,self.flux.shape,self.wcs,self.flux,self.column,catalogue_path,maglim,scale,strapsize,ref=self.ref,col_offset=self._col_offset) + mask, cat = Cat_mask(self.ra,self.dec,self.flux.shape,self.wcs,self.flux,self.column,catalogue_path,maglim,scale,strapsize,col_offset=self._col_offset) # Generate sky background as the inverse of mask sky = ((mask & 1)+1 == 1) * 1. @@ -1409,7 +1434,7 @@ def get_ref(self,start = None, stop = None): start = int(start) stop = int(stop) - ind = self.tpf.quality[start:stop] == 0 + ind = self.quality[start:stop] == 0 d = deepcopy(data[start:stop])[ind] summed = np.nanmedian(d,axis=(1,2)) summed[summed <=0] = 1e5 @@ -1508,7 +1533,7 @@ def centroids_shifts_starfind(self,plot=None,savename=None): meds = np.nanmedian(shifts,axis = 2) meds[~np.isfinite(meds)] = 0 - smooth = Smooth_motion(meds,self.tpf) + smooth = Smooth_motion(meds,self.mjd,self.flux_raw) nans = np.nansum(f,axis=(1,2)) ==0 smooth[nans] = np.nan self.shift = meds #smooth @@ -1579,7 +1604,7 @@ def fit_shift(self,smooth=True,plot=None,savename=None): shifts[i,:] = difference_shifts(f[i],m) sraw = deepcopy(shifts) if smooth: - shifts = Smooth_motion(shifts,self.tpf) + shifts = Smooth_motion(shifts,self.mjd,self.flux_raw) if self.shift is not None: self.shift += shifts @@ -2826,8 +2851,7 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, _times['alignment'] = time.perf_counter() - _t print('alignment done') - rawcube = self.rawflux if hasattr(self,'rawflux') else self.tpf.flux.value - _bad_frames = np.nansum(rawcube,axis=(1,2))==0 + _bad_frames = np.nansum(self.flux_raw,axis=(1,2))==0 if not self.diff: if self.align: @@ -2846,8 +2870,7 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, # reseting to do diffim _t = time.perf_counter() - self.flux = rawcube - self.flux = self.flux / self.qe + self.flux = self.flux_raw / self.qe if self.align: self.shift_images() @@ -2857,7 +2880,6 @@ def reduce(self, aper = None, align = None, parallel = None, calibrate=None, self._flux_aligned = deepcopy(self.flux) if test_seed is not None: self.flux += test_seed - rawcube = self.rawflux if hasattr(self,'rawflux') else self.tpf.flux.value self.flux[_bad_frames] = np.nan # subtract reference if self._ref_type.lower() == 'single': From ebb9e393dfa7381318e499e6e076ae271d799af8 Mon Sep 17 00:00:00 2001 From: Ryan Ridden Date: Thu, 6 Aug 2026 10:28:22 +1200 Subject: [PATCH 22/22] Update sector_mjd.csv with revised end times and extend through sector 107 mjd_end values for sectors 56-102 were corrected (previously truncated to the start of the next sector's gap rather than its true end), and new rows were added for sectors 103-107. --- tessreduce/sector_mjd.csv | 100 ++++++++++++++++++++------------------ 1 file changed, 52 insertions(+), 48 deletions(-) diff --git a/tessreduce/sector_mjd.csv b/tessreduce/sector_mjd.csv index d24941b..2ca36ed 100755 --- a/tessreduce/sector_mjd.csv +++ b/tessreduce/sector_mjd.csv @@ -54,51 +54,55 @@ Sector,mjd_start,mjd_end 53,59743.49652777778,59768.48611111111 54,59769.399305555555,59795.635416666664 55,59796.600694444445,59823.770833333336 -56,59824.756944444445,59838.020833333336 -57,59852.854166666664,59866.75347222222 -58,59881.82986111111,59895.5625 -59,59909.76388888889,59922.854166666664 -60,59936.40277777778,59949.274305555555 -61,59962.29861111111,59974.70486111111 -62,59987.94097222222,60000.430555555555 -63,60013.868055555555,60026.899305555555 -64,60040.614583333336,60054.586805555555 -65,60068.239583333336,60082.98263888889 -66,60097.17361111111,60111.48263888889 -67,60126.14236111111,60139.756944444445 -68,60154.11111111111,60167.77777777778 -69,60181.854166666664,60194.5 -70,60207.854166666664,60220.631944444445 -71,60233.53472222222,60246.51388888889 -72,60259.68402777778,60272.506944444445 -73,60285.29861111111,60298.447916666664 -74,60312.364583333336,60325.354166666664 -75,60339.28472222222,60353.0625 -76,60367.197916666664,60380.89236111111 -77,60394.98611111111,60408.774305555555 -78,60423.26388888889,60437.166666666664 -79,60452.038194444445,60465.10763888889 -80,60479.38888888889,60492.850694444445 -81,60506.052083333336,60519.07986111111 -82,60532.89236111111,60544.96875 -83,60558.927083333336,60571.21527777778 -84,60584.09027777778,60596.805555555555 -85,60610.055555555555,60622.739583333336 -86,60635.760416666664,60649.114583333336 -87,60662.541666666664,60675.947916666664 -88,60689.65277777778,60703.395833333336 -89,60717.63888888889,60732.45138888889 -90,60746.65972222222,60760.239583333336 -91,60774.791666666664,60788.76388888889 -92,60802.47222222222,60815.98611111111 -93,60829.35763888889,60842.63888888889 -94,60855.76388888889,60868.23611111111 -95,60881.833333333336,60894.17361111111 -96,60907.274305555555,60919.989583333336 -97,60933.15972222222,60946.37152777778 -98,60988.0,61002.302083333336 -99,61045.614583333336,61059.23611111111 -100,61073.493055555555,61086.78125 -101,61100.180555555555,61113.04513888889 -102,61126.520833333336,61139.024305555555 -103,61151.53125,61164.729166666664 +56,59824.756944444445,59852.645833333336 +57,59852.854166666664,59881.62152777778 +58,59881.82986111111,59909.555555555555 +59,59909.76388888889,59936.194444444445 +60,59936.40277777778,59962.09027777778 +61,59962.29861111111,59987.73263888889 +62,59987.94097222222,60013.65972222222 +63,60013.868055555555,60040.40625 +64,60040.614583333336,60068.03125 +65,60068.239583333336,60096.96527777778 +66,60097.17361111111,60125.93402777778 +67,60126.14236111111,60153.90277777778 +68,60154.11111111111,60181.645833333336 +69,60181.854166666664,60207.645833333336 +70,60207.854166666664,60233.34375 +71,60233.53472222222,60259.475694444445 +72,60259.68402777778,60285.09027777778 +73,60285.29861111111,60312.15625 +74,60312.364583333336,60339.07638888889 +75,60339.28472222222,60366.989583333336 +76,60367.197916666664,60394.77777777778 +77,60394.98611111111,60423.055555555555 +78,60423.26388888889,60451.82986111111 +79,60452.038194444445,60479.180555555555 +80,60479.38888888889,60505.84375 +81,60506.052083333336,60532.68402777778 +82,60532.89236111111,60558.75347222222 +83,60558.927083333336,60583.881944444445 +84,60584.09027777778,60609.84722222222 +85,60610.055555555555,60635.552083333336 +86,60635.760416666664,60662.333333333336 +87,60662.541666666664,60689.444444444445 +88,60689.65277777778,60717.430555555555 +89,60717.63888888889,60746.45138888889 +90,60746.65972222222,60774.583333333336 +91,60774.791666666664,60802.26388888889 +92,60802.47222222222,60829.149305555555 +93,60829.35763888889,60855.555555555555 +94,60855.76388888889,60881.625 +95,60881.833333333336,60907.06597222222 +96,60907.274305555555,60932.95138888889 +97,60933.15972222222,60987.791666666664 +98,60988.0,61045.40625 +99,61045.614583333336,61073.28472222222 +100,61073.493055555555,61099.97222222222 +101,61100.180555555555,61126.3125 +102,61126.520833333336,61164.729166666664 +103,61164.9375,61177.87847222222 +104,61178.086805555555,61204.708333333336 +105,61204.916666666664,61232.72222222222 +106,61232.930555555555,61261.57638888889 +107,61261.78472222222,61290.28125