Coverage for python / lsst / multiprofit / fitting / fit_psf.py: 15%
312 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-17 08:58 +0000
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/>.
22__all__ = [
23 "CatalogExposurePsfABC",
24 "CatalogPsfFitterConfig",
25 "CatalogPsfFitterConfigData",
26 "CatalogPsfFitter",
27]
29from abc import abstractmethod
30from functools import cached_property
31import logging
32import math
33import time
34from typing import Any, ClassVar, Mapping, Type
36import astropy
37import lsst.gauss2d as g2
38import lsst.gauss2d.fit as g2f
39import lsst.pex.config as pexConfig
40from lsst.utils.logging import PeriodicLogger
41import numpy as np
42import pydantic
44from ..componentconfig import (
45 Fluxes,
46 FluxFractionParameterConfig,
47 FluxParameterConfig,
48 GaussianComponentConfig,
49 ParameterConfig,
50)
51from ..errors import PsfRebuildFitFlagError
52from ..modeller import FitInputsDummy, LinearGaussians, Modeller, make_psf_model_null
53from ..sourceconfig import ComponentConfigs, ComponentGroupConfig, SourceConfig
54from ..utils import frozen_arbitrary_allowed_config, get_params_uniq
55from .fit_catalog import CatalogExposureABC, CatalogFitterConfig, ColumnInfo
57SourceRow = astropy.table.Row | Mapping[str, Any]
60class CatalogPsfFitterConfig(CatalogFitterConfig):
61 """Configuration for MultiProFit PSF image fitter."""
63 model = pexConfig.ConfigField[SourceConfig](
64 default=SourceConfig(
65 component_groups={
66 "": ComponentGroupConfig(
67 components_gauss={
68 "gauss1": GaussianComponentConfig(
69 size_x=ParameterConfig(value_initial=1.5),
70 size_y=ParameterConfig(value_initial=1.5),
71 fluxfrac=FluxFractionParameterConfig(value_initial=0.5),
72 flux=FluxParameterConfig(value_initial=1.0, fixed=True),
73 ),
74 "gauss2": GaussianComponentConfig(
75 size_x=ParameterConfig(value_initial=3.0),
76 size_y=ParameterConfig(value_initial=3.0),
77 fluxfrac=FluxFractionParameterConfig(value_initial=1.0, fixed=True),
78 ),
79 },
80 is_fractional=True,
81 )
82 }
83 ),
84 doc="PSF model configuration",
85 )
86 prior_axrat_mean = pexConfig.Field[float](default=0.95, doc="Mean for axis ratio prior")
87 sigma_min = pexConfig.Field[float](
88 default=0.8,
89 doc="Minimum sigma in pixels for PSF components. Must be >=0.8 to avoid undersampling.",
90 check=lambda x: x >= 0.8,
91 )
93 def make_psf_model(
94 self,
95 component_group_fluxes: list[list[Fluxes]] | None = None,
96 ) -> [g2f.PsfModel, list[g2f.Prior]]:
97 """Make a PsfModel object for a given source.
99 Parameters
100 ----------
101 component_group_fluxes
102 Initial fluxes for each constituent ComponentGroup.
104 Returns
105 -------
106 psf_model
107 The rebuilt PSF model.
109 Notes
110 -----
111 This function does not initialize the PSF model.
112 """
113 if component_group_fluxes is None:
114 channels = (g2f.Channel.NONE,)
115 component_group_fluxes = [
116 component_group.get_fluxes_default(
117 channels=channels,
118 component_configs=component_group.get_component_configs(),
119 is_fractional=component_group.is_fractional,
120 )
121 for component_group in self.model.component_groups.values()
122 ]
124 psf_model, _ = self.model.make_psf_model(component_group_fluxes=component_group_fluxes)
125 return psf_model
127 def schema_configurable(self) -> list[ColumnInfo]:
128 columns = []
129 if self.config_fit.eval_residual:
130 columns.append(ColumnInfo(key="n_eval_jac", dtype="i4"))
131 return columns
133 def schema(
134 self,
135 bands: list[str] | None = None,
136 ) -> list[ColumnInfo]:
137 """Return the schema as an ordered list of columns.
139 Parameters
140 ----------
141 bands
142 The bands to add band-dependent columns for.
143 """
144 if bands is not None:
145 if len(bands) != 1:
146 raise ValueError("CatalogPsfFitter must have exactly one band")
147 schema = super().schema(bands)
148 parameters = CatalogPsfFitterConfigData(config=self).parameters
149 unit_size = "pix"
150 units = {
151 g2f.ReffXParameterD: unit_size,
152 g2f.ReffYParameterD: unit_size,
153 g2f.SizeXParameterD: unit_size,
154 g2f.SizeYParameterD: unit_size,
155 }
156 schema.extend(
157 [
158 ColumnInfo(key=key, dtype="f8", unit=units.get(type(param)))
159 for key, param in parameters.items()
160 ]
161 )
162 schema.extend(self.schema_configurable())
164 return schema
166 def setDefaults(self) -> None:
167 self.prefix_column = "mpf_psf_"
168 self.compute_errors = "NONE"
171class CatalogPsfFitterConfigData(pydantic.BaseModel):
172 """A PSF fit configuration that can initialize models and images thereof.
174 This class relies on cached properties being computed once, mostly shortly
175 after initialization. Therefore, it and the config field must be frozen to
176 ensure that the model remains unchanged.
177 """
179 model_config: ClassVar[pydantic.ConfigDict] = frozen_arbitrary_allowed_config
181 config: CatalogPsfFitterConfig = pydantic.Field(title="A CatalogPsfFitterConfig to be frozen")
183 @pydantic.field_validator("config")
184 @classmethod
185 def validate_config(cls, v: CatalogPsfFitterConfig) -> CatalogPsfFitterConfig:
186 v.validate()
187 return v
189 @cached_property
190 def components(self) -> dict[str, g2f.Component]:
191 """Return the components of the PSF model by name."""
192 components = self.psf_model.components
193 names = self.component_configs.keys()
194 if len(components) != len(names):
195 raise RuntimeError(f"{len(components)=} != {len(names)=}")
196 components_names = {name: component for name, component in zip(names, components)}
197 return components_names
199 @cached_property
200 def component_configs(self) -> ComponentConfigs:
201 """Return the config for each component in the PSF model."""
202 return self.config.model.get_component_configs()
204 @cached_property
205 def componentgroup_configs(self) -> dict[str, ComponentGroupConfig]:
206 """Return the config for each component group in the PSF model
207 by name.
208 """
209 return {k: v for k, v in self.config.model.component_groups.items()}
211 def init_psf_model(
212 self,
213 params: SourceRow,
214 ) -> None:
215 """Initialize the PSF model for a single source.
217 Parameters
218 ----------
219 params : astropy.table.Row | typing.Mapping[str, typing.Any]
220 A mapping with parameter values for the best-fit PSF model at the
221 centroid of a single source.
222 """
223 # TODO: Improve _flag checking (add a total _flag column)
224 for flag in (col for col in params.keys() if col.endswith("_flag")):
225 if params[flag]:
226 raise PsfRebuildFitFlagError(f"Failed to rebuild PSF; {flag} set")
228 for name, param in self.parameters.items():
229 param.value = params[f"{self.config.prefix_column}{name}"]
231 def model_post_init(self, __context: Any) -> None:
232 self.config.freeze()
233 n_component_configs = len(self.component_configs)
234 n_components = len(self.psf_model.components)
235 if n_components != n_component_configs:
236 raise AssertionError(f"{n_components=} != {n_component_configs=}")
238 @cached_property
239 def parameters(self) -> dict[str, g2f.ParameterD]:
240 """Return the free parameters for the PSF model by name."""
241 parameters = {}
242 config = self.config
243 has_prefix_group = config.model.has_prefix_group()
244 components = self.psf_model.components
245 idx_comp_first = 0
247 label_cen = config.get_key_cen()
248 label_rho = config.get_key_rho()
249 label_flux = config.get_key_flux("")
250 label_fluxfrac = f"{label_flux}frac"
251 suffix_x, suffix_y = config.get_suffix_x(), config.get_suffix_y()
253 # Iterate over each component group
254 for name_group, config_group in self.componentgroup_configs.items():
255 prefix_group = f"{name_group}_" if has_prefix_group else ""
256 is_fractional = config_group.is_fractional
257 # Does this group have independent centroids?
258 # Should generally be False if not is_fractional, since there is
259 # little reason to make a group otherwise.
260 multicen = len(config_group.centroids) > 1
261 configs_comp = config_group.get_component_configs()
262 idx_last = len(configs_comp) - 1
263 n_params_flux_frac = 0
265 # Iterate over each component in the group
266 for idx_comp_group, (name_comp, config_comp) in enumerate(configs_comp.items()):
267 # The last component needs special handling if is_fractional
268 is_last = idx_comp_group == idx_last
269 component = components[idx_comp_first + idx_comp_group]
270 prefix_comp = f"{prefix_group}{name_comp}"
271 key_size = config.get_prefixed_label(
272 config.get_key_size(config_comp.get_size_label()),
273 prefix_comp,
274 )
275 key_rho = config.get_prefixed_label(label_rho, prefix_comp)
277 # Give the centroid parameters an appropriate prefix
278 if multicen or (idx_comp_group == 0):
279 prefix_cen = prefix_comp if multicen else prefix_group
280 # Avoid redundant -underscores if there's nothing to prefix
281 # or an existing prefix starting with an underscore
282 key_cen = config.get_prefixed_label(label_cen, prefix_cen)
283 parameters[f"{key_cen}{suffix_x}"] = component.centroid.x_param
284 parameters[f"{key_cen}{suffix_y}"] = component.centroid.y_param
285 # Add each free shape parameter
286 if not config_comp.size_x.fixed:
287 parameters[f"{key_size}{suffix_x}"] = component.ellipse.size_x_param
288 if not config_comp.size_y.fixed:
289 parameters[f"{key_size}{suffix_y}"] = component.ellipse.size_y_param
290 if not config_comp.rho.fixed:
291 parameters[key_rho] = component.ellipse.rho_param
293 # TODO: return this to component.integralmodel
294 # when binding for g2f.FractionalIntegralModel is fixed
295 params_flux = get_params_uniq(component, fixed=False, nonlinear=False)
296 # Fractional models should have one fixed flux parameter
297 # (with value=1), and hence an empty params_flux
298 has_params_flux = not config_comp.flux.fixed and (
299 (not is_fractional) or (idx_comp_group == 0)
300 )
301 n_params_flux_expect = int(has_params_flux)
302 if len(params_flux) != n_params_flux_expect:
303 raise RuntimeError(
304 f"{params_flux=} has len={len(params_flux)} but expected {n_params_flux_expect}"
305 )
306 if has_params_flux:
307 parameters[f"{prefix_comp}{label_flux}"] = params_flux[0]
308 # TODO: return this to component.integralmodel
309 # when binding for g2f.FractionalIntegralModel is fixed
310 params_fluxfrac = [
311 param
312 for param in get_params_uniq(component, fixed=False, linear=False)
313 if isinstance(param, g2f.ProperFractionParameterD)
314 ]
315 if is_fractional:
316 # The last flux fraction must be fixed at 1.0, since it
317 # by definition contains all of the remaining flux
318 if is_last:
319 if (config_comp.fluxfrac.value_initial != 1.0) or (not config_comp.fluxfrac.fixed):
320 raise ValueError(
321 f"{config_comp=} {is_last=} and must be fixed with value_initial==1.0"
322 )
323 else:
324 if not config_comp.fluxfrac.fixed:
325 parameters[f"{prefix_comp}{label_fluxfrac}"] = params_fluxfrac[-1]
326 n_params_flux_frac += 1
327 if len(params_fluxfrac) != n_params_flux_frac:
328 raise RuntimeError(
329 f"{config_comp=} has {params_fluxfrac=} but expected {n_params_flux_frac=}"
330 )
331 else:
332 if len(params_fluxfrac) > 0:
333 raise RuntimeError(f"{config_group=} has {params_fluxfrac=} but {is_fractional=}")
335 return parameters
337 @cached_property
338 def psf_model(self) -> g2f.PsfModel:
339 """Return the PSF model."""
340 psf_model = self.config.make_psf_model()
341 return psf_model
343 @cached_property
344 def psf_model_gaussians(self) -> g2.Gaussians:
345 """Return the Gaussians comprising the PSF model."""
346 gaussians = self.psf_model.gaussians()
347 return gaussians
350class CatalogExposurePsfABC(CatalogExposureABC):
351 """A CatalogExposure for PSF fitting."""
353 @abstractmethod
354 def get_psf_image(
355 self,
356 source: SourceRow,
357 ) -> np.ndarray:
358 """Get a PSF image for a specific source.
360 Parameters
361 ----------
362 source
363 The source row/dict.
365 Returns
366 -------
367 psf
368 The image of the PSF.
370 Notes
371 -----
372 The PSF image should be normalized, and centered in a 2D array of odd
373 dimensions on both sides.
374 """
377class CatalogPsfFitter:
378 """Fit a Gaussian mixture model to a pixelated PSF image.
380 Parameters
381 ----------
382 modeller : `multiprofit.Modeller`
383 A Modeller instance to use for fitting.
384 errors_expected : dict[Type[Exception], str]
385 A dictionary keyed by an Exception type, with a string value of the
386 flag column key to assign if this Exception is raised.
388 Notes
389 -----
390 Any exceptions raised and not in errors_expected will be logged in a
391 generic unknown_flag failure column.
392 """
394 def __init__(
395 self,
396 modeller: Modeller | None = None,
397 errors_expected: dict[Type[Exception], str] | None = None,
398 ):
399 if modeller is None:
400 modeller = Modeller()
401 if errors_expected is None:
402 errors_expected = {}
403 self.errors_expected = errors_expected
404 self.modeller = modeller
406 @staticmethod
407 def _get_data_default(img_psf: np.ndarray, gain: float = 1e5) -> g2f.DataD:
408 """Return a default-initialized data for a given PSF image.
410 Parameters
411 ----------
412 img_psf
413 A normalized image of the PSF, ideally with noise.
414 gain
415 The multiplicative factor needed to change the image to counts.
417 Returns
418 -------
419 data
420 A DataD instance than can be used to fit a PSF model.
421 """
422 # TODO: Try to improve these ad-hoc settings after DM-49008 merges
423 # e.g. estimate PSF model noise from S/N of PSF stars?
424 # If more than 1% of pixels are negative, try to estimate noise
425 if np.sum(img_psf < 0) / img_psf.size > 0.01:
426 background = np.std(img_psf[img_psf < 2 * np.abs(np.min(img_psf))])
427 # otherwise, try to estimate noise from the lowest-values pixels
428 # There might not actually be much noise; in that case, the
429 # background value is not going to matter much anyway
430 else:
431 background = np.max((1e-10, np.std(img_psf[img_psf <= np.percentile(img_psf, 5.0)])))
432 min_psf = np.min(img_psf)
433 # ensure there are no negative values before sqrt
434 if not (background > -min_psf):
435 background = -1.1 * min_psf
436 img_sig_inv = np.sqrt(gain / (img_psf + background))
437 return g2f.DataD(
438 [
439 g2f.ObservationD(
440 channel=g2f.Channel.NONE,
441 image=g2.ImageD(img_psf),
442 sigma_inv=g2.ImageD(img_sig_inv),
443 mask_inv=g2.ImageB(np.ones_like(img_psf)),
444 )
445 ]
446 )
448 def _get_data(self, img_psf: np.ndarray, gain: float = 1e5) -> g2f.DataD:
449 """Build a Model-able gauss2d.fit.Data from a normalized PSF image.
451 Parameters
452 ----------
453 img_psf
454 A normalized PSF image array.
455 gain
456 The number of counts in the image, used as a multiplicative
457 factor for the inverse variance.
459 Returns
460 -------
461 data
462 A Data object that can be passed to a Model(ler).
463 """
464 return self._get_data_default(img_psf=img_psf, gain=gain)
466 @staticmethod
467 def _get_logger() -> logging.Logger:
468 """Return a suitably-named and configured logger."""
469 logger = logging.getLogger(__name__)
470 return logger
472 @abstractmethod
473 def check_source(self, source: SourceRow, config: CatalogPsfFitterConfig) -> None:
474 """Check whether a source can have its PSF model fit.
476 Parameters
477 ----------
478 source
479 The source row to check.
480 config
481 The fitter config.
483 Notes
484 -----
485 Derived classes may use the source row as they deem fit. For example,
486 if the source has poor quality flags, a fitter may choose not to fit
487 the PSF model if it will not end up being used anyway.
488 """
490 def fit(
491 self,
492 catexp: CatalogExposurePsfABC,
493 config_data: CatalogPsfFitterConfigData | None = None,
494 logger: logging.Logger | None = None,
495 **kwargs: Any,
496 ) -> astropy.table.Table:
497 """Fit PSF models for a catalog with MultiProFit.
499 Each source has its PSF fit with a configureable Gaussian mixture PSF
500 model, given a pixellated PSF image from the CatalogExposure.
502 Parameters
503 ----------
504 catexp
505 An exposure to fit a model PSF at the position of all
506 sources in the corresponding catalog.
507 config_data
508 Configuration settings for fitting and output.
509 logger
510 The logger. Defaults to calling `_getlogger`.
511 **kwargs
512 Additional keyword arguments to pass to self.modeller.
514 Returns
515 -------
516 catalog
517 A table with fit parameters for the PSF model at the location
518 of each source.
519 """
520 if config_data is None:
521 config_data = CatalogPsfFitterConfigData(config=CatalogPsfFitterConfig())
522 if logger is None:
523 logger = CatalogPsfFitter._get_logger()
524 config = config_data.config
525 if config.compute_errors != "NONE":
526 raise ValueError("CatalogPsfFitter doesn't support computing errors")
528 errors_expected = set(self.errors_expected.values())
529 n_errors_expected = len(errors_expected)
530 if n_errors_expected != len(self.errors_expected):
531 raise ValueError(f"{self.errors_expected=} has duplicate values; they must be unique")
532 if n_errors_expected != len(config.flag_errors):
533 raise ValueError(f"len({self.errors_expected=}) != len({config.flag_errors=})")
535 priors = []
536 sigmas = [
537 np.linalg.norm((comp.size_x.value_initial, comp.size_y.value_initial))
538 for comp in config_data.component_configs.values()
539 ]
541 psf_model = config_data.psf_model
542 model_source = g2f.Source(psf_model.components)
544 for idx, (comp, config_comp) in enumerate(
545 zip(psf_model.components, config_data.component_configs.values())
546 ):
547 prior = config_comp.make_shape_prior(comp.ellipse)
548 if prior:
549 if prior_size := prior.prior_size:
550 prior_size.mean = sigmas[idx]
551 if prior_axrat := prior.prior_axrat:
552 prior_axrat.mean = config.prior_axrat_mean
553 priors.append(prior)
555 params = config_data.parameters
556 flux_total = tuple(get_params_uniq(psf_model, nonlinear=False, channel=g2f.Channel.NONE))
557 if len(flux_total) != 1:
558 raise RuntimeError(f"len({flux_total=}) != 1; PSF model is badly-formed")
559 flux_total = flux_total[0]
560 gaussians_linear = None
561 if config.fit_linear_init:
562 # The total flux must be freed first or else LinearGaussians.make
563 # will fail to find the required number of free linear params
564 flux_total.fixed = False
565 gaussians_linear = LinearGaussians.make(model_source, is_psf=True)
566 flux_total.fixed = True
568 # TODO: Remove isinstance when channel filtering is fixed: DM-45665
569 fluxfracs = tuple(
570 param
571 for param in get_params_uniq(model_source, linear=False, channel=g2f.Channel.NONE, fixed=False)
572 if isinstance(param, g2f.ProperFractionParameterD)
573 )
574 # We're fitting the PSF, so make a single Gaussian
575 model_psf = make_psf_model_null()
576 # Set the size to the minimum sigma to avoid undersampling
577 ellipse = model_psf.components[0].ellipse
578 ellipse.sigma_x_param.value = config.sigma_min
579 ellipse.sigma_y_param.value = config.sigma_min
580 sigma_min_sq = config.sigma_min**2
582 catalog = catexp.get_catalog()
583 n_rows = len(catalog)
584 range_idx = range(n_rows)
586 results, columns = config.make_catalog(n_rows)
587 prefix = config.prefix_column
588 columns_param = {}
589 for key, param in params.items():
590 is_sigma = isinstance(param, g2f.SigmaXParameterD) or isinstance(param, g2f.SigmaYParameterD)
591 columns_param[f"{prefix}{key}"] = param, is_sigma
592 if is_sigma:
593 param.value = math.sqrt(max(param.value**2 - sigma_min_sq, 0.1))
595 # dummy size for first iteration
596 size, size_new = 0, 0
597 fitInputs = FitInputsDummy()
598 time_init_all = time.process_time()
599 logger_periodic = PeriodicLogger(logger)
600 n_skipfail = 0
602 for idx in range_idx:
603 time_init = time.process_time()
604 row = results[idx]
605 source = catalog[idx]
606 id_source = source[config.column_id]
607 row[config.column_id] = id_source
608 time_final = time_init
610 try:
611 self.check_source(source, config=config)
612 img_psf = catexp.get_psf_image(source)
613 data = self._get_data(img_psf)
614 model = g2f.ModelD(data=data, psfmodels=[model_psf], sources=[model_source], priors=priors)
615 self.initialize_model(model=model, config_data=config_data)
617 # Caches the jacobian residual if the kernel size is unchanged
618 if img_psf.size != size:
619 fitInputs = None
620 size = int(img_psf.size)
621 # Some algorithms might not even use fitInputs
622 elif fitInputs is not None:
623 fitInputs = fitInputs if not fitInputs.validate_for_model(model) else None
625 if config.fit_linear_init:
626 result = self.modeller.fit_gaussians_linear(gaussians_linear, data[0])
627 result = list(result.values())[0]
628 # Re-normalize fluxes (hopefully close already)
629 result = np.clip(
630 result * np.array([x[1].value for x in gaussians_linear.gaussians_free]),
631 1e-2,
632 0.99,
633 )
634 result /= np.sum(result)
635 for idx_param, param in enumerate(fluxfracs):
636 param.value = result[idx_param]
637 # fluxfracs are the fraction of the remaining flux
638 # this renormalizes to the remaining parameters
639 result /= np.sum(result[idx_param + 1 :])
641 result_full = self.modeller.fit_model(model, fitinputs=fitInputs, **kwargs)
642 fitInputs = result_full.inputs
643 results[f"{prefix}n_iter"][idx] = result_full.n_eval_func
644 results[f"{prefix}time_eval"][idx] = result_full.time_eval
645 results[f"{prefix}time_fit"][idx] = result_full.time_run
646 results[f"{prefix}chisq_reduced"][idx] = result_full.chisq_best / size
647 if config.config_fit.eval_residual:
648 results[f"{prefix}n_eval_jac"][idx] = result_full.n_eval_jac
650 for (key, (param, is_sigma)), value in zip(columns_param.items(), result_full.params_best):
651 param.value_transformed = value
652 value = param.value
653 if is_sigma:
654 value = math.sqrt(sigma_min_sq**2 + value**2)
655 results[key][idx] = value
657 time_final = time.process_time()
658 results[f"{prefix}time_full"][idx] = time_final - time_init
659 except Exception as e:
660 n_skipfail += 1
661 size = 0 if fitInputs is None else size_new
662 column = self.errors_expected.get(e.__class__, "")
663 if column:
664 row[f"{prefix}{column}"] = True
665 logger.debug(
666 "id_source=%i (idx=%i/%i) PSF fit failed with known exception: %s",
667 id_source,
668 idx,
669 n_rows,
670 e,
671 )
672 else:
673 row[f"{prefix}unknown_flag"] = True
674 logger.info(
675 "id_source=%i (idx=%i/%i) PSF fit failed with unexpected exception: %s",
676 id_source,
677 idx,
678 n_rows,
679 e,
680 exc_info=1,
681 )
682 logger_periodic.log(
683 "Fit idx=%i/%i PSFs (%i skipped/failed) in %.2f",
684 idx,
685 n_rows,
686 n_skipfail,
687 time_final - time_init_all,
688 )
690 n_unknown = np.sum(row[f"{prefix}unknown_flag"])
691 if n_unknown > 0:
692 logger.warning("%i/%i PSF fits failed with unexpected exceptions", n_unknown, n_rows)
694 return results
696 def initialize_model(
697 self,
698 model: g2f.ModelD,
699 config_data: CatalogPsfFitterConfigData,
700 limits_x: g2f.LimitsD | None = None,
701 limits_y: g2f.LimitsD | None = None,
702 ) -> None:
703 """Initialize a Model for a single source row.
705 Parameters
706 ----------
707 model
708 The model object to initialize.
709 config_data
710 The fitter config with cached data.
711 limits_x
712 Hard limits for the source's x centroid. Default (0, n_columns).
713 limits_y
714 Hard limits for the source's y centroid. Default (0, n_rows).
715 """
716 n_rows, n_cols = model.data[0].image.data.shape
717 cen_x, cen_y = n_cols / 2.0, n_rows / 2.0
718 centroids = set()
719 if limits_x is None:
720 limits_x = g2f.LimitsD(0, n_cols)
721 if limits_y is None:
722 limits_y = g2f.LimitsD(0, n_rows)
724 for component, config_comp in zip(
725 config_data.components.values(), config_data.component_configs.values()
726 ):
727 centroid = component.centroid
728 if centroid not in centroids:
729 centroid.x_param.value = cen_x
730 centroid.x_param.limits = limits_x
731 centroid.y_param.value = cen_y
732 centroid.y_param.limits = limits_y
733 centroids.add(centroid)
734 ellipse = component.ellipse
735 ellipse.size_x_param.limits = limits_x
736 ellipse.size_x = config_comp.size_x.value_initial
737 ellipse.size_y_param.limits = limits_y
738 ellipse.size_y = config_comp.size_y.value_initial
739 ellipse.rho = config_comp.rho.value_initial