From 2e42f7fd4a4b88a423fc10c05758a1cf5486e830 Mon Sep 17 00:00:00 2001 From: BrunoSanchez Date: Wed, 30 Sep 2026 14:19:20 -0700 Subject: [PATCH 1/3] Update yaml for new fakes creation task --- pipelines/CreateInjectionCatalogs.yaml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pipelines/CreateInjectionCatalogs.yaml b/pipelines/CreateInjectionCatalogs.yaml index 2a7db788..9c8d31d1 100644 --- a/pipelines/CreateInjectionCatalogs.yaml +++ b/pipelines/CreateInjectionCatalogs.yaml @@ -16,11 +16,13 @@ tasks: templateFakeFraction: 0.20 doAddHostedFakes: True fracHostedFakes: 0.1 - minHostedFakes: 20 doAddVariableFakes: True variableFakeFraction: 0.1 variableFakeMean: 0.0 variableFakeStd: 0.5 + doAddBlendedFakes: False + fracBlendedFakes: 0.1 + fracHostedBlendedFakes: 0.2 magMin: 20 magMax: 25 subsets: From 771096014a00c5bf76ebc12cc03e6b668f3bef6c Mon Sep 17 00:00:00 2001 From: BrunoSanchez Date: Wed, 30 Sep 2026 14:25:40 -0700 Subject: [PATCH 2/3] Refactor the fake catalog creation --- python/lsst/ap/pipe/createApFakes.py | 894 ++++++++++++++++++--------- 1 file changed, 603 insertions(+), 291 deletions(-) diff --git a/python/lsst/ap/pipe/createApFakes.py b/python/lsst/ap/pipe/createApFakes.py index fc6328f3..e240dc13 100644 --- a/python/lsst/ap/pipe/createApFakes.py +++ b/python/lsst/ap/pipe/createApFakes.py @@ -34,14 +34,18 @@ from lsst.source.injection import generate_injection_catalog +from dataclasses import dataclass from deprecated.sphinx import deprecated -__all__ = ["CreateRandomApFakesTask", - "CreateRandomApFakesConfig", - "CreateRandomApFakesConnections", - "CreateVisitDetectorFakesTask", - "CreateVisitDetectorFakesConfig", - "CreateVisitDetectorFakesConnections"] + +__all__ = [ + "CreateRandomApFakesTask", + "CreateRandomApFakesConfig", + "CreateRandomApFakesConnections", + "CreateVisitDetectorFakesTask", + "CreateVisitDetectorFakesConfig", + "CreateVisitDetectorFakesConnections", +] class CreateRandomApFakesConnections(PipelineTaskConnections, @@ -328,7 +332,7 @@ class CreateVisitDetectorFakesConfig( dtype=float, default=0.25, min=0, - max=1, + max=0.5001, ) doAddVariableFakes = pexConfig.Field( doc="Whether to add variable fakes to the visit detector.", @@ -341,7 +345,7 @@ class CreateVisitDetectorFakesConfig( dtype=float, default=0.1, min=0, - max=1, + max=0.2501, ) variableFakeMean = pexConfig.RangeField( doc="Mean magnitude variation for variable fakes.", @@ -353,7 +357,7 @@ class CreateVisitDetectorFakesConfig( variableFakeStd = pexConfig.RangeField( doc="Standard deviation of magnitude variation for variable fakes.", dtype=float, - default=0.5, + default=0.1, min=0, max=1, ) @@ -390,13 +394,71 @@ 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.1, + min=0, + max=0.2501, + ) + fracHostedBlendedFakes = pexConfig.RangeField( + doc="Fraction of blended fakes that are hosted by stars.", + dtype=float, + default=0.1, + min=0, + max=0.5001, + ) + 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=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, + ) + + +@dataclass +class FakeGenerationPlan: + n_total: int = 0 + n_science: int = 0 + n_template: int = 0 + n_variable: int = 0 + n_blended: int = 0 + + n_hosted: int = 0 + n_hosted_science: int = 0 + n_hosted_template: int = 0 + n_hosted_variable: int = 0 + n_hosted_blended: int = 0 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 - sources into visit detector images. - """ + """Create visit-detector fakes in explicit planning and assembly stages.""" _DefaultName = "createVisitDetectorFakes" ConfigClass = CreateVisitDetectorFakesConfig @@ -404,315 +466,565 @@ class CreateVisitDetectorFakesTask(PipelineTask): def __init__(self, **kwargs): super().__init__(**kwargs) self.log = logging.getLogger(__name__) + self._table_dtypes = [ + ("x", " 0: - self.log.info(f"Generating random visit fakes with nRandomFakes={self.config.nRandomFakes}.") - n_random_fakes = self.config.nRandomFakes + plan = self._make_full_plan(visit_image) + if plan.n_total <= 0: + raise RuntimeError( + "No fake sources will be generated." + ) + # Build populations one by one + # Bear in mind that fractions describe primary sources/pairs, not catalog rows. + populations = [] + + # Random science fakes + if plan.n_science > 0 and self.config.doAddRandomVisitFakes: + random_science_cat = self._generate_seed_fake_catalog(plan.n_science, visit_image, wcs, rng) + random_science_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=plan.n_science) + random_science_cat["isVisitSource"] = np.ones(plan.n_science, dtype=bool) + random_science_cat["isTemplateSource"] = np.zeros(plan.n_science, dtype=bool) + random_science_cat["injection_id"] = self._make_unique_injection_ids( + len(random_science_cat), used_ids=None) + populations.append(random_science_cat) + + # Hosted science + hosted_science_cat, science_hosts = self._generate_hosted_fake_catalog( + plan.n_hosted_science, sourceCat, wcs, rng, host_type='extended') + if science_hosts is None: + self.log.warning( + "Hosted science fake generation requested, but no valid extended hosts were selected." + ) else: - self.log.info( - f"Generating random visit fakes with randomFakeDensity={self.config.randomFakeDensity}.") - n_random_fakes = self.get_n_fakes_from_density( - visit_image=visit_image, - density=self.config.randomFakeDensity + hosted_science_cat["isVisitSource"] = np.ones(plan.n_hosted_science, dtype=bool) + hosted_science_cat["isTemplateSource"] = np.zeros(plan.n_hosted_science, dtype=bool) + hosted_science_cat["hosted_fake"] = np.ones(plan.n_hosted_science, dtype=bool) + hosted_science_cat["injection_id"] = self._make_unique_injection_ids( + len(hosted_science_cat), used_ids=[cat["injection_id"] for cat in populations]) + populations.append(hosted_science_cat) + + # Template fakes + if plan.n_template > 0 and self.config.doAddRandomTemplateFakes: + random_template_cat = self._generate_seed_fake_catalog(plan.n_template, visit_image, wcs, rng) + random_template_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=plan.n_template) + random_template_cat["isVisitSource"] = np.zeros(plan.n_template, dtype=bool) + random_template_cat["isTemplateSource"] = np.ones(plan.n_template, dtype=bool) + random_template_cat["injection_id"] = self._make_unique_injection_ids( + len(random_template_cat), used_ids=[cat["injection_id"] for cat in populations]) + populations.append(random_template_cat) + + # Hosted template fakes + hosted_template_cat, template_hosts = self._generate_hosted_fake_catalog( + plan.n_hosted_template, sourceCat, wcs, rng, host_type='extended') + if template_hosts is None: + self.log.warning( + "Hosted template fake generation requested, but no valid extended hosts were selected." ) - self.log.info(f"Calculated n_random_fakes={n_random_fakes}.") - - # draw random x-y coordinates - x_ssi = rng.uniform(xmin, xmax, size=n_random_fakes) - y_ssi = rng.uniform(ymin, ymax, size=n_random_fakes) - mags = rng.uniform(self.config.magMin, max_mag, size=n_random_fakes) - ra_ssi, dec_ssi = wcs.pixelToSkyArray(x_ssi, y_ssi, degrees=True) - - random_catalog = Table() - random_catalog["x"] = x_ssi - random_catalog["y"] = y_ssi - random_catalog["mag"] = mags - random_catalog["ra"] = ra_ssi - random_catalog["dec"] = dec_ssi - random_catalog["source_type"] = "Star" - random_catalog["isVisitSource"] = True - random_catalog["isTemplateSource"] = False - catalog_set.append(random_catalog) - # Ignore now the possibility of _just_ template fakes - if self.config.doAddModelFakes: - # Generate model fakes - self.log.info("Not implemented yet model fakes.") - # Placeholder for actual model fake generation logic - pass - - if self.config.doAddHostedFakes: - # Generate hosted fakes - self.log.info("Generating hosted fakes.") - # Select hosts that look like extended sources. - hostcatalog = photoCalib.calibrateCatalog(sourceCat).asAstropy() - hostcatalog = self.select_hosts(hostcatalog) - n_hosts = len(hostcatalog) - if n_hosts == 0: - self.log.warning("Hosted fake generation requested, but no valid hosts were selected.") else: - requested_n_fakes = max( - int(self.config.fracHostedFakes * n_hosts), self.config.minHostedFakes + hosted_template_cat["isVisitSource"] = np.zeros(plan.n_hosted_template, dtype=bool) + hosted_template_cat["isTemplateSource"] = np.ones(plan.n_hosted_template, dtype=bool) + hosted_template_cat["hosted_fake"] = np.ones(plan.n_hosted_template, dtype=bool) + hosted_template_cat["injection_id"] = self._make_unique_injection_ids( + len(hosted_template_cat), used_ids=[cat["injection_id"] for cat in populations]) + populations.append(hosted_template_cat) + + # Variable fakes + if plan.n_variable > 0 and self.config.doAddVariableFakes: + random_variable_cat = self._generate_seed_fake_catalog(plan.n_variable, visit_image, wcs, rng) + random_variable_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=plan.n_variable) + random_variable_cat["isVisitSource"] = np.ones(plan.n_variable, dtype=bool) + random_variable_cat["isTemplateSource"] = np.zeros(plan.n_variable, dtype=bool) + random_variable_cat["isVariable"] = np.ones(plan.n_variable, dtype=bool) + random_variable_cat["injection_id"] = self._make_unique_injection_ids( + len(random_variable_cat), used_ids=[cat["injection_id"] for cat in populations]) + # make the twins for variable sources + twins = self._generate_variable_twins(random_variable_cat, rng) + twins["injection_id"] = self._make_unique_injection_ids( + len(twins), used_ids=[cat["injection_id"] for cat in populations + [random_variable_cat]] + ) + # reciprocate the injection_id for the twins so they can be linked + random_variable_cat["twin_id"] = twins["injection_id"] + populations.append(random_variable_cat) + populations.append(twins) + + # Hosted variable fakes + hosted_variable_cat, variable_hosts = self._generate_hosted_fake_catalog( + plan.n_hosted_variable, sourceCat, wcs, rng, host_type='extended') + if variable_hosts is None: + self.log.warning( + "Hosted variable fake generation requested, but no valid extended hosts were selected." ) - n_fakes = min(requested_n_fakes, n_hosts) - if n_fakes < requested_n_fakes: - self.log.warning( - "Reducing hosted fake count from %d to %d because only %d hosts are available.", - requested_n_fakes, - n_fakes, - n_hosts, - ) - - idx = rng.choice(n_hosts, size=n_fakes, replace=False) - hostcat = hostcatalog[idx] - - x_hosts = hostcat['slot_Centroid_x'] - y_hosts = hostcat['slot_Centroid_y'] - mag_hosts = hostcat['slot_ModelFlux_mag'] - # the units below are pixels and radians - pa, a, b = self.get_PA_and_axes( - hostcat['slot_Shape_xx'], - hostcat['slot_Shape_xy'], - hostcat['slot_Shape_yy'] + else: + hosted_variable_cat["isVisitSource"] = np.ones(plan.n_hosted_variable, dtype=bool) + hosted_variable_cat["isTemplateSource"] = np.zeros(plan.n_hosted_variable, dtype=bool) + hosted_variable_cat["isVariable"] = np.ones(plan.n_hosted_variable, dtype=bool) + hosted_variable_cat["hosted_fake"] = np.ones(plan.n_hosted_variable, dtype=bool) + hosted_variable_cat["injection_id"] = self._make_unique_injection_ids( + len(hosted_variable_cat), used_ids=[cat["injection_id"] for cat in populations]) + # make the twins for hosted variable sources + hosted_twins = self._generate_variable_twins(hosted_variable_cat, rng) + hosted_twins["injection_id"] = self._make_unique_injection_ids( + len(hosted_twins), + used_ids=[cat["injection_id"] for cat in populations + [hosted_variable_cat]] ) - # random radius and angle for the fake around the host - theta = rng.uniform(0, 2 * np.pi, size=n_fakes) - angle = np.sqrt((a*np.cos(theta))**2 + (b*np.sin(theta))**2) - radii = angle * np.sqrt(rng.uniform(0, 6, size=n_fakes)) - - # Polar -> Cartesian wrt the host in the PA coordinate system - xs = radii * np.cos(theta) - ys = radii * np.sin(theta) - - # Retrieve the right position removing the galaxy orientation PA - x_rots = xs * np.cos(pa) - ys * np.sin(pa) - y_rots = xs * np.sin(pa) + ys * np.cos(pa) - - x_ssi = x_hosts + x_rots - y_ssi = y_hosts + y_rots - - # retrieving the global ra dec position of the injection - ra_ssi, dec_ssi = wcs.pixelToSkyArray(x_ssi, y_ssi, degrees=True) - delta_ra = (ra_ssi - np.rad2deg(hostcat['coord_ra'])) * 3600. - delta_dec = (dec_ssi - np.rad2deg(hostcat['coord_dec'])) * 3600. - - delta_mag = rng.normal(loc=1, scale=1, size=n_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_ModelFlux_flux'] - hosted_fakes["host_mag"] = hostcat['slot_ModelFlux_mag'] - hosted_fakes["host_ra"] = np.rad2deg(hostcat['coord_ra']) - hosted_fakes["host_dec"] = np.rad2deg(hostcat['coord_dec']) - hosted_fakes["delta_ra"] = delta_ra - hosted_fakes["delta_dec"] = delta_dec - hosted_fakes["delta_mag"] = delta_mag - hosted_fakes["host_a"] = a - hosted_fakes["host_b"] = b - hosted_fakes["host_pa"] = pa - hosted_fakes["source_type"] = "Star" - hosted_fakes["hosted_fake"] = True - hosted_fakes["isVisitSource"] = True - hosted_fakes["isTemplateSource"] = False - - catalog_set.append(hosted_fakes) - - if not catalog_set: - raise RuntimeError( - "No fake sources were generated. Enable at least one fakes mode or provide usable hosts." + hosted_variable_cat["twin_id"] = hosted_twins["injection_id"] + populations.append(hosted_variable_cat) + populations.append(hosted_twins) + + # Blended fakes + if plan.n_blended > 0 and self.config.doAddBlendedFakes: + n_sources = plan.n_blended // 2 + blended_cat = self._generate_seed_fake_catalog(n_sources, visit_image, wcs, rng) + blended_cat["injection_id"] = self._make_unique_injection_ids( + len(blended_cat), used_ids=[cat["injection_id"] for cat in populations]) + blended_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=n_sources) + blended_cat["isVisitSource"] = rng.choice([True, False], size=n_sources) + blended_cat["isTemplateSource"] = ~blended_cat["isVisitSource"] + blended_cat["isBlended"] = np.ones(n_sources, dtype=bool) + # make the twins for blends + blended_twins = self._generate_blended_twins(blended_cat, rng, wcs) + blended_twins["isVisitSource"] = rng.choice([True, False], size=n_sources) + blended_twins["isTemplateSource"] = ~blended_twins["isVisitSource"] + + blended_twins["injection_id"] = self._make_unique_injection_ids( + len(blended_twins), used_ids=[cat["injection_id"] for cat in populations + [blended_cat]] ) + # reciprocate the injection_id for the twins so they can be linked + blended_cat["blend_twin_id"] = blended_twins["injection_id"] + populations.append(blended_cat) + populations.append(blended_twins) + + # Hosted blended fakes are in stars so we don't add twins. + n_sources = plan.n_hosted_blended + hosted_blended_cat, blend_hosts = self._generate_hosted_fake_catalog( + n_sources, sourceCat, wcs, rng, host_type='star') + if blend_hosts is None: + self.log.warning( + "Hosted blended fake generation requested, but no valid star hosts were selected." + ) + else: + hosted_blended_cat, blend_hosts = self._cap_hosted_blended_count( + hosted_blended_cat, blend_hosts + ) + n_sources = len(hosted_blended_cat) + hosted_blended_cat["injection_id"] = self._make_unique_injection_ids( + len(hosted_blended_cat), used_ids=[cat["injection_id"] for cat in populations]) + hosted_blended_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=n_sources) + hosted_blended_cat["isVisitSource"] = rng.choice([True, False], size=n_sources) + hosted_blended_cat["isTemplateSource"] = ~hosted_blended_cat["isVisitSource"] + hosted_blended_cat["isBlended"] = np.ones(n_sources, dtype=bool) + hosted_blended_cat["hosted_fake"] = np.ones(n_sources, dtype=bool) + populations.append(hosted_blended_cat) + + # Merge all populations into a single catalog + if len(populations) == 0: + raise RuntimeError("No fakes were generated for this visit detector.") + + populations = [self._set_population_defaults(population) for population in populations] + catalog = vstack(populations) if populations else Table(dtype=self._table_dtypes) + + self._validate_catalog(catalog) + self._normalize_catalog(catalog, visit_id, detector_id) + return Struct(outputCat=catalog) - catalog = vstack(catalog_set) - catalog['injection_id'] = self._make_unique_injection_ids(len(catalog)) + def _make_full_plan(self, visit_image): + """Plan exactly the number of fakes to inject, and their types, + based on the configuration and the visit image. + This is prior to knowing available hosts for hosted fakes, so the plan + may be adjusted later.""" - if self.config.doAddRandomTemplateFakes: - is_tmplt_fake = rng.random(len(catalog)) < self.config.templateFakeFraction - catalog["isTemplateSource"] = is_tmplt_fake - catalog["isVisitSource"] = ~is_tmplt_fake + n_total = 0 + plan = FakeGenerationPlan() + if self.config.nRandomFakes > 0: + n_total = self.config.nRandomFakes else: - catalog["isVisitSource"] = True - catalog["isTemplateSource"] = False + n_total = self.get_n_fakes_from_density( + visit_image, self.config.randomFakeDensity + ) + + plan.n_total = n_total + if self.config.doAddRandomTemplateFakes: + plan.n_template = int(n_total * self.config.templateFakeFraction) if self.config.doAddVariableFakes: - # Generate variable fakes by duplicating some fakes and adding the counterpart - # either science or template with a magnitude offset drawn from a normal - # distribution with mean and std defined in the config. - self.log.info("Generating variable fakes.") - n_variable_fakes = int(len(catalog) * self.config.variableFakeFraction) - idx = rng.choice(len(catalog), size=n_variable_fakes, replace=False) - variable_fakes = catalog[idx].copy() - variable_fakes["mag_offset"] = rng.normal( - loc=self.config.variableFakeMean, - scale=self.config.variableFakeStd, - size=n_variable_fakes + plan.n_variable = int(n_total * self.config.variableFakeFraction) + if self.config.doAddBlendedFakes: + plan.n_blended = int(n_total * self.config.fracBlendedFakes) + # we will generate n_blends pairs of blended sources, + # so the total number of blended sources is 2*n_blends + n_blends = plan.n_blended//2 + plan.n_blended = n_blends * 2 + if self.config.doAddRandomVisitFakes: + n_science = n_total - plan.n_template - plan.n_variable - plan.n_blended + if n_science < 0: + self.log.warning( + "Requested number of fakes exceeds the total planned. Adjusting n_science to 0." + ) + n_science = 0 + plan.n_science = n_science + + if self.config.doAddHostedFakes: + plan.n_hosted = int(n_total * self.config.fracHostedFakes) + plan.n_hosted_science = int(plan.n_science * self.config.fracHostedFakes) + plan.n_science = plan.n_science - plan.n_hosted_science + plan.n_hosted_template = int(plan.n_template * self.config.fracHostedFakes) + plan.n_template = plan.n_template - plan.n_hosted_template + plan.n_hosted_variable = int(plan.n_variable * self.config.fracHostedFakes) + plan.n_variable = plan.n_variable - plan.n_hosted_variable + + if self.config.fracHostedBlendedFakes > 0: + plan.n_hosted_blended = int(plan.n_blended * self.config.fracHostedBlendedFakes) + plan.n_hosted_blended = plan.n_hosted_blended + plan.n_blended = plan.n_blended - plan.n_hosted_blended + + return plan + + def _set_population_defaults(self, catalog): + """Ensure population flags and relationship IDs are concrete values.""" + defaults = { + "isVariable": False, + "isBlended": False, + "twin_id": 0, + "blend_twin_id": 0, + } + for name, default in defaults.items(): + if name not in catalog.colnames: + catalog[name] = np.full(len(catalog), default) + else: + values = np.ma.asarray(catalog[name]) + catalog[name] = np.ma.filled(values, default) + return catalog + + def _generate_seed_fake_catalog(self, n_fakes, visit_image, wcs, rng): + """Generate a seed catalog of random fakes with positions and magnitudes.""" + if n_fakes <= 0: + return Table(dtype=self._table_dtypes, masked=False) + + bbox = visit_image.getBBox() + x = rng.uniform(bbox.getMinX(), bbox.getMaxX(), size=n_fakes) + y = rng.uniform(bbox.getMinY(), bbox.getMaxY(), size=n_fakes) + ra, dec = wcs.pixelToSkyArray(x, y, degrees=True) + + zero_table = Table(dtype=self._table_dtypes, masked=False) + catalog = Table(masked=False) + catalog["x"] = x + catalog["y"] = y + catalog["ra"] = ra + catalog["dec"] = dec + catalog["source_type"] = "Star" + catalog = vstack([zero_table, catalog]) + + return catalog + + def _build_hosted_population(self, hosts, wcs, rng, host_type='extended'): + """Build a population of hosted fakes around selected hosts.""" + if host_type == 'extended': + pa, a, b = self.get_PA_and_axes( + hosts["slot_Shape_xx"], hosts["slot_Shape_xy"], hosts["slot_Shape_yy"] ) - variable_fakes["mag"] += variable_fakes["mag_offset"] - # we flip the source, so for example if it was a visit, we trasnform it into a template - # with the idea of having duplicate injections, in the same location - variable_fakes["isVisitSource"] = ~variable_fakes["isVisitSource"] - variable_fakes["isTemplateSource"] = ~variable_fakes["isTemplateSource"] - - variable_fakes["twin_id"] = variable_fakes["injection_id"] - variable_fakes["injection_id"] = self._make_unique_injection_ids( - len(variable_fakes), - used_ids=catalog["injection_id"], + flux_type = "slot_ModelFlux_flux" + mag_type = "slot_ModelFlux_mag" + else: + pa, a, b = 0, 1, 1 # For stars, we can treat them as circular with no orientation + flux_type = "slot_PsfFlux_flux" + mag_type = "slot_PsfFlux_mag" + + n_fakes = len(hosts) + theta = rng.uniform(0, 2 * np.pi, size=n_fakes) + radii = np.sqrt((a * np.cos(theta)) ** 2 + (b * np.sin(theta)) ** 2) + radii *= np.sqrt(rng.uniform(0, 6, size=n_fakes)) + x = hosts["slot_Centroid_x"] + radii * np.cos(theta) * np.cos(pa) - radii * np.sin(theta) * np.sin(pa) + y = hosts["slot_Centroid_y"] + radii * np.cos(theta) * np.sin(pa) + radii * np.sin(theta) * np.cos(pa) + ra, dec = wcs.pixelToSkyArray(x, y, degrees=True) + + zero_table = Table(dtype=self._table_dtypes, masked=False) + catalog = Table(masked=False) + catalog["x"] = x + catalog["y"] = y + # magnitudes are drawn from a normal distribution around host magnitude + catalog["mag"] = hosts[mag_type] + rng.normal(1, 1, size=n_fakes) + catalog["ra"] = ra + catalog["dec"] = dec + catalog["host_id"] = hosts["id"] + catalog["host_flux"] = hosts[flux_type] + catalog["host_mag"] = hosts[mag_type] + catalog["host_ra"] = np.rad2deg(hosts["coord_ra"]) + catalog["host_dec"] = np.rad2deg(hosts["coord_dec"]) + catalog["delta_ra"] = (ra - catalog["host_ra"]) * 3600 + catalog["delta_dec"] = (dec - catalog["host_dec"]) * 3600 + catalog["delta_mag"] = catalog["mag"] - catalog["host_mag"] + catalog["host_a"] = a + catalog["host_b"] = b + catalog["host_pa"] = pa + catalog["source_type"] = "Star" + catalog["host_type"] = host_type + catalog["hosted_fake"] = np.ones(len(catalog), dtype=bool) + catalog = vstack([zero_table, catalog]) + + return catalog, hosts + + def _generate_hosted_fake_catalog(self, n_fakes, sourceCat, wcs, rng, host_type='extended'): + """Generate a catalog of hosted fakes around selected host galaxies or stars.""" + if n_fakes <= 0: + return Table(dtype=self._table_dtypes, masked=False), None + + host_catalog = self.select_hosts(sourceCat, host_type=host_type) + if len(host_catalog) == 0: + self.log.warning("No valid hosts were selected for hosted fakes.") + return Table(dtype=self._table_dtypes, masked=False), None + + # we allow more than one fake per host since we can deal with blending now + if n_fakes > len(host_catalog): + self.log.warning( + "Requested %d hosted fakes, but only %d valid hosts were found. " + "Some hosts will have multiple fakes.", + n_fakes, + len(host_catalog), ) - # create column of isVariable flag - catalog["isVariable"] = np.where(np.isin(np.arange(len(catalog)), idx), True, False) - variable_fakes["isVariable"] = True - - catalog = vstack([catalog, variable_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, - ) + idx = rng.choice(len(host_catalog), size=n_fakes, replace=True) + else: + idx = rng.choice(len(host_catalog), size=n_fakes, replace=False) + + selected_hosts = host_catalog[idx] + return self._build_hosted_population(selected_hosts, wcs, rng, host_type=host_type) + + def _generate_variable_twins(self, variable_catalog, rng): + """Generate twin sources for variable fakes with magnitude offsets.""" + n_twins = len(variable_catalog) + if n_twins <= 0: + return Table(dtype=self._table_dtypes, masked=False) + # magnitudes are drawn from a normal distribution around the variable source's magnitude + twins = variable_catalog.copy() + twins["mag_offset"] = rng.normal( + loc=self.config.variableFakeMean, + scale=self.config.variableFakeStd, + size=n_twins + ) + twins["mag"] += twins["mag_offset"] + twins["isVisitSource"] = ~twins["isVisitSource"] + twins["isTemplateSource"] = ~twins["isTemplateSource"] + twins["twin_id"] = twins["injection_id"] + return twins + + def _generate_blended_twins(self, blended_catalog, rng, wcs): + """Generate twin sources for blended fakes with magnitude and positional offsets.""" + n_twins = len(blended_catalog) + if n_twins <= 0: + return Table(dtype=self._table_dtypes, masked=False) + + twins = blended_catalog.copy() + twins["mag_offset"] = rng.normal( + loc=0.0, + scale=self.config.blendedFakeMagOffset, + size=n_twins + ) + # we are reusing the delta_ra and delta_dec columns to store the offsets for the twins + # original host information will be in the entry for the parent blend source + # stored in the blend_parent_id column + delta_ra, delta_dec = self._draw_offset_components_arcsec(rng, n_twins) + twins["delta_ra"] = delta_ra + twins["delta_dec"] = delta_dec + twins["mag"] += twins["mag_offset"] - catalog["visit"] = visitId - catalog["detector"] = detId + twins["ra"] += twins["delta_ra"] * np.cos(np.deg2rad(twins["dec"])) / 3600 + twins["dec"] += twins["delta_dec"] / 3600 + twins["x"], twins["y"] = wcs.skyToPixelArray(twins["ra"], twins["dec"], degrees=True) - return Struct(outputCat=catalog) + twins["blend_twin_id"] = twins["injection_id"] + return twins - def select_hosts(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. - """ + def _validate_catalog(self, catalog): + """Validate IDs and reciprocal source relationships.""" + injection_ids = np.asarray(catalog["injection_id"], dtype=np.int64) + if len(np.unique(injection_ids)) != len(injection_ids): + raise RuntimeError("Fake injection IDs are not unique.") - # Avoid calibration stars or psf stars; remove flagged sources sky_sources - skySourceCut = ~sourceCat['sky_source'] + id_set = set(injection_ids) - 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'] + variable = np.asarray(catalog["isVariable"], dtype=bool) + blended = np.asarray(catalog["isBlended"], dtype=bool) - extendednessCut = sourceCat['base_ClassificationSizeExtendedness_value'] > 0.9 - extendednessCut &= sourceCat['base_ClassificationExtendedness_value'] == 1 + if np.any(variable & blended): + raise RuntimeError("Variable and blended fake populations overlap.") - snrCut = sourceCat['slot_ModelFlux_flux']/sourceCat['slot_ModelFlux_fluxErr'] > 15 + twin_ids = np.asarray(catalog["twin_id"], dtype=np.int64) + variable_twin_ids = twin_ids[variable] - hostCat = sourceCat[ - skySourceCut & flagCut & extendednessCut & snrCut].copy() - return hostCat + if np.any(variable_twin_ids <= 0): + raise RuntimeError("Variable fake twin links are incomplete.") - def get_PA_and_axes(self, Ixx, Ixy, Iyy): - ''' - Calculates the orientation and extent of an object based on its second moments. - - Parameters: - Ixx (float): Second moment of the object along the x-axis. Often in degree² - Ixy (float): Second moment of the object along the x and y-axes. Often in degree² - Iyy (float): Second moment of the object along the y-axis. Often in degree² - - Returns: - tuple: A tuple containing: - - theta (float): The orientation angle of the object in radians. - - a (float): The semi-major axis length of the object. - - b (float): The semi-minor axis length of the object. - ''' - # Calculate position angle (orientation) - theta = 0.5 * np.arctan2(2 * Ixy, Ixx - Iyy) + if not set(variable_twin_ids).issubset(id_set): + raise RuntimeError("Variable twin IDs do not refer to catalog rows.") + + variable_rows_by_id = { + int(injection_id): index + for index, injection_id in enumerate(injection_ids) + if variable[index] + } - # Calculate eigenvalues of the moment matrix - term1 = (Ixx + Iyy) / 2 - term2 = np.sqrt(((Ixx - Iyy) / 2) ** 2 + Ixy ** 2) - lambda1 = term1 + term2 - lambda2 = term1 - term2 + for index in np.flatnonzero(variable): + twin_index = variable_rows_by_id.get(int(twin_ids[index])) + if twin_index is None: + raise RuntimeError("Variable twin ID does not identify a variable row.") + if twin_ids[twin_index] != injection_ids[index]: + raise RuntimeError("Variable twin links are not reciprocal.") - a = np.sqrt(lambda1) - b = np.sqrt(lambda2) + # removing hosted blends from the check + non_hosted_blends = blended & ~np.asarray(catalog["hosted_fake"], dtype=bool) + blend_twin_ids = np.asarray(catalog["blend_twin_id"], dtype=np.int64) + blended_twin_ids = blend_twin_ids[non_hosted_blends] - return theta, a, b + if np.any(blended_twin_ids <= 0): + raise RuntimeError("Blended fake twin links are incomplete.") + + if not set(blended_twin_ids).issubset(id_set): + raise RuntimeError("Blend twin IDs do not refer to catalog rows.") + + rows_by_id = { + int(injection_id): index + for index, injection_id in enumerate(injection_ids) + } + + for index in np.flatnonzero(blend_twin_ids): + twin_index = rows_by_id[int(blend_twin_ids[index])] + if blend_twin_ids[twin_index] != injection_ids[index]: + raise RuntimeError("Blend twin links are not reciprocal.") + + def _normalize_catalog(self, catalog, visit_id, detector_id): + """Normalize catalog columns to remove units and add visit/detector identifiers.""" + catalog["visit"] = visit_id + catalog["detector"] = detector_id + for name in ("ra", "dec", "delta_ra", "delta_dec", "host_ra", "host_dec"): + if name in catalog.colnames and catalog[name].unit is not None: + catalog[name] = catalog[name].value + + def select_hosts(self, sourceCat, host_type='extended'): + """Select hosts from the source catalog.""" + sky = ~sourceCat["sky_source"] + flags = ~sourceCat["base_ClassificationSizeExtendedness_flag"] + flags &= ~sourceCat["base_ClassificationExtendedness_flag"] + flags &= ~sourceCat["slot_Shape_flag"] + flags &= ~sourceCat["slot_Centroid_flag"] + flags &= ~sourceCat["base_PixelFlags_flag"] + + if host_type == 'extended': + extended = sourceCat["base_ClassificationSizeExtendedness_value"] > 0.9 + extended &= sourceCat["base_ClassificationExtendedness_value"] == 1 + snr = sourceCat["slot_ModelFlux_flux"] / sourceCat["slot_ModelFlux_fluxErr"] > 15 + return sourceCat[sky & flags & extended & snr].copy() + elif host_type == 'star': + stellar = sourceCat["base_ClassificationSizeExtendedness_value"] < 0.9 + stellar &= sourceCat["base_ClassificationExtendedness_value"] == 0 + snr = sourceCat["slot_PsfFlux_flux"] / sourceCat["slot_PsfFlux_fluxErr"] > 30 + return sourceCat[sky & flags & stellar & snr].copy() + else: + raise ValueError("host_type must be 'extended' or 'star'.") + + def get_PA_and_axes(self, Ixx, Ixy, Iyy): + """Compute the position angle and semi-major/minor axes from second moments.""" + theta = 0.5 * np.arctan2(2 * Ixy, Ixx - Iyy) + term = np.sqrt(((Ixx - Iyy) / 2) ** 2 + Ixy ** 2) + return theta, np.sqrt((Ixx + Iyy) / 2 + term), np.sqrt((Ixx + Iyy) / 2 - term) + + def _cap_hosted_blended_count(self, hosted_fakes, hosts): + """Cap the number of hosted blended fakes based on configuration limits.""" + if len(hosted_fakes) == 0 or len(hosts) == 0: + return hosted_fakes, hosts + + # check if we have multiple ids for the same host and cap the number of fakes per host + unique_hosts, counts = np.unique(hosted_fakes["host_id"], return_counts=True) + capped_hosts = [] + for host_id, count in zip(unique_hosts, counts): + if count > self.config.maxHostedBlendedFakesPerHost: + self.log.warning( + "Capping hosted blended fakes for host_id %d from %d to %d.", + host_id, count, self.config.maxHostedBlendedFakesPerHost + ) + capped_hosts.append(host_id) + row_ids = np.where(hosted_fakes["host_id"] == host_id) + row_ids_to_remove = row_ids[0][self.config.maxHostedBlendedFakesPerHost:] + hosted_fakes.remove_rows(row_ids_to_remove) + if len(hosted_fakes) > self.config.maxHostedBlendedFakesTotal: + self.log.warning( + "Capping total hosted blended fakes from %d to %d.", + len(hosted_fakes), self.config.maxHostedBlendedFakesTotal + ) + hosted_fakes = hosted_fakes[:self.config.maxHostedBlendedFakesTotal] + return hosted_fakes, hosts + + def _draw_offset_components_arcsec(self, rng, n_points): + """Draw random offsets in arcseconds for blended fakes. + The offsets are drawn from a uniform distribution in radius and angle.""" + if n_points <= 0: + return np.zeros(0), np.zeros(0) + rmin = max(0.0, self.config.blendedFakeMinOffset) + rmax = max(rmin, self.config.blendedFakeMaxOffset) + radius = np.sqrt(rng.uniform(rmin**2, rmax**2, n_points)) + theta = rng.uniform(0, 2 * np.pi, n_points) + return radius * np.cos(theta), radius * np.sin(theta) def get_n_fakes_from_density(self, visit_image, density): - """Calculate the area of the injection limits in square degrees based on the RA and Dec limits.""" - image_area = visit_image.getConvexPolygon().getBoundingBox().getArea() - image_area *= (180 / np.pi) ** 2 - number = np.round(density * image_area).astype(int) - return number + area = visit_image.getConvexPolygon().getBoundingBox().getArea() + return np.round(density * area * (180 / np.pi) ** 2).astype(int) + + def _make_unique_injection_ids(self, n_ids, used_ids=None): + """Generate unique 24-bit IDs for a catalog considering IDs already in use.""" + new_ids = set() + if used_ids is not None: + # compile used ids from the existing catalog to avoid collisions + if isinstance(used_ids, list): + used_ids = np.concatenate(used_ids) + used_ids = set(used_ids) + else: + used_ids = set() + # check that the new ids are not in the used set already + while len(new_ids) < n_ids: + candidate_id = uuid.uuid4().int & ((1 << 24) - 1) + if candidate_id not in used_ids.union(new_ids): + new_ids.add(candidate_id) + + return np.asarray(list(new_ids), dtype=np.int64) From e697c9ccc54040345474764ee087b930b9518158 Mon Sep 17 00:00:00 2001 From: BrunoSanchez Date: Wed, 30 Sep 2026 14:26:15 -0700 Subject: [PATCH 3/3] Fix tests for new fake catalog creation tool --- python/lsst/ap/pipe/createApFakes.py | 2 +- tests/test_createApFakes.py | 180 +++++++++++++++++++-------- 2 files changed, 127 insertions(+), 55 deletions(-) diff --git a/python/lsst/ap/pipe/createApFakes.py b/python/lsst/ap/pipe/createApFakes.py index e240dc13..0e24865d 100644 --- a/python/lsst/ap/pipe/createApFakes.py +++ b/python/lsst/ap/pipe/createApFakes.py @@ -584,7 +584,7 @@ def run(self, sourceCat, visit_image): # make the twins for variable sources twins = self._generate_variable_twins(random_variable_cat, rng) twins["injection_id"] = self._make_unique_injection_ids( - len(twins), used_ids=[cat["injection_id"] for cat in populations + [random_variable_cat]] + len(twins), used_ids=[cat["injection_id"] for cat in populations + [random_variable_cat]] ) # reciprocate the injection_id for the twins so they can be linked random_variable_cat["twin_id"] = twins["injection_id"] diff --git a/tests/test_createApFakes.py b/tests/test_createApFakes.py index 5c019d17..7ad0428b 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): @@ -297,57 +350,6 @@ def testEmptyCatalogRaises(self): with self.assertRaises(RuntimeError): task.run(self.source_cat, self.visit_image) - # ------------------------------------------------------------------ - # Hosted fakes: sparse-host cap (fewer hosts than minHostedFakes) - # ------------------------------------------------------------------ - def testHostedFakesSparseHostCap(self): - """When n_hosts < minHostedFakes the task should clamp, not crash.""" - cfg = CreateVisitDetectorFakesConfig() - cfg.doAddRandomVisitFakes = False - cfg.doAddHostedFakes = True - cfg.doAddRandomTemplateFakes = False - cfg.doAddVariableFakes = False - cfg.doAddModelFakes = False - cfg.fracHostedFakes = 0.99 # near-100 % but within [0, 1) range - cfg.minHostedFakes = 50 # request 50 but we'll only supply 5 hosts - 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), - }) - - # Patch photoCalib so calibrateCatalog().asAstropy() returns our table - img = _make_mock_visit_image() - img.getPhotoCalib().calibrateCatalog.return_value.asAstropy.return_value = host_table - - result = task.run(self.source_cat, img) - cat = result.outputCat - - # Number of hosted fakes should be clamped to n_hosts, not minHostedFakes - self.assertEqual(len(cat), n_hosts) - self.assertIn("hosted_fake", cat.colnames) - self.assertTrue(np.all(cat["hosted_fake"])) - self.assertEqual(len(np.unique(cat["injection_id"])), n_hosts) - # ------------------------------------------------------------------ # Hosted fakes: zero valid hosts — warning issued, no crash when # random fakes are also enabled as fallback @@ -355,8 +357,9 @@ def testHostedFakesSparseHostCap(self): def testHostedFakesNoHostsWarning(self): task = self._make_task( doAddRandomVisitFakes=True, - nRandomFakes=10, + nRandomFakes=20, doAddHostedFakes=True, + fracHostedFakes=0.1, ) empty_table = Table({ "slot_Centroid_x": np.array([]), @@ -388,7 +391,76 @@ def testHostedFakesNoHostsWarning(self): self.assertTrue(any("no valid hosts" in msg.lower() for msg in cm.output)) # Random fakes still produced - self.assertEqual(len(result.outputCat), 10) + self.assertEqual(len(result.outputCat), 18) + + def testHostedBlendedFakesPerHostCap(self): + task = self._make_task( + doAddRandomVisitFakes=True, + nRandomFakes=100, + doAddBlendedFakes=True, + fracBlendedFakes=0.25, + fracHostedBlendedFakes=0.5, + maxHostedBlendedFakesPerHost=2, + maxHostedBlendedFakesTotal=-1, + ) + star_hosts = _make_calibrated_source_table(5, 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"] & cat["hosted_fake"]] + hosted_mask = ~np.ma.getmaskarray(hosted_blended["host_flux"]) + hosted_star_blended = hosted_blended[hosted_mask] + + # Requested 20, capped to n_hosts * per_host = 5 * 2 = 6 + self.assertLessEqual(len(hosted_star_blended), 8) + host_ids, counts = np.unique(hosted_star_blended["host_id"], return_counts=True) + self.assertLessEqual(len(host_ids), 5) + self.assertTrue(np.all(counts[host_ids > 0] <= 2)) + + def testHostedBlendedFakesTotalCap(self): + task = self._make_task( + doAddRandomVisitFakes=True, + nRandomFakes=300, + doAddBlendedFakes=True, + fracBlendedFakes=0.25, + fracHostedBlendedFakes=0.5, + 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"] & cat["hosted_fake"]] + 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 = True + cfg.doAddRandomTemplateFakes = False + cfg.doAddHostedFakes = True + cfg.doAddVariableFakes = False + cfg.doAddModelFakes = False + cfg.fracHostedFakes = 0.5 + 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)) + hosted_cat = cat[cat['hosted_fake']] + self.assertTrue(np.all(np.asarray(hosted_cat["host_ra"], dtype=float) > 2.0 * np.pi)) class MemoryTester(lsst.utils.tests.MemoryTestCase):