Coverage for python/lsst/multiprofit/fitting/fit_psf.py: 82%
312 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-06-23 08:24 +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: 113 ↛ 124line 113 didn't jump to line 124 because the condition on line 113 was always true
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: 129 ↛ 131line 129 didn't jump to line 131 because the condition on line 129 was always true
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: 144 ↛ 145line 144 didn't jump to line 145 because the condition on line 144 was never true
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): 194 ↛ 195line 194 didn't jump to line 195 because the condition on line 194 was never true
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]: 225 ↛ 226line 225 didn't jump to line 226 because the condition on line 225 was never true
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: 286 ↛ 288line 286 didn't jump to line 288 because the condition on line 286 was always true
287 parameters[f"{key_size}{suffix_x}"] = component.ellipse.size_x_param
288 if not config_comp.size_y.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_y}"] = component.ellipse.size_y_param
290 if not config_comp.rho.fixed: 290 ↛ 295line 290 didn't jump to line 295 because the condition on line 290 was always true
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: 302 ↛ 303line 302 didn't jump to line 303 because the condition on line 302 was never true
303 raise RuntimeError(
304 f"{params_flux=} has len={len(params_flux)} but expected {n_params_flux_expect}"
305 )
306 if has_params_flux: 306 ↛ 307line 306 didn't jump to line 307 because the condition on line 306 was never true
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: 315 ↛ 332line 315 didn't jump to line 332 because the condition on line 315 was always true
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): 319 ↛ 320line 319 didn't jump to line 320 because the condition on line 319 was never true
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: 324 ↛ 327line 324 didn't jump to line 327 because the condition on line 324 was always true
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: 327 ↛ 328line 327 didn't jump to line 328 because the condition on line 327 was never true
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: 399 ↛ 401line 399 didn't jump to line 401 because the condition on line 399 was always true
400 modeller = Modeller()
401 if errors_expected is None: 401 ↛ 403line 401 didn't jump to line 403 because the condition on line 401 was always true
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: 425 ↛ 431line 425 didn't jump to line 431 because the condition on line 425 was always true
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): 434 ↛ 436line 434 didn't jump to line 436 because the condition on line 434 was always true
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: 520 ↛ 521line 520 didn't jump to line 521 because the condition on line 520 was never true
521 config_data = CatalogPsfFitterConfigData(config=CatalogPsfFitterConfig())
522 if logger is None: 522 ↛ 524line 522 didn't jump to line 524 because the condition on line 522 was always true
523 logger = CatalogPsfFitter._get_logger()
524 config = config_data.config
525 if config.compute_errors != "NONE": 525 ↛ 526line 525 didn't jump to line 526 because the condition on line 525 was never true
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): 530 ↛ 531line 530 didn't jump to line 531 because the condition on line 530 was never true
531 raise ValueError(f"{self.errors_expected=} has duplicate values; they must be unique")
532 if n_errors_expected != len(config.flag_errors): 532 ↛ 533line 532 didn't jump to line 533 because the condition on line 532 was never true
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: 548 ↛ 549line 548 didn't jump to line 549 because the condition on line 548 was never true
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: 557 ↛ 558line 557 didn't jump to line 558 because the condition on line 557 was never true
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: 561 ↛ 569line 561 didn't jump to line 569 because the condition on line 561 was always true
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: 622 ↛ 625line 622 didn't jump to line 625 because the condition on line 622 was always true
623 fitInputs = fitInputs if not fitInputs.validate_for_model(model) else None
625 if config.fit_linear_init: 625 ↛ 641line 625 didn't jump to line 641 because the condition on line 625 was always true
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: 647 ↛ 650line 647 didn't jump to line 650 because the condition on line 647 was always true
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: 691 ↛ 692line 691 didn't jump to line 692 because the condition on line 691 was never true
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: 719 ↛ 721line 719 didn't jump to line 721 because the condition on line 719 was always true
720 limits_x = g2f.LimitsD(0, n_cols)
721 if limits_y is None: 721 ↛ 724line 721 didn't jump to line 724 because the condition on line 721 was always true
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