Merge dev: HPC/parallelism fixes, PSF photometry, and tpf-dependency removal

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,

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

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()

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)
 

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