diff --git a/tessreduce/adaptive_background.py b/tessreduce/adaptive_background.py index b636cf2..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)(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)( + 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)(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 751c806..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)( + 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/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/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..26b1ad9 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] @@ -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 @@ -305,7 +307,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 +324,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 +333,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 +380,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 +400,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 +414,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 +480,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/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/helpers.py b/tessreduce/helpers.py index 4e93ba5..c911f85 100644 --- a/tessreduce/helpers.py +++ b/tessreduce/helpers.py @@ -166,29 +166,37 @@ 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 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)( - delayed(_background2d_frame)(frame, box_size, filter_size, sc, estimator, mask) +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)) def Smooth_bkg(data, gauss_smooth=0, interpolate=False, extrapolate=True): @@ -247,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 @@ -376,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 ------- @@ -403,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) @@ -415,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) @@ -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[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, backend='loky', verbose=0): 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, 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: + 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 — 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) - 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 (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 - 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[1:] >= max_size) | (b_obj_size[1:] <= 9))[0] + 1 + 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, backend='loky', verbose=0): 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, backend=backend, verbose=verbose)( + delayed(_strap_fit_col)(frame[:, i], norm[:, i]) + for i in sind) + for col, q in zip(sind, results): + qe[:, col] = q return qe @@ -1391,11 +1405,201 @@ 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 + 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) 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, + sigma_clip=sc, bkg_estimator=MedianBackground(), + mask=transient_mask, fill_value=0.0) + corr = b.background + except Exception: + 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 + 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, 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. @@ -1454,227 +1658,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 - - 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'] - 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) - 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 - circ = np.sqrt((xx - cx)**2 + (yy - cy)**2) <= true_r - sep_mask |= circ - - 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) - 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 - 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) - 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 - - return fixed, excess, sharp_mask, bad_bkg_mask - - results = Parallel(n_jobs=n_jobs, prefer='threads')(delayed(_process)(i) for i in range(T)) + def _bkgmask_i(i): + if bkgmask_arr is None: + return None + return bkgmask_arr[i] if bkgmask_arr.ndim == 3 else bkgmask_arr + + 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, + _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) - 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 - return corr - residuals = flux - bkg_fixed - corrections = Parallel(n_jobs=n_jobs, prefer='threads')( - delayed(_fit_residual)( + 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) + 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) @@ -1776,7 +1794,25 @@ 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): + _, _, 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, 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 + 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, 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 @@ -1803,25 +1839,13 @@ 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 + 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=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)) + return np.array(results) diff --git a/tessreduce/lastpercent.py b/tessreduce/lastpercent.py index 2a9486c..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)(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)( + 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/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/rescale_straps.py b/tessreduce/rescale_straps.py index 85dfc08..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)(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/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/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 diff --git a/tessreduce/sep_aligner.py b/tessreduce/sep_aligner.py index 4b58475..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='loky')( + 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 2246b25..c79dacc 100644 --- a/tessreduce/tessreduce.py +++ b/tessreduce/tessreduce.py @@ -22,6 +22,35 @@ import multiprocessing from joblib import Parallel, delayed +def _available_cores(verbose=1): + """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: + n = int(slurm) + 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 >= 3: + print(f'[tessreduce] _available_cores: sched_getaffinity → {n} cores') + return n + except AttributeError: + pass + n = multiprocessing.cpu_count() + if verbose >= 3: + print(f'[tessreduce] _available_cores: cpu_count fallback → {n} cores') + return n + from .catalog_tools import * from .calibration_tools import * from .ground_tools import ground @@ -42,12 +71,28 @@ # 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 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 @@ -79,35 +124,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 = -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, 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) return bkg @@ -116,8 +143,9 @@ def _fit_frame(residual): 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,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, @@ -165,6 +193,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 @@ -180,7 +212,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 none, 1 is verbose. The default is 1. + 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. @@ -206,9 +240,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 isinstance(num_cores, str): - self.num_cores = multiprocessing.cpu_count() + if num_cores == -1 or isinstance(num_cores, str): + self.num_cores = _available_cores(verbose=verbose) else: self.num_cores = num_cores self._assign_phot_method(phot_method) @@ -223,6 +258,15 @@ 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 + + # SLURM environment diagnostics + 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} + print(f'[tessreduce] SLURM env: {_slurm_env}') + print(f'[tessreduce] num_cores resolved to: {self.num_cores}') # Offline Paths if catalogue_path is None: @@ -230,7 +274,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 @@ -240,16 +283,26 @@ 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 - # Calculated + # 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.flux_raw = None + self.quality = None + self.column = None + self.row = None self.delta_kernel = None self.ref = None self.ref_ind = None - self.wcs = None self.qe = None self.lc = None self.sky = None @@ -316,6 +369,51 @@ 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 + 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: + 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.quality = np.zeros(len(self.flux),dtype=int) + self.flux_raw = deepcopy(self.flux) + + # Retrieve TPF elif self.check_coord(): if self.verbose>0: @@ -328,6 +426,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 (>=3).""" + return 1 if self.verbose >= 3 else 0 + def check_coord(self): """ Checks if target coordinate / name input is valid. @@ -360,7 +463,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', @@ -477,12 +580,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: - 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: @@ -498,6 +598,21 @@ 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 + 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): """ @@ -537,9 +652,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,col_offset=self._col_offset) # Generate sky background as the inverse of mask sky = ((mask & 1)+1 == 1) * 1. @@ -599,8 +714,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 @@ -613,16 +728,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}') @@ -634,7 +749,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, 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) @@ -654,8 +769,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,:] @@ -760,8 +876,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() - bkg_smth = Parallel(n_jobs=self.num_cores)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*m) + 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 >= 2: + print(f'smooth background done ({_times["initial smooth background"]:.2f}s)') if rerun_negative: _t = time.perf_counter() if self._use_error_image: @@ -782,9 +902,12 @@ 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) + 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 >= 2: + print(f'smooth background rerun done ({_times["negative over-subtraction rerun"]:.2f}s)') if rerun_diff: _t = time.perf_counter() @@ -803,13 +926,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 @@ -835,10 +962,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) - bkg_smth = Parallel(n_jobs=self.num_cores)(delayed(Smooth_bkg)(frame,0,interpolate) for frame in flux*new_mask) + 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) + 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() @@ -860,30 +991,38 @@ 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) - _bkg_median = np.array([np.nanmedian(self.bkg[i]) for i in range(len(self.bkg))]) + _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) _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 >= 2: + 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 + if self.verbose >= 2: + print(f'fixing background anomalies done ({_times["anomaly fixing"]:.2f}s)') _t = time.perf_counter() + 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, - data_path=self._vector_path,n_jobs=self.num_cores) + 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 self.bkg = smoothed _times['adaptive temporal smoothing'] = time.perf_counter() - _t + 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: @@ -897,7 +1036,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) @@ -978,7 +1117,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, 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) @@ -1002,7 +1141,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, 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) @@ -1028,7 +1167,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, 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) @@ -1297,7 +1436,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 @@ -1369,8 +1508,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) @@ -1384,7 +1523,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, backend=self.backend, verbose=self._joblib_verbose)( delayed(Calculate_shifts)(frame,mx,my,finder) for frame in f) shifts = np.array(shifts) else: @@ -1396,13 +1535,13 @@ 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 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)) @@ -1456,7 +1595,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, 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) @@ -1467,7 +1606,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 @@ -1475,7 +1614,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)) @@ -1493,7 +1632,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 @@ -1524,20 +1663,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, 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) + 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, 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) + 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): """ @@ -1622,7 +1772,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) @@ -1737,14 +1887,14 @@ 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(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]) @@ -1817,8 +1967,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]) @@ -1880,8 +2030,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) @@ -2171,8 +2321,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)) @@ -2188,9 +2338,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 @@ -2238,8 +2388,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)(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 = [] @@ -2251,7 +2401,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, 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 = [] @@ -2288,8 +2438,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: @@ -2297,7 +2447,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: @@ -2337,14 +2487,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, 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)(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 @@ -2408,7 +2558,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, 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 = [] @@ -2454,9 +2604,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, 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)(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])] @@ -2472,9 +2622,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, 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)(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])] @@ -2487,6 +2637,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. @@ -2498,8 +2800,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) @@ -2519,7 +2821,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, backend=self.backend, verbose=self._joblib_verbose)(delayed(parallel_delta_diff)(frame,self.ref,mask,size) for frame in flux)) else: d = [] kernel = [] @@ -2612,7 +2914,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: @@ -2622,32 +2924,34 @@ 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 + if self.verbose > 0: + 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 + if self.verbose > 0: + print(f'background pass 1 done ({_times["background (pass 1)"]:.2f}s)') if np.isnan(self.bkg).all(): - # check to see if the background worked raise ValueError('bkg all nans') # flux = strip_units(self.flux) @@ -2663,7 +2967,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') @@ -2675,7 +2979,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 @@ -2697,18 +3001,23 @@ 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 - + if self.verbose > 0: + print(f'alignment done ({_times["alignment"]:.2f}s)') + + _bad_frames = np.nansum(self.flux_raw,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') + print(f'images shifted ({_times["image shifting"]:.2f}s)') #if self.kernel_match: # self.kernel_matching(diff=False) @@ -2717,8 +3026,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 = self.flux / self.qe + + self.flux = self.flux_raw / self.qe if self.align: self.shift_images() @@ -2728,13 +3037,13 @@ 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 + self.flux[_bad_frames] = np.nan # subtract reference if self._ref_type.lower() == 'single': 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] @@ -2760,27 +3069,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 - if self.verbose > 0: - print('remade mask') - # background - if self.verbose > 0: - print('background') + print(f'diff setup done ({_times["difference imaging setup"]:.2f}s)') + self.bkg_orig = deepcopy(self.bkg) + 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._grad_bkg_clip() self.flux -= self.bkg _times['background (pass 2)'] = time.perf_counter() - _t + if self.verbose > 0: + print(f'background pass 2 done ({_times["background (pass 2)"]:.2f}s)') 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 + 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: @@ -2793,14 +3104,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: @@ -2817,6 +3133,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): """ @@ -2859,7 +3184,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]) @@ -3317,13 +3642,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: @@ -3455,7 +3780,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.') @@ -3463,7 +3788,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 @@ -3652,7 +3977,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') @@ -3675,7 +4000,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: @@ -3694,6 +4019,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()