Coverage for python / lsst / meas / extensions / multiprofit / fit_coadd_multiband.py: 21%
498 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-15 07:25 -0300
1# This file is part of meas_extensions_multiprofit.
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/>.
22__all__ = (
23 "PsfFitSuccessActionBase",
24 "PsfComponentsActionBase",
25 "SourceTablePsfFitSuccessAction",
26 "SourceTablePsfComponentsAction",
27 "MagnitudeDependentSizePriorConfig",
28 "ModelInitializer",
29 "MakeInitializerActionBase",
30 "BasicModelInitializer",
31 "CachedBasicModelInitializer",
32 "InitialInputData",
33 "MakeBasicInitializerAction",
34 "MakeCachedBasicInitializerAction",
35 "MultiProFitSourceConfig",
36 "CatalogExposurePsfs",
37 "MultiProFitSourceFitter",
38 "MultiProFitSourceTask",
39)
41from abc import ABC, abstractmethod
42from functools import cached_property
43import logging
44import math
45from typing import Any, ClassVar, Iterable, Mapping, Sequence
47from astropy.table import Table
48import astropy.units as u
49import lsst.afw.geom
50from lsst.daf.butler.formatters.parquet import astropy_to_arrow
51import lsst.gauss2d as g2
52import lsst.gauss2d.fit as g2f
53from lsst.multiprofit.errors import NoDataError, PsfRebuildFitFlagError
54from lsst.multiprofit.fitting.fit_psf import CatalogPsfFitterConfig, CatalogPsfFitterConfigData
55from lsst.multiprofit.fitting.fit_source import (
56 CatalogExposureSourcesABC,
57 CatalogSourceFitterABC,
58 CatalogSourceFitterConfig,
59 CatalogSourceFitterConfigData,
60)
61from lsst.multiprofit.modeller import Model
62from lsst.multiprofit.utils import frozen_arbitrary_allowed_config, get_params_uniq, set_config_from_dict
63import lsst.pex.config as pexConfig
64from lsst.pex.config.configurableActions import ConfigurableAction, ConfigurableActionField
65import lsst.pipe.base as pipeBase
66import lsst.pipe.tasks.fit_coadd_multiband as fitMB
67import lsst.utils.timer as utilsTimer
68import numpy as np
69import pydantic
71from .errors import IsParentError, NotPrimaryError
72from .input_config import InputConfig
73from .utils import get_spanned_image
75_LOG = logging.getLogger(__name__)
76TWO_SQRT_PI = 2 * math.sqrt(np.pi)
79class PsfFitSuccessActionBase(ConfigurableAction):
80 """Base action to return whether a source had a succesful PSF fit."""
82 def get_schema(self) -> list[str]:
83 """Return the list of columns required to call this action."""
84 raise NotImplementedError("This method must be overloaded in subclasses")
86 def __call__(self, source: Mapping[str, Any], *args: Any, **kwargs: Any) -> bool:
87 raise NotImplementedError("This method must be overloaded in subclasses")
90class PsfComponentsActionBase(ConfigurableAction):
91 """Base action to return a list of Gaussians from a source mapping.
93 This base class should be used as a sentinel when using a MultiProFit PSF
94 fit table, and only needs to be specialized for external PSF fitters.
95 """
97 def get_schema(self) -> list[str]:
98 """Return the list of columns required to call this action."""
99 raise NotImplementedError("This method must be overloaded in subclasses")
101 def __call__(self, source: Mapping[str, Any], *args: Any, **kwargs: Any) -> list[g2.Gaussian]:
102 raise NotImplementedError("This method must be overloaded in subclasses")
105class SourceTablePsfFitSuccessAction(PsfFitSuccessActionBase):
106 """Action to return PSF fit status from a SourceTable row."""
108 flag_format = pexConfig.Field[str](
109 doc="Format for the flag field; flag_prefix, flag_suffix and flag_sub are substituted",
110 default="{flag_prefix}{flag_suffix}{flag_sub}",
111 )
112 flag_prefix = pexConfig.Field[str](
113 doc="Prefix for the key for the summed flag field",
114 default="modelfit_DoubleShapeletPsfApprox",
115 )
116 flag_suffix = pexConfig.Field[str](
117 doc="Suffix for all flag fields",
118 default="_flag",
119 )
120 flags_sub = pexConfig.ListField[str](
121 doc="Suffixes for specific flag fields that must not be true",
122 default=["_invalidPointForPsf", "_invalidMoments", "_maxIterations"],
123 )
125 def _format(self, flag_sub: str) -> str:
126 return self.flag_format.format(
127 flag_prefix=self.flag_prefix,
128 flag_sub=flag_sub,
129 flag_suffix=self.flag_suffix,
130 )
132 def get_schema(self) -> Iterable[str]:
133 for flag_sub in self.flags_sub:
134 yield self._format(flag_sub=flag_sub)
136 def __call__(self, source: Mapping[str, Any], *args: Any, **kwargs: Any) -> bool:
137 good = True
138 for flag_sub in self.flags_sub:
139 good &= not source[self._format(flag_sub=flag_sub)]
140 return good
143class SourceTablePsfComponentsAction(PsfComponentsActionBase):
144 """Action to return PSF components from a SourceTable.
146 This is anticipated to be a deepCoadd_meas with PSF fit parameters from a
147 measurement plugin returning covariance matrix terms.
148 """
150 action_source = ConfigurableActionField[PsfFitSuccessActionBase](
151 doc="Action to return whether the PSF fit was successful for a single source row",
152 default=SourceTablePsfFitSuccessAction,
153 )
154 format = pexConfig.Field[str](
155 doc="Format for the field names, where {idx_comp} is the index of the component and {moment}"
156 "is the name of the moment (xx, xy or yy, integral)",
157 default="modelfit_DoubleShapeletPsfApprox_{idx_comp}_{moment}",
158 )
159 name_moment_xx = pexConfig.Field[str](doc="Name of the xx (2nd x-axis) moment", default="xx")
160 name_moment_xy = pexConfig.Field[str](doc="Name of the xy (covariance term) moment", default="xy")
161 name_moment_yy = pexConfig.Field[str](doc="Name of the yy (2nd y-axis) moment", default="yy")
162 name_moment_integral = pexConfig.Field[str](doc="Name of the integral (zeroth) moment", default="0")
163 n_components = pexConfig.Field[int](
164 doc="Number of Gaussian components",
165 default=2,
166 check=lambda x: x >= 2,
167 )
169 @staticmethod
170 def get_integral(moment_zero) -> float:
171 """Get the total integrated flux from a zeroth moment value.
173 The zeroth moment is simply the integrated flux divided by a
174 constant value of 2*sqrt(pi).
176 Parameters
177 ----------
178 moment_zero
179 The zeroth moment value.
181 Returns
182 -------
183 integral
184 The total integrated weight (flux).
185 """
186 return moment_zero * TWO_SQRT_PI
188 def get_schema(self) -> list[str]:
189 names_moments = (
190 self.name_moment_xx,
191 self.name_moment_yy,
192 self.name_moment_xy,
193 self.name_moment_integral,
194 )
195 columns = [
196 column
197 for idx_comp in range(self.n_components)
198 for column in (
199 self.format.format(name_moment=name_moment, idx_comp=idx_comp)
200 for name_moment in names_moments
201 )
202 ] + self.action_source.get_schema()
203 return columns
205 def __call__(self, source: Mapping[str, Any], *args: Any, **kwargs: Any) -> list[g2.Gaussian]:
206 if not self.action_source(source):
207 raise PsfRebuildFitFlagError(
208 f"PSF fit failed due to action based on schema: {self.action_source.get_schema()}"
209 )
210 gaussians = [None] * self.n_components
211 for idx_comp in range(self.n_components):
212 gaussian = g2.Gaussian(
213 ellipse=g2.Ellipse(
214 g2.Covariance(
215 sigma_x_sq=source[self.format.format(moment=self.name_moment_xx, idx_comp=idx_comp)],
216 sigma_y_sq=source[self.format.format(moment=self.name_moment_yy, idx_comp=idx_comp)],
217 cov_xy=source[self.format.format(moment=self.name_moment_xy, idx_comp=idx_comp)],
218 )
219 ),
220 integral=g2.GaussianIntegralValue(
221 value=self.get_integral(
222 source[self.format.format(moment=self.name_moment_integral, idx_comp=idx_comp)]
223 )
224 ),
225 )
226 gaussians[idx_comp] = gaussian
227 return gaussians
230class MagnitudeDependentSizePriorConfig(pexConfig.Config):
231 """Configuration for a magnitude-dependent size prior.
233 Defaults are for ugrizy total mag and log10(r_eff/arcsec).
234 """
236 intercept_mag = pexConfig.Field[float](
237 doc="The magnitude at which no adjustment is applied",
238 default=18.0,
239 )
240 slope_median_per_mag = pexConfig.Field[float](
241 doc="The slope in the median size, in dex per mag",
242 default=-0.15,
243 )
244 slope_stddev_per_mag = pexConfig.Field[float](
245 doc="The slope in the standard deviation of the size, in dex per mag",
246 default=0.0,
247 )
250class ModelInitializer(ABC, pydantic.BaseModel):
251 """An interface for a configurable model initializer based on priors
252 and optional external data.
253 """
255 model_config: ClassVar[pydantic.ConfigDict] = frozen_arbitrary_allowed_config
257 inputs: dict[str, Any] = pydantic.Field(
258 title="Additional external inputs used in initialization",
259 default_factory=dict,
260 )
261 priors_shape_mag: dict = pydantic.Field(
262 title="Magnitude-dependent shape prior configurations",
263 default_factory=dict,
264 )
266 @abstractmethod
267 def initialize_model(
268 self,
269 model: Model,
270 source: Mapping[str, Any],
271 catexps: list[CatalogExposureSourcesABC],
272 config_data: CatalogSourceFitterConfigData,
273 values_init: Mapping[g2f.ParameterD, float] | None = None,
274 **kwargs,
275 ):
276 """Initialize a MultiProFit model for a single object corresponding
277 to a row in a catalog.
279 Parameters
280 ----------
281 model
282 The model to initialize parameter values for.
283 source
284 A mapping with fields expected to be populated in the
285 corresponding source catalog for initialization.
286 catexps
287 Per-band catalog-exposure pairs.
288 config_data
289 Fitter configuration and data.
290 values_init
291 Default initial values for parameters.
292 **kwargs
293 Additional keyword arguments for any purpose.
294 """
295 raise NotImplementedError(f"{self.__name__} must implement initialize_model")
298class MakeInitializerActionBase(ConfigurableAction):
299 """An interface for an action that creates an initializer."""
301 def __call__(
302 self,
303 catalog_multi: Sequence,
304 catexps: list[fitMB.CatalogExposureInputs],
305 config_data: CatalogSourceFitterConfigData,
306 **kwargs,
307 ) -> ModelInitializer:
308 """Make a ModelInitializer object that can initialize model
309 parameter values for a given object in a catalog.
311 Parameters
312 ----------
313 catalog_multi
314 The multiband catalog with one row per object to fit.
315 catexps
316 Per-band catalog-exposure pairs.
317 config_data
318 Fitter configuration and data.
319 **kwargs
320 Additional arguments to pass to add to ModelInitializer.inputs.
322 Returns
323 -------
324 initializer
325 The configured ModelInitializer.
326 """
327 raise NotImplementedError(f"{self.__name__} must implement __call__")
330class BasicModelInitializerConfig(pexConfig.Config):
331 """Configuration for a BasicModelInitializer."""
333 psf_factor_shrink = pexConfig.Field[float](
334 doc="Multiplicative factor to shrink PSF sizes by for deconvolution",
335 default=0.9,
336 check=lambda x: 0.0 <= x < 1.0,
337 )
338 psf_factor_minimum = pexConfig.Field[float](
339 doc="Factor to multiply the PSF size by for a minimum initialize size",
340 default=0.5,
341 check=lambda x: x >= 0,
342 )
343 size_minimum = pexConfig.Field[float](
344 doc="Absolute minimum initial size in pixels",
345 default=0.5,
346 check=lambda x: x >= 0,
347 )
348 rho_abs_max = pexConfig.Field[float](
349 doc="Maximum absolute initial value of rho",
350 default=0.8,
351 check=lambda x: x >= 0,
352 )
355class BasicModelInitializer(ModelInitializer):
356 """A generic model initializer that should work on most kinds of models
357 with a single source.
358 """
360 config: BasicModelInitializerConfig = pydantic.Field(title="A BasicModelInitializerConfig to be frozen")
362 def _get_params_init(self, model_sources: tuple[g2f.Source]) -> tuple[g2f.ParameterD]:
363 """Return an ordered set of free parameters from a model's sources.
365 Parameters
366 ----------
367 model_sources
368 The sources in the model.
370 Returns
371 -------
372 params_init
373 The parameter objects for sources in the model.
375 Notes
376 -----
377 Only free and/or centroid parameters are returned (centroids are
378 always needed even if they are fixed).
379 """
380 # TODO: There ought to be a better way to not get the PSF centroids
381 # (those are part of model.data's fixed parameters)
382 params_init = (
383 tuple(
384 (
385 param
386 for param in get_params_uniq(model_sources[0])
387 if param.free
388 or (
389 isinstance(param, g2f.CentroidXParameterD)
390 or isinstance(param, g2f.CentroidYParameterD)
391 )
392 )
393 )
394 if (len(model_sources) == 1)
395 else tuple(
396 {
397 param: None
398 for source in model_sources
399 for param in get_params_uniq(source)
400 if param.free
401 or (
402 isinstance(param, g2f.CentroidXParameterD)
403 or isinstance(param, g2f.CentroidYParameterD)
404 )
405 }.keys()
406 )
407 )
408 return params_init
410 def _get_priors_type(
411 self,
412 priors: tuple[g2f.Prior],
413 ) -> tuple[tuple[g2f.GaussianPrior], tuple[g2f.ShapePrior]]:
414 """Return the list of priors of known type, by type.
416 Parameters
417 ----------
418 priors
419 A list of priors of any type, typically from a model.
421 Returns
422 -------
423 priors_gauss
424 A list of all of the Gaussian priors, in the order they occurred.
425 priors_shape
426 A list of all of the shape priors, in the order they occurred.
427 """
428 priors_gauss: list[g2f.GaussianPrior] = []
429 priors_shape: list[g2f.ShapePrior] = []
430 for prior in priors:
431 if isinstance(prior, g2f.GaussianPrior):
432 priors_gauss.append(prior)
433 elif isinstance(prior, g2f.ShapePrior):
434 priors_shape.append(prior)
435 return tuple(priors_gauss), tuple(priors_shape)
437 def get_centroid_and_shape(
438 self,
439 source: Mapping[str, Any],
440 catexps: list[CatalogExposureSourcesABC],
441 config_data: CatalogSourceFitterConfigData,
442 values_init: Mapping[g2f.ParameterD, float] | None = None,
443 ) -> tuple[tuple[float, float], tuple[float, float, float]]:
444 """Get the centroid and shape for a source.
446 Parameters
447 ----------
448 source
449 A mapping with fields expected to be populated in the
450 corresponding source catalog for initialization.
451 catexps
452 A list of (source and psf) catalog-exposure pairs.
453 config_data
454 Configuration settings and data for fitting and output.
455 values_init
456 Initial parameter values from the model configuration.
458 Returns
459 -------
460 centroid
461 The x- and y-axis centroid values.
462 sig_x, sig_y, rho
463 The x- and y-axis Gaussian sigma and rho values defining the
464 estimated elliptical shape of the source.
465 """
466 centroid = source["slot_Centroid_x"], source["slot_Centroid_y"]
467 # Attempt partial deconvolution of observed moments
468 psf_factor_shrink = self.config.psf_factor_shrink**2
469 psf_factor_minimum = self.config.psf_factor_minimum**2
470 rho_min, rho_max = -self.config.rho_abs_max, self.config.rho_abs_max
471 psf_xx = source["base_SdssShape_psf_xx"]
472 psf_yy = source["base_SdssShape_psf_yy"]
473 sig_x, sig_y = (
474 math.sqrt(
475 np.nanmax(
476 (
477 source[f"slot_Shape_{suffix}"] - moment_sq * psf_factor_shrink,
478 moment_sq * psf_factor_minimum,
479 self.config.size_minimum,
480 )
481 )
482 )
483 for suffix, moment_sq in (("xx", psf_xx), ("yy", psf_yy))
484 )
485 psf_xy = source["base_SdssShape_psf_xy"]
486 sig_xy = sig_x * sig_y
487 if not (sig_xy > 0):
488 rho = 0
489 else:
490 rho = np.clip((source["slot_Shape_xy"] - psf_xy * psf_factor_shrink) / sig_xy, rho_min, rho_max)
491 shape = sig_x, sig_y, rho
492 return centroid, shape
494 def get_params_init(self, model: Model) -> tuple[g2f.ParameterD]:
495 """Return the free and/or centroid parameters for a model.
497 Parameters
498 ----------
499 model
500 The model to return parameters for.
502 Returns
503 -------
504 parameters
505 The ordered list of parameters for the model.
506 """
507 return self._get_params_init(model_sources=model.sources)
509 def get_priors_type(self, model: Model) -> tuple[tuple[g2f.GaussianPrior], tuple[g2f.ShapePrior]]:
510 """Return the list of priors of known type, by type.
512 Parameters
513 ----------
514 model
515 The model to return priors for.
517 Returns
518 -------
519 priors_gauss
520 A list of all of the Gaussian priors, in the order they occurred.
521 priors_shape
522 A list of all of the shape priors, in the order they occurred.
523 """
524 return self._get_priors_type(model.priors)
526 def initialize_model(
527 self,
528 model: Model,
529 source: Mapping[str, Any],
530 catexps: list[CatalogExposureSourcesABC],
531 config_data: CatalogSourceFitterConfigData,
532 values_init: Mapping[g2f.ParameterD, float] | None = None,
533 **kwargs,
534 ):
535 if values_init is None:
536 values_init = {}
537 set_flux_limits = kwargs.pop("set_flux_limits", True)
538 flux_init_min = kwargs.pop("value_init_min", 1e-10)
539 flux_limit_min = kwargs.pop("flux_limit_min", 1e-12)
540 if kwargs:
541 raise ValueError(f"Unexpected {kwargs=}")
542 centroid_pixel_offset = config_data.config.centroid_pixel_offset
543 (cen_x, cen_y), (sig_x, sig_y, rho) = self.get_centroid_and_shape(
544 source,
545 catexps,
546 config_data,
547 values_init=values_init,
548 )
549 # If we couldn't get a shape at all, make it small and roundish
550 if not np.isfinite(rho):
551 # Note rho=0 (circular) is generally disfavoured by shape priors
552 # However, setting it to a non-zero value seems to make scipy
553 # fail to move off initial conditions, as do sizes below 2 pixels
554 sig_x, sig_y, rho = 2.0, 2.0, 0.0
556 # Make restrictive centroid limits (intersection, not union)
557 x_min, y_min, x_max, y_max = -np.inf, -np.inf, np.inf, np.inf
559 fluxes_init = {}
560 fluxes_limits = {}
562 # This is the maximum number of potential observations
563 # They might not all have made it into the data
564 n_catexps = len(catexps)
565 n_components = len(model.sources[0].components)
567 # If not true, some bands must have no data to fit
568 if len(catexps) != len(model.data):
569 catexps_obs = []
570 for catexp in catexps:
571 fluxes_init[catexp.channel] = flux_init_min
572 fluxes_limits[catexp.channel] = (0, np.inf)
573 # No associated catalog means we can't fit (and should be
574 # because there's no exposure for this band in this patch)
575 if len(catexp.get_catalog()) > 0:
576 catexps_obs.append(catexp)
577 else:
578 catexps_obs = catexps
580 for idx_obs, observation in enumerate(model.data):
581 coordsys = observation.image.coordsys
582 catexp = catexps_obs[idx_obs]
583 band = catexp.band
585 x_min = max(x_min, coordsys.x_min)
586 y_min = max(y_min, coordsys.y_min)
587 x_max = min(x_max, coordsys.x_min + float(observation.image.n_cols))
588 y_max = min(y_max, coordsys.y_min + float(observation.image.n_rows))
590 flux_total = np.nansum(observation.image.data[observation.mask_inv.data])
592 column_ref = f"merge_measurement_{band}"
593 if column_ref in source.schema.getNames() and source[column_ref]:
594 row = source
595 else:
596 row = catexp.catalog.find(source["id"])
598 if not row["base_SdssShape_flag"]:
599 flux_init = row["base_SdssShape_instFlux"]
600 else:
601 flux_init = row["slot_GaussianFlux_instFlux"]
602 if not (flux_init > 0):
603 flux_init = row["slot_PsfFlux_instFlux"]
605 calib = catexp.exposure.photoCalib
606 flux_init = calib.instFluxToNanojansky(flux_init) if (flux_init > 0) else max(flux_total, 1.0)
607 if set_flux_limits:
608 flux_max = 10 * max((flux_init, flux_total))
609 flux_min = min(flux_limit_min, flux_max / 1000)
610 else:
611 flux_min, flux_max = 0, np.inf
612 if not (flux_init > flux_min):
613 flux_upper = flux_max if (flux_max < np.inf) else 10.0 * flux_min
614 flux_init = flux_min + 0.01 * (flux_upper - flux_min)
615 fluxes_init[observation.channel] = flux_init / n_components
616 fluxes_limits[observation.channel] = (flux_min, flux_max)
618 if not np.isfinite(cen_x):
619 cen_x = observation.image.n_cols / 2.0
620 else:
621 cen_x -= centroid_pixel_offset
622 if not np.isfinite(cen_y):
623 # TODO: Add bbox coords or remove
624 cen_y = observation.image.n_rows / 2.0
625 else:
626 cen_y -= centroid_pixel_offset
628 # An R_eff larger than the box size is problematic. This should also
629 # stop unreasonable size proposals; a log10 transform isn't enough.
630 # TODO: Try logit for r_eff?
631 size_major = g2.EllipseMajor(g2.Ellipse(sigma_x=sig_x, sigma_y=sig_y, rho=rho)).r_major
632 limits_size = max(5.0 * size_major, 2.0 * np.hypot(x_max - x_min, y_max - y_min))
633 limits_xy = (1e-5, limits_size)
634 params_limits_init = {
635 g2f.CentroidXParameterD: (cen_x, (x_min, x_max)),
636 g2f.CentroidYParameterD: (cen_y, (y_min, y_max)),
637 g2f.ReffXParameterD: (sig_x, limits_xy),
638 g2f.ReffYParameterD: (sig_y, limits_xy),
639 g2f.SigmaXParameterD: (sig_x, limits_xy),
640 g2f.SigmaYParameterD: (sig_y, limits_xy),
641 g2f.RhoParameterD: (rho, None),
642 # TODO: get guess from configs?
643 g2f.SersicMixComponentIndexParameterD: (1.0, None),
644 }
646 fluxes_init_tuple = tuple(fluxes_init.values())
647 fluxes_limits_tuple = tuple(fluxes_limits.values())
648 idx_obs = 0
649 for param in self.params_init:
650 if param.linear:
651 value_init = fluxes_init_tuple[idx_obs]
652 limits_new = fluxes_limits_tuple[idx_obs]
653 idx_obs += 1
654 if idx_obs == n_catexps:
655 idx_obs = 0
656 else:
657 type_param = type(param)
658 value_init, limits_new = params_limits_init.get(type_param, (values_init.get(param), None))
659 if limits_new:
660 param.limits = g2f.LimitsD(limits_new[0], limits_new[1])
661 if value_init is not None:
662 param.value = np.clip(value_init, param.limits.min, param.limits.max)
664 priors_shape_mag = self.priors_shape_mag
665 has_priors_mag = len(priors_shape_mag) > 0
666 if has_priors_mag:
667 mag_total = u.nJy.to(u.ABmag, np.nansum(fluxes_init_tuple))
669 # TODO: Add centroid prior
670 priors_gauss, priors_shape = self.get_priors_type(model)
671 for prior in priors_shape:
672 if has_priors_mag and ((prior_adjustments := priors_shape_mag.get(prior)) is not None):
673 mag_dep_prior, prior_shape_new = prior_adjustments
674 prior_size_new = prior_shape_new.prior_size
675 # the size-apparent mag relation probably flattens
676 # for very bright/faint objects - maybe not so
677 # sharply, but clipping a broad mag range ought to be fine
678 prior.prior_size.mean_parameter.value = prior_size_new.mean_parameter.value * 10 ** (
679 mag_dep_prior.slope_median_per_mag
680 * np.clip(
681 mag_total - mag_dep_prior.intercept_mag,
682 -12.5,
683 12.5,
684 )
685 )
686 # it's uncertain how the intrinsic scatter behaves
687 # educated guess is it doesn't change much, also
688 # one runs out of bright galaxies to measure it anyway
689 prior.prior_size.stddev_parameter.value = prior_size_new.stddev_parameter.value * 10 ** (
690 mag_dep_prior.slope_stddev_per_mag
691 * np.clip(
692 mag_total - mag_dep_prior.intercept_mag,
693 -12.5,
694 12.5,
695 )
696 )
697 else:
698 prior.prior_size.mean_parameter.value = size_major
701class CachedBasicModelInitializer(BasicModelInitializer):
702 """A basic initializer with a cached list of model sources and priors."""
704 priors: tuple[g2f.Prior, ...] = pydantic.Field(title="The gauss2d_fit model priors")
705 sources: tuple[g2f.Source, ...] = pydantic.Field(title="The gauss2d_fit model sources")
707 @cached_property
708 def params_init(self) -> tuple[g2f.ParameterD]:
709 """Return a cached reference to the result of _get_params_init."""
710 return self._get_params_init(model_sources=self.sources)
712 @cached_property
713 def priors_type(self) -> tuple[tuple[g2f.GaussianPrior], tuple[g2f.ShapePrior]]:
714 """Return a cached reference to the result of _get_priors_type."""
715 return self._get_priors_type(self.priors)
717 def get_params_init(self, model: Model) -> tuple[g2f.ParameterD]:
718 assert tuple(model.sources) == self.sources
719 return self.params_init
721 def get_priors_type(self, model: Model) -> tuple[tuple[g2f.GaussianPrior], tuple[g2f.ShapePrior]]:
722 assert tuple(model.priors) == self.priors
723 return self.priors_type
726class InitialInputData(pydantic.BaseModel):
727 """A configurable wrapper to retrieve formatted columns from a catalog.
729 This provides a common interface to typical MultiProFit table outputs.
730 """
732 model_config: ClassVar[pydantic.ConfigDict] = frozen_arbitrary_allowed_config
734 column_id: str | None = pydantic.Field(
735 title="Override for id column specified in config_input",
736 default=None,
737 )
738 config_input: InputConfig = pydantic.Field(title="Configuration for the data table")
739 data: Table = pydantic.Field(title="The data table")
740 name_model: str = pydantic.Field(title="The name of the model in columns")
741 prefix_column: str = pydantic.Field(title="The prefix for all fitted column names")
742 size_column: str = pydantic.Field(title="The name of the size column", default="reff")
744 def get_column_id(self):
745 """Return the name of the object ID column."""
746 return self.column_id or self.config_input.column_id
748 def get_column(self, name_column: str, data=None):
749 """Get the values from a column.
751 Parameters
752 ----------
753 name_column
754 The name of the column to retrieve.
755 data
756 The catalog to retrieve the column from. Default is self.data.
758 Returns
759 -------
760 values
761 The column values.
762 """
763 if data is None:
764 data = self.data
765 return data[f"{self.prefix_column}{name_column}"]
767 def model_post_init(self, __context: Any) -> None:
768 # Initialize a mapping of the row number for a given object ID value
769 # This is implemented in afw catalogs but not most other tabular types
770 id_index = {idnum: idx for idx, idnum in enumerate(self.data[self.get_column_id()])}
771 object.__setattr__(self, "id_index", id_index)
774class MakeBasicInitializerAction(MakeInitializerActionBase):
775 """An action to construct an initializer for a single-component,
776 single-source model.
777 """
779 config = pexConfig.ConfigField[BasicModelInitializerConfig](
780 doc="Configuration for the initializer to be constructed",
781 )
783 def _make_initializer(
784 self,
785 catalog_multi: Sequence,
786 catexps: list[fitMB.CatalogExposureInputs],
787 config_data: CatalogSourceFitterConfigData,
788 ) -> ModelInitializer:
789 return BasicModelInitializer(config=self.config)
791 def __call__(
792 self,
793 catalog_multi: Sequence,
794 catexps: list[fitMB.CatalogExposureInputs],
795 config_data: CatalogSourceFitterConfigData,
796 **kwargs,
797 ) -> ModelInitializer:
798 initializer = self._make_initializer(
799 catalog_multi=catalog_multi,
800 catexps=catexps,
801 config_data=config_data,
802 )
803 for name, (config_input, data) in kwargs.items():
804 if not isinstance(data, Table) and hasattr(data, "meta"):
805 _LOG.warning(
806 f"Ignoring extra input {name=} because it is of type {type(data)} and is either not an"
807 f" astropy.table.Table or missing a 'meta' attr"
808 )
809 config_data = data.meta["config"]
810 prefix_column = config_data["prefix_column"]
811 config_source = next(iter(config_data["config_model"]["sources"].values()))
812 config_group = next(iter(config_source["component_groups"].values()))
813 is_sersic = len(config_group["components_sersic"]) > 0
814 name_model = next(
815 iter(config_group["components_sersic"] if is_sersic else config_group["components_gaussian"])
816 )
817 initializer.inputs[name] = InitialInputData(
818 column_id=config_data.get("column_id"),
819 config_input=config_input,
820 data=data,
821 name_model=name_model,
822 prefix_column=prefix_column,
823 size_column="reff" if is_sersic else "sig",
824 )
825 return initializer
828class MakeCachedBasicInitializerAction(MakeBasicInitializerAction):
829 """A MakeBasicInitializerAction that caches references to the source
830 and prior objects of the model.
832 This is solely a performance optimization and should be favored over
833 MakeBasicInitializerAction unless the caching is shown to be slower.
834 """
836 def _make_initializer(
837 self,
838 catalog_multi: Sequence,
839 catexps: list[fitMB.CatalogExposureInputs],
840 config_data: CatalogSourceFitterConfigData,
841 ) -> ModelInitializer:
842 sources, priors = config_data.sources_priors
843 return CachedBasicModelInitializer(config=self.config, priors=priors, sources=sources)
846class MultiProFitSourceConfig(CatalogSourceFitterConfig, fitMB.CoaddMultibandFitSubConfig):
847 """Configuration for the MultiProFit profile fitter."""
849 action_initializer = ConfigurableActionField[MakeInitializerActionBase](
850 doc="The action to return an initializer",
851 default=MakeCachedBasicInitializerAction,
852 )
853 action_psf = ConfigurableActionField[PsfComponentsActionBase](
854 doc="The action to return PSF component values from catalogs, if implemented",
855 default=None,
856 )
857 columns_copy = pexConfig.DictField[str, str](
858 doc="Mapping of input/output column names to copy from the input"
859 "multiband catalog to the output fit catalog.",
860 default={},
861 dictCheck=lambda x: len(set(x.values())) == len(x.values()),
862 )
863 mask_names_zero = pexConfig.ListField[str](
864 doc="Mask bits to mask out",
865 default=["BAD", "EDGE", "SAT", "NO_DATA"],
866 )
867 psf_sigma_subtract = pexConfig.Field[float](
868 doc="PSF x/y sigma value to subtract in quadrature from best-fit values",
869 default=0.1,
870 check=lambda x: np.isfinite(x) and (x >= 0),
871 )
872 prefix_column = pexConfig.Field[str](default="mpf_", doc="Column name prefix")
873 size_priors = pexConfig.ConfigDictField[str, MagnitudeDependentSizePriorConfig](
874 doc="Per-component magnitude-dependent size prior configurations."
875 " Will be added to component with existing configs.",
876 default={},
877 )
879 def bands_read_only(self) -> set[str]:
880 # TODO: Re-implement determination of prior-only bands once
881 # data-driven priors are re-implemented (DM-4xxxx)
882 return set()
884 def requires_psf(self):
885 """Return whether the PSF action is not None."""
886 return type(self.action_psf) is PsfComponentsActionBase
888 def setDefaults(self):
889 super().setDefaults()
890 self.flag_errors = {
891 IsParentError.column_name(): "IsParentError",
892 NoDataError.column_name(): "NoDataError",
893 NotPrimaryError.column_name(): "NotPrimaryError",
894 PsfRebuildFitFlagError.column_name(): "PsfRebuildFitFlagError",
895 }
896 self.centroid_pixel_offset = -0.5
897 self.naming_scheme = "lsst"
898 self.prefix_column = ""
899 self.suffix_error = "Err"
902@pydantic.dataclasses.dataclass(frozen=True, kw_only=True, config=fitMB.CatalogExposureConfig)
903class CatalogExposurePsfs(fitMB.CatalogExposureInputs, CatalogExposureSourcesABC):
904 """Input data from lsst pipelines, parsed for MultiProFit."""
906 channel: g2f.Channel = pydantic.Field(title="Channel for the image's band")
907 config_fit: MultiProFitSourceConfig = pydantic.Field(title="Config for fitting options")
909 @cached_property
910 def _psf_flux_params(self) -> tuple[list[g2f.ParameterD], bool]:
911 psf_model = self.psf_model_data.psf_model
912 n_comps = len(psf_model.components)
913 params_flux = [None] * n_comps
914 is_frac = [False] * n_comps
915 for idx_comp, comp in enumerate(psf_model.components):
916 # TODO: Change to comp.integralmodel when DM-44344 is fixed
917 # integralmodels will still need to be handled differently
918 params_all = get_params_uniq(comp)
919 params_frac = [param for param in params_all if isinstance(param, g2f.ProperFractionParameterD)]
920 if params_frac:
921 is_last = idx_comp == (n_comps - 1)
922 if len(params_frac) != (idx_comp + 1 - is_last):
923 raise RuntimeError(
924 f"Got unexpected {params_frac=} for"
925 f" {self.psf_model_data.psf_model.components[idx_comp]=} ({idx_comp=});"
926 f" len should be idx_comp+1"
927 )
928 params_flux[idx_comp] = None if is_last else params_frac[idx_comp]
929 is_frac[idx_comp] = True
930 else:
931 params_integral = [param for param in params_all if isinstance(param, g2f.IntegralParameterD)]
932 if len(params_integral != 1):
933 raise RuntimeError(
934 f"Got unexpected {params_integral=} != 1 for"
935 f" {self.psf_model_data.psf_model.components[idx_comp]=} ({idx_comp=})"
936 )
937 params_flux[idx_comp] = params_integral[0]
938 is_frac_any = any(is_frac)
939 if is_frac_any and not all(is_frac):
940 # TODO: This should work by iterating through componentgroups
941 # But that's not trivial or supported now
942 raise RuntimeError("Got PSF model with a mix of fractional and linear models; cannot initialize")
944 return params_flux, is_frac_any
946 def get_psf_model(self, source):
947 psf_model = self.psf_model_data.psf_model
948 # PsfComponentsActionBase is an abstract class, so check if the action
949 # is a subclass that needs to be called
950 if not self.config_fit.requires_psf():
951 try:
952 gaussians = self.config_fit.action_psf(source)
953 except PsfRebuildFitFlagError:
954 return None
955 n_comps = len(psf_model.components)
956 fluxes = [0.0] * n_comps
957 params_flux, is_frac = self._psf_flux_params
958 for idx_comp, (comp, gaussian) in enumerate(zip(psf_model.components, gaussians)):
959 ellipse_out = comp.ellipse
960 ellipse_in = gaussian.ellipse
961 ellipse_out.sigma_x = ellipse_in.sigma_x
962 ellipse_out.sigma_y = ellipse_in.sigma_y
963 ellipse_out.rho = ellipse_in.rho
964 fluxes[idx_comp] = gaussian.integral.value
965 # Apparently negative fluxes are possible. Not much can be done to
966 # fix that but set them to a tiny value (zero might work)
967 fluxes = np.clip(fluxes, 1e-3, np.inf)
968 flux_total = sum(fluxes)
969 if is_frac:
970 flux_remaining = 1.0
971 for flux, param_frac in zip(fluxes, params_flux[:-1]):
972 flux_component = flux / flux_total
973 param_frac.value = flux_component / flux_remaining
974 flux_remaining -= flux_component
975 else:
976 for flux, param_flux in zip(fluxes, params_flux):
977 param_flux.value = flux / flux_total
978 else:
979 # TODO: this should probably use .index or something
980 match = np.argwhere(
981 self.table_psf_fits[self.psf_model_data.config.column_id] == source[self.config_fit.column_id]
982 )[0][0]
983 psf_model = self.psf_model_data.psf_model
984 self.psf_model_data.init_psf_model(self.table_psf_fits[match])
986 sigma_subtract = self.config_fit.psf_sigma_subtract
987 if sigma_subtract > 0:
988 sigma_subtract_sq = sigma_subtract * sigma_subtract
989 # 1/10 of PSF sigma should suffice as a minimum size
990 sigma_min_sq = sigma_subtract_sq / 100.0
991 for param in self.psf_model_data.parameters.values():
992 if isinstance(
993 param,
994 g2f.SigmaXParameterD | g2f.SigmaYParameterD | g2f.ReffXParameterD | g2f.ReffYParameterD,
995 ):
996 param.value = math.sqrt(max(param.value**2 - sigma_subtract_sq, sigma_min_sq))
997 return psf_model
999 def get_source_observation(self, source, **kwargs) -> g2f.ObservationD | None:
1000 if not kwargs.get("skip_flags"):
1001 if (not source["detect_isPrimary"]) or source["merge_peak_sky"]:
1002 raise NotPrimaryError(f"source {source[self.config_fit.column_id]} has invalid flags for fit")
1003 footprint = source.getFootprint()
1004 bbox = footprint.getBBox()
1005 if not (bbox.getArea() > 0):
1006 return None
1007 bitmask = 0
1008 mask = self.exposure.mask[bbox]
1009 spans = footprint.spans.asArray()
1010 for bitname in self.config_fit.mask_names_zero:
1011 bitval = mask.getPlaneBitMask(bitname)
1012 bitmask |= bitval
1013 mask = ((mask.array & bitmask) != 0) & (spans != 0)
1014 mask = ~mask
1016 is_deblended_child = source["parent"] != 0
1018 img, _, sigma_inv = get_spanned_image(
1019 exposure=self.exposure,
1020 footprint=footprint if is_deblended_child else None,
1021 bbox=bbox,
1022 spans=spans,
1023 get_sig_inv=True,
1024 )
1025 x_min_bbox, y_min_bbox = bbox.beginX, bbox.beginY
1026 # Crop to tighter box for deblended model if edges are unusable
1027 # ... this rarely ever seems to happen though
1028 if is_deblended_child:
1029 coords = np.argwhere(np.isfinite(img) & (sigma_inv > 0) & np.isfinite(sigma_inv))
1030 if len(coords) == 0:
1031 return None
1032 x_min, y_min = coords.min(axis=0)
1033 x_max, y_max = coords.max(axis=0)
1034 x_max += 1
1035 y_max += 1
1037 if (x_min > 0) or (y_min > 0) or (x_max < img.shape[0]) or (y_max < img.shape[1]):
1038 # Ensure the nominal centroid is still inside the box
1039 # ... although it's a bad sign if that row/column is all bad
1040 x_cen = source["slot_Centroid_x"] - x_min_bbox
1041 y_cen = source["slot_Centroid_y"] - y_min_bbox
1042 x_min = min(x_min, int(np.floor(x_cen)))
1043 x_max = max(x_max, int(np.ceil(x_cen)))
1044 y_min = min(y_min, int(np.floor(y_cen)))
1045 y_max = max(y_max, int(np.ceil(y_cen)))
1046 x_min_bbox += x_min
1047 y_min_bbox += y_min
1048 img = img[x_min:x_max, y_min:y_max]
1049 sigma_inv = sigma_inv[x_min:x_max, y_min:y_max]
1050 mask = mask[x_min:x_max, y_min:y_max]
1052 sigma_inv[~mask] = 0
1054 coordsys = g2.CoordinateSystem(1.0, 1.0, x_min_bbox, y_min_bbox)
1056 obs = g2f.ObservationD(
1057 image=g2.ImageD(img, coordsys),
1058 sigma_inv=g2.ImageD(sigma_inv, coordsys),
1059 mask_inv=g2.ImageB(mask, coordsys),
1060 channel=self.channel,
1061 )
1062 return obs
1064 def __post_init__(self):
1065 # TODO: Can/should this be the derived type (MultiProFitPsfConfig)?
1066 config = CatalogPsfFitterConfig()
1067 config_dict = self.table_psf_fits.meta.get("config")
1068 if config_dict:
1069 set_config_from_dict(config, config_dict)
1070 else:
1071 # TODO: How should this be set?
1072 # If using external PSF fits, it needs to be configured normally
1073 pass
1074 config_data = CatalogPsfFitterConfigData(config=config)
1075 object.__setattr__(self, "psf_model_data", config_data)
1078class MultiProFitSourceFitter(CatalogSourceFitterABC):
1079 """A MultiProFit source fitter.
1081 Parameters
1082 ----------
1083 wcs
1084 A WCS solution that applies to all exposures.
1085 errors_expected
1086 A dictionary of exceptions that are expected to sometimes be raised
1087 during processing (e.g. for missing data) keyed by the name of the
1088 flag column used to record the failure.
1089 add_missing_errors
1090 Whether to add all of the standard MultiProFit errors with default
1091 column names to errors_expected, if not already present.
1092 **kwargs
1093 Keyword arguments to pass to the superclass constructor.
1094 """
1096 initializer: ModelInitializer = pydantic.Field(
1097 title="The model parameter initializer",
1098 default_factory=lambda: BasicModelInitializer(),
1099 )
1100 wcs: lsst.afw.geom.SkyWcs = pydantic.Field(
1101 title="The WCS object to use to convert pixel coordinates to RA/dec",
1102 )
1104 def __init__(
1105 self,
1106 wcs: lsst.afw.geom.SkyWcs,
1107 errors_expected: dict[str, Exception] | None = None,
1108 add_missing_errors: bool = True,
1109 **kwargs: Any,
1110 ):
1111 if errors_expected is None:
1112 errors_expected = {}
1113 if add_missing_errors:
1114 for error_catalog in (IsParentError, NoDataError, NotPrimaryError, PsfRebuildFitFlagError):
1115 if error_catalog not in errors_expected:
1116 errors_expected[error_catalog] = error_catalog.column_name()
1117 super().__init__(wcs=wcs, errors_expected=errors_expected, **kwargs)
1119 def copy_centroid_errors(
1120 self,
1121 columns_cenx_err_copy: tuple[str],
1122 columns_ceny_err_copy: tuple[str],
1123 results: Table,
1124 catalog_multi: Sequence,
1125 catexps: list[CatalogExposureSourcesABC],
1126 config_data: CatalogSourceFitterConfigData,
1127 ):
1128 for column in columns_cenx_err_copy:
1129 results[column] = catalog_multi["slot_Centroid_xErr"]
1130 for column in columns_ceny_err_copy:
1131 results[column] = catalog_multi["slot_Centroid_yErr"]
1133 def get_model_radec(self, source: Mapping[str, Any], cen_x: float, cen_y: float):
1134 # no extra conversions are needed here - cen_x, cen_y are in catalog
1135 # coordinates already
1136 ra, dec = self.wcs.pixelToSky(cen_x, cen_y)
1137 return ra.asDegrees(), dec.asDegrees()
1139 def initialize_model(
1140 self,
1141 model: g2f.ModelD,
1142 source: Mapping[str, Any],
1143 catexps: list[CatalogExposureSourcesABC],
1144 config_data: CatalogSourceFitterConfigData,
1145 values_init: Mapping[g2f.ParameterD, float] | None = None,
1146 **kwargs,
1147 ):
1148 self.initializer.initialize_model(
1149 model=model,
1150 source=source,
1151 catexps=catexps,
1152 config_data=config_data,
1153 values_init=values_init,
1154 **kwargs,
1155 )
1157 def make_CatalogExposurePsfs(
1158 self,
1159 catexp: fitMB.CatalogExposureInputs,
1160 config: MultiProFitSourceConfig,
1161 ) -> CatalogExposurePsfs:
1162 """Make a CatalogExposurePsfs from a list of inputs and a fit config.
1164 Parameters
1165 ----------
1166 catexp
1167 The input catalog-exposure pairs.
1168 config
1169 The MultiProFit source fitting config.
1171 Returns
1172 -------
1173 catexp_psf
1174 The resulting CatalogExposurePsfs.
1175 """
1176 catexp_psf = CatalogExposurePsfs(
1177 # dataclasses.asdict(catexp)_makes a recursive deep copy.
1178 # That must be avoided.
1179 **{key: getattr(catexp, key) for key in catexp.__dataclass_fields__.keys()},
1180 channel=g2f.Channel.get(catexp.band),
1181 config_fit=config,
1182 )
1183 return catexp_psf
1185 def validate_fit_inputs(
1186 self,
1187 catalog_multi: Sequence,
1188 catexps: list[CatalogExposurePsfs],
1189 config_data: CatalogSourceFitterConfigData = None,
1190 logger: logging.Logger = None,
1191 **kwargs: Any,
1192 ) -> None:
1193 errors = []
1194 for idx, catexp in enumerate(catexps):
1195 if not isinstance(catexp, CatalogExposurePsfs):
1196 errors.append(f"catexps[{idx=} {type(catexp)=} !isinstance(CatalogExposurePsfs)")
1197 # Pre-validate the model
1198 config_sources = config_data.config.config_model.sources
1199 model_sources, priors = config_data.sources_priors
1200 priors_shape = [prior for prior in priors if isinstance(prior, g2f.ShapePrior)]
1202 if len(config_sources.keys()) > 1:
1203 errors.append(f"model config has multiple sources: {list(config_sources.keys())}")
1204 elif len(priors_shape) > 0:
1205 idx_prior_found = 0
1206 name_source, config_source = next(iter(config_sources.items()))
1207 source = model_sources[0]
1208 config_groups = config_source.component_groups
1209 if len(config_groups.keys()) > 1:
1210 errors.append(f"model {name_source=} has multiple groups: {list(config_source.keys())}")
1211 else:
1212 name_group, config_group = next(iter(config_groups.items()))
1213 for idx_comp, (name_comp, config_comp) in enumerate(
1214 config_group.get_component_configs().items()
1215 ):
1216 ellipse = source.components[idx_comp].ellipse
1217 # component.ellipse returns a const ref and must be copied
1218 # The ellipse classes might need copy constructors
1219 ellipse_copy = type(ellipse)(
1220 # No kwargs here, since they are unfortunately not
1221 # standardized (e.g. Gaussian is sigma_x not size_x)
1222 # but the arg order is
1223 ellipse.size_x,
1224 ellipse.size_y,
1225 ellipse.rho,
1226 )
1227 prior_shape_new = config_comp.make_shape_prior(ellipse_copy)
1228 if prior_shape_new is not None:
1229 if idx_prior_found == len(priors_shape):
1230 errors.append(
1231 f"Could not validate prior for {name_source=} {name_group=} {name_comp=}"
1232 )
1233 break
1234 prior_shape_old = priors_shape[idx_prior_found]
1235 ll_new, ll_old = (
1236 prior.evaluate().loglike for prior in (prior_shape_new, prior_shape_old)
1237 )
1238 # The necessary tolerance for this check is uncertain
1239 if not np.isclose(ll_new, ll_old):
1240 logger.warning(
1241 f"shape prior for {name_comp=} got inconsistent {ll_new=} vs {ll_old}"
1242 )
1243 if (prior_shape_mod := config_data.config.size_priors.get(name_comp)) is not None:
1244 self.initializer.priors_shape_mag[prior_shape_old] = (
1245 prior_shape_mod,
1246 prior_shape_new,
1247 )
1249 if errors:
1250 raise RuntimeError("\n".join(errors))
1253class MultiProFitSourceTask(fitMB.CoaddMultibandFitSubTask):
1254 """Run MultiProFit on Exposure/SourceCatalog pairs in multiple bands.
1256 This task uses MultiProFit to fit a single model to all sources in a coadd,
1257 using a previously-fit PSF model for each exposure. The task may also use
1258 prior measurements from single- or merged multiband catalogs for
1259 initialization.
1260 """
1262 ConfigClass: ClassVar = MultiProFitSourceConfig
1263 _DefaultName: ClassVar = "multiProFitSource"
1265 def make_default_fitter(
1266 self,
1267 catalog_multi: Sequence,
1268 catexps: list[fitMB.CatalogExposureInputs],
1269 config_data: CatalogSourceFitterConfigData,
1270 **kwargs,
1271 ) -> MultiProFitSourceFitter:
1272 """Make a default MultiProFitSourceFitter.
1274 Parameters
1275 ----------
1276 catalog_multi
1277 A multi-band, indexable source catalog.
1278 catexps
1279 Catalog-exposure-PSF model tuples to fit source models for.
1280 config_data
1281 Configuration and data for the initalizer.
1282 **kwargs
1283 Additional keyword arguments to pass to
1284 self.config.action_initializer.
1286 Returns
1287 -------
1288 fitter
1289 A MultiProFitSourceFitter using the first catexp's wcs.
1290 """
1291 initializer = self.config.action_initializer(
1292 catalog_multi=catalog_multi, catexps=catexps, config_data=config_data, **kwargs
1293 )
1294 # Look for the first WCS - they ought to be identical
1295 # If they are not, the patch coadd data model must have changed
1296 wcs = None
1297 for catexp in catexps:
1298 if catexp.exposure is not None:
1299 wcs = catexp.exposure.wcs
1300 break
1301 if wcs is None:
1302 raise RuntimeError(f"Could not find valid wcs in any of {catexps=}")
1303 fitter = MultiProFitSourceFitter(wcs=wcs, initializer=initializer)
1304 return fitter
1306 @utilsTimer.timeMethod
1307 def run(
1308 self,
1309 catalog_multi: Sequence,
1310 catexps: list[fitMB.CatalogExposureInputs],
1311 fitter: MultiProFitSourceFitter | None = None,
1312 **kwargs,
1313 ) -> pipeBase.Struct:
1314 """Run the MultiProFit source fit task on catalog-exposure pairs.
1316 Parameters
1317 ----------
1318 catalog_multi
1319 A multi-band, indexable source catalog.
1320 catexps
1321 Catalog-exposure-PSF model tuples to fit source models for.
1322 fitter
1323 The fitter instance to use. Default-initialized if not provided.
1324 **kwargs
1325 Additional keyword arguments to pass to self.fit.
1327 Returns
1328 -------
1329 catalog : `astropy.Table`
1330 A table with fit parameters for the PSF model at the location
1331 of each source.
1332 """
1333 n_catexps = len(catexps)
1334 if n_catexps == 0:
1335 raise ValueError("Must provide at least one catexp")
1336 catexps_conv: list[CatalogExposurePsfs] = [None] * n_catexps
1337 channels = [g2f.Channel.get(catexp.band) for catexp in catexps]
1338 config_data = CatalogSourceFitterConfigData(channels=channels, config=self.config)
1339 if fitter is None:
1340 inputs_init = kwargs.get("inputs_init")
1341 if inputs_init:
1342 del kwargs["inputs_init"]
1343 else:
1344 inputs_init = {}
1345 fitter = self.make_default_fitter(
1346 catalog_multi=catalog_multi, catexps=catexps, config_data=config_data, **inputs_init
1347 )
1348 for idx, catexp in enumerate(catexps):
1349 if not isinstance(catexp, CatalogExposurePsfs):
1350 catexp = fitter.make_CatalogExposurePsfs(catexp, config=self.config)
1351 catexps_conv[idx] = catexp
1352 catalog = fitter.fit(
1353 catalog_multi=catalog_multi, catexps=catexps_conv, config_data=config_data, **kwargs
1354 )
1355 for name_in, name_out in self.config.columns_copy.items():
1356 catalog[name_out] = catalog_multi[name_in]
1357 catalog[name_out].description = catalog_multi.schema.find(name_in).field.getDoc()
1358 return pipeBase.Struct(output=astropy_to_arrow(catalog))