Coverage for python/lsst/ap/pipe/createApFakes.py: 90%
260 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-04 17:46 +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 deprecated.sphinx import deprecated
39__all__ = ["CreateRandomApFakesTask",
40 "CreateRandomApFakesConfig",
41 "CreateRandomApFakesConnections",
42 "CreateVisitDetectorFakesTask",
43 "CreateVisitDetectorFakesConfig",
44 "CreateVisitDetectorFakesConnections"]
47class CreateRandomApFakesConnections(PipelineTaskConnections,
48 dimensions=("tract", "skymap")):
49 skyMap = connTypes.Input(
50 doc="Input definition of geometry/bbox and projection/wcs for "
51 "template exposures",
52 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
53 dimensions=("skymap",),
54 storageClass="SkyMap",
55 )
56 fakeCat = connTypes.Output(
57 doc="Catalog of fake sources to draw inputs from.",
58 name="fakeSourceCat",
59 storageClass="DataFrame",
60 dimensions=("tract", "skymap")
61 )
64@deprecated(
65 reason="This task will be removed in v28.0 as it is replaced by `source_injection` tasks.",
66 version="v28.0",
67 category=FutureWarning,
68)
69class CreateRandomApFakesConfig(
70 InsertFakesConfig,
71 pipelineConnections=CreateRandomApFakesConnections):
72 """Config for CreateRandomApFakesTask. Copy from the InsertFakesConfig to
73 assert that columns created with in this task match that those expected in
74 the InsertFakes and related tasks.
75 """
76 fakeDensity = pexConfig.RangeField(
77 doc="Goal density of random fake sources per square degree. Default "
78 "value is roughly the density per square degree for ~10k sources "
79 "visit.",
80 dtype=float,
81 default=1000,
82 min=0,
83 )
84 filterSet = pexConfig.ListField(
85 doc="Set of Abstract filter names to produce magnitude columns for.",
86 dtype=str,
87 default=["u", "g", "r", "i", "z", "y"],
88 )
89 fraction = pexConfig.RangeField(
90 doc="Fraction of the created source that should be inserted into both "
91 "the visit and template images. Values less than 1 will result in "
92 "(1 - fraction) / 2 inserted into only visit or the template.",
93 dtype=float,
94 default=1/3,
95 min=0,
96 max=1,
97 )
98 magMin = pexConfig.RangeField(
99 doc="Minimum magnitude the mag distribution. All magnitudes requested "
100 "are set to the same value.",
101 dtype=float,
102 default=20,
103 min=1,
104 max=40,
105 )
106 magMax = pexConfig.RangeField(
107 doc="Maximum magnitude the mag distribution. All magnitudes requested "
108 "are set to the same value.",
109 dtype=float,
110 default=30,
111 min=1,
112 max=40,
113 )
114 visitSourceFlagCol = pexConfig.Field(
115 doc="Name of the column flagging objects for insertion into the visit "
116 "image.",
117 dtype=str,
118 default="isVisitSource"
119 )
120 templateSourceFlagCol = pexConfig.Field(
121 doc="Name of the column flagging objects for insertion into the "
122 "template image.",
123 dtype=str,
124 default="isTemplateSource"
125 )
128@deprecated(
129 reason="This task will be removed in v28.0 as it is replaced by `source_injection` tasks.",
130 version="v28.0",
131 category=FutureWarning,
132)
133class CreateRandomApFakesTask(PipelineTask):
134 """Create and store a set of spatially uniform star fakes over the sphere
135 for use in AP processing. Additionally assign random magnitudes to said
136 fakes and assign them to be inserted into either a visit exposure or
137 template exposure.
138 """
140 _DefaultName = "createApFakes"
141 ConfigClass = CreateRandomApFakesConfig
143 def runQuantum(self, butlerQC, inputRefs, outputRefs):
144 inputs = butlerQC.get(inputRefs)
145 inputs["tractId"] = butlerQC.quantum.dataId["tract"]
147 outputs = self.run(**inputs)
148 butlerQC.put(outputs, outputRefs)
150 def run(self, tractId, skyMap):
151 """Create a set of uniform random points that covers a tract.
153 Parameters
154 ----------
155 tractId : `int`
156 Tract id to produce randoms over.
157 skyMap : `lsst.skymap.SkyMap`
158 Skymap to produce randoms over.
160 Returns
161 -------
162 randoms : `pandas.DataFrame`
163 Catalog of random points covering the given tract. Follows the
164 columns and format expected in `lsst.pipe.tasks.InsertFakes`.
165 """
167 tract = skyMap.generateTract(tractId)
168 tractArea = tract.getOuterSkyPolygon().getBoundingBox().getArea()
169 tractArea *= (180 / np.pi) ** 2
170 tractWcs = tract.getWcs()
171 vertexList = tract.getVertexList()
172 vertexRas = [vertex.getRa().asDegrees() for vertex in vertexList]
173 vertexDecs = [vertex.getDec().asDegrees() for vertex in vertexList]
175 catalog = generate_injection_catalog(
176 ra_lim=sorted([np.min(vertexRas), np.max(vertexRas)]),
177 dec_lim=sorted([np.min(vertexDecs), np.max(vertexDecs)]),
178 mag_lim=(self.config.magMin, self.config.magMax),
179 density=self.config.fakeDensity,
180 source_type="Star",
181 seed=str(tractId),
182 wcs=tractWcs
183 )
185 nFakes = len(catalog)
187 self.log.info(
188 f"Creating {nFakes} star fakes over tractId={tractId} with "
189 f" RA in ({sorted([np.min(vertexRas), np.max(vertexRas)])} "
190 f" Dec in ({sorted([np.min(vertexDecs), np.max(vertexDecs)])}), "
191 f"area={tractArea:.4f} deg^2 and "
192 f"magnitude range: [{self.config.magMin, self.config.magMax}]")
194 onesColumn = np.ones(nFakes, dtype="float")
195 zerosColumn = np.zeros(nFakes, dtype="float")
196 # Concatenate the data and add dummy values for the unused variables.
197 # Set all data to PSF like objects.
198 mags = np.asarray(catalog["mag"], dtype=float)
199 randData = {
200 "fakeId": [uuid.uuid4().int & (1 << 64) - 1 for n in range(nFakes)],
201 self.config.ra_col: np.asarray(catalog["ra"], dtype=float),
202 self.config.dec_col: np.asarray(catalog["dec"], dtype=float),
203 **self.createVisitCoaddSubdivision(nFakes),
204 **self.createMagnitudeColumns(mags),
205 self.config.disk_semimajor_col: onesColumn,
206 self.config.bulge_semimajor_col: onesColumn,
207 self.config.disk_n_col: onesColumn,
208 self.config.bulge_n_col: onesColumn,
209 self.config.disk_axis_ratio_col: onesColumn,
210 self.config.bulge_axis_ratio_col: onesColumn,
211 self.config.disk_pa_col: zerosColumn,
212 self.config.bulge_pa_col: onesColumn,
213 self.config.sourceType: np.asarray(catalog["source_type"], dtype=str),
214 "source_type": np.asarray(catalog["source_type"], dtype=str),
215 "injection_id": np.asarray(catalog["injection_id"], dtype=np.int64)
216 }
218 return Struct(fakeCat=pd.DataFrame(data=randData))
220 def createVisitCoaddSubdivision(self, nFakes):
221 """Assign a given fake either a visit image or coadd or both based on
222 the ``faction`` config value.
224 Parameters
225 ----------
226 nFakes : `int`
227 Number of fakes to create.
229 Returns
230 -------
231 output : `dict`[`str`, `numpy.ndarray`]
232 Dictionary of boolean arrays specifying which image to put a
233 given fake into.
234 """
235 nBoth = int(self.config.fraction * nFakes)
236 nOnly = int((1 - self.config.fraction) / 2 * nFakes)
237 isVisitSource = np.zeros(nFakes, dtype=bool)
238 isTemplateSource = np.zeros(nFakes, dtype=bool)
239 if nBoth > 0: 239 ↛ 242line 239 didn't jump to line 242 because the condition on line 239 was always true
240 isVisitSource[:nBoth] = True
241 isTemplateSource[:nBoth] = True
242 if nOnly > 0: 242 ↛ 246line 242 didn't jump to line 246 because the condition on line 242 was always true
243 isVisitSource[nBoth:(nBoth + nOnly)] = True
244 isTemplateSource[(nBoth + nOnly):] = True
246 return {self.config.visitSourceFlagCol: isVisitSource,
247 self.config.templateSourceFlagCol: isTemplateSource}
249 def createMagnitudeColumns(self, mags):
250 """Create magnitude columns from a 1D magnitude array.
252 Parameters
253 ----------
254 mags : `numpy.ndarray`
255 Magnitudes to copy to all configured filter-band columns.
257 Returns
258 -------
259 randMags : `dict`[`str`, `numpy.ndarray`]
260 Dictionary containing per-band magnitudes plus a ``mag`` column
261 compatible with ``source_injection`` catalogs.
262 """
263 randMags = {}
264 for fil in self.config.filterSet:
265 randMags[self.config.mag_col % fil] = mags
266 randMags["mag"] = mags
267 return randMags
270class CreateVisitDetectorFakesConnections(
271 PipelineTaskConnections,
272 defaultTemplates={"coaddName": "deep"},
273 dimensions=("instrument",
274 "visit",
275 "detector")):
277 sourceCat = connTypes.Input(
278 doc="Catalog of sources detected on the calibrated exposure; ",
279 name="single_visit_star_reprocessed_footprints",
280 storageClass="SourceCatalog",
281 dimensions=["instrument", "visit", "detector"],
282 )
283 visit_image = connTypes.Input(
284 doc="Calibrated exposure to inject synthetic sources into.",
285 name="preliminary_visit_image",
286 storageClass="ExposureF",
287 dimensions=["instrument", "visit", "detector"],
288 )
289 outputCat = connTypes.Output(
290 doc="Catalog of fake sources to draw inputs from.",
291 name="VisitDetectorFakeSourceCat",
292 storageClass="ArrowAstropy",
293 dimensions=["instrument", "visit", "detector"],
294 )
297class CreateVisitDetectorFakesConfig(
298 PipelineTaskConfig,
299 pipelineConnections=CreateVisitDetectorFakesConnections
300):
301 """Config for CreateVisitDetectorFakesTask."""
302 randomFakeDensity = pexConfig.RangeField(
303 doc="Goal density of visit detector fake sources per square degree.",
304 dtype=float,
305 default=1000,
306 min=1,
307 )
308 nRandomFakes = pexConfig.RangeField(
309 doc="Number of random fakes to add to the visit detector. Overrides "
310 "the randomFakeDensity if set to a positive value.",
311 dtype=int,
312 default=-1,
313 min=-1,
314 )
315 doAddRandomVisitFakes = pexConfig.Field(
316 doc="Whether to add random positive fakes to the visit detector.",
317 dtype=bool,
318 default=True,
319 )
320 doAddRandomTemplateFakes = pexConfig.Field(
321 doc="Whether to add random template fakes to the visit detector (negatives).",
322 dtype=bool,
323 default=True,
324 )
325 templateFakeFraction = pexConfig.RangeField(
326 doc="Fraction of random fakes that should be added to the template image."
327 " The rest will be added to the visit image.",
328 dtype=float,
329 default=0.25,
330 min=0,
331 max=1,
332 )
333 doAddVariableFakes = pexConfig.Field(
334 doc="Whether to add variable fakes to the visit detector.",
335 dtype=bool,
336 default=False,
337 )
338 variableFakeFraction = pexConfig.RangeField(
339 doc="Fraction of variable fakes that should be added to the template image."
340 " The rest will be added to the visit image.",
341 dtype=float,
342 default=0.1,
343 min=0,
344 max=1,
345 )
346 variableFakeMean = pexConfig.RangeField(
347 doc="Mean magnitude variation for variable fakes.",
348 dtype=float,
349 default=0.0,
350 min=-1,
351 max=1,
352 )
353 variableFakeStd = pexConfig.RangeField(
354 doc="Standard deviation of magnitude variation for variable fakes.",
355 dtype=float,
356 default=0.5,
357 min=0,
358 max=1,
359 )
360 doAddHostedFakes = pexConfig.Field(
361 doc="Whether to add hosted fakes to the visit detector.",
362 dtype=bool,
363 default=False,
364 )
365 fracHostedFakes = pexConfig.RangeField(
366 doc="Fraction of hosts with fakes to add to the visit detector.",
367 dtype=float,
368 default=0.1,
369 min=0,
370 max=1,
371 )
372 minHostedFakes = pexConfig.RangeField(
373 doc="Minimum number of hosted fakes to add to the visit detector.",
374 dtype=int,
375 default=20,
376 min=1,
377 )
378 doAddModelFakes = pexConfig.Field(
379 doc="Whether to add model fakes to the visit detector.",
380 dtype=bool,
381 default=False,
382 )
383 magMin = pexConfig.Field(
384 doc="Minimum magnitude for the fake sources.",
385 dtype=float,
386 default=20,
387 )
388 magMax = pexConfig.Field(
389 doc="Maximum magnitude for the fake sources.",
390 dtype=float,
391 default=26,
392 )
395class CreateVisitDetectorFakesTask(PipelineTask):
396 """Create and store a set of visit detector fakes for use in AP processing.
397 This task creates a catalog of fake sources that can be used to inject
398 sources into visit detector images.
399 """
401 _DefaultName = "createVisitDetectorFakes"
402 ConfigClass = CreateVisitDetectorFakesConfig
404 def __init__(self, **kwargs):
405 super().__init__(**kwargs)
406 self.log = logging.getLogger(__name__)
408 def runQuantum(self, butlerQC, inputRefs, outputRefs):
409 inputs = butlerQC.get(inputRefs)
410 outputs = self.run(**inputs)
411 butlerQC.put(outputs, outputRefs)
413 def _make_unique_injection_ids(self, n_ids, used_ids=None):
414 """Generate collision-free injection IDs within the 24-bit ID space."""
415 used = set() if used_ids is None else {int(value) for value in used_ids}
416 injection_ids = []
418 while len(injection_ids) < n_ids:
419 candidate = uuid.uuid4().int & ((1 << 24) - 1)
420 if candidate in used: 420 ↛ 421line 420 didn't jump to line 421 because the condition on line 420 was never true
421 continue
422 used.add(candidate)
423 injection_ids.append(candidate)
425 return np.asarray(injection_ids, dtype=np.int64)
427 def run(self, sourceCat, visit_image):
428 """Create a set of visit detector fakes.
430 Parameters
431 ----------
432 sourceCat : `lsst.afw.table.SourceCatalog`
433 Catalog of sources detected on the calibrated exposure.
434 visit_image : `lsst.afw.image.Exposure`
435 Visit image to inject synthetic sources into.
437 Returns
438 -------
439 outputCat : `astropy.table.Table`
440 Catalog of fake sources to draw inputs from.
441 """
442 # Use the visit+detector ids as the random seed.
443 visitId = visit_image.getInfo().getVisitInfo().id
444 detId = visit_image.detector.getId()
445 rng = np.random.default_rng([visitId, detId])
447 # set of catalogs to concatenate at the end
448 catalog_set = []
450 photoCalib = visit_image.getPhotoCalib()
451 bbox = visit_image.getBBox()
452 xmin, xmax = bbox.getMinX(), bbox.getMaxX()
453 ymin, ymax = bbox.getMinY(), bbox.getMaxY()
454 visit_stats = visit_image.getInfo().getSummaryStats()
456 wcs = visit_image.getWcs()
457 magLim = visit_stats.magLim
458 max_mag = np.min([magLim+1, self.config.magMax])
460 if self.config.doAddRandomVisitFakes:
461 # Generate random visit fakes
462 self.log.info("Generating random visit fakes.")
463 if self.config.nRandomFakes > 0: 463 ↛ 467line 463 didn't jump to line 467 because the condition on line 463 was always true
464 self.log.info(f"Generating random visit fakes with nRandomFakes={self.config.nRandomFakes}.")
465 n_random_fakes = self.config.nRandomFakes
466 else:
467 self.log.info(
468 f"Generating random visit fakes with randomFakeDensity={self.config.randomFakeDensity}.")
469 n_random_fakes = self.get_n_fakes_from_density(
470 visit_image=visit_image,
471 density=self.config.randomFakeDensity
472 )
473 self.log.info(f"Calculated n_random_fakes={n_random_fakes}.")
475 # draw random x-y coordinates
476 x_ssi = rng.uniform(xmin, xmax, size=n_random_fakes)
477 y_ssi = rng.uniform(ymin, ymax, size=n_random_fakes)
478 mags = rng.uniform(self.config.magMin, max_mag, size=n_random_fakes)
479 ra_ssi, dec_ssi = wcs.pixelToSkyArray(x_ssi, y_ssi, degrees=True)
481 random_catalog = Table()
482 random_catalog["x"] = x_ssi
483 random_catalog["y"] = y_ssi
484 random_catalog["mag"] = mags
485 random_catalog["ra"] = ra_ssi
486 random_catalog["dec"] = dec_ssi
487 random_catalog["source_type"] = "Star"
488 random_catalog["isVisitSource"] = True
489 random_catalog["isTemplateSource"] = False
490 catalog_set.append(random_catalog)
491 # Ignore now the possibility of _just_ template fakes
492 if self.config.doAddModelFakes: 492 ↛ 494line 492 didn't jump to line 494 because the condition on line 492 was never true
493 # Generate model fakes
494 self.log.info("Not implemented yet model fakes.")
495 # Placeholder for actual model fake generation logic
496 pass
498 if self.config.doAddHostedFakes:
499 # Generate hosted fakes
500 self.log.info("Generating hosted fakes.")
501 # Select hosts that look like extended sources.
502 hostcatalog = photoCalib.calibrateCatalog(sourceCat).asAstropy()
503 hostcatalog = self.select_hosts(hostcatalog)
504 n_hosts = len(hostcatalog)
505 if n_hosts == 0:
506 self.log.warning("Hosted fake generation requested, but no valid hosts were selected.")
507 else:
508 requested_n_fakes = max(
509 int(self.config.fracHostedFakes * n_hosts), self.config.minHostedFakes
510 )
511 n_fakes = min(requested_n_fakes, n_hosts)
512 if n_fakes < requested_n_fakes: 512 ↛ 520line 512 didn't jump to line 520 because the condition on line 512 was always true
513 self.log.warning(
514 "Reducing hosted fake count from %d to %d because only %d hosts are available.",
515 requested_n_fakes,
516 n_fakes,
517 n_hosts,
518 )
520 idx = rng.choice(n_hosts, size=n_fakes, replace=False)
521 hostcat = hostcatalog[idx]
523 x_hosts = hostcat['slot_Centroid_x']
524 y_hosts = hostcat['slot_Centroid_y']
525 mag_hosts = hostcat['slot_ModelFlux_mag']
526 # the units below are pixels and radians
527 pa, a, b = self.get_PA_and_axes(
528 hostcat['slot_Shape_xx'],
529 hostcat['slot_Shape_xy'],
530 hostcat['slot_Shape_yy']
531 )
532 # random radius and angle for the fake around the host
533 theta = rng.uniform(0, 2 * np.pi, size=n_fakes)
534 angle = np.sqrt((a*np.cos(theta))**2 + (b*np.sin(theta))**2)
535 radii = angle * np.sqrt(rng.uniform(0, 6, size=n_fakes))
537 # Polar -> Cartesian wrt the host in the PA coordinate system
538 xs = radii * np.cos(theta)
539 ys = radii * np.sin(theta)
541 # Retrieve the right position removing the galaxy orientation PA
542 x_rots = xs * np.cos(pa) - ys * np.sin(pa)
543 y_rots = xs * np.sin(pa) + ys * np.cos(pa)
545 x_ssi = x_hosts + x_rots
546 y_ssi = y_hosts + y_rots
548 # retrieving the global ra dec position of the injection
549 ra_ssi, dec_ssi = wcs.pixelToSkyArray(x_ssi, y_ssi, degrees=True)
550 delta_ra = (ra_ssi - np.rad2deg(hostcat['coord_ra'])) * 3600.
551 delta_dec = (dec_ssi - np.rad2deg(hostcat['coord_dec'])) * 3600.
553 delta_mag = rng.normal(loc=1, scale=1, size=n_fakes)
554 mags = mag_hosts + delta_mag
556 # Create the table of hosted fakes
557 hosted_fakes = Table()
558 hosted_fakes["x"] = x_ssi
559 hosted_fakes["y"] = y_ssi
560 hosted_fakes["mag"] = mags
561 hosted_fakes["ra"] = ra_ssi
562 hosted_fakes["dec"] = dec_ssi
563 hosted_fakes["host_id"] = hostcat['id']
564 hosted_fakes["host_flux"] = hostcat['slot_ModelFlux_flux']
565 hosted_fakes["host_mag"] = hostcat['slot_ModelFlux_mag']
566 hosted_fakes["host_ra"] = np.rad2deg(hostcat['coord_ra'])
567 hosted_fakes["host_dec"] = np.rad2deg(hostcat['coord_dec'])
568 hosted_fakes["delta_ra"] = delta_ra
569 hosted_fakes["delta_dec"] = delta_dec
570 hosted_fakes["delta_mag"] = delta_mag
571 hosted_fakes["host_a"] = a
572 hosted_fakes["host_b"] = b
573 hosted_fakes["host_pa"] = pa
574 hosted_fakes["source_type"] = "Star"
575 hosted_fakes["hosted_fake"] = True
576 hosted_fakes["isVisitSource"] = True
577 hosted_fakes["isTemplateSource"] = False
579 catalog_set.append(hosted_fakes)
581 if not catalog_set:
582 raise RuntimeError(
583 "No fake sources were generated. Enable at least one fakes mode or provide usable hosts."
584 )
586 catalog = vstack(catalog_set)
587 catalog['injection_id'] = self._make_unique_injection_ids(len(catalog))
589 if self.config.doAddRandomTemplateFakes:
590 is_tmplt_fake = rng.random(len(catalog)) < self.config.templateFakeFraction
591 catalog["isTemplateSource"] = is_tmplt_fake
592 catalog["isVisitSource"] = ~is_tmplt_fake
593 else:
594 catalog["isVisitSource"] = True
595 catalog["isTemplateSource"] = False
597 if self.config.doAddVariableFakes:
598 # Generate variable fakes by duplicating some fakes and adding the counterpart
599 # either science or template with a magnitude offset drawn from a normal
600 # distribution with mean and std defined in the config.
601 self.log.info("Generating variable fakes.")
602 n_variable_fakes = int(len(catalog) * self.config.variableFakeFraction)
603 idx = rng.choice(len(catalog), size=n_variable_fakes, replace=False)
604 variable_fakes = catalog[idx].copy()
605 variable_fakes["mag_offset"] = rng.normal(
606 loc=self.config.variableFakeMean,
607 scale=self.config.variableFakeStd,
608 size=n_variable_fakes
609 )
610 variable_fakes["mag"] += variable_fakes["mag_offset"]
611 # we flip the source, so for example if it was a visit, we trasnform it into a template
612 # with the idea of having duplicate injections, in the same location
613 variable_fakes["isVisitSource"] = ~variable_fakes["isVisitSource"]
614 variable_fakes["isTemplateSource"] = ~variable_fakes["isTemplateSource"]
616 variable_fakes["twin_id"] = variable_fakes["injection_id"]
617 variable_fakes["injection_id"] = self._make_unique_injection_ids(
618 len(variable_fakes),
619 used_ids=catalog["injection_id"],
620 )
621 # create column of isVariable flag
622 catalog["isVariable"] = np.where(np.isin(np.arange(len(catalog)), idx), True, False)
623 variable_fakes["isVariable"] = True
625 catalog = vstack([catalog, variable_fakes])
627 if len(catalog) > len(np.unique(catalog["injection_id"])): 627 ↛ 628line 627 didn't jump to line 628 because the condition on line 627 was never true
628 self.log.warning("Duplicate injection IDs detected after catalog assembly; reassigning them.")
629 old_injection_ids = np.asarray(catalog["injection_id"], dtype=np.int64)
630 new_injection_ids = self._make_unique_injection_ids(len(catalog))
631 # re-assign fresh injection ids
632 catalog["injection_id"] = new_injection_ids
633 if "twin_id" in catalog.colnames:
634 id_map = {old_id: new_id for old_id, new_id in zip(old_injection_ids, new_injection_ids)}
635 catalog["twin_id"] = np.asarray(
636 [id_map.get(int(twin_id), int(twin_id)) for twin_id in catalog["twin_id"]],
637 dtype=np.int64,
638 )
640 catalog["visit"] = visitId
641 catalog["detector"] = detId
643 return Struct(outputCat=catalog)
645 def select_hosts(self, sourceCat):
646 """
647 Selects host sources from a given source catalog based on a series of classification and flux cuts.
648 The selection criteria are:
649 - The 'base_ClassificationSizeExtendedness_flag' and
650 'base_ClassificationExtendedness_flag' must both be False.
651 - The 'base_ClassificationSizeExtendedness_value' must be greater than 0.9.
652 - The 'base_ClassificationExtendedness_value' must be equal to 1.
653 - The 'base_PsfFlux_flux' must be greater than 0.
654 Parameters
655 ----------
656 sourceCat : SourceCatalog
657 The source catalog containing the columns required for selection.
658 *args, **kwargs
659 Additional arguments (not used).
660 Returns
661 -------
662 hostCat : ArrowAstropy
663 A deep copy of the subset of the source catalog that passes all selection criteria.
664 """
666 # Avoid calibration stars or psf stars; remove flagged sources sky_sources
667 skySourceCut = ~sourceCat['sky_source']
669 flagCut = ~sourceCat['base_ClassificationSizeExtendedness_flag']
670 flagCut &= ~sourceCat['base_ClassificationExtendedness_flag']
671 flagCut &= ~sourceCat['slot_Shape_flag']
672 flagCut &= ~sourceCat['slot_Centroid_flag']
673 flagCut &= ~sourceCat['base_PixelFlags_flag']
675 extendednessCut = sourceCat['base_ClassificationSizeExtendedness_value'] > 0.9
676 extendednessCut &= sourceCat['base_ClassificationExtendedness_value'] == 1
678 snrCut = sourceCat['slot_ModelFlux_flux']/sourceCat['slot_ModelFlux_fluxErr'] > 15
680 hostCat = sourceCat[
681 skySourceCut & flagCut & extendednessCut & snrCut].copy()
682 return hostCat
684 def get_PA_and_axes(self, Ixx, Ixy, Iyy):
685 '''
686 Calculates the orientation and extent of an object based on its second moments.
688 Parameters:
689 Ixx (float): Second moment of the object along the x-axis. Often in degree²
690 Ixy (float): Second moment of the object along the x and y-axes. Often in degree²
691 Iyy (float): Second moment of the object along the y-axis. Often in degree²
693 Returns:
694 tuple: A tuple containing:
695 - theta (float): The orientation angle of the object in radians.
696 - a (float): The semi-major axis length of the object.
697 - b (float): The semi-minor axis length of the object.
698 '''
699 # Calculate position angle (orientation)
700 theta = 0.5 * np.arctan2(2 * Ixy, Ixx - Iyy)
702 # Calculate eigenvalues of the moment matrix
703 term1 = (Ixx + Iyy) / 2
704 term2 = np.sqrt(((Ixx - Iyy) / 2) ** 2 + Ixy ** 2)
705 lambda1 = term1 + term2
706 lambda2 = term1 - term2
708 a = np.sqrt(lambda1)
709 b = np.sqrt(lambda2)
711 return theta, a, b
713 def get_n_fakes_from_density(self, visit_image, density):
714 """Calculate the area of the injection limits in square degrees based on the RA and Dec limits."""
715 image_area = visit_image.getConvexPolygon().getBoundingBox().getArea()
716 image_area *= (180 / np.pi) ** 2
717 number = np.round(density * image_area).astype(int)
718 return number