From 5a11b2d08859643dff92a919d90561e7cdfa2e10 Mon Sep 17 00:00:00 2001 From: BrunoSanchez Date: Mon, 13 Jul 2026 02:46:00 -0700 Subject: [PATCH 1/3] Saving work on new blended fakes --- python/lsst/ap/pipe/createApFakes.py | 188 ++++++++++++++++++++++++++- 1 file changed, 187 insertions(+), 1 deletion(-) diff --git a/python/lsst/ap/pipe/createApFakes.py b/python/lsst/ap/pipe/createApFakes.py index fc6328f3..e2fe4400 100644 --- a/python/lsst/ap/pipe/createApFakes.py +++ b/python/lsst/ap/pipe/createApFakes.py @@ -390,7 +390,40 @@ class CreateVisitDetectorFakesConfig( dtype=float, default=26, ) - + doAddBlendedFakes = pexConfig.Field( + doc="Whether to add blended fakes to the visit detector.", + dtype=bool, + default=False, + ) + fracBlendedFakes = pexConfig.RangeField( + doc="Fraction of blended fakes to add to the visit detector.", + dtype=float, + default=0.5, + min=0, + max=1, + ) + fracHostedBlendedFakes = pexConfig.RangeField( + doc="Fraction of blended fakes that are hosted by stars.", + dtype=float, + default=0.5, + min=0, + max=1, + ) + blendedFakeMagOffset = pexConfig.Field( + doc="Standard deviation of magnitude offset for blended fakes.", + dtype=float, + default=0.5, + ) + blendedFakeMaxOffset = pexConfig.Field( + doc="Maximum positional offset for blended fakes in arcseconds.", + dtype=float, + default=5.0, + ) + blendedFakeMinOffset = pexConfig.Field( + doc="Minimum positional offset for blended fakes in arcseconds.", + dtype=float, + default=0.2, + ) class CreateVisitDetectorFakesTask(PipelineTask): """Create and store a set of visit detector fakes for use in AP processing. @@ -623,6 +656,119 @@ def run(self, sourceCat, visit_image): variable_fakes["isVariable"] = True catalog = vstack([catalog, variable_fakes]) + catalog["isBlended"] = False + + if self.config.doAddBlendedFakes: + self.log.info("Generating blended fakes.") + n_blended_fakes_total = int(len(catalog) * self.config.fracBlendedFakes) + n_star_hosted_blended_fakes = int(n_blended_fakes_total * self.config.fracHostedBlendedFakes) + idx = rng.choice( + len(catalog), size=n_blended_fakes_total - n_star_hosted_blended_fakes, replace=False) + # striaghtforward blended fakes + blended_fakes = catalog[idx].copy() + n_blended_fakes = len(blended_fakes) + # blended fakes will be copies of randomly chosen fakes, with a small + # magnitude offset and a small positional offset, and will be flagged as blended. + # The original fakes will also be flagged as blended. + # The idea is to have two fakes in the same location, but with different + # magnitudes, to test the deblending capabilities of the pipeline. + # We work as if the original fake was a "host"of the blend + blended_fakes["delta_mag"] = rng.normal( + loc=0.0, + scale=self.config.blendedFakeMagOffset, + size=n_blended_fakes + ) + blended_fakes["delta_ra"] = rng.uniform( + low=-self.config.blendedFakeMaxOffset, + high=self.config.blendedFakeMaxOffset, + size=n_blended_fakes + ) + blended_fakes["delta_ra"] *= rng.choice([-1, 1], size=n_blended_fakes) + + blended_fakes["delta_dec"] = np.sqrt( + self.config.blendedFakeMaxOffset**2 - blended_fakes["delta_ra"]**2 + ) + blended_fakes["delta_dec"] *= rng.choice([-1, 1], size=n_blended_fakes) + + blended_fakes["host_mag"] = blended_fakes["mag"] + blended_fakes["mag"] += blended_fakes["delta_mag"] + blended_fakes["host_ra"] = blended_fakes["ra"] + blended_fakes["host_dec"] = blended_fakes["dec"] + + blended_fakes["ra"] += blended_fakes["delta_ra"] / 3600.0 + blended_fakes["dec"] += blended_fakes["delta_dec"] / 3600 + blended_fakes["x"], blended_fakes["y"] = wcs.skyToPixelArray( + np.deg2rad(blended_fakes["ra"]), np.deg2rad(blended_fakes["dec"]) + ) + + blended_fakes["source_type"] = "Star" + blended_fakes["host_id"] = blended_fakes["injection_id"] + blended_fakes["twin_id"] = blended_fakes["injection_id"] + blended_fakes["isBlended"] = True + blended_fakes["hosted_fake"] = True + + if self.config.fracHostedBlendedFakes > 0: + hostcatalog = photoCalib.calibrateCatalog(sourceCat).asAstropy() + star_hosts = self.select_host_stars(hostcatalog) + # if len(star_hosts) is less than the blended fakes, then use replacement + idx = rng.choice(len(star_hosts), size=n_star_hosted_blended_fakes, replace=True) + hostcat = star_hosts[idx] + + x_hosts = hostcat['slot_Centroid_x'] + y_hosts = hostcat['slot_Centroid_y'] + ra_hosts = np.rad2deg(hostcat['coord_ra']) + dec_hosts = np.rad2deg(hostcat['coord_dec']) + mag_hosts = hostcat['slot_PsfFlux_mag'] + + # retrieving the global ra dec position of the injection + delta_ra = rng.uniform( + low=-self.config.blendedFakeMaxOffset, + high=self.config.blendedFakeMaxOffset, + size=n_star_hosted_blended_fakes + ) + delta_ra *= rng.choice([-1, 1], size=n_star_hosted_blended_fakes) + + delta_dec = np.sqrt( + self.config.blendedFakeMaxOffset**2 - delta_ra**2 + ) + delta_dec *= rng.choice([-1, 1], size=n_star_hosted_blended_fakes) + + ra_ssi = ra_hosts + delta_ra / 3600.0 + dec_ssi = dec_hosts + delta_dec / 3600.0 + + x_ssi, y_ssi = wcs.skyToPixelArray(np.deg2rad(ra_ssi), np.deg2rad(dec_ssi)) + + delta_mag = rng.normal(loc=1, scale=1, size=n_star_hosted_blended_fakes) + mags = mag_hosts + delta_mag + + # Create the table of hosted fakes + hosted_fakes = Table() + hosted_fakes["x"] = x_ssi + hosted_fakes["y"] = y_ssi + hosted_fakes["mag"] = mags + hosted_fakes["ra"] = ra_ssi + hosted_fakes["dec"] = dec_ssi + hosted_fakes["host_id"] = hostcat['id'] + hosted_fakes["host_flux"] = hostcat['slot_PsfFlux_flux'] + hosted_fakes["host_mag"] = hostcat['slot_PsfFlux_mag'] + hosted_fakes["host_ra"] = ra_hosts + hosted_fakes["host_dec"] = dec_hosts + hosted_fakes["delta_ra"] = delta_ra + hosted_fakes["delta_dec"] = delta_dec + hosted_fakes["delta_mag"] = delta_mag + hosted_fakes["source_type"] = "Star" + hosted_fakes["hosted_fake"] = True + hosted_fakes["isVisitSource"] = True + hosted_fakes["isTemplateSource"] = False + hosted_fakes["isBlended"] = True + + blended_fakes = vstack([blended_fakes, hosted_fakes]) + + blended_fakes["injection_id"] = self._make_unique_injection_ids( + len(blended_fakes), + used_ids=catalog["injection_id"], + ) + catalog = vstack([catalog, blended_fakes]) if len(catalog) > len(np.unique(catalog["injection_id"])): self.log.warning("Duplicate injection IDs detected after catalog assembly; reassigning them.") @@ -681,6 +827,46 @@ def select_hosts(self, sourceCat): skySourceCut & flagCut & extendednessCut & snrCut].copy() return hostCat + def select_host_stars(self, sourceCat): + """ + Selects host sources from a given source catalog based on a series of classification and flux cuts. + The selection criteria are: + - The 'base_ClassificationSizeExtendedness_flag' and + 'base_ClassificationExtendedness_flag' must both be False. + - The 'base_ClassificationSizeExtendedness_value' must be greater than 0.9. + - The 'base_ClassificationExtendedness_value' must be equal to 1. + - The 'base_PsfFlux_flux' must be greater than 0. + Parameters + ---------- + sourceCat : SourceCatalog + The source catalog containing the columns required for selection. + *args, **kwargs + Additional arguments (not used). + Returns + ------- + hostCat : ArrowAstropy + A deep copy of the subset of the source catalog that passes all selection criteria. + """ + + # Avoid calibration stars or psf stars; remove flagged sources sky_sources + skySourceCut = ~sourceCat['sky_source'] + + flagCut = ~sourceCat['base_ClassificationSizeExtendedness_flag'] + flagCut &= ~sourceCat['base_ClassificationExtendedness_flag'] + flagCut &= ~sourceCat['slot_Shape_flag'] + flagCut &= ~sourceCat['slot_Centroid_flag'] + flagCut &= ~sourceCat['base_PixelFlags_flag'] + + extendednessCut = sourceCat['base_ClassificationSizeExtendedness_value'] < 0.9 + extendednessCut &= sourceCat['base_ClassificationExtendedness_value'] != 1 + + snrCut = sourceCat['slot_PsfFlux_flux']/sourceCat['slot_PsfFlux_fluxErr'] > 30 + + hostCat = sourceCat[ + skySourceCut & flagCut & extendednessCut & snrCut].copy() + return hostCat + + def get_PA_and_axes(self, Ixx, Ixy, Iyy): ''' Calculates the orientation and extent of an object based on its second moments. From 457b37d9e3ce60ca47c871c7c293fe13c9eb5108 Mon Sep 17 00:00:00 2001 From: BrunoSanchez Date: Wed, 15 Jul 2026 03:04:41 -0700 Subject: [PATCH 2/3] Save work on fakes --- python/lsst/ap/pipe/createApFakes.py | 168 ++++++++++++++++----------- 1 file changed, 101 insertions(+), 67 deletions(-) diff --git a/python/lsst/ap/pipe/createApFakes.py b/python/lsst/ap/pipe/createApFakes.py index e2fe4400..a816cb62 100644 --- a/python/lsst/ap/pipe/createApFakes.py +++ b/python/lsst/ap/pipe/createApFakes.py @@ -417,13 +417,25 @@ class CreateVisitDetectorFakesConfig( blendedFakeMaxOffset = pexConfig.Field( doc="Maximum positional offset for blended fakes in arcseconds.", dtype=float, - default=5.0, + default=3.0, ) blendedFakeMinOffset = pexConfig.Field( doc="Minimum positional offset for blended fakes in arcseconds.", dtype=float, default=0.2, ) + maxHostedBlendedFakesPerHost = pexConfig.RangeField( + doc="Maximum number of hosted blended fakes assigned to the same host in one detector.", + dtype=int, + default=2, + min=1, + ) + maxHostedBlendedFakesTotal = pexConfig.RangeField( + doc="Hard cap on number of hosted blended fakes per detector. Set to -1 to disable.", + dtype=int, + default=-1, + min=-1, + ) class CreateVisitDetectorFakesTask(PipelineTask): """Create and store a set of visit detector fakes for use in AP processing. @@ -711,80 +723,102 @@ def run(self, sourceCat, visit_image): hostcatalog = photoCalib.calibrateCatalog(sourceCat).asAstropy() star_hosts = self.select_host_stars(hostcatalog) # if len(star_hosts) is less than the blended fakes, then use replacement - idx = rng.choice(len(star_hosts), size=n_star_hosted_blended_fakes, replace=True) - hostcat = star_hosts[idx] - - x_hosts = hostcat['slot_Centroid_x'] - y_hosts = hostcat['slot_Centroid_y'] - ra_hosts = np.rad2deg(hostcat['coord_ra']) - dec_hosts = np.rad2deg(hostcat['coord_dec']) - mag_hosts = hostcat['slot_PsfFlux_mag'] - - # retrieving the global ra dec position of the injection - delta_ra = rng.uniform( - low=-self.config.blendedFakeMaxOffset, - high=self.config.blendedFakeMaxOffset, - size=n_star_hosted_blended_fakes + if len(star_hosts) == 0: + self.log.warning( + "Hosted blended fake generation requested, but no valid star hosts were selected." ) - delta_ra *= rng.choice([-1, 1], size=n_star_hosted_blended_fakes) - - delta_dec = np.sqrt( - self.config.blendedFakeMaxOffset**2 - delta_ra**2 - ) - delta_dec *= rng.choice([-1, 1], size=n_star_hosted_blended_fakes) - - ra_ssi = ra_hosts + delta_ra / 3600.0 - dec_ssi = dec_hosts + delta_dec / 3600.0 + n_star_hosted_blended_fakes = 0 + else: + requested = n_star_hosted_blended_fakes + n_star_hosted_blended_fakes = self._cap_hosted_blended_count( + requested, len(star_hosts) + ) + if n_star_hosted_blended_fakes < requested: + self.log.warning( + "Reducing hosted blended fakes from %d to %d to respect host guard rails.", + requested, + n_star_hosted_blended_fakes, + ) + + if n_star_hosted_blended_fakes > 0: + # Build a finite host pool so each host appears at most maxHostedBlendedFakesPerHost times. + host_pool = np.repeat( + np.arange(len(star_hosts), dtype=int), + self.config.maxHostedBlendedFakesPerHost, + ) + idx = rng.choice(host_pool, size=n_star_hosted_blended_fakes, replace=False) + hostcat = star_hosts[idx] - x_ssi, y_ssi = wcs.skyToPixelArray(np.deg2rad(ra_ssi), np.deg2rad(dec_ssi)) + x_hosts = hostcat['slot_Centroid_x'] + y_hosts = hostcat['slot_Centroid_y'] + ra_hosts = np.rad2deg(hostcat['coord_ra']) + dec_hosts = np.rad2deg(hostcat['coord_dec']) + mag_hosts = hostcat['slot_PsfFlux_mag'] - delta_mag = rng.normal(loc=1, scale=1, size=n_star_hosted_blended_fakes) - mags = mag_hosts + delta_mag - - # Create the table of hosted fakes - hosted_fakes = Table() - hosted_fakes["x"] = x_ssi - hosted_fakes["y"] = y_ssi - hosted_fakes["mag"] = mags - hosted_fakes["ra"] = ra_ssi - hosted_fakes["dec"] = dec_ssi - hosted_fakes["host_id"] = hostcat['id'] - hosted_fakes["host_flux"] = hostcat['slot_PsfFlux_flux'] - hosted_fakes["host_mag"] = hostcat['slot_PsfFlux_mag'] - hosted_fakes["host_ra"] = ra_hosts - hosted_fakes["host_dec"] = dec_hosts - hosted_fakes["delta_ra"] = delta_ra - hosted_fakes["delta_dec"] = delta_dec - hosted_fakes["delta_mag"] = delta_mag - hosted_fakes["source_type"] = "Star" - hosted_fakes["hosted_fake"] = True - hosted_fakes["isVisitSource"] = True - hosted_fakes["isTemplateSource"] = False - hosted_fakes["isBlended"] = True - - blended_fakes = vstack([blended_fakes, hosted_fakes]) + delta_ra, delta_dec = self._draw_offset_components_arcsec( + rng, n_star_hosted_blended_fakes + ) - blended_fakes["injection_id"] = self._make_unique_injection_ids( - len(blended_fakes), - used_ids=catalog["injection_id"], - ) - catalog = vstack([catalog, blended_fakes]) - - if len(catalog) > len(np.unique(catalog["injection_id"])): - self.log.warning("Duplicate injection IDs detected after catalog assembly; reassigning them.") - old_injection_ids = np.asarray(catalog["injection_id"], dtype=np.int64) - new_injection_ids = self._make_unique_injection_ids(len(catalog)) - # re-assign fresh injection ids - catalog["injection_id"] = new_injection_ids - if "twin_id" in catalog.colnames: - id_map = {old_id: new_id for old_id, new_id in zip(old_injection_ids, new_injection_ids)} - catalog["twin_id"] = np.asarray( - [id_map.get(int(twin_id), int(twin_id)) for twin_id in catalog["twin_id"]], - dtype=np.int64, + ra_ssi = ra_hosts + delta_ra / 3600.0 + dec_ssi = dec_hosts + delta_dec / 3600.0 + + x_ssi, y_ssi = wcs.skyToPixelArray(ra_ssi, dec_ssi, degrees=True) + + delta_mag = rng.normal(loc=1, scale=1, size=n_star_hosted_blended_fakes) + mags = mag_hosts + delta_mag + + # Create the table of hosted fakes + hosted_fakes = Table() + hosted_fakes["x"] = x_ssi + hosted_fakes["y"] = y_ssi + hosted_fakes["mag"] = mags + hosted_fakes["ra"] = ra_ssi + hosted_fakes["dec"] = dec_ssi + hosted_fakes["host_id"] = hostcat['id'] + hosted_fakes["host_flux"] = hostcat['slot_PsfFlux_flux'] + hosted_fakes["host_mag"] = hostcat['slot_PsfFlux_mag'] + hosted_fakes["host_ra"] = ra_hosts + hosted_fakes["host_dec"] = dec_hosts + hosted_fakes["delta_ra"] = delta_ra + hosted_fakes["delta_dec"] = delta_dec + hosted_fakes["delta_mag"] = delta_mag + hosted_fakes["source_type"] = "Star" + hosted_fakes["hosted_fake"] = True + hosted_fakes["isVisitSource"] = True + hosted_fakes["isTemplateSource"] = False + hosted_fakes["isBlended"] = True + + blended_fakes = vstack([blended_fakes, hosted_fakes]) + + blended_fakes["injection_id"] = self._make_unique_injection_ids( + len(blended_fakes), + used_ids=catalog["injection_id"], ) + catalog = vstack([catalog, blended_fakes]) + + if len(catalog) > len(np.unique(catalog["injection_id"])): + self.log.warning("Duplicate injection IDs detected after catalog assembly; reassigning them.") + old_injection_ids = np.asarray(catalog["injection_id"], dtype=np.int64) + new_injection_ids = self._make_unique_injection_ids(len(catalog)) + # re-assign fresh injection ids + catalog["injection_id"] = new_injection_ids + if "twin_id" in catalog.colnames: + id_map = {old_id: new_id for old_id, new_id in zip(old_injection_ids, new_injection_ids)} + catalog["twin_id"] = np.asarray( + [id_map.get(int(twin_id), int(twin_id)) for twin_id in catalog["twin_id"]], + dtype=np.int64, + ) catalog["visit"] = visitId catalog["detector"] = detId + if catalog["ra"].unit is not None: + catalog["ra"] = catalog["ra"].value + catalog["dec"] = catalog["dec"].value + catalog["delta_ra"] = catalog["delta_ra"].value + catalog["delta_dec"] = catalog["delta_dec"].value + catalog["host_ra"] = catalog["host_ra"].value + catalog["host_dec"] = catalog["host_dec"].value + return Struct(outputCat=catalog) From e9a13508c3dcfcdf6a774b1c5775ea7627278351 Mon Sep 17 00:00:00 2001 From: BrunoSanchez Date: Thu, 3 Sep 2026 08:11:17 -0700 Subject: [PATCH 3/3] Add blended fakes option to creatApFakes --- python/lsst/ap/pipe/createApFakes.py | 78 +++++++++++--- tests/test_createApFakes.py | 145 +++++++++++++++++++++++---- 2 files changed, 186 insertions(+), 37 deletions(-) diff --git a/python/lsst/ap/pipe/createApFakes.py b/python/lsst/ap/pipe/createApFakes.py index a816cb62..7dc499ff 100644 --- a/python/lsst/ap/pipe/createApFakes.py +++ b/python/lsst/ap/pipe/createApFakes.py @@ -437,6 +437,7 @@ class CreateVisitDetectorFakesConfig( min=-1, ) + class CreateVisitDetectorFakesTask(PipelineTask): """Create and store a set of visit detector fakes for use in AP processing. This task creates a catalog of fake sources that can be used to inject @@ -449,6 +450,39 @@ class CreateVisitDetectorFakesTask(PipelineTask): def __init__(self, **kwargs): super().__init__(**kwargs) self.log = logging.getLogger(__name__) + self._table_dtypes = ( + [ + ('x', ' 0: + cap = min(cap, self.config.maxHostedBlendedFakesTotal) + return min(requested, cap) + def run(self, sourceCat, visit_image): """Create a set of visit detector fakes. @@ -484,6 +539,8 @@ def run(self, sourceCat, visit_image): outputCat : `astropy.table.Table` Catalog of fake sources to draw inputs from. """ + zero_table = Table(dtype=self._table_dtypes) + # Use the visit+detector ids as the random seed. visitId = visit_image.getInfo().getVisitInfo().id detId = visit_image.detector.getId() @@ -690,17 +747,9 @@ def run(self, sourceCat, visit_image): scale=self.config.blendedFakeMagOffset, size=n_blended_fakes ) - blended_fakes["delta_ra"] = rng.uniform( - low=-self.config.blendedFakeMaxOffset, - high=self.config.blendedFakeMaxOffset, - size=n_blended_fakes - ) - blended_fakes["delta_ra"] *= rng.choice([-1, 1], size=n_blended_fakes) - - blended_fakes["delta_dec"] = np.sqrt( - self.config.blendedFakeMaxOffset**2 - blended_fakes["delta_ra"]**2 - ) - blended_fakes["delta_dec"] *= rng.choice([-1, 1], size=n_blended_fakes) + delta_ra, delta_dec = self._draw_offset_components_arcsec(rng, n_blended_fakes) + blended_fakes["delta_ra"] = delta_ra + blended_fakes["delta_dec"] = delta_dec blended_fakes["host_mag"] = blended_fakes["mag"] blended_fakes["mag"] += blended_fakes["delta_mag"] @@ -741,7 +790,8 @@ def run(self, sourceCat, visit_image): ) if n_star_hosted_blended_fakes > 0: - # Build a finite host pool so each host appears at most maxHostedBlendedFakesPerHost times. + # Build a finite host pool so each host appears at most + # maxHostedBlendedFakesPerHost times. host_pool = np.repeat( np.arange(len(star_hosts), dtype=int), self.config.maxHostedBlendedFakesPerHost, @@ -819,8 +869,7 @@ def run(self, sourceCat, visit_image): catalog["host_ra"] = catalog["host_ra"].value catalog["host_dec"] = catalog["host_dec"].value - - return Struct(outputCat=catalog) + return Struct(outputCat=vstack([zero_table, catalog])) def select_hosts(self, sourceCat): """ @@ -900,7 +949,6 @@ def select_host_stars(self, sourceCat): skySourceCut & flagCut & extendednessCut & snrCut].copy() return hostCat - def get_PA_and_axes(self, Ixx, Ixy, Iyy): ''' Calculates the orientation and extent of an object based on its second moments. diff --git a/tests/test_createApFakes.py b/tests/test_createApFakes.py index 5c019d17..a9078350 100644 --- a/tests/test_createApFakes.py +++ b/tests/test_createApFakes.py @@ -30,6 +30,7 @@ import lsst.daf.butler.tests as butlerTests import lsst.geom as geom from astropy.table import Table +from astropy import units as u from lsst.pipe.base import testUtils import lsst.skymap as skyMap import lsst.utils.tests @@ -162,7 +163,8 @@ def testVisitCoaddSubdivision(self): def _make_mock_visit_image(visitId=2024111100094, detId=3, xmin=0, xmax=4096, ymin=0, ymax=4096, - magLim=25.0, ra_center=10.0, dec_center=-1.0): + magLim=25.0, ra_center=10.0, dec_center=-1.0, + return_quantity_angles=False): """Build a minimal MagicMock that satisfies CreateVisitDetectorFakesTask.run.""" img = MagicMock() @@ -195,8 +197,22 @@ def _make_mock_visit_image(visitId=2024111100094, detId=3, def _pix_to_sky(xs, ys, degrees=True): ra = ra_center + xs * 1e-4 dec = dec_center + ys * 1e-4 + if return_quantity_angles: + return np.asarray(ra) * u.deg, np.asarray(dec) * u.deg return np.asarray(ra), np.asarray(dec) + + def _sky_to_pix(ra, dec, degrees=False): + ra_arr = np.asarray(ra) + dec_arr = np.asarray(dec) + if not degrees: + ra_arr = np.rad2deg(ra_arr) + dec_arr = np.rad2deg(dec_arr) + xs = (ra_arr - ra_center) / 1e-4 + ys = (dec_arr - dec_center) / 1e-4 + return np.asarray(xs), np.asarray(ys) + wcs.pixelToSkyArray.side_effect = _pix_to_sky + wcs.skyToPixelArray.side_effect = _sky_to_pix img.getWcs.return_value = wcs # photoCalib — not used in non-hosted paths, but must exist @@ -205,6 +221,43 @@ def _pix_to_sky(xs, ys, degrees=True): return img +def _make_calibrated_source_table(n_sources, host_kind="galaxy", ra_deg=10.0, dec_deg=-1.0): + """Build a calibrated source table that can pass host selection cuts.""" + if host_kind == "galaxy": + size_ext = np.ones(n_sources) + ext = np.ones(n_sources, dtype=int) + elif host_kind == "star": + size_ext = np.full(n_sources, 0.2) + ext = np.zeros(n_sources, dtype=int) + else: + raise ValueError(f"Unknown host_kind={host_kind}") + + return Table({ + "slot_Centroid_x": np.linspace(1000.0, 3000.0, n_sources), + "slot_Centroid_y": np.linspace(1100.0, 3100.0, n_sources), + "slot_ModelFlux_mag": np.full(n_sources, 20.0), + "slot_ModelFlux_flux": np.full(n_sources, 1e4), + "slot_ModelFlux_fluxErr": np.full(n_sources, 100.0), + "slot_PsfFlux_mag": np.full(n_sources, 20.0), + "slot_PsfFlux_flux": np.full(n_sources, 1e4), + "slot_PsfFlux_fluxErr": np.full(n_sources, 100.0), + "slot_Shape_xx": np.full(n_sources, 4.0), + "slot_Shape_xy": np.zeros(n_sources), + "slot_Shape_yy": np.full(n_sources, 4.0), + "id": np.arange(n_sources, dtype=np.int64), + "coord_ra": np.deg2rad(np.full(n_sources, ra_deg)), + "coord_dec": np.deg2rad(np.full(n_sources, dec_deg)), + "sky_source": np.zeros(n_sources, dtype=bool), + "base_ClassificationSizeExtendedness_flag": np.zeros(n_sources, dtype=bool), + "base_ClassificationExtendedness_flag": np.zeros(n_sources, dtype=bool), + "slot_Shape_flag": np.zeros(n_sources, dtype=bool), + "slot_Centroid_flag": np.zeros(n_sources, dtype=bool), + "base_PixelFlags_flag": np.zeros(n_sources, dtype=bool), + "base_ClassificationSizeExtendedness_value": size_ext, + "base_ClassificationExtendedness_value": ext, + }) + + class TestCreateVisitDetectorFakesTask(lsst.utils.tests.TestCase): def setUp(self): @@ -313,27 +366,7 @@ def testHostedFakesSparseHostCap(self): task = CreateVisitDetectorFakesTask(config=cfg) n_hosts = 5 - host_table = Table({ - "slot_Centroid_x": np.full(n_hosts, 2048.0), - "slot_Centroid_y": np.full(n_hosts, 2048.0), - "slot_ModelFlux_mag": np.full(n_hosts, 20.0), - "slot_ModelFlux_flux": np.full(n_hosts, 1e4), - "slot_ModelFlux_fluxErr": np.full(n_hosts, 100.0), - "slot_Shape_xx": np.full(n_hosts, 4.0), - "slot_Shape_xy": np.zeros(n_hosts), - "slot_Shape_yy": np.full(n_hosts, 4.0), - "id": np.arange(n_hosts, dtype=np.int64), - "coord_ra": np.deg2rad(np.full(n_hosts, 10.0)), - "coord_dec": np.deg2rad(np.full(n_hosts, -1.0)), - "sky_source": np.zeros(n_hosts, dtype=bool), - "base_ClassificationSizeExtendedness_flag": np.zeros(n_hosts, dtype=bool), - "base_ClassificationExtendedness_flag": np.zeros(n_hosts, dtype=bool), - "slot_Shape_flag": np.zeros(n_hosts, dtype=bool), - "slot_Centroid_flag": np.zeros(n_hosts, dtype=bool), - "base_PixelFlags_flag": np.zeros(n_hosts, dtype=bool), - "base_ClassificationSizeExtendedness_value": np.ones(n_hosts), - "base_ClassificationExtendedness_value": np.ones(n_hosts, dtype=int), - }) + host_table = _make_calibrated_source_table(n_hosts, host_kind="galaxy") # Patch photoCalib so calibrateCatalog().asAstropy() returns our table img = _make_mock_visit_image() @@ -390,6 +423,74 @@ def testHostedFakesNoHostsWarning(self): # Random fakes still produced self.assertEqual(len(result.outputCat), 10) + def testHostedBlendedFakesPerHostCap(self): + task = self._make_task( + doAddRandomVisitFakes=True, + nRandomFakes=20, + doAddBlendedFakes=True, + fracBlendedFakes=0.999, + fracHostedBlendedFakes=0.999, + maxHostedBlendedFakesPerHost=2, + maxHostedBlendedFakesTotal=-1, + ) + star_hosts = _make_calibrated_source_table(3, host_kind="star") + self.visit_image.getPhotoCalib().calibrateCatalog.return_value.asAstropy.return_value = star_hosts + + cat = task.run(self.source_cat, self.visit_image).outputCat + hosted_blended = cat[cat["isBlended"]] + hosted_mask = ~np.ma.getmaskarray(hosted_blended["host_flux"]) + hosted_star_blended = hosted_blended[hosted_mask] + + # Requested 20, capped to n_hosts * per_host = 3 * 2 = 6. + self.assertEqual(len(hosted_star_blended), 6) + host_ids, counts = np.unique(hosted_star_blended["host_id"], return_counts=True) + self.assertEqual(len(host_ids), 3) + self.assertTrue(np.all(counts <= 2)) + + def testHostedBlendedFakesTotalCap(self): + task = self._make_task( + doAddRandomVisitFakes=True, + nRandomFakes=20, + doAddBlendedFakes=True, + fracBlendedFakes=0.999, + fracHostedBlendedFakes=0.999, + maxHostedBlendedFakesPerHost=10, + maxHostedBlendedFakesTotal=4, + ) + star_hosts = _make_calibrated_source_table(10, host_kind="star") + self.visit_image.getPhotoCalib().calibrateCatalog.return_value.asAstropy.return_value = star_hosts + + cat = task.run(self.source_cat, self.visit_image).outputCat + hosted_blended = cat[cat["isBlended"]] + hosted_mask = ~np.ma.getmaskarray(hosted_blended["host_flux"]) + hosted_star_blended = hosted_blended[hosted_mask] + + self.assertEqual(len(hosted_star_blended), 4) + + def testAngleColumnsAreUnitless(self): + cfg = CreateVisitDetectorFakesConfig() + cfg.doAddRandomVisitFakes = False + cfg.doAddRandomTemplateFakes = False + cfg.doAddHostedFakes = True + cfg.doAddVariableFakes = False + cfg.doAddModelFakes = False + cfg.fracHostedFakes = 0.999 + cfg.minHostedFakes = 5 + task = CreateVisitDetectorFakesTask(config=cfg) + + img = _make_mock_visit_image() + host_table = _make_calibrated_source_table(5, host_kind="galaxy") + img.getPhotoCalib().calibrateCatalog.return_value.asAstropy.return_value = host_table + + cat = task.run(self.source_cat, img).outputCat + for col in ("ra", "dec", "delta_ra", "delta_dec", "host_ra", "host_dec"): + self.assertIn(col, cat.colnames) + self.assertIsNone(cat[col].unit) + + # Sanity check that RA values are in degrees-like scale, not radians. + self.assertTrue(np.all(np.asarray(cat["ra"], dtype=float) > 2.0 * np.pi)) + self.assertTrue(np.all(np.asarray(cat["host_ra"], dtype=float) > 2.0 * np.pi)) + class MemoryTester(lsst.utils.tests.MemoryTestCase): pass