Coverage for python/lsst/multiprofit/fitting/fit_source.py: 10%
490 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-07 02:01 -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 results[f"{prefix}delta_lnL_fit_ps"][idx] = loglike_final[0] - loglike_psfmodel[0]
913 if compute_errors:
914 errors = []
915 model_eval = model
916 errors_iter = None
917 for param in params_free_missing:
918 param.fixed = True
920 if config.compute_errors_from_jacobian:
921 try:
922 errors_iter = np.sqrt(
923 self.modeller.compute_variances(
924 model_eval,
925 transformed=False,
926 use_diag_only=config.compute_errors_no_covar,
927 )
928 )
929 errors.append((errors_iter, np.sum(~(errors_iter > 0))))
930 except Exception:
931 pass
932 # If computing errors from the Jacobian didn't work, or if
933 # it was disabled in the config, try the Hessian
934 if errors_iter is None:
935 img_data_old = []
936 if errors_hessian_bestfit:
937 # Model sans prior
938 model_eval = g2f.ModelD(
939 data=model.data, psfmodels=model.psfmodels, sources=model.sources
940 )
941 model_eval.setup_evaluators(evaluatormode=g2f.EvaluatorMode.image)
942 model_eval.evaluate()
943 # Compute the errors by setting the data to the
944 # best-fit model (a quasi-parametric bootstrap
945 # with one iteration)
946 for obs, output in zip(model_eval.data, model_eval.outputs):
947 img_data_old.append(obs.image.data.copy())
948 img = obs.image.data
949 img.flat = output.data.flat
950 # To make this a real bootstrap, could do this
951 # (but would need to iterate):
952 # + rng.standard_normal(img.size)*(
953 # obs.sigma_inv.data.flat)
955 # Try without forcing all of the Hessian terms to be
956 # negative first. At the optimum they should be, but
957 # in practice the best-fit values are always at least
958 # a little off and so the sign is equally likely to be
959 # positive as negative.
960 for return_negative in (False, True):
961 kwargs_err = kwargs_err_default[return_negative]
962 if errors and errors[-1][1] == 0:
963 break
964 try:
965 errors_iter = np.sqrt(
966 self.modeller.compute_variances(
967 model_eval, transformed=False, **kwargs_err
968 )
969 )
970 errors.append((errors_iter, np.sum(~(errors_iter > 0))))
971 except Exception:
972 try:
973 errors_iter = np.sqrt(
974 self.modeller.compute_variances(
975 model_eval,
976 transformed=False,
977 use_svd=True,
978 **kwargs_err,
979 )
980 )
981 errors.append((errors_iter, np.sum(~(errors_iter > 0))))
982 except Exception:
983 pass
984 # Return the data to its original noisy values
985 # (it was replaced by the model earlier)
986 if errors_hessian_bestfit:
987 for obs, img_datum_old in zip(model.data, img_data_old):
988 obs.image.data.flat = img_datum_old.flat
989 # Save and optionally plot the errors
990 if errors:
991 idx_min = np.argmax([err[1] for err in errors])
992 errors = errors[idx_min][0]
993 if plot:
994 errors_plot = np.clip(errors, 0, 1000)
995 errors_plot[~np.isfinite(errors_plot)] = 0
996 from ..plotting import ErrorValues, plot_loglike
998 try:
999 plot_loglike(model, errors={"err": ErrorValues(values=errors_plot)})
1000 except Exception:
1001 for param in params:
1002 param.fixed = False
1004 if params_free_missing:
1005 columns_err_fitted = [
1006 column
1007 for column, param in zip(columns_err, params.values())
1008 if param not in params_free_missing
1009 ]
1010 else:
1011 columns_err_fitted = columns_err
1013 for value, column_err in zip(errors, columns_err_fitted):
1014 results[column_err][idx] = value
1016 for param in params_free_missing:
1017 param.fixed = False
1019 # Convert the x/y errors to ra/dec errors
1020 if config.convert_cen_xy_to_radec:
1021 for (
1022 key_ra_err,
1023 key_dec_err,
1024 key_cen_x,
1025 key_cen_y,
1026 key_cen_x_err,
1027 key_cen_y_err,
1028 ) in columns_params_radec_err:
1029 cen_x, cen_y = results[key_cen_x][idx], results[key_cen_y][idx]
1030 # TODO: improve this in DM-45682
1031 # For one, it won't work right at limits:
1032 # RA=359.99... or dec=+89.99...
1033 # Could also consider dividing by sqrt(2)
1034 # ...but that factor would multiply out later
1035 radec_err = np.array(
1036 self.get_model_radec(
1037 source_multi,
1038 cen_x + results[key_cen_x_err][idx],
1039 cen_y + results[key_cen_y_err][idx],
1040 )
1041 )
1042 radec_err -= radec
1043 results[key_ra_err][idx], results[key_dec_err][idx] = np.abs(radec_err)
1045 results[f"{prefix}chisq_reduced"][idx] = result_full.chisq_best / size
1046 time_final = time.process_time()
1047 results[f"{prefix}time_full"][idx] = time_final - time_init
1048 except Exception as e:
1049 n_skipfail += 1
1050 size = 0 if fitInputs is None else size_new
1051 column = self.errors_expected.get(e.__class__, "")
1052 if column:
1053 row[f"{prefix}{column}"] = True
1054 logger.debug(
1055 "id_source=%i (idx=%i/%i) fit failed with known exception: %s",
1056 id_source,
1057 idx,
1058 n_rows,
1059 e,
1060 )
1061 else:
1062 row[f"{prefix}unknown_flag"] = True
1063 logger.info(
1064 "id_source=%i (idx=%i/%i) fit failed with unexpected exception: %s",
1065 id_source,
1066 idx,
1067 n_rows,
1068 e,
1069 exc_info=1,
1070 )
1071 logger_periodic.log(
1072 "Fit idx=%i/%i sources (%i skipped/failed) in %.2f",
1073 idx,
1074 n_rows,
1075 n_skipfail,
1076 time_final - time_init_all,
1077 )
1079 n_unknown = np.sum(row[f"{prefix}unknown_flag"])
1080 if n_unknown > 0:
1081 logger.warning("%i/%i source fits failed with unexpected exceptions", n_unknown, n_rows)
1083 return results
1085 def get_channels(
1086 self,
1087 catexps: list[CatalogExposureSourcesABC],
1088 ) -> dict[str, g2f.Channel]:
1089 channels = {}
1090 for catexp in catexps:
1091 try:
1092 channel = catexp.channel
1093 except AttributeError:
1094 band = catexp.band
1095 if callable(band):
1096 band = band()
1097 channel = g2f.Channel.get(band)
1098 if channel not in channels:
1099 channels[channel.name] = channel
1100 return channels
1102 def get_model(
1103 self,
1104 idx_row: int,
1105 catalog_multi: Sequence,
1106 catexps: list[CatalogExposureSourcesABC],
1107 config_data: CatalogSourceFitterConfigData | None = None,
1108 results: astropy.table.Table | None = None,
1109 **kwargs: Any,
1110 ) -> g2f.ModelD:
1111 """Reconstruct the model for a single row of a fit catalog.
1113 Parameters
1114 ----------
1115 idx_row
1116 The index of the row in the catalog.
1117 catalog_multi
1118 The multi-band catalog originally used for initialization.
1119 catexps
1120 The catalog-exposure pairs to reconstruct the model for.
1121 config_data
1122 The configuration used to generate sources.
1123 Default-initialized if None.
1124 results
1125 The corresponding best-fit parameter catalog to initialize
1126 parameter values from. If None, the model params will be set by
1127 `self.initialize_model`, as they would be when calling `self.fit`.
1128 **kwargs
1129 Additional keyword arguments to pass to initialize_model. Not
1130 used during fitting.
1132 Returns
1133 -------
1134 model
1135 The reconstructed model.
1136 """
1137 channels = self.get_channels(catexps)
1138 if config_data is None:
1139 config_data = CatalogSourceFitterConfigData(
1140 config=CatalogSourceFitterConfig(),
1141 channels=list(channels.values()),
1142 )
1143 config = config_data.config
1145 if not idx_row >= 0:
1146 raise ValueError(f"{idx_row=} !>=0")
1147 if not len(catalog_multi) > idx_row:
1148 raise ValueError(f"{len(catalog_multi)=} !> {idx_row=}")
1149 if (results is not None) and not (len(results) > idx_row):
1150 raise ValueError(f"{len(results)=} !> {idx_row=}")
1152 model_sources, priors = config_data.sources_priors
1153 source_multi = catalog_multi[idx_row]
1155 data, psf_models = config.make_model_data(
1156 idx_row=idx_row,
1157 catexps=catexps,
1158 )
1159 model = g2f.ModelD(data=data, psfmodels=psf_models, sources=model_sources, priors=priors)
1160 self.initialize_model(model, source_multi, catexps, **kwargs)
1162 if results is not None:
1163 row = results[idx_row]
1164 for column, param in config_data.parameters.items():
1165 param.value = row[f"{config.prefix_column}{column}"]
1167 return model
1169 def get_model_radec(self, source: Mapping[str, Any], cen_x: float, cen_y: float) -> tuple[float, float]:
1170 """Return right ascension and declination values for a source.
1172 Implementing this method is necessary only when fitting data with
1173 accompanying WCS.
1175 Parameters
1176 ----------
1177 source
1178 A mapping with fields expected to be populated in the
1179 corresponding source catalog.
1180 cen_x
1181 The x-axis centroid in pixel coordinates.
1182 cen_y
1183 The y-axis centroid in pixel coordinates.
1185 Returns
1186 -------
1187 ra, dec
1188 The right ascension and declination.
1189 """
1190 return np.nan, np.nan
1192 @abstractmethod
1193 def initialize_model(
1194 self,
1195 model: g2f.ModelD,
1196 source: Mapping[str, Any],
1197 catexps: list[CatalogExposureSourcesABC],
1198 config_data: CatalogSourceFitterConfigData,
1199 values_init: Mapping[g2f.ParameterD, float] | None = None,
1200 **kwargs: Any,
1201 ) -> None:
1202 """Initialize a Model for a single source row.
1204 Parameters
1205 ----------
1206 model
1207 The model object to initialize.
1208 source
1209 A mapping with fields expected to be populated in the
1210 corresponding source catalog for initialization.
1211 catexps
1212 A list of (source and psf) catalog-exposure pairs.
1213 config_data
1214 Configuration settings and data for fitting and output.
1215 values_init
1216 Initial parameter values from the model configuration.
1217 **kwargs
1218 Additional keyword arguments that cannot be required for fitting.
1219 """
1221 @abstractmethod
1222 def validate_fit_inputs(
1223 self,
1224 catalog_multi: Sequence,
1225 catexps: list[CatalogExposureSourcesABC],
1226 config_data: CatalogSourceFitterConfigData = None,
1227 logger: logging.Logger = None,
1228 **kwargs: Any,
1229 ) -> None:
1230 """Validate inputs to self.fit.
1232 This method is called before any fitting is done. It may be used for
1233 any purpose, including checking that the inputs are a particular
1234 subclass of the base classes.
1236 Parameters
1237 ----------
1238 catalog_multi
1239 A multi-band source catalog to fit a model to.
1240 catexps
1241 A list of (source and psf) catalog-exposure pairs.
1242 config_data
1243 Configuration settings and data for fitting and output.
1244 logger
1245 The logger. Defaults to calling `_getlogger`.
1246 **kwargs
1247 Additional keyword arguments to pass to self.modeller.
1248 """
1249 pass