From 6d7b28d5644d87b49cd1452350c3298de6436963 Mon Sep 17 00:00:00 2001 From: Liam Gray Date: Mon, 13 May 2024 16:28:05 -0700 Subject: [PATCH 1/8] fix(containers): reverse_map["stack"] would not be created when using "copy_from" --- draco/core/containers.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/draco/core/containers.py b/draco/core/containers.py index ccd5dea70..6452d4a5f 100644 --- a/draco/core/containers.py +++ b/draco/core/containers.py @@ -806,6 +806,11 @@ def __init__(self, *args, **kwargs): # Call initializer from `ContainerBase` super().__init__(*args, **kwargs) + # `axes_from` can be provided via the `copy_from` argument, so + # need to update kwargs to account + if ("axes_from" not in kwargs) and ("copy_from" in kwargs): + kwargs["axes_from"] = kwargs.pop("copy_from") + reverse_map_stack = None # Create reverse map if "reverse_map_stack" in kwargs: From f8eb3d5db5b016c017e1bb50809b6f8e1305b939 Mon Sep 17 00:00:00 2001 From: Liam Gray Date: Wed, 15 May 2024 12:11:29 -0700 Subject: [PATCH 2/8] feat(containers): add an 'effective_ra' dataset to SiderealStream --- draco/core/containers.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/draco/core/containers.py b/draco/core/containers.py index 6452d4a5f..7104cc893 100644 --- a/draco/core/containers.py +++ b/draco/core/containers.py @@ -1239,6 +1239,17 @@ class SiderealStream( "compression_opts": COMPRESSION_OPTS, "chunks": (64, 128, 128), }, + "effective_ra": { + "axes": ["freq", "stack", "ra"], + "dtype": np.float32, + "initialise": False, + "distributed": True, + "distributed_axis": "freq", + "compression": COMPRESSION, + "compression_opts": COMPRESSION_OPTS, + "chunks": (64, 128, 128), + "truncate": True, + }, } @property @@ -1256,6 +1267,14 @@ def _mean(self): """Get the vis dataset.""" return self.datasets["vis"] + @property + def effective_ra(self): + """Get the effective_ra dataset if it exists, None otherwise.""" + if "effective_ra" in self.datasets: + return self.datasets["effective_ra"] + + raise KeyError("Dataset 'effective_ra' not initialised.") + class SystemSensitivity(FreqContainer, TODContainer): """A container for holding the total system sensitivity. From 8fff0616e9a6afc15d3f20419e3aeb1cd37e2695 Mon Sep 17 00:00:00 2001 From: Liam Gray Date: Wed, 15 May 2024 12:11:57 -0700 Subject: [PATCH 3/8] feat(regrid): add a method to calculate a proportional mapping rebinning matrix --- draco/util/regrid.py | 137 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 137 insertions(+) diff --git a/draco/util/regrid.py b/draco/util/regrid.py index b704eb6ad..2c4733521 100644 --- a/draco/util/regrid.py +++ b/draco/util/regrid.py @@ -4,8 +4,11 @@ `_. """ +from typing import List, Union + import numpy as np import scipy.linalg as la +import scipy.sparse as ss from ..util import _fast_tools @@ -155,3 +158,137 @@ def lanczos_inverse_matrix(x, y, a=5, cond=1e-1): """ lz_forward = lanczos_forward_matrix(x, y, a) return la.pinv(lz_forward, rcond=cond) + + +def rebin_matrix(tra: np.ndarray, ra: np.ndarray, width_t: float = 0) -> np.ndarray: + """Construct a matrix to rebin the samples. + + Parameters + ---------- + tra + The samples we have in the time stream. + ra + The target samples we want in the sidereal stream. + width_t + The width of a time sample. Set to zero to do a nearest bin assignment. + + Returns + ------- + R + A matrix to perform the rebinning. + """ + R = np.zeros((ra.shape[0], tra.shape[0])) + + inds = np.searchsorted(ra, tra) + + # Estimate the typical width of an RA bin + width_ra = np.median(np.abs(np.diff(ra))) + + lowest_ra = ra[0] - width_ra / 2 + highest_ra = ra[-1] + width_ra / 2 + + # NOTE: this is a bit of a hack to avoid zero division by zeros, but should have + # the required effect of giving 1 if the time sample is inside a bin, and zero + # otherwise. + if width_t == 0: + width_t = 1e-10 + + # NOTE: this can probably be done more efficiently, but we typically only need + # to do this once per day + for ii, (jj, t) in enumerate(zip(inds, tra)): + + lower_edge = t - width_t / 2.0 + upper_edge = t + width_t / 2.0 + + # If we are in here we have weight to assign to the sample below + if upper_edge > lowest_ra and jj < len(ra): + ra_edge = ra[jj] - width_ra / 2 + wh = np.clip((upper_edge - ra_edge) / width_t, 0.0, 1.0) + R[jj, ii] = wh + + if lower_edge < highest_ra and jj > 0: + ra_edge = ra[jj - 1] + width_ra / 2 + wl = np.clip((ra_edge - lower_edge) / width_t, 0.0, 1.0) + R[jj - 1, ii] = wl + + return R + + +def taylor_coeff( + x: np.ndarray, + N: int, + M: int, + Ni: np.ndarray, + Si: float, + period: Union[float, None] = None, + xc: Union[np.ndarray, None] = None, +) -> List[ss.csr_array]: + """Return a set of sparse matrices that estimates expansion coefficients. + + Parameters + ---------- + x + Positions of each element. + N + Number of positions each side to estimate from. + M + The number of terms in the expansion. + Ni + The weight for each position. The inverse noise variance if interpreted as a + Wiener filter. + Si + A regulariser. The inverse signal variance if interpreted as a Wiener filter. + period + If set, assume the axis is periodic with this period. + xc + An optional parameter giving the location to return the coefficients at. If not + set, just use the locations of each individual sample. + + + Returns + ------- + matrices + A set of sparse matrices that will estimate the coefficents at each location. + Each matrix is for a different coefficient. + """ + nx = x.shape[0] + + ind = np.arange(nx)[:, np.newaxis] + np.arange(-N, N + 1)[np.newaxis, :] + + xc = x if xc is None else xc + + # If periodic then just wrap back around + if period is not None: + ind = ind % nx + xf = x[ind] - xc[:, np.newaxis] + x = ((x + period / 2) % period) - period / 2 + Na = Ni[ind] + + # If not then set the weights for out of bounds entries to zero + else: + mask = (ind < 0) | (ind >= nx) + ind = np.where(mask, 0, ind) + xf = x[ind] - x[:, np.newaxis] + Na = Ni[ind] + Na[mask] = 0.0 + + # Create the required matrices at each location + X = np.stack([xf**m for m in range(M)], axis=2) + XhNi = (X * Na[:, :, np.newaxis]).transpose(0, 2, 1) + XhNiX = XhNi @ X + + # Calculate the covariance part of the filter + Ci = np.identity(M) * Si + XhNiX + C = np.zeros_like(Ci) + for i in range(nx): + C[i] = la.inv(Ci[i]) + + W = C @ XhNi + + # Create the indptr array designating the start and end of the indices for each row + # in the csr_array + indptr = (2 * N + 1) * np.arange(nx + 1, dtype=int) + return [ + ss.csr_array((W[:, i].ravel(), ind.ravel(), indptr), shape=(nx, nx)) + for i in range(M) + ] From ebf992d7992b634b9107e1122c00a318670ca858 Mon Sep 17 00:00:00 2001 From: Liam Gray Date: Wed, 15 May 2024 12:16:57 -0700 Subject: [PATCH 4/8] feat(sidereal): add a new SiderealRebinner task --- draco/analysis/sidereal.py | 96 +++++++++++++++++++++++++++++++++++++- 1 file changed, 95 insertions(+), 1 deletion(-) diff --git a/draco/analysis/sidereal.py b/draco/analysis/sidereal.py index a4e727439..79b97eca5 100644 --- a/draco/analysis/sidereal.py +++ b/draco/analysis/sidereal.py @@ -9,13 +9,15 @@ :class:`SiderealStacker` if you want to combine the different days. """ +from typing import Union + import numpy as np import scipy.linalg as la from caput import config, mpiarray, tod from cora.util import units from ..core import containers, io, task -from ..util import tools +from ..util import regrid, tools from .transform import Regridder @@ -477,6 +479,98 @@ def _regrid(self, vis, weight, lsd): return interp_grid, interp_vis, interp_weight +class SiderealRebinner(SiderealRegridder): + """Regrid a sidereal day of data using a binning method. + + Assign a fraction of each time sample to the nearest RA bin based + on the propotion of the time bin that overlaps the RA bin. + + Tracks the weighted effective centre of RA bin so that a centring + correction can be applied afterwards. A correction option is + implemented in `RebinGradientCorrection`. + """ + + def process(self, data): + """Rebin the sidereal day. + + Parameters + ---------- + data : containers.TimeStream + Timestream data for the day (must have a `LSD` attribute). + + Returns + ------- + sdata : containers.SiderealStream + The regularly gridded sidereal timestream. + """ + import scipy.sparse as ss + + self.log.info(f"Rebinning LSD:{data.attrs['lsd']:.0f}") + + # Redistribute if needed too + data.redistribute("freq") + + # Convert data timestamps into LSDs + timestamp_lsd = self.observer.unix_to_lsd(data.time) + + # Fetch which LSD this is to set bounds + self.start = data.attrs["lsd"] + self.end = self.start + 1 + + # Create the output container + sdata = containers.SiderealStream(axes_from=data, ra=self.samples) + # Initialize the `effective_ra` dataset + sdata.add_dataset("effective_ra") + sdata.redistribute("freq") + sdata.attrs["lsd"] = self.start + sdata.attrs["tag"] = f"lsd_{self.start:.0f}" + + # Get view of data + weight = data.weight[:].local_array + vis_data = data.vis[:].local_array + + # Get the median time sample width + width_t = np.median(np.abs(np.diff(timestamp_lsd))) + # Create the regular grid of RA samples + target_lsd = np.linspace(self.start, self.end, self.samples, endpoint=False) + # Create the rebinning matrix + R = regrid.rebin_matrix(timestamp_lsd, target_lsd, width_t=width_t) + Rt = ss.csr_array(R.T) + # The square is used to calculate rebinned weights + Rtsq = Rt.power(2) + + # dereference arrays before loop + sera = sdata.effective_ra[:].local_array + ssv = sdata.vis[:].local_array + ssw = sdata.weight[:].local_array + + # Create a repeating view of the gridded ra axis. This is + # used to fill the effective ra of masked samples + grid_ra = np.broadcast_to(sdata.ra, (*sera.shape[1:],)) + + for fi in range(vis_data.shape[0]): + # Normalisation for rebinned datasets + norm = tools.invert_no_zero(weight[fi] @ Rt) + + # Weighted rebin of the visibilities + ssv[fi] = norm * ((vis_data[fi] * weight[fi]) @ Rt) + + # Weighted rebin of the effective RA + effective_lsd = norm * ((timestamp_lsd[np.newaxis] * weight[fi]) @ Rt) + sera[fi] = 360 * (effective_lsd - self.start) + + # Rebin the weights + rvar = weight[fi] @ Rtsq + ssw[fi] = tools.invert_no_zero(norm**2 * rvar) + + # Correct the effective ra where weights are zero. This + # is required to avoid discontinuities + mask = ssw[fi] == 0.0 + sera[fi][mask] = grid_ra[mask] + + return sdata + + class SiderealStacker(task.SingleTask): """Take in a set of sidereal days, and stack them up. From a68782d022bc11af4e7903432c6ad95b600dc3e5 Mon Sep 17 00:00:00 2001 From: Liam Gray Date: Wed, 22 May 2024 11:52:37 -0700 Subject: [PATCH 5/8] feat(regrid): add a function to calculate a periodic gradient --- draco/util/regrid.py | 50 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/draco/util/regrid.py b/draco/util/regrid.py index 2c4733521..32aaffc04 100644 --- a/draco/util/regrid.py +++ b/draco/util/regrid.py @@ -214,6 +214,56 @@ def rebin_matrix(tra: np.ndarray, ra: np.ndarray, width_t: float = 0) -> np.ndar return R +def grad_1d( + x: np.ndarray, si: np.ndarray, mask: np.ndarray, period: Union[float, None] = None +) -> np.ndarray: + """Gradient with boundary samples wrapped. + + Parameters + ---------- + x + Data to calculate the gradient for. + si + Positions of the samples. Must be monotonically increasing. + mask + Boolean mask, True where a sample is flagged. + period + Period of `samples`. Default is None, which produces a + non-periodic gradient. + + Returns + ------- + gradient + Gradient of `x`. Gradient is set to zero where any sample in + the calculation was flagged. + mask + Boolean mask corresponding to samples for which a proper + gradient could not be calculated. True where a sample + is flagged. + """ + if period is not None: + # Wrap each array, accounting for the periodicity + # in sample positions + x = np.concatenate(([x[-1]], x, [x[0]])) + mask = np.concatenate(([mask[-1]], mask, [mask[0]])) + # Account for the possibility of multiple periods + shift = np.ceil(si[-1] / period) * period + si = np.concatenate(([si[-1] - shift], si, [si[0] + shift])) + # Return with wrapped samples removed + sel = slice(1, -1) + else: + # Calculate the gradient using `np.gradient` first order + # one-sided difference at the boundaries + sel = slice(None) + + # Extend the flagged values such that any gradient which + # includes a flagged sample is set to 0. This effectively + # involves masking any sample where an adjacent sample is masked + mask |= np.concatenate(([False], mask[:-1])) | np.concatenate((mask[1:], [False])) + + return (~mask * np.gradient(x, si))[sel], mask[sel] + + def taylor_coeff( x: np.ndarray, N: int, From 59e36c313149f407926729f5e1bb0bed1088eaa5 Mon Sep 17 00:00:00 2001 From: Liam Gray Date: Wed, 15 May 2024 12:18:16 -0700 Subject: [PATCH 6/8] feat(sidereal): add a task to apply a linear gradient correction to rebinned data --- draco/analysis/sidereal.py | 108 ++++++++++++++++++++++++++++++++++++- 1 file changed, 106 insertions(+), 2 deletions(-) diff --git a/draco/analysis/sidereal.py b/draco/analysis/sidereal.py index 79b97eca5..3dda53b1e 100644 --- a/draco/analysis/sidereal.py +++ b/draco/analysis/sidereal.py @@ -9,8 +9,6 @@ :class:`SiderealStacker` if you want to combine the different days. """ -from typing import Union - import numpy as np import scipy.linalg as la from caput import config, mpiarray, tod @@ -571,6 +569,112 @@ def process(self, data): return sdata +class RebinGradientCorrection(task.SingleTask): + """Apply a linear gradient correction to shift RA samples to bin centres. + + Requires a sidereal day with full sidereal coverage to calculate a local + gradient for each RA bin. The dataset value at the RA bin centre is + interpolated based on the local gradient and difference between bin centre + and effective bin centre. + + If the rebinned dataset has full sidereal coverage, it can be used to + create the gradient. + """ + + def setup(self, sstream_ref: containers.SiderealStream) -> None: + """Provide the dataset to use in the gradient calculation. + + This dataset must have complete sidereal coverage. + + Parameters + ---------- + sstream_ref + Reference SiderealStream + """ + self.sstream_ref = sstream_ref + + def process(self, sstream: containers.SiderealStream) -> containers.SiderealStream: + """Apply the gradient correction to the input dataset. + + Parameters + ---------- + sstream + sidereal day to apply a correction to + + Returns + ------- + sstream + Input sidereal day with gradient correction applied, if needed + """ + self.sstream_ref.redistribute("freq") + sstream.redistribute("freq") + + # Allows a normal sidereal stream to pass through this task. + # Helpful for creating generic configs + try: + era = sstream.effective_ra[:].local_array + except KeyError: + self.log.info( + f"Dataset of type ({type(sstream)}) does not have an effective " + "ra dataset. No correction will be applied." + ) + return sstream + + try: + # If the reference dataset has an effective ra dataset, use this + # when calculating the gradient. This could be true if the reference + # and target datasets are the same + ref_ra = self.sstream_ref.effective_ra[:].local_array + except KeyError: + # Use fixed ra, which should be regularly sampled + ref_ra = self.sstream_ref.ra + + vis = sstream.vis[:].local_array + weight = sstream.weight[:].local_array + + ref_vis = self.sstream_ref.vis[:].local_array + ref_weight = self.sstream_ref.weight[:].local_array + + # Iterate over frequencies and baselines for memory + for fi in range(vis.shape[0]): + # Skip if entirely masked already + if not np.any(weight[fi]): + continue + + for vi in range(vis.shape[1]): + # Skip if entire baseline is masked + if not np.any(weight[fi, vi]): + continue + + # Depends on whether the effective ra has baseline dependence + rra = ref_ra[fi, vi] if ref_ra.ndim > 1 else ref_ra + # Calculate the vis gradient at the reference RA points + mask = ref_weight[fi, vi] == 0.0 + grad, mask = regrid.grad_1d(ref_vis[fi, vi], rra, mask, period=360.0) + + # Apply the correction to estimate the sample value at the + # RA bin centre + vis[fi, vi] -= grad * (era[fi, vi] - sstream.ra) + # Zero any weights that could not be corrected + weight[fi, vi] *= (~mask).astype(weight.dtype) + + # Copy to a new container without the effective_ra dataset + # This sort of functionality should probably be properly + # implemented somewhere else + out = containers.SiderealStream(attrs_from=sstream, axes_from=sstream) + # Iterate over datasets + for name, dset in sstream.datasets.items(): + if name == "effective_ra": + continue + + if name not in out.datasets: + out.add_dataset(name) + + out.datasets[name][:] = dset[:] + + return out + + class SiderealStacker(task.SingleTask): """Take in a set of sidereal days, and stack them up. From 745e1aec056d5a1f9689a05754771a4b4ff55b0f Mon Sep 17 00:00:00 2001 From: Liam Gray Date: Wed, 15 May 2024 12:20:09 -0700 Subject: [PATCH 7/8] feat(sidereal): update stacker to track rebinned effective ra --- draco/analysis/sidereal.py | 87 ++++++++++++++++++++++++++------------ 1 file changed, 61 insertions(+), 26 deletions(-) diff --git a/draco/analysis/sidereal.py b/draco/analysis/sidereal.py index 3dda53b1e..6b891315a 100644 --- a/draco/analysis/sidereal.py +++ b/draco/analysis/sidereal.py @@ -710,6 +710,13 @@ def process(self, sdata): sdata : containers.SiderealStream Individual sidereal day to add to stack. """ + # Check that the input container is of the correct type + if (self.stack is not None) and not isinstance(sdata, type(self.stack)): + raise TypeError( + f"type(sdata) (={type(sdata)}) does not match " + f"type(stack) (={type(self.stack)})." + ) + sdata.redistribute("freq") # Get the LSD (or CSD) label out of the input's attributes. @@ -748,13 +755,22 @@ def process(self, sdata): # Keep track of the sum of squared weights # to perform Bessel's correction at the end. if self.with_sample_variance: - self.sum_coeff_sq = np.zeros_like(self.stack.weight[:].view(np.ndarray)) + self.sum_coeff_sq = mpiarray.zeros( + self.stack.weight[:].local_shape, + axis=0, + comm=self.stack.comm, + dtype=np.float32, + ) # Accumulate self.log.info(f"Adding to stack LSD(s): {input_lsd!s}") self.lsd_list += input_lsd + # Extract weight dataset + weight = sdata.weight[:] + + # Determine if the input sidereal stream is itself a stack over days if "nsample" in sdata.datasets: # The input sidereal stream is already a stack # over multiple sidereal days. Use the nsample @@ -762,34 +778,22 @@ def process(self, sdata): count = sdata.nsample[:] else: # The input sidereal stream contains a single - # sidereal day. Use a boolean array that + # sidereal day. Use a boolean array that # indicates a non-zero weight dataset as # the weight for the uniform case. - dtype = self.stack.nsample.dtype - count = (sdata.weight[:] > 0.0).astype(dtype) - - # Determine the weights to be used in the average. - if self.weight == "uniform": - coeff = count.astype(np.float32) - # Accumulate the variances in the stack.weight dataset. - self.stack.weight[:] += (coeff**2) * tools.invert_no_zero(sdata.weight[:]) - else: - coeff = sdata.weight[:] - # We are using inverse variance weights. In this case, - # we accumulate the inverse variances in the stack.weight - # dataset. Do that directly to avoid an unneccessary - # division in the more general expression above. - self.stack.weight[:] += sdata.weight[:] + count = (weight > 0.0).astype(self.stack.nsample.dtype) # Accumulate the total number of samples. self.stack.nsample[:] += count - # Below we will need to normalize by the current sum of coefficients. - # Can be found in the stack.nsample dataset for uniform case or - # the stack.weight dataset for inverse variance case. if self.weight == "uniform": - sum_coeff = self.stack.nsample[:].astype(np.float32) + coeff = count.astype(np.float32) + self.stack.weight[:] += (coeff**2) * tools.invert_no_zero(weight) + sum_coeff = self.stack.nsample[:] + else: + coeff = weight + self.stack.weight[:] += weight sum_coeff = self.stack.weight[:] # Calculate weighted difference between the new data and the current mean. @@ -798,6 +802,12 @@ def process(self, sdata): # Update the mean. self.stack.vis[:] += delta_before * tools.invert_no_zero(sum_coeff) + if "effective_ra" in self.stack.datasets: + # Accumulate the weighted average effective RA + delta_ra = coeff * (sdata.effective_ra[:] - self.stack.effective_ra[:]) + # Update the mean + self.stack.effective_ra[:] += delta_ra * tools.invert_no_zero(sum_coeff) + # The calculations below are only required if the sample variance was requested if self.with_sample_variance: # Accumulate the sum of squared coefficients. @@ -822,17 +832,23 @@ def process_finish(self): self.stack.attrs["tag"] = self.tag self.stack.attrs["lsd"] = np.array(self.lsd_list) - # For uniform weighting, normalize the accumulated variances by the total - # number of samples squared and then invert to finalize stack.weight. if self.weight == "uniform": + # For uniform weighting, normalize the accumulated variances by the total + # number of samples squared and then invert to finalize stack.weight. norm = self.stack.nsample[:].astype(np.float32) self.stack.weight[:] = tools.invert_no_zero(self.stack.weight[:]) * norm**2 + else: + # For inverse variance weighting without rebinning, + # additional normalization is not required. + norm = None + # We need to normalize the sample variance by the sum of coefficients. # Can be found in the stack.nsample dataset for uniform case # or the stack.weight dataset for inverse variance case. if self.with_sample_variance: - if self.weight != "uniform": + + if norm is None: norm = self.stack.weight[:] # Perform Bessel's correction. In the case of @@ -882,6 +898,13 @@ def process(self, sdata): sdata : containers.SiderealStream Individual sidereal day to stack up. """ + # Check that the input container is of the correct type + if (self.stack is not None) and not isinstance(sdata, type(self.stack)): + raise TypeError( + f"type(sdata) (={type(sdata)}) does not match " + f"type(stack) (={type(self.stack)})." + ) + sdata.redistribute("freq") if self.stack is None: @@ -889,8 +912,10 @@ def process(self, sdata): self.stack = containers.empty_like(sdata) self.stack.redistribute("freq") - self.stack.vis[:] = 0.0 - self.stack.weight[:] = 0.0 + + # Initialise all datasets to zero + for ds in self.stack.datasets.values(): + ds[:] = 0 self.count = 0 self.Ni_s = mpiarray.zeros( @@ -939,6 +964,16 @@ def process(self, sdata): # We need to keep the projection vector until the end self.Vm.append(v) + if "effective_ra" in self.stack.datasets: + # Track the effective ra bin centres. We're using the averaged + # weights here, so this generally isn't optimal + delta = Ni_d * (sdata.effective_ra[:] - self.stack.effective_ra[:]) + # Update the mean effective ra using the mean of the normalized weights + sum_weight = tools.invert_no_zero(self.stack.weight[:]) * self.Ni_s**2 + self.stack.effective_ra[:] += delta * tools.invert_no_zero( + sum_weight.mean(axis=1) + ) + # Get the LSD label out of the data (resort to using a CSD if it's # present). If there's no label just use a place holder and stack # anyway. From f800c5a2a6d610f39975886d5480ca5f073da0e2 Mon Sep 17 00:00:00 2001 From: Liam Gray Date: Wed, 15 May 2024 12:15:54 -0700 Subject: [PATCH 8/8] feat(flagging): restrict uncorrected rebinned data from blending --- draco/analysis/flagging.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/draco/analysis/flagging.py b/draco/analysis/flagging.py index 87939921c..00e64896c 100644 --- a/draco/analysis/flagging.py +++ b/draco/analysis/flagging.py @@ -2298,10 +2298,16 @@ def process(self, data): The modified data. This is the same object as the input, and it has been modified in place. """ - if type(self.data_stack) != type(data): + if "effective_ra" in data.datasets: + raise TypeError( + "Blending uncorrected rebinned data not supported. " + "Please apply a correction such as `sidereal.RebinGradientCorrection." + ) + + if not isinstance(data, type(self.data_stack)): raise TypeError( f"type(data) (={type(data)}) must match" - f"type(data_stack) (={type(self.type)}" + f"type(data_stack) (={type(self.data_stack)}" ) # Try and get both the stack and the incoming data to have the same