Coverage for python / lsst / meas / extensions / multiprofit / fit_coadd_multiband.py: 21%
501 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-21 01:13 -0800
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, params: Mapping[str, Any]) -> g2f.PsfModel | None:
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(params)
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] == params[self.config_fit.column_id]
982 )[0][0]
983 psf_model = self.psf_model_data.psf_model
984 try:
985 self.psf_model_data.init_psf_model(self.table_psf_fits[match])
986 except PsfRebuildFitFlagError:
987 return None
989 sigma_subtract = self.config_fit.psf_sigma_subtract
990 if sigma_subtract > 0:
991 sigma_subtract_sq = sigma_subtract * sigma_subtract
992 # 1/10 of PSF sigma should suffice as a minimum size
993 sigma_min_sq = sigma_subtract_sq / 100.0
994 for param in self.psf_model_data.parameters.values():
995 if isinstance(
996 param,
997 g2f.SigmaXParameterD | g2f.SigmaYParameterD | g2f.ReffXParameterD | g2f.ReffYParameterD,
998 ):
999 param.value = math.sqrt(max(param.value**2 - sigma_subtract_sq, sigma_min_sq))
1000 return psf_model
1002 def get_source_observation(self, source, **kwargs) -> g2f.ObservationD | None:
1003 if not kwargs.get("skip_flags"):
1004 if (not source["detect_isPrimary"]) or source["merge_peak_sky"]:
1005 raise NotPrimaryError(f"source {source[self.config_fit.column_id]} has invalid flags for fit")
1006 footprint = source.getFootprint()
1007 bbox = footprint.getBBox()
1008 if not (bbox.getArea() > 0):
1009 return None
1010 bitmask = 0
1011 mask = self.exposure.mask[bbox]
1012 spans = footprint.spans.asArray()
1013 for bitname in self.config_fit.mask_names_zero:
1014 bitval = mask.getPlaneBitMask(bitname)
1015 bitmask |= bitval
1016 mask = ((mask.array & bitmask) != 0) & (spans != 0)
1017 mask = ~mask
1019 is_deblended_child = source["parent"] != 0
1021 img, _, sigma_inv = get_spanned_image(
1022 exposure=self.exposure,
1023 footprint=footprint if is_deblended_child else None,
1024 bbox=bbox,
1025 spans=spans,
1026 get_sig_inv=True,
1027 )
1028 x_min_bbox, y_min_bbox = bbox.beginX, bbox.beginY
1029 # Crop to tighter box for deblended model if edges are unusable
1030 # ... this rarely ever seems to happen though
1031 if is_deblended_child:
1032 coords = np.argwhere(np.isfinite(img) & (sigma_inv > 0) & np.isfinite(sigma_inv))
1033 if len(coords) == 0:
1034 return None
1035 x_min, y_min = coords.min(axis=0)
1036 x_max, y_max = coords.max(axis=0)
1037 x_max += 1
1038 y_max += 1
1040 if (x_min > 0) or (y_min > 0) or (x_max < img.shape[0]) or (y_max < img.shape[1]):
1041 # Ensure the nominal centroid is still inside the box
1042 # ... although it's a bad sign if that row/column is all bad
1043 x_cen = source["slot_Centroid_x"] - x_min_bbox
1044 y_cen = source["slot_Centroid_y"] - y_min_bbox
1045 x_min = min(x_min, int(np.floor(x_cen)))
1046 x_max = max(x_max, int(np.ceil(x_cen)))
1047 y_min = min(y_min, int(np.floor(y_cen)))
1048 y_max = max(y_max, int(np.ceil(y_cen)))
1049 x_min_bbox += x_min
1050 y_min_bbox += y_min
1051 img = img[x_min:x_max, y_min:y_max]
1052 sigma_inv = sigma_inv[x_min:x_max, y_min:y_max]
1053 mask = mask[x_min:x_max, y_min:y_max]
1055 sigma_inv[~mask] = 0
1057 coordsys = g2.CoordinateSystem(1.0, 1.0, x_min_bbox, y_min_bbox)
1059 obs = g2f.ObservationD(
1060 image=g2.ImageD(img, coordsys),
1061 sigma_inv=g2.ImageD(sigma_inv, coordsys),
1062 mask_inv=g2.ImageB(mask, coordsys),
1063 channel=self.channel,
1064 )
1065 return obs
1067 def __post_init__(self):
1068 # TODO: Can/should this be the derived type (MultiProFitPsfConfig)?
1069 config = CatalogPsfFitterConfig()
1070 config_dict = self.table_psf_fits.meta.get("config")
1071 if config_dict:
1072 set_config_from_dict(config, config_dict)
1073 else:
1074 # TODO: How should this be set?
1075 # If using external PSF fits, it needs to be configured normally
1076 pass
1077 config_data = CatalogPsfFitterConfigData(config=config)
1078 object.__setattr__(self, "psf_model_data", config_data)
1081class MultiProFitSourceFitter(CatalogSourceFitterABC):
1082 """A MultiProFit source fitter.
1084 Parameters
1085 ----------
1086 wcs
1087 A WCS solution that applies to all exposures.
1088 errors_expected
1089 A dictionary of exceptions that are expected to sometimes be raised
1090 during processing (e.g. for missing data) keyed by the name of the
1091 flag column used to record the failure.
1092 add_missing_errors
1093 Whether to add all of the standard MultiProFit errors with default
1094 column names to errors_expected, if not already present.
1095 **kwargs
1096 Keyword arguments to pass to the superclass constructor.
1097 """
1099 initializer: ModelInitializer = pydantic.Field(
1100 title="The model parameter initializer",
1101 default_factory=lambda: BasicModelInitializer(),
1102 )
1103 wcs: lsst.afw.geom.SkyWcs = pydantic.Field(
1104 title="The WCS object to use to convert pixel coordinates to RA/dec",
1105 )
1107 def __init__(
1108 self,
1109 wcs: lsst.afw.geom.SkyWcs,
1110 errors_expected: dict[str, Exception] | None = None,
1111 add_missing_errors: bool = True,
1112 **kwargs: Any,
1113 ):
1114 if errors_expected is None:
1115 errors_expected = {}
1116 if add_missing_errors:
1117 for error_catalog in (IsParentError, NoDataError, NotPrimaryError, PsfRebuildFitFlagError):
1118 if error_catalog not in errors_expected:
1119 errors_expected[error_catalog] = error_catalog.column_name()
1120 super().__init__(wcs=wcs, errors_expected=errors_expected, **kwargs)
1122 def copy_centroid_errors(
1123 self,
1124 columns_cenx_err_copy: tuple[str],
1125 columns_ceny_err_copy: tuple[str],
1126 results: Table,
1127 catalog_multi: Sequence,
1128 catexps: list[CatalogExposureSourcesABC],
1129 config_data: CatalogSourceFitterConfigData,
1130 ):
1131 for column in columns_cenx_err_copy:
1132 results[column] = catalog_multi["slot_Centroid_xErr"]
1133 for column in columns_ceny_err_copy:
1134 results[column] = catalog_multi["slot_Centroid_yErr"]
1136 def get_model_radec(self, source: Mapping[str, Any], cen_x: float, cen_y: float):
1137 # no extra conversions are needed here - cen_x, cen_y are in catalog
1138 # coordinates already
1139 ra, dec = self.wcs.pixelToSky(cen_x, cen_y)
1140 return ra.asDegrees(), dec.asDegrees()
1142 def initialize_model(
1143 self,
1144 model: g2f.ModelD,
1145 source: Mapping[str, Any],
1146 catexps: list[CatalogExposureSourcesABC],
1147 config_data: CatalogSourceFitterConfigData,
1148 values_init: Mapping[g2f.ParameterD, float] | None = None,
1149 **kwargs,
1150 ):
1151 self.initializer.initialize_model(
1152 model=model,
1153 source=source,
1154 catexps=catexps,
1155 config_data=config_data,
1156 values_init=values_init,
1157 **kwargs,
1158 )
1160 def make_CatalogExposurePsfs(
1161 self,
1162 catexp: fitMB.CatalogExposureInputs,
1163 config: MultiProFitSourceConfig,
1164 ) -> CatalogExposurePsfs:
1165 """Make a CatalogExposurePsfs from a list of inputs and a fit config.
1167 Parameters
1168 ----------
1169 catexp
1170 The input catalog-exposure pairs.
1171 config
1172 The MultiProFit source fitting config.
1174 Returns
1175 -------
1176 catexp_psf
1177 The resulting CatalogExposurePsfs.
1178 """
1179 catexp_psf = CatalogExposurePsfs(
1180 # dataclasses.asdict(catexp)_makes a recursive deep copy.
1181 # That must be avoided.
1182 **{key: getattr(catexp, key) for key in catexp.__dataclass_fields__.keys()},
1183 channel=g2f.Channel.get(catexp.band),
1184 config_fit=config,
1185 )
1186 return catexp_psf
1188 def validate_fit_inputs(
1189 self,
1190 catalog_multi: Sequence,
1191 catexps: list[CatalogExposurePsfs],
1192 config_data: CatalogSourceFitterConfigData = None,
1193 logger: logging.Logger = None,
1194 **kwargs: Any,
1195 ) -> None:
1196 errors = []
1197 for idx, catexp in enumerate(catexps):
1198 if not isinstance(catexp, CatalogExposurePsfs):
1199 errors.append(f"catexps[{idx=} {type(catexp)=} !isinstance(CatalogExposurePsfs)")
1200 # Pre-validate the model
1201 config_sources = config_data.config.config_model.sources
1202 model_sources, priors = config_data.sources_priors
1203 priors_shape = [prior for prior in priors if isinstance(prior, g2f.ShapePrior)]
1205 if len(config_sources.keys()) > 1:
1206 errors.append(f"model config has multiple sources: {list(config_sources.keys())}")
1207 elif len(priors_shape) > 0:
1208 idx_prior_found = 0
1209 name_source, config_source = next(iter(config_sources.items()))
1210 source = model_sources[0]
1211 config_groups = config_source.component_groups
1212 if len(config_groups.keys()) > 1:
1213 errors.append(f"model {name_source=} has multiple groups: {list(config_source.keys())}")
1214 else:
1215 name_group, config_group = next(iter(config_groups.items()))
1216 for idx_comp, (name_comp, config_comp) in enumerate(
1217 config_group.get_component_configs().items()
1218 ):
1219 ellipse = source.components[idx_comp].ellipse
1220 # component.ellipse returns a const ref and must be copied
1221 # The ellipse classes might need copy constructors
1222 ellipse_copy = type(ellipse)(
1223 # No kwargs here, since they are unfortunately not
1224 # standardized (e.g. Gaussian is sigma_x not size_x)
1225 # but the arg order is
1226 ellipse.size_x,
1227 ellipse.size_y,
1228 ellipse.rho,
1229 )
1230 prior_shape_new = config_comp.make_shape_prior(ellipse_copy)
1231 if prior_shape_new is not None:
1232 if idx_prior_found == len(priors_shape):
1233 errors.append(
1234 f"Could not validate prior for {name_source=} {name_group=} {name_comp=}"
1235 )
1236 break
1237 prior_shape_old = priors_shape[idx_prior_found]
1238 ll_new, ll_old = (
1239 prior.evaluate().loglike for prior in (prior_shape_new, prior_shape_old)
1240 )
1241 # The necessary tolerance for this check is uncertain
1242 if not np.isclose(ll_new, ll_old):
1243 logger.warning(
1244 f"shape prior for {name_comp=} got inconsistent {ll_new=} vs {ll_old}"
1245 )
1246 if (prior_shape_mod := config_data.config.size_priors.get(name_comp)) is not None:
1247 self.initializer.priors_shape_mag[prior_shape_old] = (
1248 prior_shape_mod,
1249 prior_shape_new,
1250 )
1252 if errors:
1253 raise RuntimeError("\n".join(errors))
1256class MultiProFitSourceTask(fitMB.CoaddMultibandFitSubTask):
1257 """Run MultiProFit on Exposure/SourceCatalog pairs in multiple bands.
1259 This task uses MultiProFit to fit a single model to all sources in a coadd,
1260 using a previously-fit PSF model for each exposure. The task may also use
1261 prior measurements from single- or merged multiband catalogs for
1262 initialization.
1263 """
1265 ConfigClass: ClassVar = MultiProFitSourceConfig
1266 _DefaultName: ClassVar = "multiProFitSource"
1268 def make_default_fitter(
1269 self,
1270 catalog_multi: Sequence,
1271 catexps: list[fitMB.CatalogExposureInputs],
1272 config_data: CatalogSourceFitterConfigData,
1273 **kwargs,
1274 ) -> MultiProFitSourceFitter:
1275 """Make a default MultiProFitSourceFitter.
1277 Parameters
1278 ----------
1279 catalog_multi
1280 A multi-band, indexable source catalog.
1281 catexps
1282 Catalog-exposure-PSF model tuples to fit source models for.
1283 config_data
1284 Configuration and data for the initalizer.
1285 **kwargs
1286 Additional keyword arguments to pass to
1287 self.config.action_initializer.
1289 Returns
1290 -------
1291 fitter
1292 A MultiProFitSourceFitter using the first catexp's wcs.
1293 """
1294 initializer = self.config.action_initializer(
1295 catalog_multi=catalog_multi, catexps=catexps, config_data=config_data, **kwargs
1296 )
1297 # Look for the first WCS - they ought to be identical
1298 # If they are not, the patch coadd data model must have changed
1299 wcs = None
1300 for catexp in catexps:
1301 if catexp.exposure is not None:
1302 wcs = catexp.exposure.wcs
1303 break
1304 if wcs is None:
1305 raise RuntimeError(f"Could not find valid wcs in any of {catexps=}")
1306 fitter = MultiProFitSourceFitter(wcs=wcs, initializer=initializer)
1307 return fitter
1309 @utilsTimer.timeMethod
1310 def run(
1311 self,
1312 catalog_multi: Sequence,
1313 catexps: list[fitMB.CatalogExposureInputs],
1314 fitter: MultiProFitSourceFitter | None = None,
1315 **kwargs,
1316 ) -> pipeBase.Struct:
1317 """Run the MultiProFit source fit task on catalog-exposure pairs.
1319 Parameters
1320 ----------
1321 catalog_multi
1322 A multi-band, indexable source catalog.
1323 catexps
1324 Catalog-exposure-PSF model tuples to fit source models for.
1325 fitter
1326 The fitter instance to use. Default-initialized if not provided.
1327 **kwargs
1328 Additional keyword arguments to pass to self.fit.
1330 Returns
1331 -------
1332 catalog : `astropy.Table`
1333 A table with fit parameters for the PSF model at the location
1334 of each source.
1335 """
1336 n_catexps = len(catexps)
1337 if n_catexps == 0:
1338 raise ValueError("Must provide at least one catexp")
1339 catexps_conv: list[CatalogExposurePsfs] = [None] * n_catexps
1340 channels = [g2f.Channel.get(catexp.band) for catexp in catexps]
1341 config_data = CatalogSourceFitterConfigData(channels=channels, config=self.config)
1342 if fitter is None:
1343 inputs_init = kwargs.get("inputs_init")
1344 if inputs_init:
1345 del kwargs["inputs_init"]
1346 else:
1347 inputs_init = {}
1348 fitter = self.make_default_fitter(
1349 catalog_multi=catalog_multi, catexps=catexps, config_data=config_data, **inputs_init
1350 )
1351 for idx, catexp in enumerate(catexps):
1352 if not isinstance(catexp, CatalogExposurePsfs):
1353 catexp = fitter.make_CatalogExposurePsfs(catexp, config=self.config)
1354 catexps_conv[idx] = catexp
1355 catalog = fitter.fit(
1356 catalog_multi=catalog_multi, catexps=catexps_conv, config_data=config_data, **kwargs
1357 )
1358 for name_in, name_out in self.config.columns_copy.items():
1359 catalog[name_out] = catalog_multi[name_in]
1360 catalog[name_out].description = catalog_multi.schema.find(name_in).field.getDoc()
1361 return pipeBase.Struct(output=astropy_to_arrow(catalog))