Coverage for python/lsst/multiprofit/modeller.py: 66%
514 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-22 09:52 +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 "FitInputsBase",
24 "FitInputsDummy",
25 "FitResult",
26 "InvalidProposalError",
27 "LinearGaussians",
28 "ModelFitConfig",
29 "Modeller",
30 "fit_methods_linear",
31 "make_image_gaussians",
32 "make_psf_model_null",
33]
35import logging
36import sys
37import time
38from abc import ABC, abstractmethod
39from collections.abc import Iterable, Sequence
40from typing import Any, ClassVar, Self, TypeAlias
42import numpy as np
43import pydantic
44import scipy.optimize as spopt
46import lsst.gauss2d as g2
47import lsst.gauss2d.fit as g2f
48import lsst.pex.config as pexConfig
50from .model_utils import make_image_gaussians, make_psf_model_null
51from .utils import arbitrary_allowed_config, frozen_arbitrary_allowed_config, get_params_uniq
53_has_py_13_plus = sys.version_info >= (3, 13, 0)
54if _has_py_13_plus: 54 ↛ 57line 54 didn't jump to line 57 because the condition on line 54 was always true
55 Model: type = g2f.ModelD | g2f.ModelF
56else:
57 Model: TypeAlias = g2f.ModelD | g2f.ModelF # noqa: UP040
59try:
60 # TODO: try importlib.util.find_spec
61 from fastnnls import fnnls
63 has_fastnnls = True
64except ImportError:
65 has_fastnnls = False
67try:
68 # TODO: try importlib.util.find_spec
69 import pygmo as pg
71 has_pygmo = True
72except ImportError:
73 has_pygmo = False
76class InvalidProposalError(ValueError):
77 """Error for an invalid parameter proposal."""
80fit_methods_linear = {
81 "scipy.optimize.nnls": {},
82 "scipy.optimize.lsq_linear": {"bounds": (1e-5, np.inf), "method": "bvls"},
83 "numpy.linalg.lstsq": {"rcond": 1e-3},
84}
85if has_fastnnls: 85 ↛ 86line 85 didn't jump to line 86 because the condition on line 85 was never true
86 fit_methods_linear["fastnnls.fnnls"] = {}
89class LinearGaussians(pydantic.BaseModel):
90 """Helper for linear least-squares fitting of Gaussian mixtures."""
92 model_config: ClassVar[pydantic.ConfigDict] = frozen_arbitrary_allowed_config
94 gaussians_fixed: g2.Gaussians = pydantic.Field(title="Fixed Gaussian components")
95 gaussians_free: tuple[tuple[g2.Gaussians, g2f.ParameterD], ...] = pydantic.Field(
96 title="Free Gaussian components"
97 )
99 @staticmethod
100 def make(
101 component_mixture: g2f.ComponentMixture,
102 channel: g2f.Channel = None,
103 is_psf: bool = False,
104 ) -> Self:
105 """Make a LinearGaussians from a ComponentMixture.
107 Parameters
108 ----------
109 component_mixture
110 A component mixture to initialize Gaussians from.
111 channel
112 The channel all Gaussians are applicable for.
113 is_psf
114 Whether the components are a smoothing kernel.
116 Returns
117 -------
118 lineargaussians
119 A LinearGaussians instance initialized with the appropriate
120 fixed/free gaussians.
121 """
122 if channel is None:
123 channel = g2f.Channel.NONE
124 components = component_mixture.components
125 if len(components) == 0: 125 ↛ 126line 125 didn't jump to line 126 because the condition on line 125 was never true
126 raise ValueError(f"Can't get linear Source from {component_mixture=} with no components")
128 gaussians_free = []
129 gaussians_fixed = []
131 for component in components:
132 gaussians: g2.Gaussians = component.gaussians(channel)
133 # TODO: Support multi-Gaussian components if sensible
134 # The challenge would be in mapping linear param values back onto
135 # non-linear IntegralModels
136 if is_psf:
137 n_g = len(gaussians)
138 if n_g != 1: 138 ↛ 139line 138 didn't jump to line 139 because the condition on line 138 was never true
139 raise ValueError(f"{component=} has {gaussians=} of len {n_g=}!=1")
140 param_fluxes = component.parameters(paramfilter=g2f.ParamFilter(nonlinear=False, channel=channel))
141 if len(param_fluxes) != 1: 141 ↛ 142line 141 didn't jump to line 142 because the condition on line 141 was never true
142 raise ValueError(f"Can't make linear source from {component=} with {len(param_fluxes)=}")
143 param_flux: g2f.ParameterD = param_fluxes[0]
144 if param_flux.fixed: 144 ↛ 145line 144 didn't jump to line 145 because the condition on line 144 was never true
145 gaussians_fixed.append(gaussians.at(0))
146 else:
147 gaussians_free.append((gaussians, param_flux))
149 return LinearGaussians(
150 gaussians_fixed=g2.Gaussians(gaussians_fixed), gaussians_free=tuple(gaussians_free)
151 )
154class FitInputsBase(ABC):
155 """Interface for inputs to a model fit."""
157 @abstractmethod
158 def validate_for_model(self, model: Model) -> list[str]:
159 """Check that this FitInputs is valid for a Model.
161 Parameters
162 ----------
163 model
164 The model to validate with.
166 Returns
167 -------
168 errors
169 A list of validation errors, if any.
170 """
173class FitInputsDummy(FitInputsBase):
174 """A dummy FitInputs that always fails to validate.
176 This class can be used to initialize a FitInputsBase that may be
177 reassigned to a non-dummy derived instance in a loop.
178 """
180 def validate_for_model(self, model: Model) -> list[str]:
181 return [
182 "This is a dummy FitInputs and will never validate",
183 ]
186class FitInputs(FitInputsBase, pydantic.BaseModel):
187 """Model fit inputs for gauss2dfit."""
189 model_config: ClassVar[pydantic.ConfigDict] = arbitrary_allowed_config
191 jacobian: np.ndarray = pydantic.Field(None, title="The full Jacobian array")
192 jacobians: list[list[g2.ImageD]] = pydantic.Field(
193 title="Jacobian arrays (views) for each observation",
194 )
195 outputs_prior: tuple[g2.ImageD, ...] = pydantic.Field(
196 title="Jacobian arrays (views) for each free parameter's prior",
197 )
198 residual: np.ndarray = pydantic.Field(title="The full residual (chi) array")
199 residuals: list[g2.ImageD] = pydantic.Field(
200 default_factory=list,
201 title="Residual (chi) arrays (views) for each observation",
202 )
203 residuals_prior: g2.ImageD = pydantic.Field(
204 title="Shared residual array for all Prior instances",
205 )
207 @classmethod
208 def get_sizes(
209 cls,
210 model: Model,
211 ) -> tuple[int, int, int, np.ndarray]:
212 """Initialize Jacobian and residual arrays for a model.
214 Parameters
215 ----------
216 model : `lsst.gauss2d.fit.Model`
217 The model to initialize arrays for.
219 Returns
220 -------
221 n_obs
222 The number of observations initialized.
223 n_params_jac
224 The number of Jacobian matrix columns, which is the number of free
225 parameters plus one validation column.
226 n_prior_residuals
227 The number of residual array values required for priors.
228 shapes
229 An ndarray containing the number of rows and columns for each
230 observation in rows.
231 """
232 priors = model.priors
233 n_prior_residuals = sum(len(p) for p in priors)
234 params_free = tuple(get_params_uniq(model, fixed=False))
235 n_params_free = len(params_free)
236 # gauss2d_fit reserves the zeroth index of the jacobian array for
237 # validation, i.e. it can be used to dump terms for fixed params
238 n_params_jac = n_params_free + 1
239 if not (n_params_jac > 1): 239 ↛ 240line 239 didn't jump to line 240 because the condition on line 239 was never true
240 raise ValueError("Can't fit model with no free parameters")
242 n_obs = len(model.data)
243 shapes = np.zeros((n_obs, 2), dtype=int)
244 ranges_params = [None] * n_obs
246 for idx_obs in range(n_obs):
247 observation = model.data[idx_obs]
248 shapes[idx_obs, :] = (observation.image.n_rows, observation.image.n_cols)
249 # Get the free parameter indices for each observation
250 params = tuple(get_params_uniq(model, fixed=False, channel=observation.channel))
251 n_params_obs = len(params)
252 ranges_params_obs = [0] * (n_params_obs + 1)
253 for idx_param in range(n_params_obs):
254 ranges_params_obs[idx_param + 1] = params_free.index(params[idx_param]) + 1
255 ranges_params[idx_obs] = ranges_params_obs
257 n_free_first = len(ranges_params[0])
258 # Ensure that there are the same number of free parameters in each obs
259 # They don't need to be the same set, but the counts should equal
260 # (this assumption may be violated by future IntegralModels - TBD)
261 assert all([len(rp) == n_free_first for rp in ranges_params[1:]])
263 return n_obs, n_params_jac, n_prior_residuals, shapes
265 @classmethod
266 def from_model(
267 cls,
268 model: Model,
269 ) -> Self:
270 """Initialize Jacobian and residual arrays for a model.
272 Parameters
273 ----------
274 model : `gauss2d.fit.Model`
275 The model to initialize arrays for.
276 """
277 n_obs, n_params_jac, n_prior_residuals, shapes = cls.get_sizes(model)
278 n_pixels_cumsum = np.cumsum(np.prod(shapes, axis=1))
279 n_pixels_total = n_pixels_cumsum[-1]
280 size_data = n_pixels_total + n_prior_residuals
281 shape_jacobian = (size_data, n_params_jac)
282 jacobian = np.zeros(shape_jacobian)
283 jacobians = [None] * n_obs
284 outputs_prior = [None] * n_params_jac
285 for idx in range(n_params_jac):
286 outputs_prior[idx] = g2.ImageD(jacobian[n_pixels_total:, idx].reshape((1, n_prior_residuals)))
288 residual = np.zeros(size_data)
289 residuals = [None] * n_obs
290 residuals_prior = g2.ImageD(residual[n_pixels_total:].reshape(1, n_prior_residuals))
292 offset = 0
293 for idx_obs in range(n_obs):
294 shape = shapes[idx_obs, :]
295 size_obs = shape[0] * shape[1]
296 end = offset + size_obs
297 jacobians_obs = [None] * n_params_jac
298 for idx_jac in range(n_params_jac):
299 jacobians_obs[idx_jac] = g2.ImageD(jacobian[offset:end, idx_jac].reshape(shape))
300 jacobians[idx_obs] = jacobians_obs
301 residuals[idx_obs] = g2.ImageD(residual[offset:end].reshape(shape))
302 offset = end
303 if offset != n_pixels_cumsum[idx_obs]: 303 ↛ 304line 303 didn't jump to line 304 because the condition on line 303 was never true
304 raise RuntimeError(f"Assigned {offset=} data points != {n_pixels_cumsum[idx_obs]=}")
305 return cls(
306 jacobian=jacobian,
307 jacobians=jacobians,
308 residual=residual,
309 residuals=residuals,
310 outputs_prior=tuple(outputs_prior),
311 residuals_prior=residuals_prior,
312 )
314 def validate_for_model(self, model: Model) -> list[str]:
315 n_obs, n_params_jac, n_prior_residuals, shapes = self.get_sizes(model)
316 n_pixels_total = np.sum(np.prod(shapes, axis=1))
317 size_data = n_pixels_total + n_prior_residuals
318 shape_jacobian = (size_data, n_params_jac)
320 errors = []
322 if self.jacobian.shape != shape_jacobian: 322 ↛ 323line 322 didn't jump to line 323 because the condition on line 322 was never true
323 errors.append(f"{self.jacobian.shape=} != {shape_jacobian=}")
325 if len(self.jacobians) != n_obs: 325 ↛ 326line 325 didn't jump to line 326 because the condition on line 325 was never true
326 errors.append(f"{len(self.jacobians)=} != {n_obs=}")
328 if len(self.residuals) != n_obs: 328 ↛ 329line 328 didn't jump to line 329 because the condition on line 328 was never true
329 errors.append(f"{len(self.residuals)=} != {n_obs=}")
331 if not errors: 331 ↛ 345line 331 didn't jump to line 345 because the condition on line 331 was always true
332 for idx_obs in range(n_obs):
333 shape_obs = shapes[idx_obs, :]
334 jacobian_obs = self.jacobians[idx_obs]
335 n_jacobian_obs = len(jacobian_obs)
336 if n_jacobian_obs != n_params_jac: 336 ↛ 337line 336 didn't jump to line 337 because the condition on line 336 was never true
337 errors.append(f"len(self.jacobians[{idx_obs}])={n_jacobian_obs} != {n_params_jac=}")
338 else:
339 for idx_jac in range(n_jacobian_obs):
340 if not all(jacobian_obs[idx_jac].shape == shape_obs): 340 ↛ 341line 340 didn't jump to line 341 because the condition on line 340 was never true
341 errors.append(f"{jacobian_obs[idx_jac].shape=} != {shape_obs=}")
342 if not all(self.residuals[idx_obs].shape == shape_obs): 342 ↛ 343line 342 didn't jump to line 343 because the condition on line 342 was never true
343 errors.append(f"{self.residuals[idx_obs].shape=} != {shape_obs=}")
345 shape_residual_prior = [1, n_prior_residuals]
346 if len(self.outputs_prior) != n_params_jac: 346 ↛ 347line 346 didn't jump to line 347 because the condition on line 346 was never true
347 errors.append(f"{len(self.outputs_prior)=} != {n_params_jac=}")
348 elif n_prior_residuals > 0:
349 for idx in range(n_params_jac):
350 if self.outputs_prior[idx].shape != shape_residual_prior: 350 ↛ 351line 350 didn't jump to line 351 because the condition on line 350 was never true
351 errors.append(f"{self.outputs_prior[idx].shape=} != {shape_residual_prior=}")
353 if n_prior_residuals > 0:
354 if self.residuals_prior.shape != shape_residual_prior: 354 ↛ 355line 354 didn't jump to line 355 because the condition on line 354 was never true
355 errors.append(f"{self.residuals_prior.shape=} != {shape_residual_prior=}")
357 return errors
360class ModelFitConfig(pexConfig.Config):
361 """Configuration for model fitting."""
363 eval_residual = pexConfig.Field[bool](
364 doc="Whether to evaluate the residual every iteration before the Jacobian, which can improve "
365 "performance if most steps do not call the Jacobian function. Must be set to True if the "
366 "optimizer does not always evaluate the residual first, before the Jacobian.",
367 default=True,
368 )
369 fit_linear_iter = pexConfig.Field[int](
370 doc="The number of iterations to wait before performing a linear fit during optimization."
371 " Default 0 disables the feature.",
372 default=0,
373 )
374 optimization_library = pexConfig.ChoiceField[str](
375 doc="The optimization library to use when fitting",
376 allowed={
377 "pygmo": "Pygmo2",
378 "scipy": "scipy.optimize",
379 },
380 default="scipy",
381 )
383 def validate(self) -> None:
384 if not self.fit_linear_iter >= 0: 384 ↛ 385line 384 didn't jump to line 385 because the condition on line 384 was never true
385 raise ValueError(f"{self.fit_linear_iter=} must be >=0")
388class FitResult(pydantic.BaseModel):
389 """Results from a Modeller fit, including metadata."""
391 model_config: ClassVar[pydantic.ConfigDict] = arbitrary_allowed_config
393 chisq_best: float = pydantic.Field(default=0, title="The chi-squared (sum) of the best-fit parameters")
394 # TODO: Why does setting default=ModelFitConfig() cause a circular import?
395 config: ModelFitConfig = pydantic.Field(None, title="The configuration for fitting")
396 inputs: FitInputs | None = pydantic.Field(None, title="The fit input arrays")
397 result: Any | None = pydantic.Field(
398 None,
399 title="The result object of the fit, directly from the optimizer",
400 )
401 params: tuple[g2f.ParameterD, ...] | None = pydantic.Field(
402 None,
403 title="The parameter instances corresponding to params_best",
404 )
405 params_best: tuple[float, ...] | None = pydantic.Field(
406 None,
407 title="The best-fit parameter array (un-transformed)",
408 )
409 params_free_missing: tuple[g2f.ParameterD, ...] | None = pydantic.Field(
410 None,
411 title="Free parameters that were fixed during fitting - usually an"
412 " IntegralParameterD for a band with missing data",
413 )
414 n_eval_resid: int = pydantic.Field(0, title="Total number of self-reported residual function evaluations")
415 n_eval_func: int = pydantic.Field(
416 0, title="Total number of optimizer-reported fitness function evaluations"
417 )
418 n_eval_jac: int = pydantic.Field(
419 0, title="Total number of optimizer-reported Jacobian function evaluations"
420 )
421 time_eval: float = pydantic.Field(0, title="Total runtime spent in model/Jacobian evaluation")
422 time_run: float = pydantic.Field(0, title="Total runtime spent in fitting, excluding initial setup")
425def set_params(params: Iterable[g2f.ParameterD], params_new: Iterable[float], model_loglike: Model):
426 """Set new parameter values from an optimizer proposal.
428 Parameters
429 ----------
430 params
431 An iterable of ParameterD instances.
432 params_new
433 An iterable of new untransformed values for params.
434 model_loglike
435 A model instance configured to compute the log-likelihood.
437 Raises
438 ------
439 InvalidProposalError
440 Raised if a new value is nan, or if a RuntimeError is raised when
441 setting the new value.
442 RuntimeError
443 Raised if the new transformed value is not finite.
444 """
445 try:
446 for param, value in zip(params, params_new, strict=True):
447 if np.isnan(value): 447 ↛ 448line 447 didn't jump to line 448 because the condition on line 447 was never true
448 raise InvalidProposalError(
449 f"optimizer for {model_loglike=} proposed non-finite {value=} for {param=}"
450 )
451 param.value_transformed = value
452 if not np.isfinite(param.value): 452 ↛ 453line 452 didn't jump to line 453 because the condition on line 452 was never true
453 raise RuntimeError(f"{param=} set to (transformed) non-finite {value=}")
454 except RuntimeError as e:
455 raise InvalidProposalError(f"optimizer for {model_loglike=} proposal generated error={e}")
458def residual_scipy(
459 params_new: np.ndarray,
460 model_jacobian: Model,
461 model_loglike: Model,
462 params: tuple[g2f.ParameterD],
463 result: FitResult,
464 jacobian: np.ndarray | None,
465 never_evaluate_jacobian: bool = False,
466 return_loglike: bool = False,
467) -> np.ndarray:
468 """Compute the residual for a scipy optimizer.
470 Parameters
471 ----------
472 params_new
473 An array of new parameter values.
474 model_jacobian
475 A model instance configured to compute the Jacobian.
476 model_loglike
477 A model instance configured to compute the log-likelihood.
478 params
479 A tuple of the free parameters. The length and order must be identical
480 to params_new.
481 result
482 A FitResult instance to update.
483 jacobian
484 The Jacobian array. Unused in this function.
485 never_evaluate_jacobian
486 If True, the jacobian will never be evaluated, taking precedence
487 over result.config.eval_residual.
488 return_loglike
489 If False, will return the negative of the residual instead of the
490 log-likelihood.
492 Returns
493 -------
494 result
495 The log-likehood if return_loglike, otherwise the negative of the
496 residual from result.inputs.residual.
498 Notes
499 -----
500 Scipy requires that this function have the same args as the jacobian
501 function (jacobian_scipy), so unused args must not be removed.
503 Scipy generally calls this function every iteration, but only conditionally
504 calls the jacobian_scipy function (e.g. if the proposal is accepted). If
505 users expect proposals to (almost) always be accepted, it is more efficient
506 to compute the Jacobian here (and skip evaluating it again when
507 jacobian_scipy is called), because evaluating model_jacobian also updates
508 the residual array, and so there is no need to evaluate model_loglike.
510 To summarize, if never_evaluate_jacobian or config_fit.eval_residual:
511 There is ALWAYS one call to model_jacobian.evaluate,
512 and ZERO calls to model_loglike.evaluate.
513 else:
514 There is always one call to model_loglike.evaluate,
515 and MAYBE one call to model_jacobian.evaluate.
516 """
517 set_params(params, params_new, model_loglike)
518 config_fit = result.config
519 fit_linear_iter = config_fit.fit_linear_iter
520 if (fit_linear_iter > 0) and ((result.n_eval_resid + 1) % fit_linear_iter == 0):
521 Modeller.fit_model_linear(model_loglike, ratio_min=1e-6)
522 time_init = time.process_time()
524 if never_evaluate_jacobian or config_fit.eval_residual: 524 ↛ 531line 524 didn't jump to line 531 because the condition on line 524 was always true
525 try:
526 loglike = model_loglike.evaluate()
527 except Exception:
528 loglike = None
529 result.n_eval_resid += 1
530 else:
531 loglike = model_jacobian.evaluate()
532 result.n_eval_jac += 1
533 result.time_eval += time.process_time() - time_init
534 return loglike if return_loglike else -result.inputs.residual
537def jacobian_scipy(
538 params_new: np.ndarray,
539 model_jacobian: Model,
540 model_loglike: Model,
541 params: tuple[g2f.ParameterD],
542 result: FitResult,
543 jacobian: np.ndarray,
544 always_evaluate_jacobian: bool = False,
545) -> np.ndarray:
546 """Compute the Jacobian for a scipy optimizer.
548 Parameters
549 ----------
550 params_new
551 An array of new parameter values. Unused here.
552 model_jacobian
553 A model instance configured to compute the Jacobian.
554 model_loglike
555 A model instance configured to compute the log-likelihood. Unused here.
556 params
557 A tuple of the free parameters. Unused here.
558 result
559 A FitResult instance to update.
560 jacobian
561 The Jacobian array. Unused in this function.
562 always_evaluate_jacobian
563 If True, the jacobian will always be evaluated, taking precedence
564 over result.config.eval_residual.
566 Returns
567 -------
568 jacobian
569 A reference to jacobian, whose values may have been updated.
571 Notes
572 -----
573 Scipy requires that this function have the same args as the residual
574 function (residual_scipy), so unused args must not be removed. kwargs are
575 for the convenience of libraries other than scipy and will not be changed
576 by scipy itself.
578 Parameter objects and new values are unused here as they will have already
579 been set by the residual funciton.
581 Scipy generally does not call this function every iteration. If it is
582 configured to skip evaluating the Jacobian, it is presumed to have been
583 updated by the residual function already.
584 """
585 if always_evaluate_jacobian or result.config.eval_residual: 585 ↛ 590line 585 didn't jump to line 590 because the condition on line 585 was always true
586 time_init = time.process_time()
587 model_jacobian.evaluate()
588 result.time_eval += time.process_time() - time_init
589 result.n_eval_jac += 1
590 return jacobian
593if has_pygmo: 593 ↛ 595line 593 didn't jump to line 595 because the condition on line 593 was never true
595 class PygmoUDP:
596 """A Pygmo User-Defined Problem for a MultiProFit model.
598 Pygmo optimizers take a class with a fitness function
599 (i.e. the negative log-likelihood, although one could use some other
600 arbitrary fitness function if it made snese to do so), with a
601 fitness function and a gradient function returning the derivative of
602 the fitness w.r.t. each free parameters.
604 Pygmo optimizers do not appear to use the full residual array or
605 Jacobian the way scipy optimizers do. The gradient of the
606 log-likelihood is cheaper to compute than the full Jacobian; however,
607 using only the gradient of the fitness may cause slower convergence.
609 Parameters
610 ----------
611 params
612 A tuple of the free parameters.
613 model_loglike
614 A model configured to compute the log-likelihood.
615 model_loglike_grad
616 A model configured to compute the gradient of the
617 log-likelihood w.r.t. each free parameter.
618 bounds_lower
619 A tuple of the lower bounds of the transformed value for each
620 free parameter in params.
621 bounds_upper
622 A tuple of the upper bounds of the transformed value for each
623 free parameter in params.
624 result
625 A result object to update and read configuration from.
626 """
628 def __init__(
629 self,
630 params: tuple[g2f.ParameterD],
631 model_loglike: Model,
632 model_loglike_grad: Model,
633 bounds_lower: tuple[float],
634 bounds_upper: tuple[float],
635 result: FitResult,
636 ):
637 self.params = params
638 self.model_loglike = model_loglike
639 self.model_loglike_grad = model_loglike_grad
640 self.bounds_lower = bounds_lower
641 self.bounds_upper = bounds_upper
642 self.result = result
644 def fitness(self, x):
645 loglike = residual_scipy(
646 x,
647 model_jacobian=self.model_loglike_grad,
648 model_loglike=self.model_loglike,
649 params=self.params,
650 result=self.result,
651 jacobian=None,
652 never_evaluate_jacobian=True,
653 return_loglike=True,
654 )
655 return [
656 -sum(loglike),
657 ]
659 def get_bounds(self):
660 return self.bounds_lower, self.bounds_upper
662 def gradient(self, x):
663 set_params(params=self.params, params_new=x, model_loglike=self.model_loglike)
664 time_init = time.process_time()
665 loglike_grad = -np.array(self.model_loglike_grad.compute_loglike_grad())
666 self.result.time_eval += time.process_time() - time_init
667 self.result.n_eval_jac += 1
668 return loglike_grad
670 def __deepcopy__(self, memo):
671 """Make a deep copy of a model with a shallow copy of the data
672 which should not be duplicated.
674 Pygmo optimizers always make at least one copy of this class, and
675 some (like particle swarm) will make many more. The input data
676 must be shallow copies, both to avoid excess memory usage and
677 because Model instances cannot be deep copied.
678 """
679 fitinputs = FitInputs.from_model(self.model_loglike)
680 model_loglike, model_loglike_grad = (
681 g2f.ModelD(
682 data=model.data,
683 psfmodels=model.psfmodels,
684 sources=model.sources,
685 priors=model.priors,
686 )
687 for model in (self.model_loglike, self.model_loglike_grad)
688 )
689 model_loglike.setup_evaluators(evaluatormode=g2f.EvaluatorMode.loglike)
690 model_loglike_grad.setup_evaluators(evaluatormode=g2f.EvaluatorMode.loglike_grad)
692 copied = self.__class__(
693 params=self.params,
694 model_loglike=model_loglike,
695 model_loglike_grad=model_loglike_grad,
696 bounds_lower=self.bounds_lower,
697 bounds_upper=self.bounds_upper,
698 result=FitResult(inputs=fitinputs, config=self.result.config),
699 )
700 memo[id(self)] = copied
701 return copied
704class Modeller:
705 """Fit lsst.gauss2d.fit Model instances using Python optimizers.
707 Parameters
708 ----------
709 logger : `logging.Logger`
710 The logger. Defaults to calling `_getlogger`.
711 """
713 def __init__(self, logger: logging.Logger | None = None) -> None:
714 if logger is None: 714 ↛ 716line 714 didn't jump to line 716 because the condition on line 714 was always true
715 logger = self._get_logger()
716 self.logger = logger
718 @staticmethod
719 def _get_logger() -> logging.Logger:
720 logger = logging.getLogger(__name__)
721 return logger
723 @staticmethod
724 def compute_variances(
725 model: Model, use_diag_only: bool = False, use_svd: bool = False, **kwargs: Any
726 ) -> np.ndarray:
727 """Compute model free parameter variances from the inverse Hessian.
729 Parameters
730 ----------
731 model
732 The model to compute parameter variances for.
733 use_diag_only
734 Whether to use diagonal terms only, i.e. ignore covariance.
735 use_svd
736 Whether to use singular value decomposition to compute the inverse
737 Hessian.
738 **kwargs
739 Additional keyword arguments to pass to model.compute_hessian.
741 Returns
742 -------
743 variances
744 The free parameter variances.
745 """
746 hessian = model.compute_hessian(**kwargs).data
747 if use_diag_only:
748 return -1 / np.diag(hessian)
749 if use_svd: 749 ↛ 750line 749 didn't jump to line 750 because the condition on line 749 was never true
750 u, s, v = np.linalg.svd(-hessian)
751 inverse = np.dot(v.transpose(), np.dot(np.diag(s**-1), u.transpose()))
752 else:
753 inverse = np.linalg.inv(-hessian)
754 return np.diag(inverse)
756 @staticmethod
757 def fit_gaussians_linear(
758 gaussians_linear: LinearGaussians,
759 observation: g2f.ObservationD,
760 psf_model: g2f.PsfModel = None,
761 fit_methods: dict[str, dict[str, Any]] | None = None,
762 plot: bool = False,
763 ) -> dict[str, FitResult]:
764 """Fit normalizations for a Gaussian mixture model.
766 Parameters
767 ----------
768 gaussians_linear
769 The Gaussian components - fixed or otherwise - to fit.
770 observation
771 The observation to fit against.
772 psf_model
773 A PSF model for the observation, if fitting sources.
774 fit_methods
775 A dictionary of fitting methods to employ, keyed by method name,
776 with a value of a dict of options (kwargs) to pass on. Default
777 is "scipy.optimize.nnls".
778 plot
779 Whether to generate fit residual/diagnostic plots.
781 Returns
782 -------
783 results
784 Fit results for each method, keyed by the fit method name.
785 """
786 if psf_model is None:
787 psf_model = make_psf_model_null()
788 if fit_methods is None:
789 fit_methods = {"scipy.optimize.nnls": fit_methods_linear["scipy.optimize.nnls"]}
790 else:
791 for fit_method in fit_methods:
792 if fit_method not in fit_methods_linear: 792 ↛ 793line 792 didn't jump to line 793 because the condition on line 792 was never true
793 raise ValueError(f"Unknown linear {fit_method=}")
794 n_params = len(gaussians_linear.gaussians_free)
795 if not (n_params > 0): 795 ↛ 796line 795 didn't jump to line 796 because the condition on line 795 was never true
796 raise ValueError(f"!({len(gaussians_linear.gaussians_free)=}>0); can't fit with no free params")
797 image = observation.image
798 shape = image.shape
799 coordsys = image.coordsys
801 mask_inv = observation.mask_inv.data
802 sigma_inv = observation.sigma_inv.data
803 bad = ~(sigma_inv > 0)
804 n_bad = np.sum(bad)
805 if n_bad > 0: 805 ↛ 806line 805 didn't jump to line 806 because the condition on line 805 was never true
806 mask_inv &= ~bad
808 sigma_inv = sigma_inv[mask_inv]
809 size = np.sum(mask_inv)
811 gaussians_psf = psf_model.gaussians(g2f.Channel.NONE)
812 if len(gaussians_linear.gaussians_fixed) > 0: 812 ↛ 813line 812 didn't jump to line 813 because the condition on line 812 was never true
813 image_fixed = make_image_gaussians(
814 gaussians_source=gaussians_linear.gaussians_fixed,
815 gaussians_kernel=gaussians_psf,
816 n_rows=shape[0],
817 n_cols=shape[1],
818 ).data
819 image_fixed = image_fixed[mask_inv]
820 else:
821 image_fixed = None
823 x = np.zeros((size, n_params))
825 params = [None] * n_params
826 for idx_param, (gaussians_free, param) in enumerate(gaussians_linear.gaussians_free):
827 image_free = make_image_gaussians(
828 gaussians_source=gaussians_free,
829 gaussians_kernel=gaussians_psf,
830 n_rows=shape[0],
831 n_cols=shape[1],
832 coordsys=coordsys,
833 ).data
834 x[:, idx_param] = ((image_free if mask_inv is None else image_free[mask_inv]) * sigma_inv).flat
835 params[idx_param] = param
837 y = observation.image.data
838 if plot: 838 ↛ 839line 838 didn't jump to line 839 because the condition on line 838 was never true
839 import matplotlib.pyplot as plt
841 plt.imshow(y, origin="lower")
842 plt.show()
843 if mask_inv is not None: 843 ↛ 845line 843 didn't jump to line 845 because the condition on line 843 was always true
844 y = y[mask_inv]
845 if image_fixed is not None: 845 ↛ 846line 845 didn't jump to line 846 because the condition on line 845 was never true
846 y -= image_fixed
847 y = (y * sigma_inv).flat
849 results = {}
851 for fit_method, kwargs in fit_methods.items():
852 kwargs = kwargs if kwargs is not None else fit_methods_linear[fit_method]
853 if fit_method == "scipy.optimize.nnls":
854 values = spopt.nnls(x, y)[0]
855 elif fit_method == "scipy.optimize.lsq_linear":
856 values = spopt.lsq_linear(x, y, **kwargs).x
857 elif fit_method == "numpy.linalg.lstsq": 857 ↛ 859line 857 didn't jump to line 859 because the condition on line 857 was always true
858 values = np.linalg.lstsq(x, y, **kwargs)[0]
859 elif fit_method == "fastnnls.fnnls":
860 y = x.T.dot(y)
861 x = x.T.dot(x)
862 values = fnnls(x, y)
863 else:
864 raise RuntimeError(f"Unknown linear {fit_method=} not caught earlier (logic error)")
865 results[fit_method] = values
866 return results
868 def fit_model(
869 self,
870 model: Model,
871 fitinputs: FitInputs | None = None,
872 printout: bool = False,
873 config: ModelFitConfig | None = None,
874 **kwargs: Any,
875 ) -> FitResult:
876 """Fit a model with a nonlinear optimizer.
878 Parameters
879 ----------
880 model
881 The model to fit.
882 fitinputs
883 An existing FitInputs with jacobian/residual arrays to reuse.
884 printout
885 Whether to print diagnostic information.
886 config
887 Configuration settings for model fitting.
888 **kwargs
889 Keyword arguments to pass to the optimizer.
891 Returns
892 -------
893 result
894 The results from running the fitter.
896 Notes
897 -----
898 The only supported fitter is scipy.optimize.least_squares.
899 """
900 if config is None:
901 config = ModelFitConfig()
902 config.validate()
904 use_pygmo = config.optimization_library == "pygmo"
905 model_loglike = (
906 g2f.ModelD(
907 data=model.data,
908 psfmodels=model.psfmodels,
909 sources=model.sources,
910 priors=model.priors,
911 )
912 if (use_pygmo or config.eval_residual)
913 else None
914 )
916 if use_pygmo: 916 ↛ 917line 916 didn't jump to line 917 because the condition on line 916 was never true
917 model.setup_evaluators(g2f.EvaluatorMode.loglike_grad, force=True)
918 model_loglike.setup_evaluators(g2f.EvaluatorMode.loglike, force=True)
919 else:
920 if fitinputs is None:
921 fitinputs = FitInputs.from_model(model)
922 else:
923 errors = fitinputs.validate_for_model(model)
924 if errors: 924 ↛ 925line 924 didn't jump to line 925 because the condition on line 924 was never true
925 newline = "\n"
926 raise ValueError(f"fitinputs validation got errors:\n{newline.join(errors)}")
927 model.setup_evaluators(
928 evaluatormode=g2f.EvaluatorMode.jacobian,
929 outputs=fitinputs.jacobians,
930 residuals=fitinputs.residuals,
931 outputs_prior=fitinputs.outputs_prior,
932 residuals_prior=fitinputs.residuals_prior,
933 print=printout,
934 force=True,
935 )
937 params_psf_free = []
938 for psfmodel in model.psfmodels:
939 params_psf_free.extend(get_params_uniq(psfmodel, fixed=False))
940 if params_psf_free: 940 ↛ 941line 940 didn't jump to line 941 because the condition on line 940 was never true
941 params_psf_free = {k: None for k in params_psf_free}
942 raise ValueError(
943 f"Model has free PSF model params: {list(params_psf_free.keys())}."
944 f" All PSF model parameters must be fixed before fitting."
945 )
947 offsets_params = dict(model.offsets_parameters())
948 params_offsets = {v: k for (k, v) in offsets_params.items()}
949 params_free = tuple(params_offsets[idx] for idx in range(1, len(offsets_params) + 1))
950 params_free_sorted_all = tuple(get_params_uniq(model, fixed=False))
951 params_free_sorted = []
952 params_free_sorted_missing = []
954 # If we were forced to drop an observation, re-generate the modeller
955 # Only integral parameters should be missing
956 for param in params_free_sorted_all:
957 if param in params_offsets.values(): 957 ↛ 960line 957 didn't jump to line 960 because the condition on line 957 was always true
958 params_free_sorted.append(param)
959 else:
960 if not isinstance(param, g2f.IntegralParameterD):
961 raise RuntimeError(f"non-integral {param=} missing from {offsets_params=}")
962 param.limits = g2f.LimitsD(param.min, param.max)
963 param.value = param.min
964 param.fixed = True
965 params_free_sorted_missing.append(param)
967 try:
968 if not use_pygmo: 968 ↛ 994line 968 didn't jump to line 994 because the condition on line 968 was always true
969 if params_free_sorted_missing: 969 ↛ 970line 969 didn't jump to line 970 because the condition on line 969 was never true
970 fitinputs = FitInputs.from_model(model)
971 params_free_sorted = tuple(params_free_sorted)
972 model.setup_evaluators(
973 evaluatormode=g2f.EvaluatorMode.jacobian,
974 outputs=fitinputs.jacobians,
975 residuals=fitinputs.residuals,
976 outputs_prior=fitinputs.outputs_prior,
977 residuals_prior=fitinputs.residuals_prior,
978 print=printout,
979 force=True,
980 )
981 else:
982 params_free_sorted = params_free_sorted_all
983 if config.eval_residual: 983 ↛ 990line 983 didn't jump to line 990 because the condition on line 983 was always true
984 model_loglike.setup_evaluators(
985 evaluatormode=g2f.EvaluatorMode.loglike,
986 residuals=fitinputs.residuals,
987 residuals_prior=fitinputs.residuals_prior,
988 )
990 jac = fitinputs.jacobian[:, 1:]
991 # Assert that this is a view, otherwise this won't work
992 assert id(jac.base) == id(fitinputs.jacobian)
994 n_params_free = len(params_free)
995 bounds = ([None] * n_params_free, [None] * n_params_free)
996 params_init = [None] * n_params_free
998 for idx, param in enumerate(params_free):
999 limits = param.limits
1000 # If the transform has more restrictive limits, use those
1001 if hasattr(param.transform, "limits"):
1002 limits_transform = param.transform.limits
1003 n_within = limits.check(limits_transform.min) + limits.check(limits_transform.min)
1004 if n_within == 2:
1005 limits = limits_transform
1006 elif n_within != 0: 1006 ↛ 1007line 1006 didn't jump to line 1007 because the condition on line 1006 was never true
1007 raise ValueError(
1008 f"{param=} {param.limits=} and {param.transform.limits=}"
1009 f" intersect; one must be a subset of the other"
1010 )
1011 bounds[0][idx] = param.transform.forward(limits.min)
1012 bounds[1][idx] = param.transform.forward(limits.max)
1013 if not limits.check(param.value): 1013 ↛ 1014line 1013 didn't jump to line 1014 because the condition on line 1013 was never true
1014 raise RuntimeError(f"{param=}.value_transformed={param.value} not within {limits=}")
1015 params_init[idx] = param.value_transformed
1017 results = FitResult(inputs=fitinputs, config=config)
1018 time_init = time.process_time()
1019 if use_pygmo: 1019 ↛ 1020line 1019 didn't jump to line 1020 because the condition on line 1019 was never true
1020 uda = pg.nlopt("lbfgs")
1021 uda.ftol_abs = 1e-4
1022 algo = pg.algorithm(uda)
1024 # pygmo seems to make proposals right at the limits
1025 # parameter limits are currently set as untransformed values
1026 # and sometimes the proposal exceeds those when transformed
1027 bounds_lower = tuple(np.nextafter(x, np.inf) for x in bounds[0])
1028 bounds_upper = tuple(np.nextafter(x, -np.inf) for x in bounds[1])
1030 udp = PygmoUDP(
1031 params=params_free,
1032 model_loglike=model_loglike,
1033 model_loglike_grad=model,
1034 bounds_lower=bounds_lower,
1035 bounds_upper=bounds_upper,
1036 result=results,
1037 )
1039 # if the initial value was at one of the bounds, reset it to
1040 # one percent of the range away from the bound
1041 for idx, value_init in enumerate(params_init):
1042 bound_lower = bounds_lower[idx]
1043 bound_upper = bounds_upper[idx]
1044 if value_init >= bound_upper:
1045 params_init[idx] = bound_lower + 0.99 * (bound_upper - bound_lower)
1046 elif value_init <= bound_lower:
1047 params_init[idx] = bound_lower + 0.01 * (bound_upper - bound_lower)
1049 problem = pg.problem(udp)
1050 pop = pg.population(prob=problem, size=0)
1051 pop.push_back(np.array(params_init))
1052 result_opt = algo.evolve(pop)
1053 x_best = result_opt.champion_x
1054 results.n_eval_func = pop.problem.get_fevals()
1055 results.n_eval_jac = pop.problem.get_gevals()
1056 results.chisq_best = 2 * result_opt.champion_f
1057 else:
1058 # The initial evaluate will fill in jac for the next line
1059 # _ll_init is assigned just for convenient debugging
1060 _ll_init = model.evaluate()
1061 x_scale_jac_clipped = np.clip(1.0 / (np.sum(jac**2, axis=0) ** 0.5), 1e-5, 1e19)
1062 result_opt = spopt.least_squares(
1063 residual_scipy,
1064 params_init,
1065 jac=jacobian_scipy,
1066 bounds=bounds,
1067 args=(model, model_loglike, params_free, results, jac),
1068 x_scale=x_scale_jac_clipped,
1069 **kwargs,
1070 )
1071 x_best = result_opt.x
1072 results.n_eval_func = result_opt.nfev
1073 results.n_eval_jac = result_opt.njev if result_opt.njev else 0
1074 results.chisq_best = 2 * result_opt.cost
1076 results.time_run = time.process_time() - time_init
1077 results.result = result_opt
1078 if params_free_sorted_missing: 1078 ↛ 1079line 1078 didn't jump to line 1079 because the condition on line 1078 was never true
1079 params_best = []
1080 for param in params_free_sorted_all:
1081 if param in params_free_sorted_missing:
1082 params_best.append(param.value)
1083 param.fixed = False
1084 else:
1085 params_best.append(x_best[offsets_params[param] - 1])
1086 results.params_best = tuple(params_best)
1087 results.params = params_free_sorted_all
1088 else:
1089 results.params_best = tuple(x_best[offsets_params[param] - 1] for param in params_free_sorted)
1090 results.params = params_free_sorted
1091 results.params_free_missing = tuple(params_free_sorted_missing)
1092 except Exception as e:
1093 # Any missing params we fixed must be set free again
1094 for param in params_free_sorted_missing:
1095 param.fixed = False
1096 raise e
1098 return results
1100 # TODO: change to staticmethod if requiring py3.10+
1101 @classmethod
1102 def fit_model_linear(
1103 cls,
1104 model: Model,
1105 idx_obs: int | Sequence[int] | None = None,
1106 ratio_min: float = 0,
1107 validate: bool = False,
1108 limits_interval_min: float = 0.01,
1109 limits_interval_max: float = 1.0,
1110 ) -> tuple[np.ndarray, np.ndarray]:
1111 """Fit a model's linear parameters (integrals).
1113 Parameters
1114 ----------
1115 model
1116 The model to fit parameters for.
1117 idx_obs
1118 An index or sequence of indices of observations to fit.
1119 The default is to fit all observations.
1120 ratio_min
1121 The minimum ratio of the previous value to set any parameter to.
1122 This can prevent setting parameters to zero.
1123 validate
1124 If True, check that the model log-likelihood improves and restore
1125 the original parameter values if not.
1126 limits_interval_min
1127 A value 0<=x<limits_interval_max<=1 specifying the lower bound to
1128 clip parameter values to, as a ratio of each parameter's limits.
1129 limits_interval_max
1130 A value 0<=limits_interval_min<x<=1 specifying the upper bound to
1131 clip parameter values to, as a ratio of each parameter's limits.
1133 Returns
1134 -------
1135 loglike_init
1136 The initial log likelihood if validate is True, otherwise None.
1137 loglike_final
1138 The post-fit log likelihood if validate is True, otherwise None.
1140 Notes
1141 -----
1142 The purpose of limits_interval is to slightly offset parameters from
1143 the extrema of their limits. This is typically most useful for
1144 integral parameters with a minimum of zero, which might otherwise be
1145 stuck at zero in a subsequent nonlinear fit.
1146 """
1147 if ( 1147 ↛ 1152line 1147 didn't jump to line 1152 because the condition on line 1147 was never true
1148 not (0 <= limits_interval_min <= 1)
1149 or not (0 <= limits_interval_max <= 1)
1150 or not (limits_interval_min < limits_interval_max)
1151 ):
1152 raise ValueError(f"Must have 0 <= {limits_interval_min} < {limits_interval_max} <= 1")
1153 n_data = len(model.data)
1154 n_sources = len(model.sources)
1155 if n_sources != 1: 1155 ↛ 1156line 1155 didn't jump to line 1156 because the condition on line 1155 was never true
1156 raise ValueError("fit_model_linear does not yet support models with >1 sources")
1157 if idx_obs is not None: 1157 ↛ 1158line 1157 didn't jump to line 1158 because the condition on line 1157 was never true
1158 if isinstance(idx_obs, int):
1159 if not ((idx_obs >= 0) and (idx_obs < n_data)):
1160 raise ValueError(f"{idx_obs=} not >=0 and < {len(model.data)=}")
1161 indices = range(idx_obs, idx_obs + 1)
1162 else:
1163 if len(set(idx_obs)) != len(idx_obs):
1164 raise ValueError(f"{idx_obs=} has duplicate values")
1165 indices = tuple(idx_obs)
1166 if not all((idx_obs >= 0) and (idx_obs < n_data) for idx_obs in indices):
1167 raise ValueError(f"idx_obs={indices} has values not >=0 and < {len(model.data)=}")
1168 else:
1169 indices = range(n_data)
1171 if validate:
1172 model.setup_evaluators(evaluatormode=g2f.EvaluatorMode.loglike)
1173 loglike_init = model.evaluate()
1174 else:
1175 loglike_init = None
1176 values_init = {}
1177 values_new = {}
1179 for idx_obs in indices:
1180 obs = model.data[idx_obs]
1181 gaussians_linear = LinearGaussians.make(model.sources[0], channel=obs.channel)
1182 result = cls.fit_gaussians_linear(gaussians_linear, obs, psf_model=model.psfmodels[idx_obs])
1183 values = list(result.values())[0]
1185 for (_, parameter), ratio in zip(gaussians_linear.gaussians_free, values):
1186 values_init[parameter] = float(parameter.value)
1187 if not (ratio >= ratio_min): 1187 ↛ 1188line 1187 didn't jump to line 1188 because the condition on line 1187 was never true
1188 ratio = ratio_min
1189 value_new = max(ratio * parameter.value, parameter.limits.min)
1190 values_new[parameter] = value_new
1192 for parameter, value in values_new.items():
1193 value_min, value_max = parameter.limits.min, parameter.limits.max
1194 min_is_inf = value_min == -np.inf
1195 max_is_inf = value_max == np.inf
1196 if min_is_inf: 1196 ↛ 1197line 1196 didn't jump to line 1197 because the condition on line 1196 was never true
1197 if max_is_inf:
1198 parameter.value = value
1199 continue
1200 if not value < value_max:
1201 value = value_max - 1e-5
1202 elif max_is_inf: 1202 ↛ 1206line 1202 didn't jump to line 1206 because the condition on line 1202 was always true
1203 if not value > value_min: 1203 ↛ 1204line 1203 didn't jump to line 1204 because the condition on line 1203 was never true
1204 value = value_min + 1e-5
1205 else:
1206 limits_interval = parameter.limits.max - parameter.limits.min
1207 value = np.clip(
1208 value,
1209 value_min + limits_interval_min * limits_interval,
1210 value_min + limits_interval_max * limits_interval,
1211 )
1212 parameter.value = value
1214 if validate:
1215 loglike_new = model.evaluate()
1216 if not (sum(loglike_new) > sum(loglike_init)): 1216 ↛ 1217line 1216 didn't jump to line 1217 because the condition on line 1216 was never true
1217 for parameter, value in values_init.items():
1218 parameter.value = value
1219 else:
1220 loglike_new = None
1221 return loglike_init, loglike_new
1223 @staticmethod
1224 def make_components_linear(
1225 component_mixture: g2f.ComponentMixture,
1226 ) -> list[g2f.GaussianComponent]:
1227 """Make a list of fixed Gaussian components from a ComponentMixture.
1229 Parameters
1230 ----------
1231 component_mixture
1232 A component mixture to create a component list for.
1234 Returns
1235 -------
1236 gaussians
1237 A list of Gaussians components with fixed parameters and values
1238 matching those in the original component mixture.
1239 """
1240 components = component_mixture.components
1241 if len(components) == 0:
1242 raise ValueError(f"Can't get linear Source from {component_mixture=} with no components")
1243 components_new = [None] * len(components)
1244 for idx, component in enumerate(components):
1245 gaussians = component.gaussians(g2f.Channel.NONE)
1246 # TODO: Support multi-Gaussian components if sensible
1247 # The challenge would be in mapping linear param values back onto
1248 # non-linear IntegralModels
1249 n_g = len(gaussians)
1250 if not n_g == 1:
1251 raise ValueError(f"{component=} has {gaussians=} of len {n_g=}!=1")
1252 gaussian = gaussians.at(0)
1253 component_new = g2f.GaussianComponent(
1254 g2f.GaussianParametricEllipse(
1255 g2f.SigmaXParameterD(gaussian.ellipse.sigma_x, fixed=True),
1256 g2f.SigmaYParameterD(gaussian.ellipse.sigma_y, fixed=True),
1257 g2f.RhoParameterD(gaussian.ellipse.rho, fixed=True),
1258 ),
1259 g2f.CentroidParameters(
1260 g2f.CentroidXParameterD(gaussian.centroid.x, fixed=True),
1261 g2f.CentroidYParameterD(gaussian.centroid.y, fixed=True),
1262 ),
1263 g2f.LinearIntegralModel(
1264 [
1265 (g2f.Channel.NONE, g2f.IntegralParameterD(gaussian.integral.value)),
1266 ]
1267 ),
1268 )
1269 components_new[idx] = component_new
1270 return components_new