Coverage for python/lsst/ap/pipe/createApFakes.py: 88%
455 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-16 09:39 +0000
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 defaultTemplates={"coaddName": "deep"},
277 dimensions=("instrument",
278 "visit",
279 "detector")):
281 sourceCat = connTypes.Input(
282 doc="Catalog of sources detected on the calibrated exposure; ",
283 name="single_visit_star_reprocessed_footprints",
284 storageClass="SourceCatalog",
285 dimensions=["instrument", "visit", "detector"],
286 )
287 visit_image = connTypes.Input(
288 doc="Calibrated exposure to inject synthetic sources into.",
289 name="preliminary_visit_image",
290 storageClass="ExposureF",
291 dimensions=["instrument", "visit", "detector"],
292 )
293 outputCat = connTypes.Output(
294 doc="Catalog of fake sources to draw inputs from.",
295 name="VisitDetectorFakeSourceCat",
296 storageClass="ArrowAstropy",
297 dimensions=["instrument", "visit", "detector"],
298 )
301class CreateVisitDetectorFakesConfig(
302 PipelineTaskConfig,
303 pipelineConnections=CreateVisitDetectorFakesConnections
304):
305 """Config for CreateVisitDetectorFakesTask."""
306 randomFakeDensity = pexConfig.RangeField(
307 doc="Goal density of visit detector fake sources per square degree.",
308 dtype=float,
309 default=1000,
310 min=1,
311 )
312 nRandomFakes = pexConfig.RangeField(
313 doc="Number of random fakes to add to the visit detector. Overrides "
314 "the randomFakeDensity if set to a positive value.",
315 dtype=int,
316 default=-1,
317 min=-1,
318 )
319 doAddRandomVisitFakes = pexConfig.Field(
320 doc="Whether to add random positive fakes to the visit detector.",
321 dtype=bool,
322 default=True,
323 )
324 doAddRandomTemplateFakes = pexConfig.Field(
325 doc="Whether to add random template fakes to the visit detector (negatives).",
326 dtype=bool,
327 default=True,
328 )
329 templateFakeFraction = pexConfig.RangeField(
330 doc="Fraction of random fakes that should be added to the template image."
331 " The rest will be added to the visit image.",
332 dtype=float,
333 default=0.25,
334 min=0,
335 max=0.5,
336 inclusiveMax=True
337 )
338 doAddVariableFakes = pexConfig.Field(
339 doc="Whether to add variable fakes to the visit detector.",
340 dtype=bool,
341 default=False,
342 )
343 variableFakeFraction = pexConfig.RangeField(
344 doc="Fraction of variable fakes that should be added to the template image."
345 " The rest will be added to the visit image.",
346 dtype=float,
347 default=0.1,
348 min=0,
349 max=0.25,
350 inclusiveMax=True
351 )
352 variableFakeMean = pexConfig.RangeField(
353 doc="Mean magnitude variation for variable fakes.",
354 dtype=float,
355 default=0.0,
356 min=-1,
357 max=1,
358 )
359 variableFakeStd = pexConfig.RangeField(
360 doc="Standard deviation of magnitude variation for variable fakes.",
361 dtype=float,
362 default=0.1,
363 min=0,
364 max=1,
365 )
366 doAddHostedFakes = pexConfig.Field(
367 doc="Whether to add hosted fakes to the visit detector.",
368 dtype=bool,
369 default=False,
370 )
371 fracHostedFakes = pexConfig.RangeField(
372 doc="Fraction of hosts with fakes to add to the visit detector.",
373 dtype=float,
374 default=0.1,
375 min=0,
376 max=1,
377 )
378 magMin = pexConfig.Field(
379 doc="Minimum magnitude for the fake sources.",
380 dtype=float,
381 default=20,
382 )
383 magMax = pexConfig.Field(
384 doc="Maximum magnitude for the fake sources.",
385 dtype=float,
386 default=26,
387 )
388 doAddBlendedFakes = pexConfig.Field(
389 doc="Whether to add blended fakes to the visit detector.",
390 dtype=bool,
391 default=False,
392 )
393 fracBlendedFakes = pexConfig.RangeField(
394 doc="Fraction of blended fakes to add to the visit detector.",
395 dtype=float,
396 default=0.1,
397 min=0,
398 max=0.25,
399 inclusiveMax=True
400 )
401 fracHostedBlendedFakes = pexConfig.RangeField(
402 doc="Fraction of blended fakes that are hosted by stars.",
403 dtype=float,
404 default=0.1,
405 min=0,
406 max=0.5,
407 inclusiveMax=True
408 )
409 blendedFakeMagOffset = pexConfig.Field(
410 doc="Standard deviation of magnitude offset for blended fakes.",
411 dtype=float,
412 default=0.5,
413 )
414 blendedFakeMaxOffset = pexConfig.Field(
415 doc="Maximum positional offset for blended fakes in arcseconds.",
416 dtype=float,
417 default=3.0,
418 )
419 blendedFakeMinOffset = pexConfig.Field(
420 doc="Minimum positional offset for blended fakes in arcseconds.",
421 dtype=float,
422 default=0.2,
423 )
424 maxHostedBlendedFakesPerHost = pexConfig.RangeField(
425 doc="Maximum number of hosted blended fakes assigned to the same host in one detector.",
426 dtype=int,
427 default=2,
428 min=1,
429 )
430 maxHostedBlendedFakesTotal = pexConfig.RangeField(
431 doc="Hard cap on number of hosted blended fakes per detector. Set to -1 to disable.",
432 dtype=int,
433 default=-1,
434 min=-1,
435 )
438@dataclass
439class FakeGenerationPlan:
440 n_total: int = 0
441 n_science: int = 0
442 n_template: int = 0
443 n_variable: int = 0
444 n_blended: int = 0
446 n_hosted: int = 0
447 n_hosted_science: int = 0
448 n_hosted_template: int = 0
449 n_hosted_variable: int = 0
450 n_hosted_blended: int = 0
453class CreateVisitDetectorFakesTask(PipelineTask):
454 """Create visit-detector fakes in explicit planning and assembly stages."""
456 _DefaultName = "createVisitDetectorFakes"
457 ConfigClass = CreateVisitDetectorFakesConfig
459 def __init__(self, **kwargs):
460 super().__init__(**kwargs)
461 self.log = logging.getLogger(__name__)
462 self._table_dtypes = [
463 ("x", "<f8"),
464 ("y", "<f8"),
465 ("mag", "<f8"),
466 ("ra", "<f8"),
467 ("dec", "<f8"),
468 ("source_type", "<U4"),
469 ("isVisitSource", "?"),
470 ("isTemplateSource", "?"),
471 ("host_id", "<i8"),
472 ("host_flux", "<f8"),
473 ("host_mag", "<f8"),
474 ("host_ra", "<f8"),
475 ("host_dec", "<f8"),
476 ("delta_ra", "<f8"),
477 ("delta_dec", "<f8"),
478 ("delta_mag", "<f8"), # referred to host_mag
479 ("host_type", "<U8"),
480 ("host_a", "<f8"),
481 ("host_b", "<f8"),
482 ("host_pa", "<f8"),
483 ("hosted_fake", "?"),
484 ("injection_id", "<i8"),
485 ("isVariable", "?"),
486 ("mag_offset", "<f8"), # referred to variable parent or blend parent
487 ("twin_id", "<i8"),
488 ("isBlended", "?"),
489 ("blend_twin_id", "<i8"),
490 ("visit", "<i8"),
491 ("detector", "<i8"),
492 ]
494 def runQuantum(self, butlerQC, inputRefs, outputRefs):
495 outputs = self.run(**butlerQC.get(inputRefs))
496 butlerQC.put(outputs, outputRefs)
498 def run(self, sourceCat, visit_image):
499 """Generate, assemble, validate, and normalize a fake catalog."""
500 visit_id = visit_image.getInfo().getVisitInfo().id
501 detector_id = visit_image.detector.getId()
502 rng = np.random.default_rng([visit_id, detector_id])
503 wcs = visit_image.getWcs()
504 photo_calib = visit_image.getPhotoCalib()
505 sourceCat = photo_calib.calibrateCatalog(sourceCat).asAstropy()
506 stats = visit_image.getInfo().getSummaryStats()
507 max_mag = min(stats.magLim + 1, self.config.magMax)
509 plan = self._make_full_plan(visit_image)
510 if plan.n_total <= 0: 510 ↛ 511line 510 didn't jump to line 511 because the condition on line 510 was never true
511 raise RuntimeError(
512 "No fake sources will be generated."
513 )
514 # Build populations one by one
515 # Bear in mind that fractions describe primary sources/pairs, not catalog rows.
516 populations = []
518 # Random science fakes
519 if plan.n_science > 0 and self.config.doAddRandomVisitFakes:
520 random_science_cat = self._generate_seed_fake_catalog(plan.n_science, visit_image, wcs, rng)
521 random_science_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=plan.n_science)
522 random_science_cat["isVisitSource"] = np.ones(plan.n_science, dtype=bool)
523 random_science_cat["isTemplateSource"] = np.zeros(plan.n_science, dtype=bool)
524 random_science_cat["injection_id"] = self._make_unique_injection_ids(
525 len(random_science_cat), used_ids=None)
526 populations.append(random_science_cat)
528 # Hosted science
529 hosted_science_cat, science_hosts = self._generate_hosted_fake_catalog(
530 plan.n_hosted_science, sourceCat, wcs, rng, host_type='extended')
531 if science_hosts is None:
532 self.log.warning(
533 "Hosted science fake generation requested, but no valid extended hosts were selected."
534 )
535 else:
536 hosted_science_cat["isVisitSource"] = np.ones(plan.n_hosted_science, dtype=bool)
537 hosted_science_cat["isTemplateSource"] = np.zeros(plan.n_hosted_science, dtype=bool)
538 hosted_science_cat["hosted_fake"] = np.ones(plan.n_hosted_science, dtype=bool)
539 hosted_science_cat["injection_id"] = self._make_unique_injection_ids(
540 len(hosted_science_cat), used_ids=[cat["injection_id"] for cat in populations])
541 populations.append(hosted_science_cat)
543 # Template fakes
544 if plan.n_template > 0 and self.config.doAddRandomTemplateFakes:
545 random_template_cat = self._generate_seed_fake_catalog(plan.n_template, visit_image, wcs, rng)
546 random_template_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=plan.n_template)
547 random_template_cat["isVisitSource"] = np.zeros(plan.n_template, dtype=bool)
548 random_template_cat["isTemplateSource"] = np.ones(plan.n_template, dtype=bool)
549 random_template_cat["injection_id"] = self._make_unique_injection_ids(
550 len(random_template_cat), used_ids=[cat["injection_id"] for cat in populations])
551 populations.append(random_template_cat)
553 # Hosted template fakes
554 hosted_template_cat, template_hosts = self._generate_hosted_fake_catalog(
555 plan.n_hosted_template, sourceCat, wcs, rng, host_type='extended')
556 if template_hosts is None: 556 ↛ 561line 556 didn't jump to line 561 because the condition on line 556 was always true
557 self.log.warning(
558 "Hosted template fake generation requested, but no valid extended hosts were selected."
559 )
560 else:
561 hosted_template_cat["isVisitSource"] = np.zeros(plan.n_hosted_template, dtype=bool)
562 hosted_template_cat["isTemplateSource"] = np.ones(plan.n_hosted_template, dtype=bool)
563 hosted_template_cat["hosted_fake"] = np.ones(plan.n_hosted_template, dtype=bool)
564 hosted_template_cat["injection_id"] = self._make_unique_injection_ids(
565 len(hosted_template_cat), used_ids=[cat["injection_id"] for cat in populations])
566 populations.append(hosted_template_cat)
568 # Variable fakes
569 if plan.n_variable > 0 and self.config.doAddVariableFakes:
570 random_variable_cat = self._generate_seed_fake_catalog(plan.n_variable, visit_image, wcs, rng)
571 random_variable_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=plan.n_variable)
572 random_variable_cat["isVisitSource"] = np.ones(plan.n_variable, dtype=bool)
573 random_variable_cat["isTemplateSource"] = np.zeros(plan.n_variable, dtype=bool)
574 random_variable_cat["isVariable"] = np.ones(plan.n_variable, dtype=bool)
575 random_variable_cat["injection_id"] = self._make_unique_injection_ids(
576 len(random_variable_cat), used_ids=[cat["injection_id"] for cat in populations])
577 # make the twins for variable sources
578 twins = self._generate_variable_twins(random_variable_cat, rng)
579 twins["injection_id"] = self._make_unique_injection_ids(
580 len(twins), used_ids=[cat["injection_id"] for cat in populations + [random_variable_cat]]
581 )
582 # reciprocate the injection_id for the twins so they can be linked
583 random_variable_cat["twin_id"] = twins["injection_id"]
584 populations.append(random_variable_cat)
585 populations.append(twins)
587 # Hosted variable fakes
588 hosted_variable_cat, variable_hosts = self._generate_hosted_fake_catalog(
589 plan.n_hosted_variable, sourceCat, wcs, rng, host_type='extended')
590 if variable_hosts is None: 590 ↛ 595line 590 didn't jump to line 595 because the condition on line 590 was always true
591 self.log.warning(
592 "Hosted variable fake generation requested, but no valid extended hosts were selected."
593 )
594 else:
595 hosted_variable_cat["isVisitSource"] = np.ones(plan.n_hosted_variable, dtype=bool)
596 hosted_variable_cat["isTemplateSource"] = np.zeros(plan.n_hosted_variable, dtype=bool)
597 hosted_variable_cat["isVariable"] = np.ones(plan.n_hosted_variable, dtype=bool)
598 hosted_variable_cat["hosted_fake"] = np.ones(plan.n_hosted_variable, dtype=bool)
599 hosted_variable_cat["injection_id"] = self._make_unique_injection_ids(
600 len(hosted_variable_cat), used_ids=[cat["injection_id"] for cat in populations])
601 # make the twins for hosted variable sources
602 hosted_twins = self._generate_variable_twins(hosted_variable_cat, rng)
603 hosted_twins["injection_id"] = self._make_unique_injection_ids(
604 len(hosted_twins),
605 used_ids=[cat["injection_id"] for cat in populations + [hosted_variable_cat]]
606 )
607 hosted_variable_cat["twin_id"] = hosted_twins["injection_id"]
608 populations.append(hosted_variable_cat)
609 populations.append(hosted_twins)
611 # Blended fakes
612 if plan.n_blended > 0 and self.config.doAddBlendedFakes:
613 n_sources = plan.n_blended // 2
614 blended_cat = self._generate_seed_fake_catalog(n_sources, visit_image, wcs, rng)
615 blended_cat["injection_id"] = self._make_unique_injection_ids(
616 len(blended_cat), used_ids=[cat["injection_id"] for cat in populations])
617 blended_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=n_sources)
618 blended_cat["isVisitSource"] = rng.choice([True, False], size=n_sources)
619 blended_cat["isTemplateSource"] = ~blended_cat["isVisitSource"]
620 blended_cat["isBlended"] = np.ones(n_sources, dtype=bool)
621 # make the twins for blends
622 blended_twins = self._generate_blended_twins(blended_cat, rng, wcs)
623 blended_twins["isVisitSource"] = rng.choice([True, False], size=n_sources)
624 blended_twins["isTemplateSource"] = ~blended_twins["isVisitSource"]
626 blended_twins["injection_id"] = self._make_unique_injection_ids(
627 len(blended_twins), used_ids=[cat["injection_id"] for cat in populations + [blended_cat]]
628 )
629 # reciprocate the injection_id for the twins so they can be linked
630 blended_cat["blend_twin_id"] = blended_twins["injection_id"]
631 populations.append(blended_cat)
632 populations.append(blended_twins)
634 # Hosted blended fakes are in stars so we don't add twins.
635 n_sources = plan.n_hosted_blended
636 hosted_blended_cat, blend_hosts = self._generate_hosted_fake_catalog(
637 n_sources, sourceCat, wcs, rng, host_type='star')
638 if blend_hosts is None: 638 ↛ 639line 638 didn't jump to line 639 because the condition on line 638 was never true
639 self.log.warning(
640 "Hosted blended fake generation requested, but no valid star hosts were selected."
641 )
642 else:
643 hosted_blended_cat, blend_hosts = self._cap_hosted_blended_count(
644 hosted_blended_cat, blend_hosts
645 )
646 n_sources = len(hosted_blended_cat)
647 hosted_blended_cat["injection_id"] = self._make_unique_injection_ids(
648 len(hosted_blended_cat), used_ids=[cat["injection_id"] for cat in populations])
649 hosted_blended_cat["mag"] = rng.uniform(self.config.magMin, max_mag, size=n_sources)
650 hosted_blended_cat["isVisitSource"] = rng.choice([True, False], size=n_sources)
651 hosted_blended_cat["isTemplateSource"] = ~hosted_blended_cat["isVisitSource"]
652 hosted_blended_cat["isBlended"] = np.ones(n_sources, dtype=bool)
653 hosted_blended_cat["hosted_fake"] = np.ones(n_sources, dtype=bool)
654 populations.append(hosted_blended_cat)
656 # Merge all populations into a single catalog
657 if len(populations) == 0:
658 raise RuntimeError("No fakes were generated for this visit detector.")
660 populations = [self._set_population_defaults(population) for population in populations]
661 catalog = vstack(populations) if populations else Table(dtype=self._table_dtypes)
663 self._validate_catalog(catalog)
664 self._normalize_catalog(catalog, visit_id, detector_id)
665 return Struct(outputCat=catalog)
667 def _make_full_plan(self, visit_image):
668 """Plan exactly the number of fakes to inject, and their types,
669 based on the configuration and the visit image.
670 This is prior to knowing available hosts for hosted fakes, so the plan
671 may be adjusted later."""
673 n_total = 0
674 plan = FakeGenerationPlan()
675 if self.config.nRandomFakes > 0:
676 n_total = self.config.nRandomFakes
677 else:
678 n_total = self.get_n_fakes_from_density(
679 visit_image, self.config.randomFakeDensity
680 )
682 plan.n_total = n_total
684 if self.config.doAddRandomTemplateFakes:
685 plan.n_template = int(n_total * self.config.templateFakeFraction)
686 if self.config.doAddVariableFakes:
687 plan.n_variable = int(n_total * self.config.variableFakeFraction)
688 if self.config.doAddBlendedFakes:
689 plan.n_blended = int(n_total * self.config.fracBlendedFakes)
690 # we will generate n_blends pairs of blended sources,
691 # so the total number of blended sources is 2*n_blends
692 n_blends = plan.n_blended//2
693 plan.n_blended = n_blends * 2
694 if self.config.doAddRandomVisitFakes:
695 n_science = n_total - plan.n_template - plan.n_variable - plan.n_blended
696 if n_science < 0: 696 ↛ 697line 696 didn't jump to line 697 because the condition on line 696 was never true
697 self.log.warning(
698 "Requested number of fakes exceeds the total planned. Adjusting n_science to 0."
699 )
700 n_science = 0
701 plan.n_science = n_science
703 if self.config.doAddHostedFakes:
704 plan.n_hosted_science = int(plan.n_science * self.config.fracHostedFakes)
705 plan.n_science = plan.n_science - plan.n_hosted_science
706 plan.n_hosted_template = int(plan.n_template * self.config.fracHostedFakes)
707 plan.n_template = plan.n_template - plan.n_hosted_template
708 plan.n_hosted_variable = int(plan.n_variable * self.config.fracHostedFakes)
709 plan.n_variable = plan.n_variable - plan.n_hosted_variable
711 if self.config.fracHostedBlendedFakes > 0: 711 ↛ 715line 711 didn't jump to line 715 because the condition on line 711 was always true
712 plan.n_hosted_blended = int(plan.n_blended * self.config.fracHostedBlendedFakes)
713 plan.n_blended = plan.n_blended - plan.n_hosted_blended
715 return plan
717 def _set_population_defaults(self, catalog):
718 """Ensure population flags and relationship IDs are concrete values."""
719 defaults = {
720 "isVariable": False,
721 "isBlended": False,
722 "twin_id": 0,
723 "blend_twin_id": 0,
724 }
725 for name, default in defaults.items():
726 if name not in catalog.colnames: 726 ↛ 727line 726 didn't jump to line 727 because the condition on line 726 was never true
727 catalog[name] = np.full(len(catalog), default)
728 else:
729 values = np.ma.asarray(catalog[name])
730 catalog[name] = np.ma.filled(values, default)
731 return catalog
733 def _generate_seed_fake_catalog(self, n_fakes, visit_image, wcs, rng):
734 """Generate a seed catalog of random fakes with positions and magnitudes."""
735 if n_fakes <= 0: 735 ↛ 736line 735 didn't jump to line 736 because the condition on line 735 was never true
736 return Table(dtype=self._table_dtypes, masked=False)
738 bbox = visit_image.getBBox()
739 x = rng.uniform(bbox.getMinX(), bbox.getMaxX(), size=n_fakes)
740 y = rng.uniform(bbox.getMinY(), bbox.getMaxY(), size=n_fakes)
741 ra, dec = wcs.pixelToSkyArray(x, y, degrees=True)
743 zero_table = Table(dtype=self._table_dtypes, masked=False)
744 catalog = Table(masked=False)
745 catalog["x"] = x
746 catalog["y"] = y
747 catalog["ra"] = ra
748 catalog["dec"] = dec
749 catalog["source_type"] = "Star"
750 catalog = vstack([zero_table, catalog])
752 return catalog
754 def _build_hosted_population(self, hosts, wcs, rng, host_type='extended'):
755 """Build a population of hosted fakes around selected hosts."""
756 if host_type == 'extended':
757 pa, a, b = self.get_PA_and_axes(
758 hosts["slot_Shape_xx"], hosts["slot_Shape_xy"], hosts["slot_Shape_yy"]
759 )
760 flux_type = "slot_ModelFlux_flux"
761 mag_type = "slot_ModelFlux_mag"
762 else:
763 pa, a, b = 0, 1, 1 # For stars, we can treat them as circular with no orientation
764 flux_type = "slot_PsfFlux_flux"
765 mag_type = "slot_PsfFlux_mag"
767 n_fakes = len(hosts)
768 theta = rng.uniform(0, 2 * np.pi, size=n_fakes)
769 radii = np.sqrt((a * np.cos(theta)) ** 2 + (b * np.sin(theta)) ** 2)
770 radii *= np.sqrt(rng.uniform(0, 6, size=n_fakes))
771 x = hosts["slot_Centroid_x"] + radii * np.cos(theta) * np.cos(pa) - radii * np.sin(theta) * np.sin(pa)
772 y = hosts["slot_Centroid_y"] + radii * np.cos(theta) * np.sin(pa) + radii * np.sin(theta) * np.cos(pa)
773 ra, dec = wcs.pixelToSkyArray(x, y, degrees=True)
775 zero_table = Table(dtype=self._table_dtypes, masked=False)
776 catalog = Table(masked=False)
777 catalog["x"] = x
778 catalog["y"] = y
779 # magnitudes are drawn from a normal distribution around host magnitude
780 catalog["mag"] = hosts[mag_type] + rng.normal(1, 1, size=n_fakes)
781 catalog["ra"] = ra
782 catalog["dec"] = dec
783 catalog["host_id"] = hosts["id"]
784 catalog["host_flux"] = hosts[flux_type]
785 catalog["host_mag"] = hosts[mag_type]
786 catalog["host_ra"] = np.rad2deg(hosts["coord_ra"])
787 catalog["host_dec"] = np.rad2deg(hosts["coord_dec"])
788 catalog["delta_ra"] = (ra - catalog["host_ra"]) * 3600
789 catalog["delta_dec"] = (dec - catalog["host_dec"]) * 3600
790 catalog["delta_mag"] = catalog["mag"] - catalog["host_mag"]
791 catalog["host_a"] = a
792 catalog["host_b"] = b
793 catalog["host_pa"] = pa
794 catalog["source_type"] = "Star"
795 catalog["host_type"] = host_type
796 catalog["hosted_fake"] = np.ones(len(catalog), dtype=bool)
797 catalog = vstack([zero_table, catalog])
799 return catalog, hosts
801 def _generate_hosted_fake_catalog(self, n_fakes, sourceCat, wcs, rng, host_type='extended'):
802 """Generate a catalog of hosted fakes around selected host galaxies or stars."""
803 if n_fakes <= 0:
804 return Table(dtype=self._table_dtypes, masked=False), None
806 host_catalog = self.select_hosts(sourceCat, host_type=host_type)
807 if len(host_catalog) == 0:
808 self.log.warning("No valid hosts were selected for hosted fakes.")
809 return Table(dtype=self._table_dtypes, masked=False), None
811 # we allow more than one fake per host since we can deal with blending now
812 if n_fakes > len(host_catalog): 812 ↛ 821line 812 didn't jump to line 821 because the condition on line 812 was always true
813 self.log.warning(
814 "Requested %d hosted fakes, but only %d valid hosts were found. "
815 "Some hosts will have multiple fakes.",
816 n_fakes,
817 len(host_catalog),
818 )
819 idx = rng.choice(len(host_catalog), size=n_fakes, replace=True)
820 else:
821 idx = rng.choice(len(host_catalog), size=n_fakes, replace=False)
823 selected_hosts = host_catalog[idx]
824 return self._build_hosted_population(selected_hosts, wcs, rng, host_type=host_type)
826 def _generate_variable_twins(self, variable_catalog, rng):
827 """Generate twin sources for variable fakes with magnitude offsets."""
828 n_twins = len(variable_catalog)
829 if n_twins <= 0: 829 ↛ 830line 829 didn't jump to line 830 because the condition on line 829 was never true
830 return Table(dtype=self._table_dtypes, masked=False)
831 # magnitudes are drawn from a normal distribution around the variable source's magnitude
832 twins = variable_catalog.copy()
833 twins["mag_offset"] = rng.normal(
834 loc=self.config.variableFakeMean,
835 scale=self.config.variableFakeStd,
836 size=n_twins
837 )
838 twins["mag"] += twins["mag_offset"]
839 twins["isVisitSource"] = ~twins["isVisitSource"]
840 twins["isTemplateSource"] = ~twins["isTemplateSource"]
841 twins["twin_id"] = twins["injection_id"]
842 return twins
844 def _generate_blended_twins(self, blended_catalog, rng, wcs):
845 """Generate twin sources for blended fakes with magnitude and positional offsets."""
846 n_twins = len(blended_catalog)
847 if n_twins <= 0: 847 ↛ 848line 847 didn't jump to line 848 because the condition on line 847 was never true
848 return Table(dtype=self._table_dtypes, masked=False)
850 twins = blended_catalog.copy()
851 twins["mag_offset"] = rng.normal(
852 loc=0.0,
853 scale=self.config.blendedFakeMagOffset,
854 size=n_twins
855 )
856 delta_ra, delta_dec = self._draw_offset_components_arcsec(rng, n_twins)
857 twins["delta_ra"] = delta_ra
858 twins["delta_dec"] = delta_dec
859 twins["mag"] += twins["mag_offset"]
861 twins["ra"] += (twins["delta_ra"] / 3600.0) / np.cos(np.deg2rad(twins["dec"]))
862 twins["dec"] += twins["delta_dec"] / 3600.0
863 twins["x"], twins["y"] = wcs.skyToPixelArray(twins["ra"], twins["dec"], degrees=True)
865 twins["blend_twin_id"] = twins["injection_id"]
866 return twins
868 def _validate_catalog(self, catalog):
869 """Validate IDs and reciprocal source relationships."""
870 injection_ids = np.asarray(catalog["injection_id"], dtype=np.int64)
871 if len(np.unique(injection_ids)) != len(injection_ids): 871 ↛ 872line 871 didn't jump to line 872 because the condition on line 871 was never true
872 raise RuntimeError("Fake injection IDs are not unique.")
874 id_set = set(injection_ids)
876 variable = np.asarray(catalog["isVariable"], dtype=bool)
877 blended = np.asarray(catalog["isBlended"], dtype=bool)
879 if np.any(variable & blended): 879 ↛ 880line 879 didn't jump to line 880 because the condition on line 879 was never true
880 raise RuntimeError("Variable and blended fake populations overlap.")
882 twin_ids = np.asarray(catalog["twin_id"], dtype=np.int64)
883 variable_twin_ids = twin_ids[variable]
885 if np.any(variable_twin_ids <= 0): 885 ↛ 886line 885 didn't jump to line 886 because the condition on line 885 was never true
886 raise RuntimeError("Variable fake twin links are incomplete.")
888 if not set(variable_twin_ids).issubset(id_set): 888 ↛ 889line 888 didn't jump to line 889 because the condition on line 888 was never true
889 raise RuntimeError("Variable twin IDs do not refer to catalog rows.")
891 variable_rows_by_id = {
892 int(injection_id): index
893 for index, injection_id in enumerate(injection_ids)
894 if variable[index]
895 }
897 for index in np.flatnonzero(variable):
898 twin_index = variable_rows_by_id.get(int(twin_ids[index]))
899 if twin_index is None: 899 ↛ 900line 899 didn't jump to line 900 because the condition on line 899 was never true
900 raise RuntimeError("Variable twin ID does not identify a variable row.")
901 if twin_ids[twin_index] != injection_ids[index]: 901 ↛ 902line 901 didn't jump to line 902 because the condition on line 901 was never true
902 raise RuntimeError("Variable twin links are not reciprocal.")
904 # removing hosted blends from the check
905 non_hosted_blends = blended & ~np.asarray(catalog["hosted_fake"], dtype=bool)
906 blend_twin_ids = np.asarray(catalog["blend_twin_id"], dtype=np.int64)
907 blended_twin_ids = blend_twin_ids[non_hosted_blends]
909 if np.any(blended_twin_ids <= 0): 909 ↛ 910line 909 didn't jump to line 910 because the condition on line 909 was never true
910 raise RuntimeError("Blended fake twin links are incomplete.")
912 if not set(blended_twin_ids).issubset(id_set): 912 ↛ 913line 912 didn't jump to line 913 because the condition on line 912 was never true
913 raise RuntimeError("Blend twin IDs do not refer to catalog rows.")
915 rows_by_id = {
916 int(injection_id): index
917 for index, injection_id in enumerate(injection_ids)
918 }
920 for index in np.flatnonzero(blend_twin_ids):
921 twin_index = rows_by_id[int(blend_twin_ids[index])]
922 if blend_twin_ids[twin_index] != injection_ids[index]: 922 ↛ 923line 922 didn't jump to line 923 because the condition on line 922 was never true
923 raise RuntimeError("Blend twin links are not reciprocal.")
925 def _normalize_catalog(self, catalog, visit_id, detector_id):
926 """Normalize catalog columns to remove units and add visit/detector identifiers."""
927 catalog["visit"] = visit_id
928 catalog["detector"] = detector_id
929 for name in ("ra", "dec", "delta_ra", "delta_dec", "host_ra", "host_dec"):
930 if name in catalog.colnames and catalog[name].unit is not None: 930 ↛ 931line 930 didn't jump to line 931 because the condition on line 930 was never true
931 catalog[name] = catalog[name].value
933 def select_hosts(self, sourceCat, host_type='extended'):
934 """Select hosts from the source catalog."""
935 sky = ~sourceCat["sky_source"]
936 flags = ~sourceCat["base_ClassificationSizeExtendedness_flag"]
937 flags &= ~sourceCat["base_ClassificationExtendedness_flag"]
938 flags &= ~sourceCat["slot_Shape_flag"]
939 flags &= ~sourceCat["slot_Centroid_flag"]
940 flags &= ~sourceCat["base_PixelFlags_flag"]
942 if host_type == 'extended':
943 extended = sourceCat["base_ClassificationSizeExtendedness_value"] > 0.9
944 extended &= sourceCat["base_ClassificationExtendedness_value"] == 1
945 snr = sourceCat["slot_ModelFlux_flux"] / sourceCat["slot_ModelFlux_fluxErr"] > 15
946 return sourceCat[sky & flags & extended & snr].copy()
947 elif host_type == 'star': 947 ↛ 953line 947 didn't jump to line 953 because the condition on line 947 was always true
948 stellar = sourceCat["base_ClassificationSizeExtendedness_value"] < 0.9
949 stellar &= sourceCat["base_ClassificationExtendedness_value"] == 0
950 snr = sourceCat["slot_PsfFlux_flux"] / sourceCat["slot_PsfFlux_fluxErr"] > 30
951 return sourceCat[sky & flags & stellar & snr].copy()
952 else:
953 raise ValueError("host_type must be 'extended' or 'star'.")
955 def get_PA_and_axes(self, Ixx, Ixy, Iyy):
956 """Compute the position angle and semi-major/minor axes from second moments."""
957 theta = 0.5 * np.arctan2(2 * Ixy, Ixx - Iyy)
958 term = np.sqrt(((Ixx - Iyy) / 2) ** 2 + Ixy ** 2)
959 return theta, np.sqrt((Ixx + Iyy) / 2 + term), np.sqrt((Ixx + Iyy) / 2 - term)
961 def _cap_hosted_blended_count(self, hosted_fakes, hosts):
962 """Cap the number of hosted blended fakes based on configuration limits."""
963 if len(hosted_fakes) == 0 or len(hosts) == 0: 963 ↛ 964line 963 didn't jump to line 964 because the condition on line 963 was never true
964 return hosted_fakes, hosts
966 # check if we have multiple ids for the same host and cap the number of fakes per host
967 unique_hosts, counts = np.unique(hosted_fakes["host_id"], return_counts=True)
968 capped_hosts = []
969 for host_id, count in zip(unique_hosts, counts):
970 if count > self.config.maxHostedBlendedFakesPerHost:
971 self.log.warning(
972 "Capping hosted blended fakes for host_id %d from %d to %d.",
973 host_id, count, self.config.maxHostedBlendedFakesPerHost
974 )
975 capped_hosts.append(host_id)
976 row_ids = np.where(hosted_fakes["host_id"] == host_id)
977 row_ids_to_remove = row_ids[0][self.config.maxHostedBlendedFakesPerHost:]
978 hosted_fakes.remove_rows(row_ids_to_remove)
979 if self.config.maxHostedBlendedFakesTotal > 0:
980 if len(hosted_fakes) > self.config.maxHostedBlendedFakesTotal: 980 ↛ 986line 980 didn't jump to line 986 because the condition on line 980 was always true
981 self.log.warning(
982 "Capping total hosted blended fakes from %d to %d.",
983 len(hosted_fakes), self.config.maxHostedBlendedFakesTotal
984 )
985 hosted_fakes = hosted_fakes[:self.config.maxHostedBlendedFakesTotal]
986 return hosted_fakes, hosts
988 def _draw_offset_components_arcsec(self, rng, n_points):
989 """Draw random offsets in arcseconds for blended fakes.
990 The offsets are drawn from a uniform distribution in radius and angle."""
991 if n_points <= 0: 991 ↛ 992line 991 didn't jump to line 992 because the condition on line 991 was never true
992 return np.zeros(0), np.zeros(0)
993 rmin = max(0.0, self.config.blendedFakeMinOffset)
994 rmax = max(rmin, self.config.blendedFakeMaxOffset)
995 radius = np.sqrt(rng.uniform(rmin**2, rmax**2, n_points))
996 theta = rng.uniform(0, 2 * np.pi, n_points)
997 return radius * np.cos(theta), radius * np.sin(theta)
999 def get_n_fakes_from_density(self, visit_image, density):
1000 area = visit_image.getConvexPolygon().getBoundingBox().getArea()
1001 return np.round(density * area * (180 / np.pi) ** 2).astype(int)
1003 def _make_unique_injection_ids(self, n_ids, used_ids=None):
1004 """Generate unique 24-bit IDs for a catalog considering IDs already in use."""
1005 new_ids = set()
1006 if used_ids:
1007 # compile used ids from the existing catalog to avoid collisions
1008 if isinstance(used_ids, list): 1008 ↛ 1014line 1008 didn't jump to line 1014 because the condition on line 1008 was always true
1009 used_ids = np.concatenate(used_ids)
1010 used_ids = set(used_ids)
1011 else:
1012 used_ids = set()
1013 # check that the new ids are not in the used set already
1014 while len(new_ids) < n_ids:
1015 candidate_id = uuid.uuid4().int & ((1 << 24) - 1)
1016 if candidate_id > 0 and candidate_id not in used_ids.union(new_ids): 1016 ↛ 1014line 1016 didn't jump to line 1014 because the condition on line 1016 was always true
1017 new_ids.add(candidate_id)
1019 return np.asarray(list(new_ids), dtype=np.int64)