From d310765b1b742173c956e188e2019722b48c8810 Mon Sep 17 00:00:00 2001 From: Christina Hedges Date: Thu, 3 Feb 2022 13:30:05 -0800 Subject: [PATCH 1/8] :art: removing inplace centroids --- src/psfmachine/machine.py | 290 +++++++++++++++++++++++--------------- 1 file changed, 179 insertions(+), 111 deletions(-) diff --git a/src/psfmachine/machine.py b/src/psfmachine/machine.py index e70e633..9f445a3 100644 --- a/src/psfmachine/machine.py +++ b/src/psfmachine/machine.py @@ -123,6 +123,8 @@ def __init__( phi: numpy.ndarray Angle between pixel and source coordinates (polar coordinates), in units of radians + rough_mask: scipy.sparce.csr_matrix + Sparce mask matrix with pixels that are close to sources, simple round mask source_mask: scipy.sparce.csr_matrix Sparce mask matrix with pixels that contains flux from sources uncontaminated_source_mask: scipy.sparce.csr_matrix @@ -183,13 +185,16 @@ def __init__( self.nt = len(self.time) self.npixels = self.flux.shape[1] + self.source_flux_estimates = np.copy(np.asarray(self.sources.phot_g_mean_flux)) # Hardcoded: sparse implementation is efficient when nsourxes * npixels < 1e7 # (JMP profile this) # https://github.com/SSDataLab/psfmachine/pull/17#issuecomment-866382898 + self.ra_centroid, self.dec_centroid = np.zeros((2, self.nt)) * u.deg if self.nsources * self.npixels < 1e7: - self._create_delta_arrays() + self._update_delta_arrays() else: - self._create_delta_sparse_arrays() + self._update_delta_sparse_arrays() + self._get_rough_source_mask() @property def shape(self): @@ -198,49 +203,100 @@ def shape(self): def __repr__(self): return f"Machine (N sources, N times, N pixels): {self.shape}" - def _create_delta_arrays(self, centroid_offset=[0, 0]): + def dx(self, frame_indices="mean"): + """Delta RA, corrected for centroid shift""" + if isinstance(frame_indices, str): + if frame_indices == "mean": + frame_indices = np.where(self.time_mask)[0] + else: + frame_indices = np.atleast_1d(frame_indices) + if not sparse.issparse(self.dra): + return self.dra - self.ra_centroid[frame_indices].mean().value + else: + ra_offset = sparse.csr_matrix( + ( + np.repeat( + self.ra_centroid[frame_indices].mean().value, + self.dra.data.shape, + ), + (self.dra.nonzero()), + ), + shape=self.dra.shape, + dtype=float, + ) + return self.dra - ra_offset + + def dy(self, frame_indices="mean"): + """Delta Dec, corrected for centroid shift""" + if isinstance(frame_indices, str): + if frame_indices == "mean": + frame_indices = np.where(self.time_mask)[0] + else: + frame_indices = np.atleast_1d(frame_indices) + if not sparse.issparse(self.ddec): + return self.ddec - self.dec_centroid[frame_indices].mean().value + else: + dec_offset = sparse.csr_matrix( + ( + np.repeat( + self.dec_centroid[frame_indices].mean().value, + self.ddec.data.shape, + ), + (self.ddec.nonzero()), + ), + shape=self.ddec.shape, + dtype=float, + ) + return self.ddec - dec_offset + + def _update_delta_arrays(self, frame_indices="mean"): """ Creates dra, ddec, r and phi numpy ndarrays . Parameters ---------- - centroid_offset : list - Centroid offset for [ra, dec] to be included in dra and ddec computation. - Default is [0, 0]. + frame_indices : list or str + "mean" takes the mean of all the centroids in "time_mask" + """ # The distance in ra & dec from each source to each pixel # when centroid offset is 0 (i.e. first time creating arrays) create delta # arrays from scratch - if centroid_offset[0] == centroid_offset[1] == 0: + + if not hasattr(self, "dra"): self.dra, self.ddec = np.asarray( [ [ - self.ra - self.sources["ra"][idx] - centroid_offset[0], - self.dec - self.sources["dec"][idx] - centroid_offset[1], + self.ra - self.sources["ra"][idx], + self.dec - self.sources["dec"][idx], ] for idx in range(len(self.sources)) ] ).transpose(1, 0, 2) self.dra = self.dra * (u.deg) self.ddec = self.ddec * (u.deg) - # when offsets are != 0 (i.e. updating dra and ddec arrays) we just substract - # the ofsets avoiding the for loop - else: - self.dra -= centroid_offset[0] * u.deg - self.ddec -= centroid_offset[1] * u.deg # convertion to polar coordinates - self.r = np.hypot(self.dra, self.ddec).to("arcsec") - self.phi = np.arctan2(self.ddec, self.dra) + self.r = ( + np.hypot( + self.dx(frame_indices), + self.dy(frame_indices), + ) + * 3600 + ) + self.phi = np.arctan2( + self.dy(frame_indices), + self.dx(frame_indices), + ) - def _create_delta_sparse_arrays(self, dist_lim=40, centroid_offset=[0, 0]): + def _update_delta_sparse_arrays(self, frame_indices="mean", dist_lim=50): """ Creates dra, ddec, r and phi arrays as sparse arrays to be used for dense data, e.g. Kepler FFIs or cluster fields. Assuming that there is no flux information further than `dist_lim` for a given source, we only keep pixels within the `dist_lim`. dra, ddec, ra, and phi are unitless because they are `sparse.csr_matrix`. But - keep same scale as '_create_delta_arrays()'. + keep same scale as '_update_delta_arrays()'. dra and ddec in deg. r in arcseconds and phi in rads Parameters @@ -251,11 +307,11 @@ def _create_delta_sparse_arrays(self, dist_lim=40, centroid_offset=[0, 0]): Centroid offset for [ra, dec] to be included in dra and ddec computation. Default is [0, 0]. """ + if frame_indices == "mean": + frame_indices = np.where(self.time_mask)[0] # If not centroid offsets or centroid correction are larget than a pixel, # then we need to compute the sparse delta arrays from scratch - if (centroid_offset[0] == centroid_offset[1] == 0) or ( - np.maximum(*np.abs(centroid_offset)) > 4 / 3600 - ): + if not hasattr(self, "dra"): # iterate over sources to only keep pixels within dist_lim # this is inefficient, could be done in a tiled manner? only for squared data dra, ddec, sparse_mask = [], [], [] @@ -264,8 +320,8 @@ def _create_delta_sparse_arrays(self, dist_lim=40, centroid_offset=[0, 0]): desc="Creating delta arrays", disable=self.quiet, ): - dra_aux = self.ra - self.sources["ra"].iloc[i] - centroid_offset[0] - ddec_aux = self.dec - self.sources["dec"].iloc[i] - centroid_offset[1] + dra_aux = self.ra - self.sources["ra"].iloc[i] + ddec_aux = self.dec - self.sources["dec"].iloc[i] box_mask = sparse.csr_matrix( (np.abs(dra_aux) <= self.sparse_dist_lim.to("deg").value) & (np.abs(ddec_aux) <= self.sparse_dist_lim.to("deg").value) @@ -282,32 +338,18 @@ def _create_delta_sparse_arrays(self, dist_lim=40, centroid_offset=[0, 0]): sparse_mask.eliminate_zeros() # if centroid correction is less than 1 pixel, then we just update dra and ddec # sparse arrays arrays and r and phi. - else: - self.dra = self.dra - sparse.csr_matrix( - ( - np.repeat(centroid_offset[0], self.dra.data.shape), - (self.dra.nonzero()), - ), - shape=self.dra.shape, - dtype=float, - ) - self.ddec = self.ddec - sparse.csr_matrix( - ( - np.repeat(centroid_offset[1], self.ddec.data.shape), - (self.ddec.nonzero()), - ), - shape=self.ddec.shape, - dtype=float, - ) - sparse_mask = self.dra.astype(bool) + + sparse_mask = self.dra.astype(bool) # convertion to polar coordinates. We can't apply np.hypot or np.arctan2 to # sparse arrays. We keep track of non-zero index, do math in numpy space, # then rebuild r, phi as sparse. nnz_inds = sparse_mask.nonzero() # convert radial dist to arcseconds - r_vals = np.hypot(self.dra.data, self.ddec.data) * 3600 - phi_vals = np.arctan2(self.ddec.data, self.dra.data) + r_vals = ( + np.hypot(self.dx(frame_indices).data, self.dy(frame_indices).data) * 3600 + ) + phi_vals = np.arctan2(self.dy(frame_indices).data, self.dx(frame_indices).data) self.r = sparse.csr_matrix( (r_vals, (nnz_inds[0], nnz_inds[1])), shape=sparse_mask.shape, @@ -321,13 +363,12 @@ def _create_delta_sparse_arrays(self, dist_lim=40, centroid_offset=[0, 0]): del r_vals, phi_vals, nnz_inds, sparse_mask return - def _get_source_mask( + def _get_rough_source_mask( self, - upper_radius_limit=28.0, + upper_radius_limit=60.0, lower_radius_limit=4.5, upper_flux_limit=2e5, - lower_flux_limit=100, - correct_centroid_offset=True, + lower_flux_limit=3, plot=False, ): """Find the pixel mask that identifies pixels with contributions from ANY NUMBER of Sources @@ -351,12 +392,6 @@ def _get_source_mask( plot: bool Whether to show diagnostic plot. Default is False """ - - if not hasattr(self, "source_flux_estimates"): - # gaia estimate flux values per pixel to be used as flux priors - self.source_flux_estimates = np.copy( - np.asarray(self.sources.phot_g_mean_flux) - ) # We will use the radius a lot, this is for readibility # don't do if sparse array if isinstance(self.r, u.quantity.Quantity): @@ -484,51 +519,17 @@ def _get_source_mask( source_radius_limit < lower_radius_limit ] = lower_radius_limit - # Here we set the radius for each source. We add two pixels, to be generous - self.radius = source_radius_limit + 2 + # Here we set the radius for each source. We add six pixels, to be generous + self.radius = source_radius_limit + 6 # This sparse mask is one where where is ANY number of sources in a pixel if not isinstance(self.r, sparse.csr_matrix): - self.source_mask = sparse.csr_matrix(self.r.value < self.radius[:, None]) + self.rough_mask = sparse.csr_matrix(self.r.value < self.radius[:, None]) else: # for a sparse matrix doing < self.radius is not efficient, it # considers all zero values in the sparse matrix and set them to True. # this is a workaround to this problem. - self.source_mask = sparse_lessthan(r, self.radius) - - self._get_uncontaminated_pixel_mask() - - # Now we can update the r and phi estimates, allowing for a slight centroid - # calculate image centroids and correct dra,ddec for offset. - if correct_centroid_offset: - self._get_centroids() - # print(self.ra_centroid_avg.to("arcsec"), self.dec_centroid_avg.to("arcsec")) - # re-estimate dra, ddec with centroid shifts, check if sparse case applies. - # Hardcoded: sparse implementation is efficient when nsourxes * npixels < 1e7 - # (JMP profile this) - # https://github.com/SSDataLab/psfmachine/pull/17#issuecomment-866382898 - if self.nsources * self.npixels < 1e7: - self._create_delta_arrays( - centroid_offset=[ - self.ra_centroid_avg.value, - self.dec_centroid_avg.value, - ] - ) - else: - self._create_delta_sparse_arrays( - centroid_offset=[ - self.ra_centroid_avg.value, - self.dec_centroid_avg.value, - ] - ) - # if centroid offset id larger than 1" then we need to recalculate the - # source mask to include/reject correct pixels. - # this 1 arcsec limit only works for Kepler/K2 - if ( - np.abs(self.ra_centroid_avg.to("arcsec").value) > 1 - or np.abs(self.dec_centroid_avg.to("arcsec").value) > 1 - ): - self._get_source_mask(correct_centroid_offset=False) + self.rough_mask = sparse_lessthan(r, self.radius) if plot: k = np.isfinite(f_temp_mask) @@ -567,6 +568,61 @@ def _get_source_mask( xlabel=("log$_{10}$ Gaia Source Flux"), ) return fig + + def _get_source_mask( + self, + upper_radius_limit=28.0, + lower_radius_limit=4.5, + upper_flux_limit=2e5, + lower_flux_limit=100, + correct_centroid_offset=True, + plot=False, + ): + """Find the pixel mask that identifies pixels with contributions from ANY NUMBER of Sources + + Fits a simple polynomial model to the log of the pixel flux values, in radial dimension and source flux, + to find the optimum circular apertures for every source. + + Parameters + ---------- + upper_radius_limit: float + The radius limit at which we assume there is no flux from a source of any brightness (arcsec) + lower_radius_limit: float + The radius limit at which we assume there is flux from a source of any brightness (arcsec) + upper_flux_limit: float + The flux at which we assume as source is saturated + lower_flux_limit: float + The flux at which we assume a source is too faint to model + correct_centroid_offset: bool + Correct the dra, ddec arrays from centroid offsets. If centroid offsets are + larger than 1 arcsec, `source_mask` will be also updated. + plot: bool + Whether to show diagnostic plot. Default is False + """ + + self._get_rough_source_mask( + upper_radius_limit=upper_radius_limit, + lower_radius_limit=lower_radius_limit, + upper_flux_limit=upper_flux_limit, + lower_flux_limit=lower_flux_limit, + plot=plot, + ) + self._get_uncontaminated_pixel_mask() + self._get_centroids() + # Now we can update the r and phi estimates, allowing for a slight centroid + # calculate image centroids and correct dra,ddec for offset. + # re-estimate dra, ddec with centroid shifts, check if sparse case applies. + # Hardcoded: sparse implementation is efficient when nsourxes * npixels < 1e7 + # (JMP profile this) + # https://github.com/SSDataLab/psfmachine/pull/17#issuecomment-866382898 + if self.nsources * self.npixels < 1e7: + self._update_delta_arrays(frame_indices="mean") + else: + self._update_delta_sparse_arrays(frame_indices="mean") + # if centroid offset id larger than 1" then we need to recalculate the + # source mask to include/reject correct pixels. + # this 1 arcsec limit only works for Kepler/K2 + return def _get_uncontaminated_pixel_mask(self): @@ -589,8 +645,12 @@ def _get_uncontaminated_pixel_mask(self): return - # CH: We're not currently using this, but it might prove useful later so I will leave for now - def _get_centroids(self): + # def dra(self): + # if not hasattr(self, 'ra_centroid'): + # return self.uncontaminated_source_mask.multiply(self.dra).data + # else: + + def _get_centroids(self, plot=False): """ Find the ra and dec centroid of the image, at each time. """ @@ -600,21 +660,28 @@ def _get_centroids(self): dra_m = self.uncontaminated_source_mask.multiply(self.dra).data ddec_m = self.uncontaminated_source_mask.multiply(self.ddec).data for t in range(self.nt): - wgts = self.uncontaminated_source_mask.multiply( - np.sqrt(np.abs(self.flux[t])) - ).data - # mask out non finite values and background pixels - k = (np.isfinite(wgts)) & ( - self.uncontaminated_source_mask.multiply(self.flux[t]).data > 100 + wgts = ( + self.uncontaminated_source_mask.multiply(self.flux[t]) + .multiply(1 / self.source_flux_estimates[:, None]) + .data + ** 2 ) + # mask out non finite values and background pixels + k = (np.isfinite(wgts)) & (wgts > 0.01) self.ra_centroid[t] = np.average(dra_m[k], weights=wgts[k]) self.dec_centroid[t] = np.average(ddec_m[k], weights=wgts[k]) - del dra_m, ddec_m + if plot: + plt.figure() + plt.scatter( + dra_m, ddec_m, c=wgts ** 0.5, s=1, vmin=0, vmax=0.2, cmap="Greys" + ) + plt.scatter(dra_m[k], ddec_m[k], c=wgts[k] ** 0.5, s=1, vmin=0, vmax=0.2) + plt.scatter(self.ra_centroid[t], self.dec_centroid[t], c="r") + plt.gca().set_aspect("equal") + self.ra_centroid *= u.deg self.dec_centroid *= u.deg - self.ra_centroid_avg = self.ra_centroid.mean() - self.dec_centroid_avg = self.dec_centroid.mean() - + del dra_m, ddec_m return def _time_bin(self, npoints=200, downsample=False): @@ -859,8 +926,8 @@ def build_time_model(self, plot=False, downsample=False): self._whitened_time = time_original # not necessary to take value from Quantity to do .multiply() dx, dy = ( - self.uncontaminated_source_mask.multiply(self.dra), - self.uncontaminated_source_mask.multiply(self.ddec), + self.uncontaminated_source_mask.multiply(self.dx()), + self.uncontaminated_source_mask.multiply(self.dy()), ) dx = dx.data * u.deg.to(u.arcsecond) dy = dy.data * u.deg.to(u.arcsecond) @@ -965,8 +1032,8 @@ def plot_time_model(self): # not necessary to take value from Quantity to do .multiply() dx, dy = ( - self.uncontaminated_source_mask.multiply(self.dra), - self.uncontaminated_source_mask.multiply(self.ddec), + self.uncontaminated_source_mask.multiply(self.dx()), + self.uncontaminated_source_mask.multiply(self.dy()), ) dx = dx.data * u.deg.to(u.arcsecond) dy = dy.data * u.deg.to(u.arcsecond) @@ -1290,8 +1357,8 @@ def plot_shape_model(self, radius=20, frame_index="mean"): ) dx, dy = ( - self.uncontaminated_source_mask.multiply(self.dra), - self.uncontaminated_source_mask.multiply(self.ddec), + self.uncontaminated_source_mask.multiply(self.dx()), + self.uncontaminated_source_mask.multiply(self.dy()), ) dx = dx.data * u.deg.to(u.arcsecond) dy = dy.data * u.deg.to(u.arcsecond) @@ -1411,9 +1478,10 @@ def fit_model(self, fit_va=False): ) # not necessary to take value from Quantity to do .multiply() + # SHOULDN'T THIS BE UNCONTAMINATED SOURCE MASK? dx, dy = ( - self.source_mask.multiply(self.dra), - self.source_mask.multiply(self.ddec), + self.source_mask.multiply(self.dx()), + self.source_mask.multiply(self.dy()), ) dx = dx.data * u.deg.to(u.arcsecond) dy = dy.data * u.deg.to(u.arcsecond) From d87fdbde4adfc69151b23cf5b548fa876ceb2b31 Mon Sep 17 00:00:00 2001 From: Christina Hedges Date: Tue, 8 Mar 2022 11:05:18 -0800 Subject: [PATCH 2/8] :art: working on centroids --- src/psfmachine/machine.py | 560 ++++++++++++++++++++------------------ src/psfmachine/utils.py | 37 ++- 2 files changed, 322 insertions(+), 275 deletions(-) diff --git a/src/psfmachine/machine.py b/src/psfmachine/machine.py index 9f445a3..d470bab 100644 --- a/src/psfmachine/machine.py +++ b/src/psfmachine/machine.py @@ -167,12 +167,13 @@ def __init__( self.rmax = rmax self.cut_r = cut_r self.sparse_dist_lim = sparse_dist_lim * u.arcsecond - self.time_corrector = "polynomial" + self.time_corrector = "centroid" # if pos_corr, then we can smooth the vector with a Gaussian kernel of size # poscorr_filter_size, if this is < 0.5 -> no smoothing, default is 12 # beacause of 6hr-CDPP self.poscorr_filter_size = 12 - self.cartesian_knot_spacing = "sqrt" + self.cartesian_knot_spacing = "linear" + # disble tqdm prgress bar when running in HPC self.quiet = False @@ -189,12 +190,21 @@ def __init__( # Hardcoded: sparse implementation is efficient when nsourxes * npixels < 1e7 # (JMP profile this) # https://github.com/SSDataLab/psfmachine/pull/17#issuecomment-866382898 - self.ra_centroid, self.dec_centroid = np.zeros((2, self.nt)) * u.deg + self.ra_centroid, self.dec_centroid = np.zeros((2)) * u.deg if self.nsources * self.npixels < 1e7: self._update_delta_arrays() else: self._update_delta_sparse_arrays() self._get_rough_source_mask() + self._get_centroid() + if self.nsources * self.npixels < 1e7: + self._update_delta_arrays() + else: + self._update_delta_sparse_arrays() + # if self.nsources * self.npixels < 1e7: + # self._update_delta_arrays() + # else: + # self._update_delta_sparse_arrays() @property def shape(self): @@ -203,20 +213,16 @@ def shape(self): def __repr__(self): return f"Machine (N sources, N times, N pixels): {self.shape}" - def dx(self, frame_indices="mean"): + @property + def dx(self): """Delta RA, corrected for centroid shift""" - if isinstance(frame_indices, str): - if frame_indices == "mean": - frame_indices = np.where(self.time_mask)[0] - else: - frame_indices = np.atleast_1d(frame_indices) if not sparse.issparse(self.dra): - return self.dra - self.ra_centroid[frame_indices].mean().value + return self.dra - self.ra_centroid.value else: ra_offset = sparse.csr_matrix( ( np.repeat( - self.ra_centroid[frame_indices].mean().value, + self.ra_centroid.value, self.dra.data.shape, ), (self.dra.nonzero()), @@ -226,20 +232,16 @@ def dx(self, frame_indices="mean"): ) return self.dra - ra_offset - def dy(self, frame_indices="mean"): + @property + def dy(self): """Delta Dec, corrected for centroid shift""" - if isinstance(frame_indices, str): - if frame_indices == "mean": - frame_indices = np.where(self.time_mask)[0] - else: - frame_indices = np.atleast_1d(frame_indices) if not sparse.issparse(self.ddec): - return self.ddec - self.dec_centroid[frame_indices].mean().value + return self.ddec - self.dec_centroid.value else: dec_offset = sparse.csr_matrix( ( np.repeat( - self.dec_centroid[frame_indices].mean().value, + self.dec_centroid.value, self.ddec.data.shape, ), (self.ddec.nonzero()), @@ -279,14 +281,14 @@ def _update_delta_arrays(self, frame_indices="mean"): # convertion to polar coordinates self.r = ( np.hypot( - self.dx(frame_indices), - self.dy(frame_indices), + self.dx, + self.dy, ) * 3600 ) self.phi = np.arctan2( - self.dy(frame_indices), - self.dx(frame_indices), + self.dy, + self.dx, ) def _update_delta_sparse_arrays(self, frame_indices="mean", dist_lim=50): @@ -346,10 +348,8 @@ def _update_delta_sparse_arrays(self, frame_indices="mean", dist_lim=50): # then rebuild r, phi as sparse. nnz_inds = sparse_mask.nonzero() # convert radial dist to arcseconds - r_vals = ( - np.hypot(self.dx(frame_indices).data, self.dy(frame_indices).data) * 3600 - ) - phi_vals = np.arctan2(self.dy(frame_indices).data, self.dx(frame_indices).data) + r_vals = np.hypot(self.dx.data, self.dy.data) * 3600 + phi_vals = np.arctan2(self.dy.data, self.dx.data) self.r = sparse.csr_matrix( (r_vals, (nnz_inds[0], nnz_inds[1])), shape=sparse_mask.shape, @@ -569,46 +569,10 @@ def _get_rough_source_mask( ) return fig - def _get_source_mask( - self, - upper_radius_limit=28.0, - lower_radius_limit=4.5, - upper_flux_limit=2e5, - lower_flux_limit=100, - correct_centroid_offset=True, - plot=False, - ): - """Find the pixel mask that identifies pixels with contributions from ANY NUMBER of Sources - - Fits a simple polynomial model to the log of the pixel flux values, in radial dimension and source flux, - to find the optimum circular apertures for every source. + def _get_source_mask(self): - Parameters - ---------- - upper_radius_limit: float - The radius limit at which we assume there is no flux from a source of any brightness (arcsec) - lower_radius_limit: float - The radius limit at which we assume there is flux from a source of any brightness (arcsec) - upper_flux_limit: float - The flux at which we assume as source is saturated - lower_flux_limit: float - The flux at which we assume a source is too faint to model - correct_centroid_offset: bool - Correct the dra, ddec arrays from centroid offsets. If centroid offsets are - larger than 1 arcsec, `source_mask` will be also updated. - plot: bool - Whether to show diagnostic plot. Default is False - """ - - self._get_rough_source_mask( - upper_radius_limit=upper_radius_limit, - lower_radius_limit=lower_radius_limit, - upper_flux_limit=upper_flux_limit, - lower_flux_limit=lower_flux_limit, - plot=plot, - ) self._get_uncontaminated_pixel_mask() - self._get_centroids() + # self._get_centroids() # Now we can update the r and phi estimates, allowing for a slight centroid # calculate image centroids and correct dra,ddec for offset. # re-estimate dra, ddec with centroid shifts, check if sparse case applies. @@ -625,63 +589,91 @@ def _get_source_mask( return - def _get_uncontaminated_pixel_mask(self): - """ - creates a mask of shape nsources x npixels where targets are not contaminated. - This mask is used to select pixels to build the PSF model. - """ - - # This could be a property, but it is a pain to calculate on the fly, perhaps with lru_cache - self.uncontaminated_source_mask = self.source_mask.multiply( - np.asarray(self.source_mask.sum(axis=0) == 1)[0] - ).tocsr() - # have to remove leaked zeros - self.uncontaminated_source_mask.eliminate_zeros() - - # # reduce to good pixels - # self.uncontaminated_pixel_mask = sparse.csr_matrix( - # self.uncontaminated_source_mask.sum(axis=0) > 0 - # ) - - return - # def dra(self): # if not hasattr(self, 'ra_centroid'): # return self.uncontaminated_source_mask.multiply(self.dra).data # else: - def _get_centroids(self, plot=False): - """ - Find the ra and dec centroid of the image, at each time. - """ - # centroids are astropy quantities - self.ra_centroid = np.zeros(self.nt) - self.dec_centroid = np.zeros(self.nt) - dra_m = self.uncontaminated_source_mask.multiply(self.dra).data - ddec_m = self.uncontaminated_source_mask.multiply(self.ddec).data - for t in range(self.nt): - wgts = ( - self.uncontaminated_source_mask.multiply(self.flux[t]) - .multiply(1 / self.source_flux_estimates[:, None]) - .data - ** 2 - ) - # mask out non finite values and background pixels - k = (np.isfinite(wgts)) & (wgts > 0.01) - self.ra_centroid[t] = np.average(dra_m[k], weights=wgts[k]) - self.dec_centroid[t] = np.average(ddec_m[k], weights=wgts[k]) - if plot: - plt.figure() - plt.scatter( - dra_m, ddec_m, c=wgts ** 0.5, s=1, vmin=0, vmax=0.2, cmap="Greys" - ) - plt.scatter(dra_m[k], ddec_m[k], c=wgts[k] ** 0.5, s=1, vmin=0, vmax=0.2) - plt.scatter(self.ra_centroid[t], self.dec_centroid[t], c="r") - plt.gca().set_aspect("equal") + # def _get_centroids(self, plot=False): + # """ + # Find the ra and dec centroid of the image, at each time. + # """ + # # centroids are astropy quantities + # self.ra_centroid = np.zeros(self.nt) + # self.dec_centroid = np.zeros(self.nt) + # dra_m = self.uncontaminated_source_mask.multiply(self.dra).data + # ddec_m = self.uncontaminated_source_mask.multiply(self.ddec).data + # for t in range(self.nt): + # wgts = ( + # self.uncontaminated_source_mask.multiply(self.flux[t]) + # .multiply(1 / self.source_flux_estimates[:, None]) + # .data + # ** 2 + # ) + # # mask out non finite values and background pixels + # k = (np.isfinite(wgts)) & (wgts > 0.01) + # self.ra_centroid[t] = np.average(dra_m[k], weights=wgts[k]) + # self.dec_centroid[t] = np.average(ddec_m[k], weights=wgts[k]) + # if plot: + # plt.figure() + # plt.scatter( + # dra_m, ddec_m, c=wgts ** 0.5, s=1, vmin=0, vmax=0.2, cmap="Greys" + # ) + # plt.scatter(dra_m[k], ddec_m[k], c=wgts[k] ** 0.5, s=1, vmin=0, vmax=0.2) + # plt.scatter(self.ra_centroid[t], self.dec_centroid[t], c="r") + # plt.gca().set_aspect("equal") + # + # self.ra_centroid *= u.deg + # self.dec_centroid *= u.deg + # del dra_m, ddec_m + # return + + def _get_centroid(self, frame_indices="mean", binsize=0.0005): + if isinstance(frame_indices, str): + if frame_indices == "mean": + frame_indices = np.where(self.time_mask)[0] + else: + frame_indices = np.atleast_1d(frame_indices) + + x, y = ( + self.rough_mask.multiply(self.dra).data, + self.rough_mask.multiply(self.ddec).data, + ) + ar, X, Y = np.histogram2d( + x, + y, + (np.arange(-0.01, 0.01, binsize), np.arange(-0.01, 0.011, binsize)), + ) + X, Y = np.meshgrid(X, Y) + j = np.zeros(np.asarray(ar.T.shape) + 1, bool) + j[:-1, :-1] = ar.T != 0 - self.ra_centroid *= u.deg - self.dec_centroid *= u.deg - del dra_m, ddec_m + c = ( + self.rough_mask.multiply(self.flux[frame_indices].mean(axis=0)) + .multiply(1 / self.source_flux_estimates[:, None]) + .data + ) + k = c < 0.8 + x_k, y_k, c_k = x[k], y[k], c[k] + ar2 = np.asarray( + [ + np.nanmedian( + c_k[ + (x_k >= X1) + & (x_k < X1 + binsize) + & (y_k >= Y1) + & (y_k < Y1 + binsize) + ] + ) + for X1, Y1 in zip(X[j], Y[j]) + ] + ) + self.ra_centroid = ( + np.average(X[j] + binsize / 2, weights=np.nan_to_num(ar2 ** 3)) * u.deg + ) + self.dec_centroid = ( + np.average(Y[j] + binsize / 2, weights=np.nan_to_num(ar2 ** 3)) * u.deg + ) return def _time_bin(self, npoints=200, downsample=False): @@ -728,8 +720,8 @@ def _time_bin(self, npoints=200, downsample=False): mpc2 = np.nanmedian(self.pos_corr2, axis=0) else: # if usig centroids need to convert to pixels - mpc1 = self.ra_centroid.to("arcsec").value / 4 - mpc2 = self.dec_centroid.to("arcsec").value / 4 + mpc1 = self.ra_centroids.to("arcsec").value / 4 + mpc2 = self.dec_centroids.to("arcsec").value / 4 # find poscorr discontinuity in each axis grads1 = np.gradient(mpc1, self.time) @@ -879,6 +871,23 @@ def _time_bin(self, npoints=200, downsample=False): return ta, tm, fm_raw, fm, fem + def _time_matrix(self, mask): + dx, dy = ( + mask.multiply(self.dx), + mask.multiply(self.dy), + ) + dx = dx.data * u.deg.to(u.arcsecond) + dy = dy.data * u.deg.to(u.arcsecond) + + A_c = _make_A_cartesian( + dx, + dy, + n_knots=self.n_time_knots, + radius=self.time_radius, + spacing=self.cartesian_knot_spacing, + ) + return A_c + def build_time_model(self, plot=False, downsample=False): """ Builds a time model that moves the PRF model to account for the scene movement @@ -924,29 +933,20 @@ def build_time_model(self, plot=False, downsample=False): ) = self._time_bin(npoints=self.n_time_points, downsample=downsample) self._whitened_time = time_original - # not necessary to take value from Quantity to do .multiply() - dx, dy = ( - self.uncontaminated_source_mask.multiply(self.dx()), - self.uncontaminated_source_mask.multiply(self.dy()), + A_c = self._time_matrix(self.uncontaminated_source_mask) + A2 = sparse.vstack( + [A_c] * time_binned.shape[0], + format="csr", ) - dx = dx.data * u.deg.to(u.arcsecond) - dy = dy.data * u.deg.to(u.arcsecond) - A_c = _make_A_cartesian( - dx, - dy, - n_knots=self.n_time_knots, - radius=self.time_radius, - spacing=self.cartesian_knot_spacing, - ) - A2 = sparse.vstack([A_c] * time_binned.shape[0], format="csr") - # Cartesian spline with time dependence if hasattr(self, "pos_corr1") and self.time_corrector in [ "pos_corr", "centroid", ]: # Cartesian spline with poscor dependence - A3 = _combine_A(A2, poscorr=[poscorr1_binned, poscorr2_binned]) + A3 = _combine_A( + A2, poscorr=[poscorr1_binned, poscorr2_binned], time=time_binned + ) else: # Cartesian spline with time dependence A3 = _combine_A(A2, time=time_binned) @@ -1030,22 +1030,11 @@ def plot_time_model(self): flux_err_binned, ) = self._time_bin(npoints=self.n_time_points) - # not necessary to take value from Quantity to do .multiply() - dx, dy = ( - self.uncontaminated_source_mask.multiply(self.dx()), - self.uncontaminated_source_mask.multiply(self.dy()), + A_c = self._time_matrix(self.uncontaminated_source_mask) + A2 = sparse.vstack( + [A_c] * time_binned.shape[0], + format="csr", ) - dx = dx.data * u.deg.to(u.arcsecond) - dy = dy.data * u.deg.to(u.arcsecond) - - A_c = _make_A_cartesian( - dx, - dy, - n_knots=self.n_time_knots, - radius=self.time_radius, - spacing=self.cartesian_knot_spacing, - ) - A2 = sparse.vstack([A_c] * time_binned.shape[0], format="csr") # Cartesian spline with time dependence # Cartesian spline with time dependence if hasattr(self, "pos_corr1") and self.time_corrector in [ @@ -1053,11 +1042,17 @@ def plot_time_model(self): "centroid", ]: # Cartesian spline with poscor dependence - A3 = _combine_A(A2, poscorr=[poscorr1_binned, poscorr2_binned]) + A3 = _combine_A( + A2, poscorr=[poscorr1_binned, poscorr2_binned], time=time_binned + ) else: # Cartesian spline with time dependence A3 = _combine_A(A2, time=time_binned) + dx, dy = ( + self.uncontaminated_source_mask.multiply(self.dx).data * 3600, + self.uncontaminated_source_mask.multiply(self.dy).data * 3600, + ) model = A3.dot(self.time_model_w).reshape(flux_binned.shape) + 1 fig, ax = plt.subplots(2, 2, figsize=(7, 6), facecolor="w") k1 = self._time_masked.reshape(flux_binned.shape)[0] @@ -1115,7 +1110,7 @@ def plot_time_model(self): return fig def build_shape_model( - self, plot=False, flux_cut_off=1, frame_index="mean", **kwargs + self, plot=False, flux_cut_off=1, frame_index="max", **kwargs ): """ Builds a sparse model matrix of shape nsources x npixels to be used when @@ -1137,39 +1132,43 @@ def build_shape_model( # Mask of shape nsources x number of pixels, one where flux from a # source exists - if not hasattr(self, "source_mask"): - self._get_source_mask(**kwargs) + if not hasattr(self, "uncontaminated_source_mask"): + mask = self.rough_mask + else: + mask = self.uncontaminated_source_mask + # self._get_source_mask(**kwargs) # Mask of shape npixels (maybe by nt) where not saturated, not faint, # not contaminated etc - self._get_uncontaminated_pixel_mask() + # self._get_uncontaminated_pixel_mask() # for iter in range(niters): flux_estimates = self.source_flux_estimates[:, None] if frame_index == "mean": f = (self.flux[self.time_mask]).mean(axis=0) - elif isinstance(frame_index, int): + elif frame_index == "max": + f = (self.flux[self.time_mask]).max(axis=0) + elif isinstance(frame_index, (int, np.int64, np.int32, np.int16)): f = self.flux[frame_index] + else: + raise ValueError(f"frame_index {frame_index} not valid") # f, fe = (self.flux[self.time_mask]).mean(axis=0), ( # (self.flux_err[self.time_mask] ** 2).sum(axis=0) ** 0.5 # ) / (self.nt) mean_f = np.log10( - self.uncontaminated_source_mask.astype(float) - .multiply(f) - .multiply(1 / flux_estimates) - .data + mask.astype(float).multiply(f).multiply(1 / flux_estimates).data ) # Actual Kepler errors cause all sorts of instability # mean_f_err = ( - # self.uncontaminated_source_mask.astype(float) + # mask.astype(float) # .multiply(fe / (f * np.log(10))) # .multiply(1 / flux_estimates) # .data # ) # We only need these weights for the wings, so we'll use poisson noise mean_f_err = ( - self.uncontaminated_source_mask.astype(float) + mask.astype(float) .multiply((f ** 0.5) / (f * np.log(10))) .multiply(1 / flux_estimates) .data @@ -1177,8 +1176,8 @@ def build_shape_model( mean_f_err.data = np.abs(mean_f_err.data) # take value from Quantity is not necessary - phi_b = self.uncontaminated_source_mask.multiply(self.phi).data - r_b = self.uncontaminated_source_mask.multiply(self.r).data + phi_b = mask.multiply(self.phi).data + r_b = mask.multiply(self.r).data # build a design matrix A with b-splines basis in radius and angle axis. A = _make_A_polar( @@ -1191,7 +1190,7 @@ def build_shape_model( n_phi_knots=self.n_phi_knots, ) prior_sigma = np.ones(A.shape[1]) * 10 - prior_mu = np.zeros(A.shape[1]) - 10 + prior_mu = np.zeros(A.shape[1]) - 100 nan_mask = np.isfinite(mean_f.ravel()) @@ -1206,7 +1205,7 @@ def build_shape_model( errors=True, ) - bad = sigma_clip(mean_f.ravel() - A.dot(psf_w), sigma=5).mask + bad = sigma_clip((mean_f.ravel() - A.dot(psf_w)), sigma=5).mask psf_w, psf_w_err = solve_linear_model( A, @@ -1223,91 +1222,103 @@ def build_shape_model( # We then build the same design matrix for all pixels with flux self._get_mean_model() + # self._update_source_mask_remove_bkg_pixels() # remove background pixels and recreate mean model - self._update_source_mask_remove_bkg_pixels( - flux_cut_off=flux_cut_off, frame_index=frame_index - ) + # self._update_source_mask_remove_bkg_pixels( + # flux_cut_off=flux_cut_off, frame_index=frame_index + # ) if plot: - return self.plot_shape_model(frame_index=frame_index) + return self.plot_shape_model(mask, frame_index=frame_index) return - def _update_source_mask_remove_bkg_pixels(self, flux_cut_off=1, frame_index="mean"): - """ - Update the `source_mask` to remove pixels that do not contribuite to the PRF - shape. - First, re-estimate the source flux usign the precomputed `mean_model`. - This re-estimation is used to remove sources with bad prediction and update - the `source_mask` by removing background pixels that do not contribuite to - the PRF shape. - Pixels with normalized flux > `flux_cut_off` are kept. - - Parameters - ---------- - flux_cut_off : float - Lower limit for the normalized flux predicted from the mean model. - frame_index : string or int - The frame index to be used, if "mean" then use the - mean value across time - """ - - # Re-estimate source flux - # ----- - prior_mu = self.source_flux_estimates - prior_sigma = ( - np.ones(self.mean_model.shape[0]) * 10 * self.source_flux_estimates - ) + # + # def _update_source_mask_remove_bkg_pixels(self, flux_cut_off=1, frame_index="max"): + # """ + # Update the `source_mask` to remove pixels that do not contribuite to the PRF + # shape. + # First, re-estimate the source flux usign the precomputed `mean_model`. + # This re-estimation is used to remove sources with bad prediction and update + # the `source_mask` by removing background pixels that do not contribuite to + # the PRF shape. + # Pixels with normalized flux > `flux_cut_off` are kept. + # + # Parameters + # ---------- + # flux_cut_off : float + # Lower limit for the normalized flux predicted from the mean model. + # frame_index : string or int + # The frame index to be used, if "mean" then use the + # mean value across time + # """ + # + # # Re-estimate source flux + # # ----- + # prior_mu = self.source_flux_estimates + # prior_sigma = ( + # np.ones(self.mean_model.shape[0]) * 10 * self.source_flux_estimates + # ) + # + # if frame_index == "mean": + # f, fe = (self.flux).mean(axis=0), ( + # (self.flux_err ** 2).sum(axis=0) ** 0.5 + # ) / (self.nt) + # if frame_index == "max": + # f, fe = (self.flux).max(axis=0), ( + # (self.flux_err ** 2).sum(axis=0) ** 0.5 + # ) / (self.nt) + # elif isinstance(frame_index, int): + # f, fe = self.flux[frame_index], self.flux_err[frame_index] + # + # X = self.mean_model.copy() + # X = X.T + # + # sigma_w_inv = X.T.dot(X.multiply(1 / fe[:, None] ** 2)).toarray() + # sigma_w_inv += np.diag(1 / (prior_sigma ** 2)) + # B = X.T.dot((f / fe ** 2)) + # B += prior_mu / (prior_sigma ** 2) + # ws = np.linalg.solve(sigma_w_inv, B) + # werrs = np.linalg.inv(sigma_w_inv).diagonal() ** 0.5 + # + # # ----- + # + # # Rebuild source mask + # ok = np.abs(ws - self.source_flux_estimates) / werrs > 3 + # ok &= ((ws / self.source_flux_estimates) < 10) & ( + # (self.source_flux_estimates / ws) < 10 + # ) + # ok &= ws > 10 + # ok &= werrs > 0 + # + # self.source_flux_estimates[ok] = ws[ok] + # + # self.source_mask = ( + # self.mean_model.multiply( + # self.mean_model.T.dot(self.source_flux_estimates) + # ).tocsr() + # > flux_cut_off + # ) + # + # # Recreate uncontaminated mask + # self._get_uncontaminated_pixel_mask() + # # self.uncontaminated_source_mask = self.uncontaminated_source_mask.multiply( + # # (self.mean_model.max(axis=1) < 1) + # # ) + # + # # Recreate mean model! + # self._get_mean_model() + + def _get_mean_model(self, relative_flux_limit=0.01, absolute_flux_limit=1): - if frame_index == "mean": - f, fe = (self.flux).mean(axis=0), ( - (self.flux_err ** 2).sum(axis=0) ** 0.5 - ) / (self.nt) - elif isinstance(frame_index, int): - f, fe = self.flux[frame_index], self.flux_err[frame_index] - - X = self.mean_model.copy() - X = X.T - - sigma_w_inv = X.T.dot(X.multiply(1 / fe[:, None] ** 2)).toarray() - sigma_w_inv += np.diag(1 / (prior_sigma ** 2)) - B = X.T.dot((f / fe ** 2)) - B += prior_mu / (prior_sigma ** 2) - ws = np.linalg.solve(sigma_w_inv, B) - werrs = np.linalg.inv(sigma_w_inv).diagonal() ** 0.5 - - # ----- - - # Rebuild source mask - ok = np.abs(ws - self.source_flux_estimates) / werrs > 3 - ok &= ((ws / self.source_flux_estimates) < 10) & ( - (self.source_flux_estimates / ws) < 10 - ) - ok &= ws > 10 - ok &= werrs > 0 - - self.source_flux_estimates[ok] = ws[ok] - - self.source_mask = ( - self.mean_model.multiply( - self.mean_model.T.dot(self.source_flux_estimates) - ).tocsr() - > flux_cut_off - ) - - # Recreate uncontaminated mask - self._get_uncontaminated_pixel_mask() - # self.uncontaminated_source_mask = self.uncontaminated_source_mask.multiply( - # (self.mean_model.max(axis=1) < 1) - # ) - - # Recreate mean model! - self._get_mean_model() + if not hasattr(self, "source_mask"): + mask = self.rough_mask + else: + mask = self.source_mask - def _get_mean_model(self): """Convenience function to make the scene model""" Ap = _make_A_polar( - self.source_mask.multiply(self.phi).data, - self.source_mask.multiply(self.r).data, + mask.multiply(self.phi).data, + mask.multiply(self.r).data, rmin=self.rmin, rmax=self.rmax, cut_r=self.cut_r, @@ -1319,11 +1330,17 @@ def _get_mean_model(self): mean_model = sparse.csr_matrix(self.r.shape) m = 10 ** Ap.dot(self.psf_w) m[~np.isfinite(m)] = 0 - mean_model[self.source_mask] = m + mean_model[mask] = m + mean_model = mean_model.multiply(mean_model > relative_flux_limit) + + mean_model = mean_model.multiply( + mean_model.multiply(self.source_flux_estimates[:, None]) + > absolute_flux_limit + ) mean_model.eliminate_zeros() self.mean_model = mean_model - def plot_shape_model(self, radius=20, frame_index="mean"): + def plot_shape_model(self, mask, radius=20, frame_index="max"): """ Diagnostic plot of shape model. @@ -1343,22 +1360,29 @@ def plot_shape_model(self, radius=20, frame_index="mean"): if frame_index == "mean": mean_f = np.log10( - self.uncontaminated_source_mask.astype(float) + mask.astype(float) .multiply(self.flux[self.time_mask].mean(axis=0)) .multiply(1 / self.source_flux_estimates[:, None]) .data ) + if frame_index == "max": + mean_f = np.log10( + mask.astype(float) + .multiply(self.flux[self.time_mask].max(axis=0)) + .multiply(1 / self.source_flux_estimates[:, None]) + .data + ) elif isinstance(frame_index, int): mean_f = np.log10( - self.uncontaminated_source_mask.astype(float) + mask.astype(float) .multiply(self.flux[frame_index]) .multiply(1 / self.source_flux_estimates[:, None]) .data ) dx, dy = ( - self.uncontaminated_source_mask.multiply(self.dx()), - self.uncontaminated_source_mask.multiply(self.dy()), + mask.multiply(self.dx), + mask.multiply(self.dy), ) dx = dx.data * u.deg.to(u.arcsecond) dy = dy.data * u.deg.to(u.arcsecond) @@ -1459,6 +1483,17 @@ def fit_model(self, fit_va=False): Fitting model accounting for velocity aberration. If `True`, then a time model has to be built previously with `build_time_model`. """ + + self._n_time_components = [ + 6 + if self.time_corrector + in [ + "pos_corr", + "centroid", + ] + else 4 + ][0] + prior_mu = self.source_flux_estimates # np.zeros(A.shape[1]) prior_sigma = ( np.ones(self.mean_model.shape[0]) @@ -1478,22 +1513,8 @@ def fit_model(self, fit_va=False): ) # not necessary to take value from Quantity to do .multiply() - # SHOULDN'T THIS BE UNCONTAMINATED SOURCE MASK? - dx, dy = ( - self.source_mask.multiply(self.dx()), - self.source_mask.multiply(self.dy()), - ) - dx = dx.data * u.deg.to(u.arcsecond) - dy = dy.data * u.deg.to(u.arcsecond) - - A_cp = _make_A_cartesian( - dx, - dy, - n_knots=self.n_time_knots, - radius=self.time_radius, - spacing=self.cartesian_knot_spacing, - ) - A_cp3 = sparse.hstack([A_cp, A_cp, A_cp, A_cp], format="csr") + A_cp = self._time_matrix(self.mean_model != 0) + A_cp3 = sparse.hstack([A_cp] * self._n_time_components, format="csr") self.ws_va = np.zeros((self.nt, self.mean_model.shape[0])) self.werrs_va = np.zeros((self.nt, self.mean_model.shape[0])) @@ -1526,18 +1547,23 @@ def fit_model(self, fit_va=False): np.array( [ 1, + self._whitened_time[tdx], + self._whitened_time[tdx] ** 2, self.pos_corr1_smooth[tdx], self.pos_corr2_smooth[tdx], self.pos_corr1_smooth[tdx] * self.pos_corr2_smooth[tdx], ] )[:, None] - * np.ones(A_cp3.shape[1] // 4) + * np.ones(A_cp3.shape[1] // self._n_time_components) ) else: # use time t_mult = np.hstack( - (self._whitened_time[tdx] ** np.arange(4))[:, None] - * np.ones(A_cp3.shape[1] // 4) + ( + self._whitened_time[tdx] + ** np.arange(self._n_time_components) + )[:, None] + * np.ones(A_cp3.shape[1] // self._n_time_components) ) X.data *= A_cp3.multiply(t_mult).dot(self.time_model_w) + 1 X = X.T @@ -1578,7 +1604,7 @@ def fit_model(self, fit_va=False): self.werrs[tdx] = np.linalg.inv(sigma_w_inv).diagonal() ** 0.5 self.model_flux[tdx] = X.dot(self.ws[tdx]) - nodata = np.asarray(self.source_mask.sum(axis=1))[:, 0] == 0 + nodata = np.asarray(mask.sum(axis=1))[:, 0] == 0 # These sources are poorly estimated nodata |= (self.mean_model.max(axis=1) > 1).toarray()[:, 0] self.ws[:, nodata] *= np.nan diff --git a/src/psfmachine/utils.py b/src/psfmachine/utils.py index 7374f72..c01524b 100644 --- a/src/psfmachine/utils.py +++ b/src/psfmachine/utils.py @@ -197,21 +197,26 @@ def _make_A_polar(phi, r, cut_r=6, rmin=1, rmax=18, n_r_knots=12, n_phi_knots=15 def _make_A_cartesian(x, y, n_knots=10, radius=3.0, spacing="sqrt"): + # Must be odd + n_time_knots = n_knots if n_knots % 2 == 1 else n_knots + 1 if spacing == "sqrt": - x_knots = np.linspace(-np.sqrt(radius), np.sqrt(radius), n_knots) + x_knots = np.linspace(-np.sqrt(radius), np.sqrt(radius), n_time_knots) x_knots = np.sign(x_knots) * x_knots ** 2 else: - x_knots = np.linspace(-radius, radius, n_knots) + x_knots = np.linspace(-radius, radius, n_time_knots) x_spline = sparse.csr_matrix( np.asarray( dmatrix( "bs(x, knots=knots, degree=3, include_intercept=True)", - {"x": list(x), "knots": x_knots}, + { + "x": list(np.hstack([x_knots.min(), x, x_knots.max()])), + "knots": x_knots, + }, ) - ) + )[1:-1] ) if spacing == "sqrt": - y_knots = np.linspace(-np.sqrt(radius), np.sqrt(radius), n_knots) + y_knots = np.linspace(-np.sqrt(radius), np.sqrt(radius), n_time_knots) y_knots = np.sign(y_knots) * y_knots ** 2 else: y_knots = np.linspace(-radius, radius, n_knots) @@ -219,9 +224,12 @@ def _make_A_cartesian(x, y, n_knots=10, radius=3.0, spacing="sqrt"): np.asarray( dmatrix( "bs(x, knots=knots, degree=3, include_intercept=True)", - {"x": list(y), "knots": y_knots}, + { + "x": list(np.hstack([y_knots.min(), y, y_knots.max()])), + "knots": y_knots, + }, ) - ) + )[1:-1] ) X = sparse.hstack( [x_spline.multiply(y_spline[:, idx]) for idx in range(y_spline.shape[1])], @@ -424,9 +432,11 @@ def _combine_A(A, poscorr=None, time=None): A2 = sparse.hstack( [ A, + A.multiply(time.ravel()[:, None]), + A.multiply(time.ravel()[:, None] ** 2), A.multiply(poscorr[0].ravel()[:, None]), A.multiply(poscorr[1].ravel()[:, None]), - A.multiply((poscorr[0] * poscorr[1]).ravel()[:, None]), + A.multiply((poscorr[0].ravel() * poscorr[1].ravel())[:, None]), ], format="csr", ) @@ -443,3 +453,14 @@ def _combine_A(A, poscorr=None, time=None): format="csr", ) return A2 + + +def _find_uncontaminated_pixels(mask): + """ + creates a mask of shape nsources x npixels where targets are not contaminated. + This mask is used to select pixels to build the PSF model. + """ + + new_mask = mask.multiply(np.asarray(mask.sum(axis=0) == 1)[0]).tocsr() + new_mask.eliminate_zeros() + return new_mask From 7b3d1a62e643088c258233d4b0c4ed41d51e8ab6 Mon Sep 17 00:00:00 2001 From: Christina Hedges Date: Thu, 17 Mar 2022 14:02:07 -0700 Subject: [PATCH 3/8] tweaks for tess --- src/psfmachine/ffi.py | 2 +- src/psfmachine/machine.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/psfmachine/ffi.py b/src/psfmachine/ffi.py index ba0183e..b1cc119 100644 --- a/src/psfmachine/ffi.py +++ b/src/psfmachine/ffi.py @@ -219,7 +219,7 @@ def from_file( ra, dec, wcs, - magnitude_limit=18 if metadata["TELESCOP"] == "Kepler" else 15, + magnitude_limit=18 if metadata["TELESCOP"] == "Kepler" else 17, epoch=time.jyear.mean(), ngrid=ngrid, dr=3, diff --git a/src/psfmachine/machine.py b/src/psfmachine/machine.py index f02dfde..4a9bae7 100644 --- a/src/psfmachine/machine.py +++ b/src/psfmachine/machine.py @@ -1062,7 +1062,6 @@ def build_shape_model( **kwargs Keyword arguments to be passed to `_get_source_mask()` """ - # gaia estimate flux values per pixel to be used as flux priors self.source_flux_estimates = np.copy(np.asarray(self.sources.phot_g_mean_flux)) From bfec53720ed7982a21e72af1e7662f6e5190f4d9 Mon Sep 17 00:00:00 2001 From: Christina Hedges Date: Fri, 18 Mar 2022 14:55:20 -0700 Subject: [PATCH 4/8] :recycle: cleaning up --- src/psfmachine/ffi.py | 147 +----------------------------- src/psfmachine/machine.py | 186 -------------------------------------- 2 files changed, 2 insertions(+), 331 deletions(-) diff --git a/src/psfmachine/ffi.py b/src/psfmachine/ffi.py index e055eb0..3221afa 100644 --- a/src/psfmachine/ffi.py +++ b/src/psfmachine/ffi.py @@ -118,16 +118,7 @@ def __init__( ).mask ).ravel() self._remove_background(mask=mask) - - # remove nan pixels - # valid_pix = np.isfinite(self.flux).sum(axis=0).astype(bool) - # self.flux = self.flux[:, valid_pix] - # self.flux_err = self.flux_err[:, valid_pix] - # self.ra = self.ra[valid_pix] - # self.dec = self.dec[valid_pix] - # self.row = self.row[valid_pix] - # self.column = self.column[valid_pix] - + # init `machine` object super().__init__( time, @@ -209,20 +200,7 @@ def from_file( row, metadata, ) = _load_file(fname, extension=extension) - # create cutouts if asked - # if cutout_size is not None: - # flux, flux_err, ra, dec, column, row = _do_image_cutout( - # flux, - # flux_err, - # ra, - # dec, - # column, - # row, - # cutout_size=cutout_size, - # cutout_origin=cutout_origin, - # ) - # hardcoded: the grid size to do the Gaia tiled query. This is different for - # cutouts and full channel. TESS and Kepler also need different grid sizes. + if metadata["TELESCOP"] == "Kepler": ngrid = (2, 2) if flux.shape[1] <= 500 else (4, 4) else: @@ -1067,127 +1045,6 @@ def _get_sources(ra, dec, wcs, img_limits=[[0, 0], [0, 0]], square=True, **kwarg return sources -# -# def _do_image_cutout( -# flux, flux_err, ra, dec, column, row, cutout_size=100, cutout_origin=[0, 0] -# ): -# """ -# Creates a cutout of the full image. Return data arrays corresponding to the cutout. -# -# Parameters -# ---------- -# flux : numpy.ndarray -# Data array with Flux values, correspond to full size image. -# flux_err : numpy.ndarray -# Data array with Flux errors values, correspond to full size image. -# ra : numpy.ndarray -# Data array with RA values, correspond to full size image. -# dec : numpy.ndarray -# Data array with Dec values, correspond to full size image. -# column : numpy.ndarray -# Data array with pixel column values, correspond to full size image. -# row : numpy.ndarray -# Data array with pixel raw values, correspond to full size image. -# cutout_size : int -# Size in pixels of the cutout, assumedto be squared. Default is 100. -# cutout_origin : tuple of ints -# Origin of the cutout following matrix indexing. Default is [0 ,0]. -# -# Returns -# ------- -# flux : numpy.ndarray -# Data array with Flux values of the cutout. -# flux_err : numpy.ndarray -# Data array with Flux errors values of the cutout. -# ra : numpy.ndarray -# Data array with RA values of the cutout. -# dec : numpy.ndarray -# Data array with Dec values of the cutout. -# column : numpy.ndarray -# Data array with pixel column values of the cutout. -# row : numpy.ndarray -# Data array with pixel raw values of the cutout. -# """ -# if (cutout_size + cutout_origin[0] <= flux.shape[1]) and ( -# cutout_size + cutout_origin[1] <= flux.shape[2] -# ): -# column = column[ -# cutout_origin[0] : cutout_origin[0] + cutout_size, -# cutout_origin[1] : cutout_origin[1] + cutout_size, -# ] -# row = row[ -# cutout_origin[0] : cutout_origin[0] + cutout_size, -# cutout_origin[1] : cutout_origin[1] + cutout_size, -# ] -# flux = flux[ -# :, -# cutout_origin[0] : cutout_origin[0] + cutout_size, -# cutout_origin[1] : cutout_origin[1] + cutout_size, -# ] -# flux_err = flux_err[ -# :, -# cutout_origin[0] : cutout_origin[0] + cutout_size, -# cutout_origin[1] : cutout_origin[1] + cutout_size, -# ] -# ra = ra[ -# cutout_origin[0] : cutout_origin[0] + cutout_size, -# cutout_origin[1] : cutout_origin[1] + cutout_size, -# ] -# dec = dec[ -# cutout_origin[0] : cutout_origin[0] + cutout_size, -# cutout_origin[1] : cutout_origin[1] + cutout_size, -# ] -# else: -# raise ValueError("Cutout size is larger than image shape ", flux.shape) -# -# return flux, flux_err, ra, dec, column, row - - -# -# def _remove_overscan(telescope, imgs): -# """ -# Removes overscan of the CCD. Return the image data with overscan columns and rows -# removed, also return 2D data arrays with pixel columns and row values. -# -# Parameters -# ---------- -# telescope : string -# Name of the telescope. -# imgs : numpy.ndarray -# Array of 2D images to. Has shape of [n_times, image_height, image_width]. -# -# Returns -# ------- -# row_2d : numpy.ndarray -# Data array with pixel row values -# col_2d : numpy.ndarray -# Data array with pixel column values -# flux_2d : numpy.ndarray -# Data array with flux values -# """ -# if telescope == "Kepler": -# # CCD overscan for Kepler -# r_min = 20 -# r_max = 1044 -# c_min = 12 -# c_max = 1112 -# elif telescope == "TESS": -# # CCD overscan for TESS -# r_min = 0 -# r_max = 2048 -# c_min = 45 -# c_max = 2093 -# else: -# raise TypeError("File is not from Kepler or TESS mission") -# # remove overscan -# row_2d, col_2d = np.mgrid[: imgs[0].shape[0], : imgs[0].shape[1]] -# col_2d = col_2d[r_min:r_max, c_min:c_max] -# row_2d = row_2d[r_min:r_max, c_min:c_max] -# flux_2d = imgs[:, r_min:r_max, c_min:c_max] -# -# return row_2d, col_2d, flux_2d - - def _compute_coordinate_offset(ra, dec, flux, sources, plot=True): """ Compute coordinate offsets if the RA Dec of objects in source catalog don't align diff --git a/src/psfmachine/machine.py b/src/psfmachine/machine.py index 4773ce5..05d562d 100644 --- a/src/psfmachine/machine.py +++ b/src/psfmachine/machine.py @@ -425,117 +425,7 @@ def _get_source_mask(self, source_flux_limit=1): self.source_mask ) return - # self.radius = ( - # rbins[k][ - # np.where( - # np.cumsum(10 ** fbins[k]) / np.nansum(10 ** fbins[k]) > 0.95 - # )[0][0] - # ] - # + 5 - # ) - # if not isinstance(self.r, sparse.csr_matrix): - # self.source_mask = sparse.csr_matrix(self.r < self.radius[:, None] * 5) - # else: - # self.source_mask = sparse_lessthan(self.r, self.radius * 5) - # - # m = self.source_mask.multiply(self.r).copy() - # # scipy sparse does -not- like that we are setting the `data` attribute - # m.data = 10 ** np.polyval(l, self.source_mask.multiply(self.r).data) - # m = m.multiply(self.source_flux_estimates[:, None]) - # self.source_mask = m > source_flux_limit - # # if not isinstance(self.r, sparse.csr_matrix): - # # self.rough_mask = sparse.csr_matrix(self.r < self.radius[:, None]) - # # else: - # # self.rough_mask = sparse_lessthan(self.r, self.radius) - # self.source_mask.eliminate_zeros() - # self.uncontaminated_source_mask = _find_uncontaminated_pixels(self.source_mask) - # def dra(self): - # if not hasattr(self, 'ra_centroid'): - # return self.uncontaminated_source_mask.multiply(self.dra).data - # else: - - # def _get_centroids(self, plot=False): - # """ - # Find the ra and dec centroid of the image, at each time. - # """ - # # centroids are astropy quantities - # self.ra_centroid = np.zeros(self.nt) - # self.dec_centroid = np.zeros(self.nt) - # dra_m = self.uncontaminated_source_mask.multiply(self.dra).data - # ddec_m = self.uncontaminated_source_mask.multiply(self.ddec).data - # for t in range(self.nt): - # wgts = ( - # self.uncontaminated_source_mask.multiply(self.flux[t]) - # .multiply(1 / self.source_flux_estimates[:, None]) - # .data - # ** 2 - # ) - # # mask out non finite values and background pixels - # k = (np.isfinite(wgts)) & (wgts > 0.01) - # self.ra_centroid[t] = np.average(dra_m[k], weights=wgts[k]) - # self.dec_centroid[t] = np.average(ddec_m[k], weights=wgts[k]) - # if plot: - # plt.figure() - # plt.scatter( - # dra_m, ddec_m, c=wgts ** 0.5, s=1, vmin=0, vmax=0.2, cmap="Greys" - # ) - # plt.scatter(dra_m[k], ddec_m[k], c=wgts[k] ** 0.5, s=1, vmin=0, vmax=0.2) - # plt.scatter(self.ra_centroid[t], self.dec_centroid[t], c="r") - # plt.gca().set_aspect("equal") - # - # self.ra_centroid *= u.deg - # self.dec_centroid *= u.deg - # del dra_m, ddec_m - # return - - # def _get_centroid(self, frame_indices="mean", binsize=0.0005): - # if isinstance(frame_indices, str): - # if frame_indices == "mean": - # frame_indices = np.where(self.time_mask)[0] - # else: - # frame_indices = np.atleast_1d(frame_indices) - # - # x, y = ( - # self.rough_mask.multiply(self.dra).data, - # self.rough_mask.multiply(self.ddec).data, - # ) - # ar, X, Y = np.histogram2d( - # x, - # y, - # (np.arange(-0.01, 0.01, binsize), np.arange(-0.01, 0.011, binsize)), - # ) - # X, Y = np.meshgrid(X, Y) - # j = np.zeros(np.asarray(ar.T.shape) + 1, bool) - # j[:-1, :-1] = ar.T != 0 - # - # c = ( - # self.rough_mask.multiply(self.flux[frame_indices].mean(axis=0)) - # .multiply(1 / self.source_flux_estimates[:, None]) - # .data - # ) - # k = c < 0.8 - # x_k, y_k, c_k = x[k], y[k], c[k] - # ar2 = np.asarray( - # [ - # np.nanmedian( - # c_k[ - # (x_k >= X1) - # & (x_k < X1 + binsize) - # & (y_k >= Y1) - # & (y_k < Y1 + binsize) - # ] - # ) - # for X1, Y1 in zip(X[j], Y[j]) - # ] - # ) - # self.ra_centroid = ( - # np.average(X[j] + binsize / 2, weights=np.nan_to_num(ar2 ** 3)) * u.deg - # ) - # self.dec_centroid = ( - # np.average(Y[j] + binsize / 2, weights=np.nan_to_num(ar2 ** 3)) * u.deg - # ) - # return def _get_centroid(self, plot=False): radius = 1.5 * self.pixel_scale.to(u.arcsecond).value if not isinstance(self.r, sparse.csr_matrix): @@ -1155,82 +1045,6 @@ def build_shape_model( return self.plot_shape_model(self.source_mask, frame_index=frame_index) return - # - # def _update_source_mask_remove_bkg_pixels(self, flux_cut_off=1, frame_index="max"): - # """ - # Update the `source_mask` to remove pixels that do not contribuite to the PRF - # shape. - # First, re-estimate the source flux usign the precomputed `mean_model`. - # This re-estimation is used to remove sources with bad prediction and update - # the `source_mask` by removing background pixels that do not contribuite to - # the PRF shape. - # Pixels with normalized flux > `flux_cut_off` are kept. - # - # Parameters - # ---------- - # flux_cut_off : float - # Lower limit for the normalized flux predicted from the mean model. - # frame_index : string or int - # The frame index to be used, if "mean" then use the - # mean value across time - # """ - # - # # Re-estimate source flux - # # ----- - # prior_mu = self.source_flux_estimates - # prior_sigma = ( - # np.ones(self.mean_model.shape[0]) * 10 * self.source_flux_estimates - # ) - # - # if frame_index == "mean": - # f, fe = (self.flux).mean(axis=0), ( - # (self.flux_err ** 2).sum(axis=0) ** 0.5 - # ) / (self.nt) - # if frame_index == "max": - # f, fe = (self.flux).max(axis=0), ( - # (self.flux_err ** 2).sum(axis=0) ** 0.5 - # ) / (self.nt) - # elif isinstance(frame_index, int): - # f, fe = self.flux[frame_index], self.flux_err[frame_index] - # - # X = self.mean_model.copy() - # X = X.T - # - # sigma_w_inv = X.T.dot(X.multiply(1 / fe[:, None] ** 2)).toarray() - # sigma_w_inv += np.diag(1 / (prior_sigma ** 2)) - # B = X.T.dot((f / fe ** 2)) - # B += prior_mu / (prior_sigma ** 2) - # ws = np.linalg.solve(sigma_w_inv, B) - # werrs = np.linalg.inv(sigma_w_inv).diagonal() ** 0.5 - # - # # ----- - # - # # Rebuild source mask - # ok = np.abs(ws - self.source_flux_estimates) / werrs > 3 - # ok &= ((ws / self.source_flux_estimates) < 10) & ( - # (self.source_flux_estimates / ws) < 10 - # ) - # ok &= ws > 10 - # ok &= werrs > 0 - # - # self.source_flux_estimates[ok] = ws[ok] - # - # self.source_mask = ( - # self.mean_model.multiply( - # self.mean_model.T.dot(self.source_flux_estimates) - # ).tocsr() - # > flux_cut_off - # ) - # - # # Recreate uncontaminated mask - # self._get_uncontaminated_pixel_mask() - # # self.uncontaminated_source_mask = self.uncontaminated_source_mask.multiply( - # # (self.mean_model.max(axis=1) < 1) - # # ) - # - # # Recreate mean model! - # self._get_mean_model() - def _get_mean_model( self, mask=None # , relative_flux_limit=0.001, absolute_flux_limit=1 ): From ca48f4532884201bdb0d574923da199c87c2fb6c Mon Sep 17 00:00:00 2001 From: Christina Hedges Date: Fri, 18 Mar 2022 15:22:37 -0700 Subject: [PATCH 5/8] :recycle: flake8 --- src/psfmachine/ffi.py | 20 +------------------- src/psfmachine/machine.py | 4 +--- src/psfmachine/utils.py | 2 +- 3 files changed, 3 insertions(+), 23 deletions(-) diff --git a/src/psfmachine/ffi.py b/src/psfmachine/ffi.py index 3221afa..8329b45 100644 --- a/src/psfmachine/ffi.py +++ b/src/psfmachine/ffi.py @@ -7,12 +7,10 @@ import matplotlib.colors as colors from astropy.io import fits -from astropy.stats import SigmaClip from astropy.time import Time from astropy.wcs import WCS import astropy.units as u from astropy.stats import sigma_clip -from photutils import Background2D, MedianBackground, BkgZoomInterpolator # from . import PACKAGEDIR from .utils import do_tiled_query, _make_A_cartesian, solve_linear_model, _load_image @@ -97,8 +95,6 @@ def __init__( self.meta = meta self.wcs = wcs - # keep 2d image for easy plotting - # self.flux_2d = flux self.image_shape = flux.shape[1:] # reshape flux and flux_err as [ntimes, npix] self.flux = flux.reshape(flux.shape[0], -1) @@ -118,7 +114,7 @@ def __init__( ).mask ).ravel() self._remove_background(mask=mask) - + # init `machine` object super().__init__( time, @@ -592,15 +588,6 @@ def _mask_pixels(self, pixel_saturation_limit=1.2e5, magnitude_bright_limit=8): ).tocsr() self.uncontaminated_source_mask.eliminate_zeros() - # self.column = self.column[good_pixels] - # self.row = self.row[good_pixels] - # self.ra = self.ra[good_pixels] - # self.dec = self.dec[good_pixels] - # self.flux = self.flux[:, good_pixels] - # self.flux_err = self.flux_err[:, good_pixels] - # - # return - def residuals(self, plot=False, zoom=False, metric="residuals"): """ Get the residuals (model - image) and compute statistics. It creates a model @@ -966,12 +953,7 @@ def _load_file(fname, extension=1, cutout_size=256, cutout_origin=[0, 0]): tdx = np.argsort(times) times = times[tdx] - # remove overscan of image - # row_2d, col_2d, flux_2d = _remove_overscan(meta["TELESCOP"], np.array(imgs)[tdx]) - # kepler FFIs have uncent maps stored in different files, so we use Poison noise as - # flux error for now. imgs = np.asarray(imgs) - # flux_err_2d = np.sqrt(np.abs(imgs)) imgs_err = np.asarray(imgs_err) # convert to RA and Dec diff --git a/src/psfmachine/machine.py b/src/psfmachine/machine.py index 05d562d..1bec617 100644 --- a/src/psfmachine/machine.py +++ b/src/psfmachine/machine.py @@ -401,9 +401,6 @@ def _get_source_mask(self, source_flux_limit=1): .multiply(1 / self.source_flux_estimates[:, None]) .data ) - sourcef = np.log10( - mask.astype(float).multiply(self.source_flux_estimates[:, None]).data - ) rbins = np.linspace(0, self.radius * 5, 50) masks = np.asarray( [ @@ -458,6 +455,7 @@ def _get_centroid(self, plot=False): marker="*", s=100, ) + return fig def _time_bin(self, npoints=200, downsample=False): """Bin the flux data down in time. If using `pos_corr`s as corrector, it will diff --git a/src/psfmachine/utils.py b/src/psfmachine/utils.py index 585c15e..67ef1f2 100644 --- a/src/psfmachine/utils.py +++ b/src/psfmachine/utils.py @@ -603,7 +603,7 @@ def _load_image( c_min += cutout_origin[1] if (r_min > r_max) | (c_min > c_max): raise ValueError("`cutout_origin` must be within the image.") - if cutout_size != None: + if cutout_size is not None: r_max = np.min([r_min + cutout_size, r_max]) c_max = np.min([c_min + cutout_size, c_max]) if return_coords: From 0bfeab0c689cc1475bf29c548fdd3d5d0d91c150 Mon Sep 17 00:00:00 2001 From: Christina Hedges Date: Fri, 18 Mar 2022 15:28:01 -0700 Subject: [PATCH 6/8] :recycle: fixing SSMachine --- src/psfmachine/superstamp.py | 75 +++++++++++++++++++++++++++++++++++- 1 file changed, 74 insertions(+), 1 deletion(-) diff --git a/src/psfmachine/superstamp.py b/src/psfmachine/superstamp.py index 1bc5791..1ada51d 100644 --- a/src/psfmachine/superstamp.py +++ b/src/psfmachine/superstamp.py @@ -15,7 +15,7 @@ from astropy.wcs import WCS import astropy.units as u -from .ffi import FFIMachine, _get_sources, _do_image_cutout +from .ffi import FFIMachine, _get_sources __all__ = ["SSMachine"] @@ -616,3 +616,76 @@ def _load_file(fname): poscorr2, meta, ) + + +def _do_image_cutout( + flux, flux_err, ra, dec, column, row, cutout_size=100, cutout_origin=[0, 0] +): + """ + Creates a cutout of the full image. Return data arrays corresponding to the cutout. + Parameters + ---------- + flux : numpy.ndarray + Data array with Flux values, correspond to full size image. + flux_err : numpy.ndarray + Data array with Flux errors values, correspond to full size image. + ra : numpy.ndarray + Data array with RA values, correspond to full size image. + dec : numpy.ndarray + Data array with Dec values, correspond to full size image. + column : numpy.ndarray + Data array with pixel column values, correspond to full size image. + row : numpy.ndarray + Data array with pixel raw values, correspond to full size image. + cutout_size : int + Size in pixels of the cutout, assumedto be squared. Default is 100. + cutout_origin : tuple of ints + Origin of the cutout following matrix indexing. Default is [0 ,0]. + Returns + ------- + flux : numpy.ndarray + Data array with Flux values of the cutout. + flux_err : numpy.ndarray + Data array with Flux errors values of the cutout. + ra : numpy.ndarray + Data array with RA values of the cutout. + dec : numpy.ndarray + Data array with Dec values of the cutout. + column : numpy.ndarray + Data array with pixel column values of the cutout. + row : numpy.ndarray + Data array with pixel raw values of the cutout. + """ + if (cutout_size + cutout_origin[0] <= flux.shape[1]) and ( + cutout_size + cutout_origin[1] <= flux.shape[2] + ): + column = column[ + cutout_origin[0] : cutout_origin[0] + cutout_size, + cutout_origin[1] : cutout_origin[1] + cutout_size, + ] + row = row[ + cutout_origin[0] : cutout_origin[0] + cutout_size, + cutout_origin[1] : cutout_origin[1] + cutout_size, + ] + flux = flux[ + :, + cutout_origin[0] : cutout_origin[0] + cutout_size, + cutout_origin[1] : cutout_origin[1] + cutout_size, + ] + flux_err = flux_err[ + :, + cutout_origin[0] : cutout_origin[0] + cutout_size, + cutout_origin[1] : cutout_origin[1] + cutout_size, + ] + ra = ra[ + cutout_origin[0] : cutout_origin[0] + cutout_size, + cutout_origin[1] : cutout_origin[1] + cutout_size, + ] + dec = dec[ + cutout_origin[0] : cutout_origin[0] + cutout_size, + cutout_origin[1] : cutout_origin[1] + cutout_size, + ] + else: + raise ValueError("Cutout size is larger than image shape ", flux.shape) + + return flux, flux_err, ra, dec, column, row From d5d271f6bcd2423c47cf91aab16ad37a65756501 Mon Sep 17 00:00:00 2001 From: Christina Hedges Date: Tue, 22 Mar 2022 11:41:04 -0700 Subject: [PATCH 7/8] :white_check_mark: making them tests work --- src/psfmachine/ffi.py | 42 ++++++++++++++------------- src/psfmachine/machine.py | 60 +++++++++++++++++++++++---------------- src/psfmachine/tpf.py | 3 +- src/psfmachine/utils.py | 11 +++++-- tests/test_ffimachine.py | 15 +++++----- tests/test_machine.py | 18 ++++++------ 6 files changed, 85 insertions(+), 64 deletions(-) diff --git a/src/psfmachine/ffi.py b/src/psfmachine/ffi.py index 8329b45..0240118 100644 --- a/src/psfmachine/ffi.py +++ b/src/psfmachine/ffi.py @@ -13,7 +13,12 @@ from astropy.stats import sigma_clip # from . import PACKAGEDIR -from .utils import do_tiled_query, _make_A_cartesian, solve_linear_model, _load_image +from .utils import ( + do_tiled_query, + _make_A_cartesian, + solve_linear_model, + _load_ffi_image, +) from .tpf import _clean_source_list from .machine import Machine @@ -185,17 +190,12 @@ def from_file( A Machine class object built from the FFI. """ # load FITS files and parse arrays - ( - wcs, - time, - flux, - flux_err, - ra, - dec, - column, - row, - metadata, - ) = _load_file(fname, extension=extension) + (wcs, time, flux, flux_err, ra, dec, column, row, metadata,) = _load_file( + fname, + extension=extension, + cutout_size=cutout_size, + cutout_origin=cutout_origin, + ) if metadata["TELESCOP"] == "Kepler": ngrid = (2, 2) if flux.shape[1] <= 500 else (4, 4) @@ -460,7 +460,7 @@ def _remove_background(self, mask=None, pixel_knot_spacing=10): self.column.max() - self.column.min() + 1, ] ) - n_knots = cutout_size // pixel_knot_spacing + n_knots = np.min([cutout_size // pixel_knot_spacing, 1]) if mask is None: k = np.ones(self.flux.shape[1]) else: @@ -473,7 +473,9 @@ def _remove_background(self, mask=None, pixel_knot_spacing=10): n_knots=n_knots, ) ws = np.linalg.solve( - X[k].T.dot(X[k]).toarray(), X[k].T.dot(self.flux[:, k].T) + X[k].T.dot(X[k]).toarray() + + np.diag(1 / (np.ones(X.shape[1]) * 1000000)), + X[k].T.dot(self.flux[:, k].T), ) bkg = X.dot(ws).T @@ -900,7 +902,7 @@ def _load_file(fname, extension=1, cutout_size=256, cutout_origin=[0, 0]): raise ValueError("FFI is not from Kepler or TESS.") if i == 0: wcs = WCS(hdr) - col_2d, row_2d, f2d = _load_image( + col_2d, row_2d, f2d = _load_ffi_image( telescopes[-1], f, extension, @@ -911,14 +913,16 @@ def _load_file(fname, extension=1, cutout_size=256, cutout_origin=[0, 0]): imgs.append(f2d) else: imgs.append( - _load_image(telescopes[-1], f, extension, cutout_size, cutout_origin) + _load_ffi_image( + telescopes[-1], f, extension, cutout_size, cutout_origin + ) ) - if telescopes[-1].lower() in ["tess", "kepler"]: + if telescopes[-1].lower() in ["tess"]: imgs_err.append( - _load_image(telescopes[-1], f, 2, cutout_size, cutout_origin) + _load_ffi_image(telescopes[-1], f, 2, cutout_size, cutout_origin) ) else: - raise ValueError("Must be Kepler or TESS") + imgs_err.append(imgs[-1] ** 0.5) # check for integrity of files, same telescope, all FFIs and same quarter/campaign if len(set(telescopes)) != 1: diff --git a/src/psfmachine/machine.py b/src/psfmachine/machine.py index 1bec617..f347d04 100644 --- a/src/psfmachine/machine.py +++ b/src/psfmachine/machine.py @@ -176,7 +176,10 @@ def __init__( self.poscorr_filter_size = 12 self.cartesian_knot_spacing = "linear" self.pixel_scale = ( - np.hypot(np.min(np.diff(self.ra)), np.min(np.diff(self.dec))) * u.deg + np.hypot( + np.min(np.abs(np.diff(self.ra))), np.min(np.abs(np.diff(self.dec))) + ) + * u.deg ).to(u.arcsecond) # disble tqdm prgress bar when running in HPC @@ -278,16 +281,13 @@ def _update_delta_numpy_arrays(self, frame_indices="mean"): for idx in range(len(self.sources)) ] ).transpose(1, 0, 2) - self.dra = self.dra * (u.deg) - self.ddec = self.ddec * (u.deg) + # self.dra = self.dra * (u.deg) + # self.ddec = self.ddec * (u.deg) # convertion to polar coordinates - self.r = ( - np.hypot( - self.dx, - self.dy, - ) - * 3600 + self.r = np.hypot( + self.dx * (u.deg.to(u.arcsecond)), + self.dy * (u.deg.to(u.arcsecond)), ) self.phi = np.arctan2( self.dy, @@ -351,7 +351,7 @@ def _update_delta_sparse_arrays(self, frame_indices="mean", dist_lim=50): # then rebuild r, phi as sparse. nnz_inds = sparse_mask.nonzero() # convert radial dist to arcseconds - r_vals = np.hypot(self.dx.data, self.dy.data) * 3600 + r_vals = np.hypot(self.dx.data, self.dy.data) * (u.deg.to(u.arcsecond)) phi_vals = np.arctan2(self.dy.data, self.dx.data) self.r = sparse.csr_matrix( (r_vals, (nnz_inds[0], nnz_inds[1])), @@ -382,8 +382,8 @@ def _get_source_mask(self, source_flux_limit=1): Upper limit on radius of circular apertures in arcsecond. """ self.radius = 2 * self.pixel_scale.to(u.arcsecond).value - if not isinstance(self.r, sparse.csr_matrix): - self.rough_mask = sparse.csr_matrix(self.r < self.radius[:, None]) + if not sparse.issparse(self.r): + self.rough_mask = sparse.csr_matrix(self.r < self.radius) else: self.rough_mask = sparse_lessthan(self.r, self.radius) self.source_mask = self.rough_mask.copy() @@ -397,7 +397,7 @@ def _get_source_mask(self, source_flux_limit=1): r = mask.multiply(self.r).data max_f = np.log10( mask.astype(float) - .multiply(np.max(self.flux[self.time_mask], axis=0)) + .multiply(np.nanmax(self.flux[self.time_mask], axis=0)) .multiply(1 / self.source_flux_estimates[:, None]) .data ) @@ -411,13 +411,22 @@ def _get_source_mask(self, source_flux_limit=1): fbins = np.asarray([np.nanpercentile(max_f[m], 20) for m in masks]) rbins = rbins[1:] - np.median(np.diff(rbins)) k = np.isfinite(fbins) + if not k.any(): + raise ValueError("Can not find source mask") l = np.polyfit(rbins[k], fbins[k], 1) - mean_model = self.r.copy() - mean_model.data = 10 ** np.polyval(l, mean_model.data) - self.source_mask = ( - mean_model.multiply(self.source_flux_estimates[:, None]) - ) > source_flux_limit + if sparse.issparse(r): + mean_model = self.r.copy() + mean_model.data = 10 ** np.polyval(l, mean_model.data) + self.source_mask = ( + mean_model.multiply(self.source_flux_estimates[:, None]) + ) > source_flux_limit + else: + mean_model = 10 ** np.polyval(l, self.r) + self.source_mask = ( + sparse.csr_matrix(mean_model * self.source_flux_estimates[:, None]) + > source_flux_limit + ) self.uncontaminated_source_mask = _find_uncontaminated_pixels( self.source_mask ) @@ -425,8 +434,8 @@ def _get_source_mask(self, source_flux_limit=1): def _get_centroid(self, plot=False): radius = 1.5 * self.pixel_scale.to(u.arcsecond).value - if not isinstance(self.r, sparse.csr_matrix): - mask = sparse.csr_matrix(self.r < radius[:, None]) + if not sparse.issparse(self.r): + mask = sparse.csr_matrix(self.r < radius) else: mask = sparse_lessthan(self.r, radius) mask = _find_uncontaminated_pixels(mask) @@ -440,10 +449,9 @@ def _get_centroid(self, plot=False): .data ) bm, nx, ny, nz, nze = threshold_bin(x, y, c, bins=20) - k = bm > 10 self.ra_centroid, self.dec_centroid = ( - np.average(nx[k], weights=(nz / nze)[k]) * u.deg, - np.average(ny[k], weights=(nz / nze)[k]) * u.deg, + np.average(nx, weights=(nz / nze)) * u.deg, + np.average(ny, weights=(nz / nze)) * u.deg, ) if plot: fig = plt.figure(figsize=(5, 5)) @@ -831,8 +839,10 @@ def plot_time_model(self): A3 = _combine_A(A2, time=time_binned) dx, dy = ( - self.uncontaminated_source_mask.multiply(self.dx).data * 3600, - self.uncontaminated_source_mask.multiply(self.dy).data * 3600, + self.uncontaminated_source_mask.multiply(self.dx).data + * (u.deg.to(u.arcsecond)), + self.uncontaminated_source_mask.multiply(self.dy).data + * (u.deg.to(u.arcsecond)), ) model = A3.dot(self.time_model_w).reshape(flux_binned.shape) + 1 fig, ax = plt.subplots(2, 2, figsize=(7, 6), facecolor="w") diff --git a/src/psfmachine/tpf.py b/src/psfmachine/tpf.py index 9da96bd..3f840dd 100644 --- a/src/psfmachine/tpf.py +++ b/src/psfmachine/tpf.py @@ -611,7 +611,6 @@ def load_shape_model(self, input=None, plot=False): # create source mask and uncontaminated pixel mask self._get_source_mask() - self._get_uncontaminated_pixel_mask() # open file hdu = fits.open(input) @@ -647,7 +646,7 @@ def load_shape_model(self, input=None, plot=False): # work in arcseconds self._get_mean_model() # remove background pixels and recreate mean model - self._update_source_mask_remove_bkg_pixels() + # self._update_source_mask_remove_bkg_pixels() if plot: return self.plot_shape_model() diff --git a/src/psfmachine/utils.py b/src/psfmachine/utils.py index 67ef1f2..736a3b4 100644 --- a/src/psfmachine/utils.py +++ b/src/psfmachine/utils.py @@ -566,7 +566,7 @@ def _find_uncontaminated_pixels(mask): return new_mask -def _load_image( +def _load_ffi_image( telescope, fname, extension, @@ -585,6 +585,7 @@ def _load_image( extension: int Extension to cut out of the image """ + f = fitsio.FITS(fname)[extension] if telescope.lower() == "kepler": # CCD overscan for Kepler r_min = 20 @@ -599,6 +600,10 @@ def _load_image( c_max = 2093 else: raise TypeError("File is not from Kepler or TESS mission") + # If the image dimension is not the FFI shape, we change the r_max and c_max + dims = f.get_dims() + if dims != [r_max, c_max]: + r_max, c_max = np.asarray(dims) r_min += cutout_origin[0] c_min += cutout_origin[1] if (r_min > r_max) | (c_min > c_max): @@ -608,5 +613,5 @@ def _load_image( c_max = np.min([c_min + cutout_size, c_max]) if return_coords: row_2d, col_2d = np.mgrid[r_min:r_max, c_min:c_max] - return col_2d, row_2d, fitsio.FITS(fname)[extension][r_min:r_max, c_min:c_max] - return fitsio.FITS(fname)[extension][r_min:r_max, c_min:c_max] + return col_2d, row_2d, f[r_min:r_max, c_min:c_max] + return f[r_min:r_max, c_min:c_max] diff --git a/tests/test_ffimachine.py b/tests/test_ffimachine.py index e870a04..d9a216b 100644 --- a/tests/test_ffimachine.py +++ b/tests/test_ffimachine.py @@ -21,13 +21,14 @@ def test_ffi_from_file(): assert isinstance(ffi, Machine) # test attributes have the right shapes assert ffi.time.shape == (1,) - assert ffi.flux.shape == (1, 33812) - assert ffi.flux_2d.shape == (1, 180, 188) - assert ffi.flux_err.shape == (1, 33812) - assert ffi.column.shape == (33812,) - assert ffi.row.shape == (33812,) - assert ffi.ra.shape == (33812,) - assert ffi.dec.shape == (33812,) + img_shape = ffi.image_shape + assert ffi.flux.shape == (1, img_shape[0] * img_shape[1]) + assert ffi.flux_2d.shape == (1, *img_shape) + assert ffi.flux_err.shape == (1, img_shape[0] * img_shape[1]) + assert ffi.column.shape == (img_shape[0] * img_shape[1],) + assert ffi.row.shape == (img_shape[0] * img_shape[1],) + assert ffi.ra.shape == (img_shape[0] * img_shape[1],) + assert ffi.dec.shape == (img_shape[0] * img_shape[1],) assert ffi.sources.shape == (259, 13) diff --git a/tests/test_machine.py b/tests/test_machine.py index 5e37968..d6232ef 100644 --- a/tests/test_machine.py +++ b/tests/test_machine.py @@ -5,6 +5,7 @@ from scipy import sparse import pytest import lightkurve as lk +import numpy as np from astropy.utils.data import get_pkg_data_filename from psfmachine import TPFMachine @@ -19,7 +20,7 @@ def test_create_delta_sparse_arrays(): machine = TPFMachine.from_TPFs(tpfs) # create numpy arrays - machine._create_delta_arrays() + # machine._create_delta_arrays() non_sparse_arr = machine.__dict__.copy() # check for main attrs shape @@ -44,12 +45,13 @@ def test_create_delta_sparse_arrays(): # manually mask numpy arrays to compare them vs sparse array dist_lim = 40 - mask = (np.abs(non_sparse_arr["dra"].value) <= dist_lim / 3600) & ( - np.abs(non_sparse_arr["ddec"].value) <= dist_lim / 3600 + mask = (np.abs(non_sparse_arr["dra"]) <= dist_lim / 3600) & ( + np.abs(non_sparse_arr["ddec"]) <= dist_lim / 3600 ) # create sparse arrays - machine._create_delta_sparse_arrays(dist_lim=dist_lim) + del machine.dra + machine._update_delta_sparse_arrays(dist_lim=dist_lim) sparse_arr = machine.__dict__.copy() assert sparse_arr["dra"].shape == non_sparse_arr["dra"].shape @@ -68,10 +70,10 @@ def test_create_delta_sparse_arrays(): assert sparse_arr["phi"].data.shape == (853,) # compare sparse array vs numpy array values - assert (non_sparse_arr["dra"][mask].value == sparse_arr["dra"].data).all() - assert (non_sparse_arr["ddec"][mask].value == sparse_arr["ddec"].data).all() - assert (non_sparse_arr["r"][mask].value == sparse_arr["r"].data).all() - assert (non_sparse_arr["phi"][mask].value == sparse_arr["phi"].data).all() + assert (non_sparse_arr["dra"][mask] == sparse_arr["dra"].data).all() + assert (non_sparse_arr["ddec"][mask] == sparse_arr["ddec"].data).all() + assert np.allclose(non_sparse_arr["r"][mask], sparse_arr["r"].data) + assert np.allclose(non_sparse_arr["phi"][mask], sparse_arr["phi"].data) @pytest.mark.remote_data From 1aacdf47debce363d95d2c805a51cc6a113469a9 Mon Sep 17 00:00:00 2001 From: Christina Hedges Date: Tue, 22 Mar 2022 12:46:04 -0700 Subject: [PATCH 8/8] :white_check_mark: flake8 --- tests/test_machine.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_machine.py b/tests/test_machine.py index d6232ef..7c084f1 100644 --- a/tests/test_machine.py +++ b/tests/test_machine.py @@ -5,7 +5,6 @@ from scipy import sparse import pytest import lightkurve as lk -import numpy as np from astropy.utils.data import get_pkg_data_filename from psfmachine import TPFMachine