Coverage for python/lsst/ap/pipe/createApFakes.py: 88%
455 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-11-05 02:05 -0800
1# This file is part of ap_pipe.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22import numpy as np
23import pandas as pd
24import uuid
26import logging
27from astropy.table import Table, vstack
29import lsst.pex.config as pexConfig
30from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct
31import lsst.pipe.base.connectionTypes as connTypes
32from lsst.pipe.tasks.insertFakes import InsertFakesConfig
33from lsst.skymap import BaseSkyMap
35from lsst.source.injection import generate_injection_catalog
37from dataclasses import dataclass
38from deprecated.sphinx import deprecated
41__all__ = [
42 "CreateRandomApFakesTask",
43 "CreateRandomApFakesConfig",
44 "CreateRandomApFakesConnections",
45 "CreateVisitDetectorFakesTask",
46 "CreateVisitDetectorFakesConfig",
47 "CreateVisitDetectorFakesConnections",
48]
51class CreateRandomApFakesConnections(PipelineTaskConnections,
52 dimensions=("tract", "skymap")):
53 skyMap = connTypes.Input(
54 doc="Input definition of geometry/bbox and projection/wcs for "
55 "template exposures",
56 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
57 dimensions=("skymap",),
58 storageClass="SkyMap",
59 )
60 fakeCat = connTypes.Output(
61 doc="Catalog of fake sources to draw inputs from.",
62 name="fakeSourceCat",
63 storageClass="DataFrame",
64 dimensions=("tract", "skymap")
65 )
68@deprecated(
69 reason="This task will be removed in v28.0 as it is replaced by `source_injection` tasks.",
70 version="v28.0",
71 category=FutureWarning,
72)
73class CreateRandomApFakesConfig(
74 InsertFakesConfig,
75 pipelineConnections=CreateRandomApFakesConnections):
76 """Config for CreateRandomApFakesTask. Copy from the InsertFakesConfig to
77 assert that columns created with in this task match that those expected in
78 the InsertFakes and related tasks.
79 """
80 fakeDensity = pexConfig.RangeField(
81 doc="Goal density of random fake sources per square degree. Default "
82 "value is roughly the density per square degree for ~10k sources "
83 "visit.",
84 dtype=float,
85 default=1000,
86 min=0,
87 )
88 filterSet = pexConfig.ListField(
89 doc="Set of Abstract filter names to produce magnitude columns for.",
90 dtype=str,
91 default=["u", "g", "r", "i", "z", "y"],
92 )
93 fraction = pexConfig.RangeField(
94 doc="Fraction of the created source that should be inserted into both "
95 "the visit and template images. Values less than 1 will result in "
96 "(1 - fraction) / 2 inserted into only visit or the template.",
97 dtype=float,
98 default=1/3,
99 min=0,
100 max=1,
101 )
102 magMin = pexConfig.RangeField(
103 doc="Minimum magnitude the mag distribution. All magnitudes requested "
104 "are set to the same value.",
105 dtype=float,
106 default=20,
107 min=1,
108 max=40,
109 )
110 magMax = pexConfig.RangeField(
111 doc="Maximum magnitude the mag distribution. All magnitudes requested "
112 "are set to the same value.",
113 dtype=float,
114 default=30,
115 min=1,
116 max=40,
117 )
118 visitSourceFlagCol = pexConfig.Field(
119 doc="Name of the column flagging objects for insertion into the visit "
120 "image.",
121 dtype=str,
122 default="isVisitSource"
123 )
124 templateSourceFlagCol = pexConfig.Field(
125 doc="Name of the column flagging objects for insertion into the "
126 "template image.",
127 dtype=str,
128 default="isTemplateSource"
129 )
132@deprecated(
133 reason="This task will be removed in v28.0 as it is replaced by `source_injection` tasks.",
134 version="v28.0",
135 category=FutureWarning,
136)
137class CreateRandomApFakesTask(PipelineTask):
138 """Create and store a set of spatially uniform star fakes over the sphere
139 for use in AP processing. Additionally assign random magnitudes to said
140 fakes and assign them to be inserted into either a visit exposure or
141 template exposure.
142 """
144 _DefaultName = "createApFakes"
145 ConfigClass = CreateRandomApFakesConfig
147 def runQuantum(self, butlerQC, inputRefs, outputRefs):
148 inputs = butlerQC.get(inputRefs)
149 inputs["tractId"] = butlerQC.quantum.dataId["tract"]
151 outputs = self.run(**inputs)
152 butlerQC.put(outputs, outputRefs)
154 def run(self, tractId, skyMap):
155 """Create a set of uniform random points that covers a tract.
157 Parameters
158 ----------
159 tractId : `int`
160 Tract id to produce randoms over.
161 skyMap : `lsst.skymap.SkyMap`
162 Skymap to produce randoms over.
164 Returns
165 -------
166 randoms : `pandas.DataFrame`
167 Catalog of random points covering the given tract. Follows the
168 columns and format expected in `lsst.pipe.tasks.InsertFakes`.
169 """
171 tract = skyMap.generateTract(tractId)
172 tractArea = tract.getOuterSkyPolygon().getBoundingBox().getArea()
173 tractArea *= (180 / np.pi) ** 2
174 tractWcs = tract.getWcs()
175 vertexList = tract.getVertexList()
176 vertexRas = [vertex.getRa().asDegrees() for vertex in vertexList]
177 vertexDecs = [vertex.getDec().asDegrees() for vertex in vertexList]
179 catalog = generate_injection_catalog(
180 ra_lim=sorted([np.min(vertexRas), np.max(vertexRas)]),
181 dec_lim=sorted([np.min(vertexDecs), np.max(vertexDecs)]),
182 mag_lim=(self.config.magMin, self.config.magMax),
183 density=self.config.fakeDensity,
184 source_type="Star",
185 seed=str(tractId),
186 wcs=tractWcs
187 )
189 nFakes = len(catalog)
191 self.log.info(
192 f"Creating {nFakes} star fakes over tractId={tractId} with "
193 f" RA in ({sorted([np.min(vertexRas), np.max(vertexRas)])} "
194 f" Dec in ({sorted([np.min(vertexDecs), np.max(vertexDecs)])}), "
195 f"area={tractArea:.4f} deg^2 and "
196 f"magnitude range: [{self.config.magMin, self.config.magMax}]")
198 onesColumn = np.ones(nFakes, dtype="float")
199 zerosColumn = np.zeros(nFakes, dtype="float")
200 # Concatenate the data and add dummy values for the unused variables.
201 # Set all data to PSF like objects.
202 mags = np.asarray(catalog["mag"], dtype=float)
203 randData = {
204 "fakeId": [uuid.uuid4().int & (1 << 64) - 1 for n in range(nFakes)],
205 self.config.ra_col: np.asarray(catalog["ra"], dtype=float),
206 self.config.dec_col: np.asarray(catalog["dec"], dtype=float),
207 **self.createVisitCoaddSubdivision(nFakes),
208 **self.createMagnitudeColumns(mags),
209 self.config.disk_semimajor_col: onesColumn,
210 self.config.bulge_semimajor_col: onesColumn,
211 self.config.disk_n_col: onesColumn,
212 self.config.bulge_n_col: onesColumn,
213 self.config.disk_axis_ratio_col: onesColumn,
214 self.config.bulge_axis_ratio_col: onesColumn,
215 self.config.disk_pa_col: zerosColumn,
216 self.config.bulge_pa_col: onesColumn,
217 self.config.sourceType: np.asarray(catalog["source_type"], dtype=str),
218 "source_type": np.asarray(catalog["source_type"], dtype=str),
219 "injection_id": np.asarray(catalog["injection_id"], dtype=np.int64)
220 }
222 return Struct(fakeCat=pd.DataFrame(data=randData))
224 def createVisitCoaddSubdivision(self, nFakes):
225 """Assign a given fake either a visit image or coadd or both based on
226 the ``faction`` config value.
228 Parameters
229 ----------
230 nFakes : `int`
231 Number of fakes to create.
233 Returns
234 -------
235 output : `dict`[`str`, `numpy.ndarray`]
236 Dictionary of boolean arrays specifying which image to put a
237 given fake into.
238 """
239 nBoth = int(self.config.fraction * nFakes)
240 nOnly = int((1 - self.config.fraction) / 2 * nFakes)
241 isVisitSource = np.zeros(nFakes, dtype=bool)
242 isTemplateSource = np.zeros(nFakes, dtype=bool)
243 if nBoth > 0: 243 ↛ 246line 243 didn't jump to line 246 because the condition on line 243 was always true
244 isVisitSource[:nBoth] = True
245 isTemplateSource[:nBoth] = True
246 if nOnly > 0: 246 ↛ 250line 246 didn't jump to line 250 because the condition on line 246 was always true
247 isVisitSource[nBoth:(nBoth + nOnly)] = True
248 isTemplateSource[(nBoth + nOnly):] = True
250 return {self.config.visitSourceFlagCol: isVisitSource,
251 self.config.templateSourceFlagCol: isTemplateSource}
253 def createMagnitudeColumns(self, mags):
254 """Create magnitude columns from a 1D magnitude array.
256 Parameters
257 ----------
258 mags : `numpy.ndarray`
259 Magnitudes to copy to all configured filter-band columns.
261 Returns
262 -------
263 randMags : `dict`[`str`, `numpy.ndarray`]
264 Dictionary containing per-band magnitudes plus a ``mag`` column
265 compatible with ``source_injection`` catalogs.
266 """
267 randMags = {}
268 for fil in self.config.filterSet:
269 randMags[self.config.mag_col % fil] = mags
270 randMags["mag"] = mags
271 return randMags
274class CreateVisitDetectorFakesConnections(
275 PipelineTaskConnections,
276 dimensions=("instrument",
277 "visit",
278 "detector")):
280 sourceCat = connTypes.Input(
281 doc="Catalog of sources detected on the calibrated exposure; ",
282 name="single_visit_star_reprocessed_footprints",
283 storageClass="SourceCatalog",
284 dimensions=["instrument", "visit", "detector"],
285 )
286 visit_image = connTypes.Input(
287 doc="Calibrated exposure to inject synthetic sources into.",
288 name="preliminary_visit_image",
289 storageClass="ExposureF",
290 dimensions=["instrument", "visit", "detector"],
291 )
292 outputCat = connTypes.Output(
293 doc="Catalog of fake sources to draw inputs from.",
294 name="VisitDetectorFakeSourceCat",
295 storageClass="ArrowAstropy",
296 dimensions=["instrument", "visit", "detector"],
297 )
300class CreateVisitDetectorFakesConfig(
301 PipelineTaskConfig,
302 pipelineConnections=CreateVisitDetectorFakesConnections
303):
304 """Config for CreateVisitDetectorFakesTask."""
305 randomFakeDensity = pexConfig.RangeField(
306 doc="Goal density of visit detector fake sources per square degree.",
307 dtype=float,
308 default=1000,
309 min=1,
310 )
311 nRandomFakes = pexConfig.RangeField(
312 doc="Number of random fakes to add to the visit detector. Overrides "
313 "the randomFakeDensity if set to a positive value.",
314 dtype=int,
315 default=-1,
316 min=-1,
317 )
318 doAddRandomVisitFakes = pexConfig.Field(
319 doc="Whether to add random positive fakes to the visit detector.",
320 dtype=bool,
321 default=True,
322 )
323 doAddRandomTemplateFakes = pexConfig.Field(
324 doc="Whether to add random template fakes to the visit detector (negatives).",
325 dtype=bool,
326 default=True,
327 )
328 templateFakeFraction = pexConfig.RangeField(
329 doc="Fraction of random fakes that should be added to the template image."
330 " The rest will be added to the visit image.",
331 dtype=float,
332 default=0.25,
333 min=0,
334 max=0.5,
335 inclusiveMax=True
336 )
337 doAddVariableFakes = pexConfig.Field(
338 doc="Whether to add variable fakes to the visit detector.",
339 dtype=bool,
340 default=False,
341 )
342 variableFakeFraction = pexConfig.RangeField(
343 doc="Fraction of variable fakes that should be added to the template image."
344 " The rest will be added to the visit image.",
345 dtype=float,
346 default=0.1,
347 min=0,
348 max=0.25,
349 inclusiveMax=True
350 )
351 variableFakeMean = pexConfig.RangeField(
352 doc="Mean magnitude variation for variable fakes.",
353 dtype=float,
354 default=0.0,
355 min=-1,
356 max=1,
357 )
358 variableFakeStd = pexConfig.RangeField(
359 doc="Standard deviation of magnitude variation for variable fakes.",
360 dtype=float,
361 default=0.1,
362 min=0,
363 max=1,
364 )
365 doAddHostedFakes = pexConfig.Field(
366 doc="Whether to add hosted fakes to the visit detector.",
367 dtype=bool,
368 default=False,
369 )
370 fracHostedFakes = pexConfig.RangeField(
371 doc="Fraction of hosts with fakes to add to the visit detector.",
372 dtype=float,
373 default=0.1,
374 min=0,
375 max=1,
376 )
377 magMin = pexConfig.Field(
378 doc="Minimum magnitude for the fake sources.",
379 dtype=float,
380 default=20,
381 )
382 magMax = pexConfig.Field(
383 doc="Maximum magnitude for the fake sources.",
384 dtype=float,
385 default=26,
386 )
387 doAddBlendedFakes = pexConfig.Field(
388 doc="Whether to add blended fakes to the visit detector.",
389 dtype=bool,
390 default=False,
391 )
392 fracBlendedFakes = pexConfig.RangeField(
393 doc="Fraction of blended fakes to add to the visit detector.",
394 dtype=float,
395 default=0.1,
396 min=0,
397 max=0.25,
398 inclusiveMax=True
399 )
400 fracHostedBlendedFakes = pexConfig.RangeField(
401 doc="Fraction of blended fakes that are hosted by stars.",
402 dtype=float,
403 default=0.1,
404 min=0,
405 max=0.5,
406 inclusiveMax=True
407 )
408 blendedFakeMagOffset = pexConfig.Field(
409 doc="Standard deviation of magnitude offset for blended fakes.",
410 dtype=float,
411 default=0.5,
412 )
413 blendedFakeMaxOffset = pexConfig.Field(
414 doc="Maximum positional offset for blended fakes in arcseconds.",
415 dtype=float,
416 default=3.0,
417 )
418 blendedFakeMinOffset = pexConfig.Field(
419 doc="Minimum positional offset for blended fakes in arcseconds.",
420 dtype=float,
421 default=0.2,
422 )
423 maxHostedBlendedFakesPerHost = pexConfig.RangeField(
424 doc="Maximum number of hosted blended fakes assigned to the same host in one detector.",
425 dtype=int,
426 default=2,
427 min=1,
428 )
429 maxHostedBlendedFakesTotal = pexConfig.RangeField(
430 doc="Hard cap on number of hosted blended fakes per detector. Set to -1 to disable.",
431 dtype=int,
432 default=-1,
433 min=-1,
434 )
437@dataclass
438class FakeGenerationPlan:
439 n_total: int = 0
440 n_science: int = 0
441 n_template: int = 0
442 n_variable: int = 0
443 n_blended: int = 0
445 n_hosted: int = 0
446 n_hosted_science: int = 0
447 n_hosted_template: int = 0
448 n_hosted_variable: int = 0
449 n_hosted_blended: int = 0
452class CreateVisitDetectorFakesTask(PipelineTask):
453 """Create visit-detector fakes in explicit planning and assembly stages."""
455 _DefaultName = "createVisitDetectorFakes"
456 ConfigClass = CreateVisitDetectorFakesConfig
458 def __init__(self, **kwargs):
459 super().__init__(**kwargs)
460 self.log = logging.getLogger(__name__)
461 self._table_dtypes = [
462 ("x", "<f8"),
463 ("y", "<f8"),
464 ("mag", "<f8"),
465 ("ra", "<f8"),
466 ("dec", "<f8"),
467 ("source_type", "<U4"),
468 ("isVisitSource", "?"),
469 ("isTemplateSource", "?"),
470 ("host_id", "<i8"),
471 ("host_flux", "<f8"),
472 ("host_mag", "<f8"),
473 ("host_ra", "<f8"),
474 ("host_dec", "<f8"),
475 ("delta_ra", "<f8"),
476 ("delta_dec", "<f8"),
477 ("delta_mag", "<f8"), # referred to host_mag
478 ("host_type", "<U8"),
479 ("host_a", "<f8"),
480 ("host_b", "<f8"),
481 ("host_pa", "<f8"),
482 ("hosted_fake", "?"),
483 ("injection_id", "<i8"),
484 ("isVariable", "?"),
485 ("mag_offset", "<f8"), # referred to variable parent or blend parent
486 ("twin_id", "<i8"),
487 ("isBlended", "?"),
488 ("blend_twin_id", "<i8"),
489 ("visit", "<i8"),
490 ("detector", "<i8"),
491 ]
493 def runQuantum(self, butlerQC, inputRefs, outputRefs):
494 outputs = self.run(**butlerQC.get(inputRefs))
495 butlerQC.put(outputs, outputRefs)
497 def run(self, sourceCat, visit_image):
498 """Generate, assemble, validate, and normalize a fake catalog."""
499 visit_id = visit_image.getInfo().getVisitInfo().id
500 detector_id = visit_image.detector.getId()
501 rng = np.random.default_rng([visit_id, detector_id])
502 wcs = visit_image.getWcs()
503 photo_calib = visit_image.getPhotoCalib()
504 sourceCat = photo_calib.calibrateCatalog(sourceCat).asAstropy()
505 stats = visit_image.getInfo().getSummaryStats()
506 max_mag = min(stats.magLim + 1, self.config.magMax)
508 plan = self._make_full_plan(visit_image)
509 if plan.n_total <= 0: 509 ↛ 510line 509 didn't jump to line 510 because the condition on line 509 was never true
510 raise RuntimeError(
511 "No fake sources will be generated."
512 )
513 # Build populations one by one
514 # Bear in mind that fractions describe primary sources/pairs, not catalog rows.
515 populations = []
517 # Random science fakes
518 if plan.n_science > 0 and self.config.doAddRandomVisitFakes:
519 random_science_cat = self._generate_seed_fake_catalog(plan.n_science, visit_image, wcs, rng)
520 random_science_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=plan.n_science)
521 random_science_cat["isVisitSource"] = np.ones(plan.n_science, dtype=bool)
522 random_science_cat["isTemplateSource"] = np.zeros(plan.n_science, dtype=bool)
523 random_science_cat["injection_id"] = self._make_unique_injection_ids(
524 len(random_science_cat), used_ids=None)
525 populations.append(random_science_cat)
527 # Hosted science
528 hosted_science_cat, science_hosts = self._generate_hosted_fake_catalog(
529 plan.n_hosted_science, sourceCat, wcs, rng, host_type='extended')
530 if science_hosts is None:
531 self.log.warning(
532 "Hosted science fake generation requested, but no valid extended hosts were selected."
533 )
534 else:
535 hosted_science_cat["isVisitSource"] = np.ones(plan.n_hosted_science, dtype=bool)
536 hosted_science_cat["isTemplateSource"] = np.zeros(plan.n_hosted_science, dtype=bool)
537 hosted_science_cat["hosted_fake"] = np.ones(plan.n_hosted_science, dtype=bool)
538 hosted_science_cat["injection_id"] = self._make_unique_injection_ids(
539 len(hosted_science_cat), used_ids=[cat["injection_id"] for cat in populations])
540 populations.append(hosted_science_cat)
542 # Template fakes
543 if plan.n_template > 0 and self.config.doAddRandomTemplateFakes:
544 random_template_cat = self._generate_seed_fake_catalog(plan.n_template, visit_image, wcs, rng)
545 random_template_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=plan.n_template)
546 random_template_cat["isVisitSource"] = np.zeros(plan.n_template, dtype=bool)
547 random_template_cat["isTemplateSource"] = np.ones(plan.n_template, dtype=bool)
548 random_template_cat["injection_id"] = self._make_unique_injection_ids(
549 len(random_template_cat), used_ids=[cat["injection_id"] for cat in populations])
550 populations.append(random_template_cat)
552 # Hosted template fakes
553 hosted_template_cat, template_hosts = self._generate_hosted_fake_catalog(
554 plan.n_hosted_template, sourceCat, wcs, rng, host_type='extended')
555 if template_hosts is None: 555 ↛ 560line 555 didn't jump to line 560 because the condition on line 555 was always true
556 self.log.warning(
557 "Hosted template fake generation requested, but no valid extended hosts were selected."
558 )
559 else:
560 hosted_template_cat["isVisitSource"] = np.zeros(plan.n_hosted_template, dtype=bool)
561 hosted_template_cat["isTemplateSource"] = np.ones(plan.n_hosted_template, dtype=bool)
562 hosted_template_cat["hosted_fake"] = np.ones(plan.n_hosted_template, dtype=bool)
563 hosted_template_cat["injection_id"] = self._make_unique_injection_ids(
564 len(hosted_template_cat), used_ids=[cat["injection_id"] for cat in populations])
565 populations.append(hosted_template_cat)
567 # Variable fakes
568 if plan.n_variable > 0 and self.config.doAddVariableFakes:
569 random_variable_cat = self._generate_seed_fake_catalog(plan.n_variable, visit_image, wcs, rng)
570 random_variable_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=plan.n_variable)
571 random_variable_cat["isVisitSource"] = np.ones(plan.n_variable, dtype=bool)
572 random_variable_cat["isTemplateSource"] = np.zeros(plan.n_variable, dtype=bool)
573 random_variable_cat["isVariable"] = np.ones(plan.n_variable, dtype=bool)
574 random_variable_cat["injection_id"] = self._make_unique_injection_ids(
575 len(random_variable_cat), used_ids=[cat["injection_id"] for cat in populations])
576 # make the twins for variable sources
577 twins = self._generate_variable_twins(random_variable_cat, rng)
578 twins["injection_id"] = self._make_unique_injection_ids(
579 len(twins), used_ids=[cat["injection_id"] for cat in populations + [random_variable_cat]]
580 )
581 # reciprocate the injection_id for the twins so they can be linked
582 random_variable_cat["twin_id"] = twins["injection_id"]
583 populations.append(random_variable_cat)
584 populations.append(twins)
586 # Hosted variable fakes
587 hosted_variable_cat, variable_hosts = self._generate_hosted_fake_catalog(
588 plan.n_hosted_variable, sourceCat, wcs, rng, host_type='extended')
589 if variable_hosts is None: 589 ↛ 594line 589 didn't jump to line 594 because the condition on line 589 was always true
590 self.log.warning(
591 "Hosted variable fake generation requested, but no valid extended hosts were selected."
592 )
593 else:
594 hosted_variable_cat["isVisitSource"] = np.ones(plan.n_hosted_variable, dtype=bool)
595 hosted_variable_cat["isTemplateSource"] = np.zeros(plan.n_hosted_variable, dtype=bool)
596 hosted_variable_cat["isVariable"] = np.ones(plan.n_hosted_variable, dtype=bool)
597 hosted_variable_cat["hosted_fake"] = np.ones(plan.n_hosted_variable, dtype=bool)
598 hosted_variable_cat["injection_id"] = self._make_unique_injection_ids(
599 len(hosted_variable_cat), used_ids=[cat["injection_id"] for cat in populations])
600 # make the twins for hosted variable sources
601 hosted_twins = self._generate_variable_twins(hosted_variable_cat, rng)
602 hosted_twins["injection_id"] = self._make_unique_injection_ids(
603 len(hosted_twins),
604 used_ids=[cat["injection_id"] for cat in populations + [hosted_variable_cat]]
605 )
606 hosted_variable_cat["twin_id"] = hosted_twins["injection_id"]
607 populations.append(hosted_variable_cat)
608 populations.append(hosted_twins)
610 # Blended fakes
611 if plan.n_blended > 0 and self.config.doAddBlendedFakes:
612 n_sources = plan.n_blended // 2
613 blended_cat = self._generate_seed_fake_catalog(n_sources, visit_image, wcs, rng)
614 blended_cat["injection_id"] = self._make_unique_injection_ids(
615 len(blended_cat), used_ids=[cat["injection_id"] for cat in populations])
616 blended_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=n_sources)
617 blended_cat["isVisitSource"] = rng.choice([True, False], size=n_sources)
618 blended_cat["isTemplateSource"] = ~blended_cat["isVisitSource"]
619 blended_cat["isBlended"] = np.ones(n_sources, dtype=bool)
620 # make the twins for blends
621 blended_twins = self._generate_blended_twins(blended_cat, rng, wcs)
622 blended_twins["isVisitSource"] = rng.choice([True, False], size=n_sources)
623 blended_twins["isTemplateSource"] = ~blended_twins["isVisitSource"]
625 blended_twins["injection_id"] = self._make_unique_injection_ids(
626 len(blended_twins), used_ids=[cat["injection_id"] for cat in populations + [blended_cat]]
627 )
628 # reciprocate the injection_id for the twins so they can be linked
629 blended_cat["blend_twin_id"] = blended_twins["injection_id"]
630 populations.append(blended_cat)
631 populations.append(blended_twins)
633 # Hosted blended fakes are in stars so we don't add twins.
634 n_sources = plan.n_hosted_blended
635 hosted_blended_cat, blend_hosts = self._generate_hosted_fake_catalog(
636 n_sources, sourceCat, wcs, rng, host_type='star')
637 if blend_hosts is None: 637 ↛ 638line 637 didn't jump to line 638 because the condition on line 637 was never true
638 self.log.warning(
639 "Hosted blended fake generation requested, but no valid star hosts were selected."
640 )
641 else:
642 hosted_blended_cat, blend_hosts = self._cap_hosted_blended_count(
643 hosted_blended_cat, blend_hosts
644 )
645 n_sources = len(hosted_blended_cat)
646 hosted_blended_cat["injection_id"] = self._make_unique_injection_ids(
647 len(hosted_blended_cat), used_ids=[cat["injection_id"] for cat in populations])
648 hosted_blended_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=n_sources)
649 hosted_blended_cat["isVisitSource"] = rng.choice([True, False], size=n_sources)
650 hosted_blended_cat["isTemplateSource"] = ~hosted_blended_cat["isVisitSource"]
651 hosted_blended_cat["isBlended"] = np.ones(n_sources, dtype=bool)
652 hosted_blended_cat["hosted_fake"] = np.ones(n_sources, dtype=bool)
653 populations.append(hosted_blended_cat)
655 # Merge all populations into a single catalog
656 if len(populations) == 0:
657 raise RuntimeError("No fakes were generated for this visit detector.")
659 populations = [self._set_population_defaults(population) for population in populations]
660 catalog = vstack(populations) if populations else Table(dtype=self._table_dtypes)
662 self._validate_catalog(catalog)
663 self._normalize_catalog(catalog, visit_id, detector_id)
664 return Struct(outputCat=catalog)
666 def _make_full_plan(self, visit_image):
667 """Plan exactly the number of fakes to inject, and their types,
668 based on the configuration and the visit image.
669 This is prior to knowing available hosts for hosted fakes, so the plan
670 may be adjusted later."""
672 n_total = 0
673 plan = FakeGenerationPlan()
674 if self.config.nRandomFakes > 0:
675 n_total = self.config.nRandomFakes
676 else:
677 n_total = self.get_n_fakes_from_density(
678 visit_image, self.config.randomFakeDensity
679 )
681 plan.n_total = n_total
683 if self.config.doAddRandomTemplateFakes:
684 plan.n_template = int(n_total * self.config.templateFakeFraction)
685 if self.config.doAddVariableFakes:
686 plan.n_variable = int(n_total * self.config.variableFakeFraction)
687 if self.config.doAddBlendedFakes:
688 plan.n_blended = int(n_total * self.config.fracBlendedFakes)
689 # we will generate n_blends pairs of blended sources,
690 # so the total number of blended sources is 2*n_blends
691 n_blends = plan.n_blended//2
692 plan.n_blended = n_blends * 2
693 if self.config.doAddRandomVisitFakes:
694 n_science = n_total - plan.n_template - plan.n_variable - plan.n_blended
695 if n_science < 0: 695 ↛ 696line 695 didn't jump to line 696 because the condition on line 695 was never true
696 self.log.warning(
697 "Requested number of fakes exceeds the total planned. Adjusting n_science to 0."
698 )
699 n_science = 0
700 plan.n_science = n_science
702 if self.config.doAddHostedFakes:
703 plan.n_hosted_science = int(plan.n_science * self.config.fracHostedFakes)
704 plan.n_science = plan.n_science - plan.n_hosted_science
705 plan.n_hosted_template = int(plan.n_template * self.config.fracHostedFakes)
706 plan.n_template = plan.n_template - plan.n_hosted_template
707 plan.n_hosted_variable = int(plan.n_variable * self.config.fracHostedFakes)
708 plan.n_variable = plan.n_variable - plan.n_hosted_variable
710 if self.config.fracHostedBlendedFakes > 0: 710 ↛ 714line 710 didn't jump to line 714 because the condition on line 710 was always true
711 plan.n_hosted_blended = int(plan.n_blended * self.config.fracHostedBlendedFakes)
712 plan.n_blended = plan.n_blended - plan.n_hosted_blended
714 return plan
716 def _set_population_defaults(self, catalog):
717 """Ensure population flags and relationship IDs are concrete values."""
718 defaults = {
719 "isVariable": False,
720 "isBlended": False,
721 "twin_id": 0,
722 "blend_twin_id": 0,
723 }
724 for name, default in defaults.items():
725 if name not in catalog.colnames: 725 ↛ 726line 725 didn't jump to line 726 because the condition on line 725 was never true
726 catalog[name] = np.full(len(catalog), default)
727 else:
728 values = np.ma.asarray(catalog[name])
729 catalog[name] = np.ma.filled(values, default)
730 return catalog
732 def _generate_seed_fake_catalog(self, n_fakes, visit_image, wcs, rng):
733 """Generate a seed catalog of random fakes with positions and magnitudes."""
734 if n_fakes <= 0: 734 ↛ 735line 734 didn't jump to line 735 because the condition on line 734 was never true
735 return Table(dtype=self._table_dtypes, masked=False)
737 bbox = visit_image.getBBox()
738 x = rng.uniform(bbox.getMinX(), bbox.getMaxX(), size=n_fakes)
739 y = rng.uniform(bbox.getMinY(), bbox.getMaxY(), size=n_fakes)
740 ra, dec = wcs.pixelToSkyArray(x, y, degrees=True)
742 zero_table = Table(dtype=self._table_dtypes, masked=False)
743 catalog = Table(masked=False)
744 catalog["x"] = x
745 catalog["y"] = y
746 catalog["ra"] = ra
747 catalog["dec"] = dec
748 catalog["source_type"] = "Star"
749 catalog = vstack([zero_table, catalog])
751 return catalog
753 def _build_hosted_population(self, hosts, wcs, rng, host_type='extended'):
754 """Build a population of hosted fakes around selected hosts."""
755 if host_type == 'extended':
756 pa, a, b = self.get_PA_and_axes(
757 hosts["slot_Shape_xx"], hosts["slot_Shape_xy"], hosts["slot_Shape_yy"]
758 )
759 flux_type = "slot_ModelFlux_flux"
760 mag_type = "slot_ModelFlux_mag"
761 else:
762 pa, a, b = 0, 1, 1 # For stars, we can treat them as circular with no orientation
763 flux_type = "slot_PsfFlux_flux"
764 mag_type = "slot_PsfFlux_mag"
766 n_fakes = len(hosts)
767 theta = rng.uniform(0, 2 * np.pi, size=n_fakes)
768 radii = np.sqrt((a * np.cos(theta)) ** 2 + (b * np.sin(theta)) ** 2)
769 radii *= np.sqrt(rng.uniform(0, 6, size=n_fakes))
770 x = hosts["slot_Centroid_x"] + radii * np.cos(theta) * np.cos(pa) - radii * np.sin(theta) * np.sin(pa)
771 y = hosts["slot_Centroid_y"] + radii * np.cos(theta) * np.sin(pa) + radii * np.sin(theta) * np.cos(pa)
772 ra, dec = wcs.pixelToSkyArray(x, y, degrees=True)
774 zero_table = Table(dtype=self._table_dtypes, masked=False)
775 catalog = Table(masked=False)
776 catalog["x"] = x
777 catalog["y"] = y
778 # magnitudes are drawn from a normal distribution around host magnitude
779 catalog["mag"] = hosts[mag_type] + rng.normal(1, 1, size=n_fakes)
780 catalog["ra"] = ra
781 catalog["dec"] = dec
782 catalog["host_id"] = hosts["id"]
783 catalog["host_flux"] = hosts[flux_type]
784 catalog["host_mag"] = hosts[mag_type]
785 catalog["host_ra"] = np.rad2deg(hosts["coord_ra"])
786 catalog["host_dec"] = np.rad2deg(hosts["coord_dec"])
787 catalog["delta_ra"] = (ra - catalog["host_ra"]) * 3600
788 catalog["delta_dec"] = (dec - catalog["host_dec"]) * 3600
789 catalog["delta_mag"] = catalog["mag"] - catalog["host_mag"]
790 catalog["host_a"] = a
791 catalog["host_b"] = b
792 catalog["host_pa"] = pa
793 catalog["source_type"] = "Star"
794 catalog["host_type"] = host_type
795 catalog["hosted_fake"] = np.ones(len(catalog), dtype=bool)
796 catalog = vstack([zero_table, catalog])
798 return catalog, hosts
800 def _generate_hosted_fake_catalog(self, n_fakes, sourceCat, wcs, rng, host_type='extended'):
801 """Generate a catalog of hosted fakes around selected host galaxies or stars."""
802 if n_fakes <= 0:
803 return Table(dtype=self._table_dtypes, masked=False), None
805 host_catalog = self.select_hosts(sourceCat, host_type=host_type)
806 if len(host_catalog) == 0:
807 self.log.warning("No valid hosts were selected for hosted fakes.")
808 return Table(dtype=self._table_dtypes, masked=False), None
810 # we allow more than one fake per host since we can deal with blending now
811 if n_fakes > len(host_catalog): 811 ↛ 820line 811 didn't jump to line 820 because the condition on line 811 was always true
812 self.log.warning(
813 "Requested %d hosted fakes, but only %d valid hosts were found. "
814 "Some hosts will have multiple fakes.",
815 n_fakes,
816 len(host_catalog),
817 )
818 idx = rng.choice(len(host_catalog), size=n_fakes, replace=True)
819 else:
820 idx = rng.choice(len(host_catalog), size=n_fakes, replace=False)
822 selected_hosts = host_catalog[idx]
823 return self._build_hosted_population(selected_hosts, wcs, rng, host_type=host_type)
825 def _generate_variable_twins(self, variable_catalog, rng):
826 """Generate twin sources for variable fakes with magnitude offsets."""
827 n_twins = len(variable_catalog)
828 if n_twins <= 0: 828 ↛ 829line 828 didn't jump to line 829 because the condition on line 828 was never true
829 return Table(dtype=self._table_dtypes, masked=False)
830 # magnitudes are drawn from a normal distribution around the variable source's magnitude
831 twins = variable_catalog.copy()
832 twins["mag_offset"] = rng.normal(
833 loc=self.config.variableFakeMean,
834 scale=self.config.variableFakeStd,
835 size=n_twins
836 )
837 twins["mag"] += twins["mag_offset"]
838 twins["isVisitSource"] = ~twins["isVisitSource"]
839 twins["isTemplateSource"] = ~twins["isTemplateSource"]
840 twins["twin_id"] = twins["injection_id"]
841 return twins
843 def _generate_blended_twins(self, blended_catalog, rng, wcs):
844 """Generate twin sources for blended fakes with magnitude and positional offsets."""
845 n_twins = len(blended_catalog)
846 if n_twins <= 0: 846 ↛ 847line 846 didn't jump to line 847 because the condition on line 846 was never true
847 return Table(dtype=self._table_dtypes, masked=False)
849 twins = blended_catalog.copy()
850 twins["mag_offset"] = rng.normal(
851 loc=0.0,
852 scale=self.config.blendedFakeMagOffset,
853 size=n_twins
854 )
855 delta_ra, delta_dec = self._draw_offset_components_arcsec(rng, n_twins)
856 twins["delta_ra"] = delta_ra
857 twins["delta_dec"] = delta_dec
858 twins["mag"] += twins["mag_offset"]
860 twins["ra"] += (twins["delta_ra"] / 3600.0) / np.cos(np.deg2rad(twins["dec"]))
861 twins["dec"] += twins["delta_dec"] / 3600.0
862 twins["x"], twins["y"] = wcs.skyToPixelArray(twins["ra"], twins["dec"], degrees=True)
864 twins["blend_twin_id"] = twins["injection_id"]
865 return twins
867 def _validate_catalog(self, catalog):
868 """Validate IDs and reciprocal source relationships."""
869 injection_ids = np.asarray(catalog["injection_id"], dtype=np.int64)
870 if len(np.unique(injection_ids)) != len(injection_ids): 870 ↛ 871line 870 didn't jump to line 871 because the condition on line 870 was never true
871 raise RuntimeError("Fake injection IDs are not unique.")
873 id_set = set(injection_ids)
875 variable = np.asarray(catalog["isVariable"], dtype=bool)
876 blended = np.asarray(catalog["isBlended"], dtype=bool)
878 if np.any(variable & blended): 878 ↛ 879line 878 didn't jump to line 879 because the condition on line 878 was never true
879 raise RuntimeError("Variable and blended fake populations overlap.")
881 twin_ids = np.asarray(catalog["twin_id"], dtype=np.int64)
882 variable_twin_ids = twin_ids[variable]
884 if np.any(variable_twin_ids <= 0): 884 ↛ 885line 884 didn't jump to line 885 because the condition on line 884 was never true
885 raise RuntimeError("Variable fake twin links are incomplete.")
887 if not set(variable_twin_ids).issubset(id_set): 887 ↛ 888line 887 didn't jump to line 888 because the condition on line 887 was never true
888 raise RuntimeError("Variable twin IDs do not refer to catalog rows.")
890 variable_rows_by_id = {
891 int(injection_id): index
892 for index, injection_id in enumerate(injection_ids)
893 if variable[index]
894 }
896 for index in np.flatnonzero(variable):
897 twin_index = variable_rows_by_id.get(int(twin_ids[index]))
898 if twin_index is None: 898 ↛ 899line 898 didn't jump to line 899 because the condition on line 898 was never true
899 raise RuntimeError("Variable twin ID does not identify a variable row.")
900 if twin_ids[twin_index] != injection_ids[index]: 900 ↛ 901line 900 didn't jump to line 901 because the condition on line 900 was never true
901 raise RuntimeError("Variable twin links are not reciprocal.")
903 # removing hosted blends from the check
904 non_hosted_blends = blended & ~np.asarray(catalog["hosted_fake"], dtype=bool)
905 blend_twin_ids = np.asarray(catalog["blend_twin_id"], dtype=np.int64)
906 blended_twin_ids = blend_twin_ids[non_hosted_blends]
908 if np.any(blended_twin_ids <= 0): 908 ↛ 909line 908 didn't jump to line 909 because the condition on line 908 was never true
909 raise RuntimeError("Blended fake twin links are incomplete.")
911 if not set(blended_twin_ids).issubset(id_set): 911 ↛ 912line 911 didn't jump to line 912 because the condition on line 911 was never true
912 raise RuntimeError("Blend twin IDs do not refer to catalog rows.")
914 rows_by_id = {
915 int(injection_id): index
916 for index, injection_id in enumerate(injection_ids)
917 }
919 for index in np.flatnonzero(blend_twin_ids):
920 twin_index = rows_by_id[int(blend_twin_ids[index])]
921 if blend_twin_ids[twin_index] != injection_ids[index]: 921 ↛ 922line 921 didn't jump to line 922 because the condition on line 921 was never true
922 raise RuntimeError("Blend twin links are not reciprocal.")
924 def _normalize_catalog(self, catalog, visit_id, detector_id):
925 """Normalize catalog columns to remove units and add visit/detector identifiers."""
926 catalog["visit"] = visit_id
927 catalog["detector"] = detector_id
928 for name in ("ra", "dec", "delta_ra", "delta_dec", "host_ra", "host_dec"):
929 if name in catalog.colnames and catalog[name].unit is not None: 929 ↛ 930line 929 didn't jump to line 930 because the condition on line 929 was never true
930 catalog[name] = catalog[name].value
932 def select_hosts(self, sourceCat, host_type='extended'):
933 """Select hosts from the source catalog."""
934 sky = ~sourceCat["sky_source"]
935 flags = ~sourceCat["base_ClassificationSizeExtendedness_flag"]
936 flags &= ~sourceCat["base_ClassificationExtendedness_flag"]
937 flags &= ~sourceCat["slot_Shape_flag"]
938 flags &= ~sourceCat["slot_Centroid_flag"]
939 flags &= ~sourceCat["base_PixelFlags_flag"]
941 if host_type == 'extended':
942 extended = sourceCat["base_ClassificationSizeExtendedness_value"] > 0.9
943 extended &= sourceCat["base_ClassificationExtendedness_value"] == 1
944 snr = sourceCat["slot_ModelFlux_flux"] / sourceCat["slot_ModelFlux_fluxErr"] > 15
945 return sourceCat[sky & flags & extended & snr].copy()
946 elif host_type == 'star': 946 ↛ 952line 946 didn't jump to line 952 because the condition on line 946 was always true
947 stellar = sourceCat["base_ClassificationSizeExtendedness_value"] < 0.9
948 stellar &= sourceCat["base_ClassificationExtendedness_value"] == 0
949 snr = sourceCat["slot_PsfFlux_flux"] / sourceCat["slot_PsfFlux_fluxErr"] > 30
950 return sourceCat[sky & flags & stellar & snr].copy()
951 else:
952 raise ValueError("host_type must be 'extended' or 'star'.")
954 def get_PA_and_axes(self, Ixx, Ixy, Iyy):
955 """Compute the position angle and semi-major/minor axes from second moments."""
956 theta = 0.5 * np.arctan2(2 * Ixy, Ixx - Iyy)
957 term = np.sqrt(((Ixx - Iyy) / 2) ** 2 + Ixy ** 2)
958 return theta, np.sqrt((Ixx + Iyy) / 2 + term), np.sqrt((Ixx + Iyy) / 2 - term)
960 def _cap_hosted_blended_count(self, hosted_fakes, hosts):
961 """Cap the number of hosted blended fakes based on configuration limits."""
962 if len(hosted_fakes) == 0 or len(hosts) == 0: 962 ↛ 963line 962 didn't jump to line 963 because the condition on line 962 was never true
963 return hosted_fakes, hosts
965 # check if we have multiple ids for the same host and cap the number of fakes per host
966 unique_hosts, counts = np.unique(hosted_fakes["host_id"], return_counts=True)
967 capped_hosts = []
968 for host_id, count in zip(unique_hosts, counts):
969 if count > self.config.maxHostedBlendedFakesPerHost:
970 self.log.warning(
971 "Capping hosted blended fakes for host_id %d from %d to %d.",
972 host_id, count, self.config.maxHostedBlendedFakesPerHost
973 )
974 capped_hosts.append(host_id)
975 row_ids = np.where(hosted_fakes["host_id"] == host_id)
976 row_ids_to_remove = row_ids[0][self.config.maxHostedBlendedFakesPerHost:]
977 hosted_fakes.remove_rows(row_ids_to_remove)
978 if self.config.maxHostedBlendedFakesTotal > 0:
979 if len(hosted_fakes) > self.config.maxHostedBlendedFakesTotal: 979 ↛ 985line 979 didn't jump to line 985 because the condition on line 979 was always true
980 self.log.warning(
981 "Capping total hosted blended fakes from %d to %d.",
982 len(hosted_fakes), self.config.maxHostedBlendedFakesTotal
983 )
984 hosted_fakes = hosted_fakes[:self.config.maxHostedBlendedFakesTotal]
985 return hosted_fakes, hosts
987 def _draw_offset_components_arcsec(self, rng, n_points):
988 """Draw random offsets in arcseconds for blended fakes.
989 The offsets are drawn from a uniform distribution in radius and angle."""
990 if n_points <= 0: 990 ↛ 991line 990 didn't jump to line 991 because the condition on line 990 was never true
991 return np.zeros(0), np.zeros(0)
992 rmin = max(0.0, self.config.blendedFakeMinOffset)
993 rmax = max(rmin, self.config.blendedFakeMaxOffset)
994 radius = np.sqrt(rng.uniform(rmin**2, rmax**2, n_points))
995 theta = rng.uniform(0, 2 * np.pi, n_points)
996 return radius * np.cos(theta), radius * np.sin(theta)
998 def get_n_fakes_from_density(self, visit_image, density):
999 area = visit_image.getConvexPolygon().getBoundingBox().getArea()
1000 return np.round(density * area * (180 / np.pi) ** 2).astype(int)
1002 def _make_unique_injection_ids(self, n_ids, used_ids=None):
1003 """Generate unique 24-bit IDs for a catalog considering IDs already in use."""
1004 new_ids = set()
1005 if used_ids:
1006 # compile used ids from the existing catalog to avoid collisions
1007 if isinstance(used_ids, list): 1007 ↛ 1013line 1007 didn't jump to line 1013 because the condition on line 1007 was always true
1008 used_ids = np.concatenate(used_ids)
1009 used_ids = set(used_ids)
1010 else:
1011 used_ids = set()
1012 # check that the new ids are not in the used set already
1013 while len(new_ids) < n_ids:
1014 candidate_id = uuid.uuid4().int & ((1 << 24) - 1)
1015 if candidate_id > 0 and candidate_id not in used_ids.union(new_ids): 1015 ↛ 1013line 1015 didn't jump to line 1013 because the condition on line 1015 was always true
1016 new_ids.add(candidate_id)
1018 return np.asarray(list(new_ids), dtype=np.int64)