Coverage for python/lsst/multiprofit/fitting/fit_source.py: 11%
447 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-09 03:36 -0700
1# This file is part of 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/>.
22from abc import ABC, abstractmethod
23from functools import cached_property
24import logging
25import time
26from typing import Any, ClassVar, Iterable, Mapping, Self, Sequence, Type
28import astropy
29from astropy.table import Table
30import astropy.units as u
31import lsst.gauss2d.fit as g2f
32import lsst.pex.config as pexConfig
33from lsst.utils.logging import PeriodicLogger
34import numpy as np
35import pydantic
37from ..componentconfig import Fluxes
38from ..errors import NoDataError
39from ..modelconfig import ModelConfig
40from ..modeller import FitInputsDummy, Modeller
41from ..utils import frozen_arbitrary_allowed_config, get_params_uniq
42from .fit_catalog import CatalogExposureABC, CatalogFitterConfig, ColumnInfo
44__all__ = [
45 "CatalogExposureSourcesABC",
46 "CatalogSourceFitterConfig",
47 "CatalogSourceFitterConfigData",
48 "CatalogSourceFitterABC",
49]
52class CatalogExposureSourcesABC(CatalogExposureABC):
53 """Interface for a CatalogExposure for source modelling."""
55 @property
56 def band(self) -> str:
57 """Return the name of the exposure's passband (e.g. 'r')."""
58 return self.channel.name
60 # Note: not named band because that's usually a string
61 @property
62 @abstractmethod
63 def channel(self) -> g2f.Channel:
64 """Return the exposure's associated channel object."""
66 @abstractmethod
67 def get_psf_model(self, params: Mapping[str, Any]) -> g2f.PsfModel | None:
68 """Get the PSF model for a given source row.
70 Parameters
71 ----------
72 params : Mapping[str, Any]
73 A mapping with parameter values for the best-fit PSF model at the
74 centroid of a single source.
76 Returns
77 -------
78 psf_model : `lsst.gauss2d.fit.PsfModel`
79 A PsfModel object initialized with the best-fit parameters, or None
80 if PSF rebuilding failed for an expected reason (i.e. the input PSF
81 fit table has a flag set).
82 """
84 @abstractmethod
85 def get_source_observation(self, source: Mapping[str, Any], **kwargs: Any) -> g2f.ObservationD | None:
86 """Get the Observation for a given source row.
88 Parameters
89 ----------
90 source : Mapping[str, Any]
91 A mapping with any values needed to retrieve an observation for a
92 single source.
93 **kwargs
94 Additional keyword arguments not used during fitting.
96 Returns
97 -------
98 observation : `lsst.gauss2d.fit.Observation`
99 An Observation object with suitable data for fitting parametric
100 models of the source, or None if the observation cannot be fit.
101 """
104class CatalogSourceFitterConfig(CatalogFitterConfig):
105 """Configuration for the MultiProFit profile fitter."""
107 centroid_pixel_offset = pexConfig.Field[float](
108 doc="Number to add to MultiProFit centroids (bottom-left corner is 0,0) to convert to catalog"
109 " coordinates (e.g. set to -0.5 if the bottom-left corner is -0.5, -0.5)",
110 default=0,
111 )
112 config_model = pexConfig.ConfigField[ModelConfig](doc="Source model configuration")
113 convert_cen_xy_to_radec = pexConfig.Field[bool](default=True, doc="Convert cen x/y params to RA/dec")
114 prior_cen_x_stddev = pexConfig.Field[float](
115 default=0, doc="Prior std. dev. on x centroid (ignored if not >0)"
116 )
117 prior_cen_y_stddev = pexConfig.Field[float](
118 default=0, doc="Prior std. dev. on y centroid (ignored if not >0)"
119 )
120 unit_flux = pexConfig.Field[str](default=None, doc="Flux unit", optional=True)
122 def make_model_data(
123 self,
124 idx_row: int,
125 catexps: list[CatalogExposureSourcesABC],
126 ) -> tuple[g2f.DataD, list[g2f.PsfModel]]:
127 """Make data and psf_models for a catalog row.
129 Parameters
130 ----------
131 idx_row
132 The index of the row in each catalog.
133 catexps
134 Catalog-exposure pairs to initialize observations from.
136 Returns
137 -------
138 data
139 The resulting data object.
140 psf_models
141 A list of psf_models, one per catexp.
143 Notes
144 -----
145 Only observations with good data and valid PSF models will be
146 returned; bad data will be excluded from the return values.
147 """
148 observations = []
149 psf_models = []
151 for catexp in catexps:
152 catalog = catexp.get_catalog()
153 # This indicates that there's no corresponding exposure
154 # (the catexp interface expects a tabular type for catalog but
155 # no interface for an exposure has been defined, yet)
156 if len(catalog) == 0:
157 continue
158 source = catalog[idx_row]
159 observation = catexp.get_source_observation(source)
160 # If the observation or PSF model is bad enough that it cannot be
161 # fit, do not add it to the data.
162 if observation is not None:
163 psf_model = catexp.get_psf_model(source)
164 if psf_model is not None:
165 observations.append(observation)
166 # PSF model parameters cannot be fit along with sources
167 for param in get_params_uniq(psf_model):
168 param.fixed = True
169 psf_models.append(psf_model)
171 data = g2f.DataD(observations)
172 return data, psf_models
174 def make_sources(
175 self,
176 channels: Iterable[g2f.Channel],
177 source_fluxes: list[list[list[Fluxes]]] | None = None,
178 ) -> tuple[list[g2f.Source], list[g2f.Prior]]:
179 """Make initialized sources given channels using `self.config_model`.
181 Parameters
182 ----------
183 channels
184 The channels to initialize fluxes for.
185 source_fluxes
186 A list of fluxes by channel for each component group in each
187 source. The default is to initialize using
188 `ComponentGroupConfig.get_fluxes_default`.
190 Returns
191 -------
192 sources
193 The list of initialized sources.
194 priors
195 The list of priors.
196 """
197 n_sources = len(self.config_model.sources)
198 if source_fluxes is None:
199 source_fluxes = [None] * n_sources
200 for idx, (config_source, component_group_fluxes) in enumerate(
201 zip(
202 self.config_model.sources.values(),
203 source_fluxes,
204 )
205 ):
206 component_group_fluxes = [
207 component_group.get_fluxes_default(
208 channels=channels,
209 component_configs=component_group.get_component_configs(),
210 is_fractional=component_group.is_fractional,
211 )
212 for component_group in config_source.component_groups.values()
213 ]
214 source_fluxes[idx] = component_group_fluxes
215 else:
216 if len(source_fluxes) != n_sources:
217 raise ValueError(f"{len(source_fluxes)=} != {len(self.config_model.sources)=}")
219 sources, priors = self.config_model.make_sources(
220 component_group_fluxes_srcs=source_fluxes,
221 )
223 has_prior_x = self.prior_cen_x_stddev > 0 and np.isfinite(self.prior_cen_x_stddev)
224 has_prior_y = self.prior_cen_y_stddev > 0 and np.isfinite(self.prior_cen_y_stddev)
225 if has_prior_x or has_prior_y:
226 for source in sources:
227 for param in get_params_uniq(source, fixed=False):
228 if has_prior_x and isinstance(param, g2f.CentroidXParameterD):
229 priors.append(g2f.GaussianPrior(param.x_param_ptr, 0, self.prior_cen_x_stddev))
230 elif has_prior_y and isinstance(param, g2f.CentroidYParameterD):
231 priors.append(g2f.GaussianPrior(param.y_param_ptr, 0, self.prior_cen_y_stddev))
233 return sources, priors
235 def schema_configurable(self) -> list[ColumnInfo]:
236 columns = []
237 if self.config_fit.eval_residual:
238 columns.append(ColumnInfo(key="n_eval_jac", dtype="i4"))
239 if self.fit_linear_final:
240 columns.append(ColumnInfo(key="delta_ll_fit_linear", dtype="f8"))
241 return columns
243 def schema(
244 self,
245 bands: list[str] | None = None,
246 ) -> list[ColumnInfo]:
247 if bands is None or not (len(bands) > 0):
248 raise ValueError("CatalogSourceFitter must provide at least one band")
249 schema = super().schema(bands)
251 parameters = CatalogSourceFitterConfigData(
252 config=self,
253 channels=tuple((g2f.Channel.get(band) for band in bands)),
254 ).parameters
255 unit_size = u.Unit("pix")
256 units = {
257 g2f.IntegralParameterD: self.unit_flux,
258 g2f.ReffXParameterD: unit_size,
259 g2f.ReffYParameterD: unit_size,
260 g2f.SizeXParameterD: unit_size,
261 g2f.SizeYParameterD: unit_size,
262 }
263 idx_start = len(schema)
264 schema.extend(
265 [
266 ColumnInfo(key=key, dtype="f8", unit=units.get(type(param)))
267 for key, param in parameters.items()
268 ]
269 )
270 if self.convert_cen_xy_to_radec:
271 label_cen = self.get_key_cen()
272 cen_underscored = label_cen.startswith("_")
273 suffix_x, suffix_y, suffix_ra, suffix_dec = (
274 f"{label_cen}{suffix}"
275 for suffix in (
276 self.get_suffix_x(),
277 self.get_suffix_y(),
278 self.get_suffix_ra(),
279 self.get_suffix_dec(),
280 )
281 )
282 suffix_ra = f"{label_cen}{self.get_suffix_ra()}"
283 suffix_dec = f"{label_cen}{self.get_suffix_dec()}"
284 for key, param in parameters.items():
285 # TODO: Update if allowing x, y <-> dec, RA mappings
286 # ... or arbitrary rotations
287 suffix_radec, suffix_xy = (
288 (suffix_ra, suffix_x)
289 if isinstance(param, g2f.CentroidXParameterD)
290 else (
291 (suffix_dec, suffix_y) if isinstance(param, g2f.CentroidYParameterD) else (None, None)
292 )
293 )
294 if suffix_radec is not None:
295 # Add whatever the corresponding prefix is, and also
296 # remove any leading underscore if there's no prefix
297 prefix, suffix = (
298 ("", suffix_radec[1:])
299 if (cen_underscored and (key == suffix_xy[1:]))
300 else (key.split(suffix_xy)[0], suffix_radec)
301 )
302 schema.append(ColumnInfo(key=f"{prefix}{suffix}", dtype="f8", unit=u.deg))
303 if self.compute_errors != "NONE":
304 suffix = self.suffix_error
305 idx_end = len(schema)
306 # Do not remove idx_end unless you like infinite recursion
307 for column in schema[idx_start:idx_end]:
308 schema.append(ColumnInfo(key=f"{column.key}{suffix}", dtype=column.dtype, unit=column.unit))
310 schema.extend(self.schema_configurable())
311 return schema
314class CatalogSourceFitterConfigData(pydantic.BaseModel):
315 """Configuration data for a fitter that can initialize lsst.gauss2d.fit
316 models and images thereof.
318 This class relies on cached properties being computed once, mostly shortly
319 after initialization. Therefore, it and the config field must be frozen to
320 ensure that the model remains unchanged.
321 """
323 model_config: ClassVar[pydantic.ConfigDict] = frozen_arbitrary_allowed_config
325 channels: list[g2f.Channel] = pydantic.Field(title="The list of channels")
326 config: CatalogSourceFitterConfig = pydantic.Field(title="A CatalogSourceFitterConfig to be frozen")
328 @pydantic.model_validator(mode="after")
329 def validate_config(self) -> Self:
330 self.config.validate()
331 return self
333 @cached_property
334 def components(self) -> tuple[g2f.Component]:
335 sources = self.sources_priors[0]
336 components = []
337 for source in sources:
338 components.extend(source.components)
339 return components
341 @cached_property
342 def parameters(self) -> dict[str, g2f.ParameterD]:
343 config = self.config
344 config_model = config.config_model
345 idx_comp_first = 0
346 has_prefix_source = config_model.has_prefix_source()
347 n_channels = len(self.channels)
348 parameters = {}
350 label_cen = config.get_key_cen()
351 label_rho = config.get_key_rho()
352 label_sersic = config.get_key_sersicindex()
353 label_x, label_y = config.get_suffix_x(), config.get_suffix_y()
355 for name_source, config_source in config_model.sources.items():
356 prefix_source = f"{name_source}_" if has_prefix_source else ""
357 has_prefix_group = config_source.has_prefix_group()
359 for name_group, config_group in config_source.component_groups.items():
360 prefix_group = f"{prefix_source}{name_group}_" if has_prefix_group else prefix_source
361 multicen = len(config_group.centroids) > 1
362 configs_comp = config_group.get_component_configs().items()
364 is_multicomp = len(configs_comp) > 1
366 for idx_comp_group, (name_comp, config_comp) in enumerate(configs_comp):
367 component = self.components[idx_comp_first + idx_comp_group]
369 key_comp = name_comp if is_multicomp else ""
370 prefix_comp = f"{prefix_group}{key_comp}"
371 key_size = config.get_prefixed_label(
372 config.get_key_size(config_comp.get_size_label()),
373 prefix_comp,
374 )
375 key_rho = config.get_prefixed_label(label_rho, prefix_comp)
377 if multicen or (idx_comp_group == 0):
378 prefix_cen = prefix_comp if multicen else prefix_group
379 # Avoid double-underscoring if there's nothing to
380 # prefix or an existing prefix
381 key_cen = config.get_prefixed_label(label_cen, prefix_cen)
382 parameters[f"{key_cen}{label_x}"] = component.centroid.x_param
383 parameters[f"{key_cen}{label_y}"] = component.centroid.y_param
384 if not config_comp.size_x.fixed:
385 parameters[f"{key_size}{label_x}"] = component.ellipse.size_x_param
386 if not config_comp.size_y.fixed:
387 parameters[f"{key_size}{label_y}"] = component.ellipse.size_y_param
388 if not config_comp.rho.fixed:
389 parameters[key_rho] = component.ellipse.rho_param
390 if not config_comp.flux.fixed:
391 # TODO: return this to component.integralmodel
392 # when binding for g2f.FractionalIntegralModel is fixed
393 params_flux = get_params_uniq(component, fixed=False, nonlinear=False)
394 if len(params_flux) != n_channels:
395 raise ValueError(f"{params_flux=} len={len(params_flux)} != {n_channels=}")
396 for channel, param_flux in zip(self.channels, params_flux):
397 key_flux = config.get_key_flux(label=prefix_comp, band=channel.name)
398 parameters[key_flux] = param_flux
399 if hasattr(config_comp, "sersic_index") and not config_comp.sersic_index.fixed:
400 parameters[config.get_prefixed_label(label_sersic, prefix_comp)] = (
401 component.sersicindex_param
402 )
404 return parameters
406 @cached_property
407 def sources_priors(self) -> tuple[tuple[g2f.Source], tuple[g2f.Prior]]:
408 sources, priors = self.config.make_sources(channels=self.channels)
409 return tuple(sources), tuple(priors)
412class CatalogSourceFitterABC(ABC, pydantic.BaseModel):
413 """Fit a Gaussian mixture source model to an image with a PSF model.
415 Notes
416 -----
417 Any exceptions raised and not in errors_expected will be logged in a
418 generic unknown_flag failure column.
419 """
421 model_config: ClassVar[pydantic.ConfigDict] = frozen_arbitrary_allowed_config
423 errors_expected: dict[Type[Exception], str] = pydantic.Field(
424 default_factory=dict,
425 title="A dictionary of Exceptions with the name of the flag column key to fill if raised.",
426 )
427 modeller: Modeller = pydantic.Field(
428 default_factory=Modeller,
429 title="A Modeller instance to use for fitting.",
430 )
432 def _get_columns_params_radec(
433 self,
434 params_radec: dict[str, tuple[g2f.CentroidXParameterD, g2f.CentroidYParameterD]],
435 compute_errors: bool,
436 config: CatalogSourceFitterConfig,
437 ) -> tuple[list[tuple[str, str, str, str]], list[tuple[str, str, str, str, str, str]]]:
438 """Get a list of the columns needed for conversion of x/y centroid
439 parameters into ra/dec.
441 Parameters
442 ----------
443 params_radec
444 Dict of tuple of x, y parameter objects by name.
445 compute_errors
446 Whether errors will be computed.
447 config
448 The configuration with column formatting parameters.
450 Returns
451 -------
452 columns_params_radec
453 Column names for RA, dec, x, and y.
454 columns_params_radec_err
455 Column names for RA_err, dec_err, x, y, x_err, y_err.
456 """
457 columns_params_radec = []
458 columns_params_radec_err = []
459 suffix_err = config.suffix_error
460 key_cen = config.get_key_cen()
461 suffix_x, suffix_y = config.get_suffix_x(), config.get_suffix_y()
462 suffix_ra, suffix_dec = config.get_suffix_ra(), config.get_suffix_dec()
464 for key_base, (param_cen_x, param_cen_y) in params_radec.items():
465 # Avoid C if there's nothing to prefix
466 # or an existing prefix ending in an underscore
467 key_base_cen = config.get_prefixed_label(f"{key_base}{key_cen}", config.prefix_column)
469 if param_cen_y is None:
470 raise RuntimeError(
471 f"Fitter failed to find corresponding cen_y param for {key_base=}; is it fixed?"
472 )
473 columns_params_radec.append(
474 (
475 f"{key_base_cen}{suffix_ra}",
476 f"{key_base_cen}{suffix_dec}",
477 f"{key_base_cen}{suffix_x}",
478 f"{key_base_cen}{suffix_y}",
479 )
480 )
481 if compute_errors:
482 columns_params_radec_err.append(
483 (
484 f"{key_base_cen}{suffix_ra}{suffix_err}",
485 f"{key_base_cen}{suffix_dec}{suffix_err}",
486 f"{key_base_cen}{suffix_x}",
487 f"{key_base_cen}{suffix_y}",
488 f"{key_base_cen}{suffix_x}{suffix_err}",
489 f"{key_base_cen}{suffix_y}{suffix_err}",
490 )
491 )
492 return columns_params_radec, columns_params_radec_err
494 @staticmethod
495 def _get_logger() -> logging.Logger:
496 logger = logging.getLogger(__name__)
498 return logger
500 def _validate_errors_expected(self, config: CatalogSourceFitterConfig) -> None:
501 """Check that self.errors_expected is set correctly.
503 Parameters
504 ----------
505 config
506 The fitting configuration.
508 Raises
509 ------
510 ValueError
511 Raised if the configuration is invalid.
512 """
513 if len(self.errors_expected) != len(config.flag_errors):
514 raise ValueError(f"{self.errors_expected=} keys not same len as {config.flag_errors=}")
515 errors_bad = {}
516 errors_recast = {}
517 for error_name, error_type in self.errors_expected.items():
518 if error_type in errors_recast:
519 errors_bad[error_name] = error_type
520 else:
521 errors_recast[error_type] = error_name
522 if errors_bad:
523 raise ValueError(f"{self.errors_expected=} keys contain duplicates from {config.flag_errors=}")
525 def copy_centroid_errors(
526 self,
527 columns_cenx_err_copy: tuple[str],
528 columns_ceny_err_copy: tuple[str],
529 results: Table,
530 catalog_multi: Sequence,
531 catexps: list[CatalogExposureSourcesABC],
532 config_data: CatalogSourceFitterConfigData,
533 ) -> None:
534 """Copy centroid errors from an input catalog.
536 This method exists to support fitting models with fixed centroids
537 derived from an input catalog. Implementers can simply copy an
538 existing column into the results catalog or use the data as needed;
539 however, there is no reasonable default implementation.
541 Parameters
542 ----------
543 columns_cenx_err_copy
544 X-axis result centroid columns to copy errors for.
545 columns_ceny_err_copy
546 Y-axis result centroid columns to copy errors for.
547 results
548 The table of fit results to copy errors into.
549 catalog_multi
550 The input multiband catalog.
551 catexps
552 The input data.
553 config_data
554 The fitter config and data.
556 Raises
557 ------
558 NotImplementedError
559 Raised if columns need to be copied but no implementation is
560 available.
561 """
562 if columns_cenx_err_copy or columns_ceny_err_copy:
563 raise NotImplementedError(
564 f"Fitter of {type(self)=} got {columns_cenx_err_copy=} and/or {columns_ceny_err_copy=}"
565 f" but has not overriden copy_centroid_errors"
566 )
568 def fit(
569 self,
570 catalog_multi: Sequence,
571 catexps: list[CatalogExposureSourcesABC],
572 config_data: CatalogSourceFitterConfigData | None = None,
573 logger: logging.Logger | None = None,
574 **kwargs: Any,
575 ) -> astropy.table.Table:
576 """Fit PSF-convolved source models with MultiProFit.
578 Each source has a single PSF-convolved model fit, given PSF model
579 parameters from a catalog, and a combination of initial source
580 model parameters and a deconvolved source image from the
581 CatalogExposureSources.
583 Parameters
584 ----------
585 catalog_multi
586 A multi-band source catalog to fit a model to.
587 catexps
588 A list of (source and psf) catalog-exposure pairs.
589 config_data
590 Configuration settings and data for fitting and output.
591 logger
592 The logger. Defaults to calling `_getlogger`.
593 **kwargs
594 Additional keyword arguments to pass to self.modeller.
596 Returns
597 -------
598 catalog : `astropy.Table`
599 A table with fit parameters for the PSF model at the location
600 of each source.
601 """
602 if config_data is None:
603 config_data = CatalogSourceFitterConfigData(
604 config=CatalogSourceFitterConfig(),
605 channels=[catexp.channel for catexp in catexps],
606 )
607 if logger is None:
608 logger = self._get_logger()
610 config = config_data.config
611 self._validate_errors_expected(config)
612 self.validate_fit_inputs(
613 catalog_multi=catalog_multi, catexps=catexps, config_data=config_data, logger=logger, **kwargs
614 )
616 model_sources, priors = config_data.sources_priors
617 # TODO: If free Observation params are ever supported, make null Data
618 # Because config_data knows nothing about the Observation(s)
619 params = config_data.parameters
620 values_init = {param: param.value for param in params.values() if param.free}
621 prefix = config.prefix_column
622 columns_param_fixed: dict[str, tuple[g2f.ParameterD, float]] = {}
623 columns_param_free: dict[str, tuple[g2f.ParameterD, float]] = {}
624 columns_param_flux: dict[str, g2f.IntegralParameterD] = {}
625 params_cen_x: dict[str, g2f.CentroidYParameterD] = {}
626 params_cen_y: dict[str, g2f.CentroidYParameterD] = {}
627 columns_err = []
629 errors_hessian: bool = config.compute_errors == "INV_HESSIAN"
630 errors_hessian_bestfit: bool = config.compute_errors == "INV_HESSIAN_BESTFIT"
631 compute_errors: bool = errors_hessian or errors_hessian_bestfit
633 columns_cenx_err_copy = []
634 columns_ceny_err_copy = []
636 suffix_err = config.suffix_error
637 key_cen = config.get_key_cen()
638 cen_underscored = key_cen.startswith("_")
639 suffix_cenx = f"{key_cen}{config.get_suffix_x()}"
640 suffix_ceny = f"{key_cen}{config.get_suffix_y()}"
642 # Add each param to appropriate and more specific pre-computed lists
643 for key, param in params.items():
644 key_full = f"{prefix}{key}"
645 is_cenx = isinstance(param, g2f.CentroidXParameterD)
646 is_ceny = isinstance(param, g2f.CentroidYParameterD)
648 # Add the corresponding error key to the appropriate list
649 if compute_errors:
650 if param.free:
651 columns_err.append(f"{key_full}{suffix_err}")
652 elif is_cenx:
653 columns_cenx_err_copy.append(f"{key_full}{suffix_err}")
654 elif is_ceny:
655 columns_ceny_err_copy.append(f"{key_full}{suffix_err}")
657 # Add this param to the appropriate dict
658 (columns_param_fixed if param.fixed else columns_param_free)[key_full] = (
659 param,
660 config_data.config.centroid_pixel_offset if (is_cenx or is_ceny) else 0,
661 )
662 if isinstance(param, g2f.IntegralParameterD):
663 columns_param_flux[key_full] = param
664 elif config.convert_cen_xy_to_radec:
665 # Infer the prefix if possible, after checking for a dropped
666 # leading underscore in case there's no prefix
667 if is_cenx:
668 prefix_cen, suffix_cen = (
669 ("", key_full)
670 if (cen_underscored and (key_full == suffix_cenx[1:]))
671 else key_full.split(suffix_cenx)
672 )
673 params_cen_x[prefix_cen] = param
674 elif is_ceny:
675 prefix_cen, suffix_cen = (
676 ("", key_full)
677 if (cen_underscored and (key_full == suffix_ceny[1:]))
678 else key_full.split(suffix_ceny)
679 )
680 params_cen_y[prefix_cen] = param
682 if config.convert_cen_xy_to_radec:
683 assert params_cen_x.keys() == params_cen_y.keys()
684 columns_params_radec, columns_params_radec_err = self._get_columns_params_radec(
685 {k: (x, params_cen_y[k]) for k, x in params_cen_x.items()},
686 compute_errors,
687 config=config,
688 )
690 # Setup the results table with correct column names
691 n_rows = len(catalog_multi)
692 channels = self.get_channels(catexps)
693 results, columns = config.make_catalog(n_rows, bands=list(channels.keys()))
695 # Copy centroid error columns into results ( if needed)
696 self.copy_centroid_errors(
697 columns_cenx_err_copy=columns_cenx_err_copy,
698 columns_ceny_err_copy=columns_ceny_err_copy,
699 results=results,
700 catalog_multi=catalog_multi,
701 catexps=catexps,
702 config_data=config_data,
703 )
705 # dummy size for first iteration
706 size, size_new = 0, 0
707 fitInputs = FitInputsDummy()
708 plot = False
710 # Configure default options for calls to compute_variances
711 # keys are for values of return_negative
712 kwargs_err_default = {
713 True: {
714 "options": g2f.HessianOptions(findiff_add=1e-3, findiff_frac=1e-3),
715 "use_diag_only": config.compute_errors_no_covar,
716 },
717 False: {"options": g2f.HessianOptions(findiff_add=1e-6, findiff_frac=1e-6)},
718 }
720 range_idx = range(n_rows)
722 # TODO: Do this check with dummy data
723 # It might not work with real data if the first row is bad
724 # data, psf_models = config.make_model_data(
725 # idx_row=range_idx[0], catexps=catexps)
726 # model = g2f.ModelD(data=data, psfmodels=psf_models,
727 # sources=model_sources, priors=priors)
728 # Remember to filter out fixed centroids from params
729 # assert list(params.values()) == get_params_uniq(model, fixed=False)
731 time_init_all = time.process_time()
732 logger_periodic = PeriodicLogger(logger)
733 n_skipfail = 0
735 for idx in range_idx:
736 time_init = time.process_time()
737 row = results[idx]
738 source_multi = catalog_multi[idx]
739 id_source = source_multi[config.column_id]
740 row[config.column_id] = id_source
741 time_final = time_init
743 try:
744 data, psf_models = config.make_model_data(idx_row=idx, catexps=catexps)
745 if data.size == 0:
746 raise NoDataError("make_model_data returned empty data")
747 model = g2f.ModelD(data=data, psfmodels=psf_models, sources=model_sources, priors=priors)
748 self.initialize_model(
749 model,
750 source_multi,
751 catexps,
752 config_data=config_data,
753 values_init=values_init,
754 )
756 # Caches the jacobian residual if the data size is unchanged
757 # Note: this will need to change with priors
758 # (data should report its own size)
759 size_new = np.sum([datum.image.size for datum in data])
760 if size_new != size:
761 fitInputs = None
762 size = size_new
763 # Some algorithms might not even use fitInputs
764 elif fitInputs is not None:
765 fitInputs = fitInputs if not fitInputs.validate_for_model(model) else None
767 # TODO: Check if flux param limits and transforms are set
768 # appropriately if config.fit_linear_init is False
769 if config.fit_linear_init:
770 self.modeller.fit_model_linear(model=model, ratio_min=0.01)
772 for observation in data:
773 observation.image.data[~np.isfinite(observation.image.data)] = 0
775 result_full = self.modeller.fit_model(
776 model, fitinputs=fitInputs, config=config.config_fit, **kwargs
777 )
778 fitInputs = result_full.inputs
779 results[f"{prefix}n_iter"][idx] = result_full.n_eval_func
780 results[f"{prefix}time_eval"][idx] = result_full.time_eval
781 results[f"{prefix}time_fit"][idx] = result_full.time_run
782 if config.config_fit.eval_residual:
783 results[f"{prefix}n_eval_jac"][idx] = result_full.n_eval_jac
785 params_free_missing = result_full.params_free_missing or tuple()
787 # Set all params to best fit values
788 # In case the optimizer doesn't
789 for (key, (param, offset)), value in zip(
790 columns_param_free.items(),
791 result_full.params_best,
792 ):
793 param.value_transformed = value
794 if param not in params_free_missing:
795 results[key][idx] = param.value + offset
797 # Also add any offset to the fixed parameters
798 # (usually centroids, if any)
799 for key, (param, offset) in columns_param_fixed.items():
800 results[key][idx] = param.value + offset
802 # Do a final linear fit
803 # If the nonlinear fit is good, the values won't change much
804 if config.fit_linear_final:
805 loglike_init, loglike_new = self.modeller.fit_model_linear(
806 model=model, ratio_min=0.01, validate=True
807 )
808 results[f"{prefix}delta_ll_fit_linear"][idx] = np.sum(loglike_new) - np.sum(loglike_init)
810 if params_free_missing:
811 columns_param_flux_fit = {
812 column: param
813 for column, param in columns_param_flux.items()
814 if param not in params_free_missing
815 }
816 else:
817 columns_param_flux_fit = columns_param_flux
819 for column, param in columns_param_flux_fit.items():
820 results[column][idx] = param.value
822 if config.convert_cen_xy_to_radec:
823 for key_ra, key_dec, key_cen_x, key_cen_y in columns_params_radec:
824 # These will have been converted back if necessary
825 cen_x, cen_y = results[key_cen_x][idx], results[key_cen_y][idx]
826 radec = self.get_model_radec(source_multi, cen_x, cen_y)
827 results[key_ra][idx], results[key_dec][idx] = radec
829 if compute_errors:
830 errors = []
831 model_eval = model
832 errors_iter = None
833 for param in params_free_missing:
834 param.fixed = True
836 if config.compute_errors_from_jacobian:
837 try:
838 errors_iter = np.sqrt(
839 self.modeller.compute_variances(
840 model_eval,
841 transformed=False,
842 use_diag_only=config.compute_errors_no_covar,
843 )
844 )
845 errors.append((errors_iter, np.sum(~(errors_iter > 0))))
846 except Exception:
847 pass
848 # If computing errors from the Jacobian didn't work, or if
849 # it was disabled in the config, try the Hessian
850 if errors_iter is None:
851 img_data_old = []
852 if errors_hessian_bestfit:
853 # Model sans prior
854 model_eval = g2f.ModelD(
855 data=model.data, psfmodels=model.psfmodels, sources=model.sources
856 )
857 model_eval.setup_evaluators(evaluatormode=g2f.EvaluatorMode.image)
858 model_eval.evaluate()
859 # Compute the errors by setting the data to the
860 # best-fit model (a quasi-parametric bootstrap
861 # with one iteration)
862 for obs, output in zip(model_eval.data, model_eval.outputs):
863 img_data_old.append(obs.image.data.copy())
864 img = obs.image.data
865 img.flat = output.data.flat
866 # To make this a real bootstrap, could do this
867 # (but would need to iterate):
868 # + rng.standard_normal(img.size)*(
869 # obs.sigma_inv.data.flat)
871 # Try without forcing all of the Hessian terms to be
872 # negative first. At the optimum they should be, but
873 # in practice the best-fit values are always at least
874 # a little off and so the sign is equally likely to be
875 # positive as negative.
876 for return_negative in (False, True):
877 kwargs_err = kwargs_err_default[return_negative]
878 if errors and errors[-1][1] == 0:
879 break
880 try:
881 errors_iter = np.sqrt(
882 self.modeller.compute_variances(
883 model_eval, transformed=False, **kwargs_err
884 )
885 )
886 errors.append((errors_iter, np.sum(~(errors_iter > 0))))
887 except Exception:
888 try:
889 errors_iter = np.sqrt(
890 self.modeller.compute_variances(
891 model_eval,
892 transformed=False,
893 use_svd=True,
894 **kwargs_err,
895 )
896 )
897 errors.append((errors_iter, np.sum(~(errors_iter > 0))))
898 except Exception:
899 pass
900 # Return the data to its original noisy values
901 # (it was replaced by the model earlier)
902 if errors_hessian_bestfit:
903 for obs, img_datum_old in zip(model.data, img_data_old):
904 obs.image.data.flat = img_datum_old.flat
905 # Save and optionally plot the errors
906 if errors:
907 idx_min = np.argmax([err[1] for err in errors])
908 errors = errors[idx_min][0]
909 if plot:
910 errors_plot = np.clip(errors, 0, 1000)
911 errors_plot[~np.isfinite(errors_plot)] = 0
912 from ..plotting import ErrorValues, plot_loglike
914 try:
915 plot_loglike(model, errors={"err": ErrorValues(values=errors_plot)})
916 except Exception:
917 for param in params:
918 param.fixed = False
920 if params_free_missing:
921 columns_err_fitted = [
922 column
923 for column, param in zip(columns_err, params.values())
924 if param not in params_free_missing
925 ]
926 else:
927 columns_err_fitted = columns_err
929 for value, column_err in zip(errors, columns_err_fitted):
930 results[column_err][idx] = value
932 for param in params_free_missing:
933 param.fixed = False
935 # Convert the x/y errors to ra/dec errors
936 if config.convert_cen_xy_to_radec:
937 for (
938 key_ra_err,
939 key_dec_err,
940 key_cen_x,
941 key_cen_y,
942 key_cen_x_err,
943 key_cen_y_err,
944 ) in columns_params_radec_err:
945 cen_x, cen_y = results[key_cen_x][idx], results[key_cen_y][idx]
946 # TODO: improve this in DM-45682
947 # For one, it won't work right at limits:
948 # RA=359.99... or dec=+89.99...
949 # Could also consider dividing by sqrt(2)
950 # ...but that factor would multiply out later
951 radec_err = np.array(
952 self.get_model_radec(
953 source_multi,
954 cen_x + results[key_cen_x_err][idx],
955 cen_y + results[key_cen_y_err][idx],
956 )
957 )
958 radec_err -= radec
959 results[key_ra_err][idx], results[key_dec_err][idx] = np.abs(radec_err)
961 results[f"{prefix}chisq_reduced"][idx] = result_full.chisq_best / size
962 time_final = time.process_time()
963 results[f"{prefix}time_full"][idx] = time_final - time_init
964 except Exception as e:
965 n_skipfail += 1
966 size = 0 if fitInputs is None else size_new
967 column = self.errors_expected.get(e.__class__, "")
968 if column:
969 row[f"{prefix}{column}"] = True
970 logger.debug(
971 "id_source=%i (idx=%i/%i) fit failed with known exception: %s",
972 id_source,
973 idx,
974 n_rows,
975 e,
976 )
977 else:
978 row[f"{prefix}unknown_flag"] = True
979 logger.info(
980 "id_source=%i (idx=%i/%i) fit failed with unexpected exception: %s",
981 id_source,
982 idx,
983 n_rows,
984 e,
985 exc_info=1,
986 )
987 logger_periodic.log(
988 "Fit idx=%i/%i sources (%i skipped/failed) in %.2f",
989 idx,
990 n_rows,
991 n_skipfail,
992 time_final - time_init_all,
993 )
995 n_unknown = np.sum(row[f"{prefix}unknown_flag"])
996 if n_unknown > 0:
997 logger.warning("%i/%i source fits failed with unexpected exceptions", n_unknown, n_rows)
999 return results
1001 def get_channels(
1002 self,
1003 catexps: list[CatalogExposureSourcesABC],
1004 ) -> dict[str, g2f.Channel]:
1005 channels = {}
1006 for catexp in catexps:
1007 try:
1008 channel = catexp.channel
1009 except AttributeError:
1010 band = catexp.band
1011 if callable(band):
1012 band = band()
1013 channel = g2f.Channel.get(band)
1014 if channel not in channels:
1015 channels[channel.name] = channel
1016 return channels
1018 def get_model(
1019 self,
1020 idx_row: int,
1021 catalog_multi: Sequence,
1022 catexps: list[CatalogExposureSourcesABC],
1023 config_data: CatalogSourceFitterConfigData | None = None,
1024 results: astropy.table.Table | None = None,
1025 **kwargs: Any,
1026 ) -> g2f.ModelD:
1027 """Reconstruct the model for a single row of a fit catalog.
1029 Parameters
1030 ----------
1031 idx_row
1032 The index of the row in the catalog.
1033 catalog_multi
1034 The multi-band catalog originally used for initialization.
1035 catexps
1036 The catalog-exposure pairs to reconstruct the model for.
1037 config_data
1038 The configuration used to generate sources.
1039 Default-initialized if None.
1040 results
1041 The corresponding best-fit parameter catalog to initialize
1042 parameter values from. If None, the model params will be set by
1043 `self.initialize_model`, as they would be when calling `self.fit`.
1044 **kwargs
1045 Additional keyword arguments to pass to initialize_model. Not
1046 used during fitting.
1048 Returns
1049 -------
1050 model
1051 The reconstructed model.
1052 """
1053 channels = self.get_channels(catexps)
1054 if config_data is None:
1055 config_data = CatalogSourceFitterConfigData(
1056 config=CatalogSourceFitterConfig(),
1057 channels=list(channels.values()),
1058 )
1059 config = config_data.config
1061 if not idx_row >= 0:
1062 raise ValueError(f"{idx_row=} !>=0")
1063 if not len(catalog_multi) > idx_row:
1064 raise ValueError(f"{len(catalog_multi)=} !> {idx_row=}")
1065 if (results is not None) and not (len(results) > idx_row):
1066 raise ValueError(f"{len(results)=} !> {idx_row=}")
1068 model_sources, priors = config_data.sources_priors
1069 source_multi = catalog_multi[idx_row]
1071 data, psf_models = config.make_model_data(
1072 idx_row=idx_row,
1073 catexps=catexps,
1074 )
1075 model = g2f.ModelD(data=data, psfmodels=psf_models, sources=model_sources, priors=priors)
1076 self.initialize_model(model, source_multi, catexps, **kwargs)
1078 if results is not None:
1079 row = results[idx_row]
1080 for column, param in config_data.parameters.items():
1081 param.value = row[f"{config.prefix_column}{column}"]
1083 return model
1085 def get_model_radec(self, source: Mapping[str, Any], cen_x: float, cen_y: float) -> tuple[float, float]:
1086 """Return right ascension and declination values for a source.
1088 Implementing this method is necessary only when fitting data with
1089 accompanying WCS.
1091 Parameters
1092 ----------
1093 source
1094 A mapping with fields expected to be populated in the
1095 corresponding source catalog.
1096 cen_x
1097 The x-axis centroid in pixel coordinates.
1098 cen_y
1099 The y-axis centroid in pixel coordinates.
1101 Returns
1102 -------
1103 ra, dec
1104 The right ascension and declination.
1105 """
1106 return np.nan, np.nan
1108 @abstractmethod
1109 def initialize_model(
1110 self,
1111 model: g2f.ModelD,
1112 source: Mapping[str, Any],
1113 catexps: list[CatalogExposureSourcesABC],
1114 config_data: CatalogSourceFitterConfigData,
1115 values_init: Mapping[g2f.ParameterD, float] | None = None,
1116 **kwargs: Any,
1117 ) -> None:
1118 """Initialize a Model for a single source row.
1120 Parameters
1121 ----------
1122 model
1123 The model object to initialize.
1124 source
1125 A mapping with fields expected to be populated in the
1126 corresponding source catalog for initialization.
1127 catexps
1128 A list of (source and psf) catalog-exposure pairs.
1129 config_data
1130 Configuration settings and data for fitting and output.
1131 values_init
1132 Initial parameter values from the model configuration.
1133 **kwargs
1134 Additional keyword arguments that cannot be required for fitting.
1135 """
1137 @abstractmethod
1138 def validate_fit_inputs(
1139 self,
1140 catalog_multi: Sequence,
1141 catexps: list[CatalogExposureSourcesABC],
1142 config_data: CatalogSourceFitterConfigData = None,
1143 logger: logging.Logger = None,
1144 **kwargs: Any,
1145 ) -> None:
1146 """Validate inputs to self.fit.
1148 This method is called before any fitting is done. It may be used for
1149 any purpose, including checking that the inputs are a particular
1150 subclass of the base classes.
1152 Parameters
1153 ----------
1154 catalog_multi
1155 A multi-band source catalog to fit a model to.
1156 catexps
1157 A list of (source and psf) catalog-exposure pairs.
1158 config_data
1159 Configuration settings and data for fitting and output.
1160 logger
1161 The logger. Defaults to calling `_getlogger`.
1162 **kwargs
1163 Additional keyword arguments to pass to self.modeller.
1164 """
1165 pass