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