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