Coverage for python / lsst / meas / extensions / multiprofit / fit_coadd_multiband.py: 21%

480 statements  

« prev     ^ index     » next       coverage.py v7.12.0, created at 2025-12-13 01:23 -0800

1# This file is part of meas_extensions_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/>. 

21 

22__all__ = ( 

23 "PsfFitSuccessActionBase", 

24 "PsfComponentsActionBase", 

25 "SourceTablePsfFitSuccessAction", 

26 "SourceTablePsfComponentsAction", 

27 "MagnitudeDependentSizePriorConfig", 

28 "ModelInitializer", 

29 "MakeInitializerActionBase", 

30 "BasicModelInitializer", 

31 "CachedBasicModelInitializer", 

32 "InitialInputData", 

33 "MakeBasicInitializerAction", 

34 "MakeCachedBasicInitializerAction", 

35 "MultiProFitSourceConfig", 

36 "CatalogExposurePsfs", 

37 "MultiProFitSourceFitter", 

38 "MultiProFitSourceTask", 

39) 

40 

41from abc import ABC, abstractmethod 

42from functools import cached_property 

43import logging 

44import math 

45from typing import Any, ClassVar, Iterable, Mapping, Sequence 

46 

47from astropy.table import Table 

48import astropy.units as u 

49import lsst.afw.geom 

50from lsst.daf.butler.formatters.parquet import astropy_to_arrow 

51import lsst.gauss2d as g2 

52import lsst.gauss2d.fit as g2f 

53from lsst.multiprofit.errors import NoDataError, PsfRebuildFitFlagError 

54from lsst.multiprofit.fitting.fit_psf import CatalogPsfFitterConfig, CatalogPsfFitterConfigData 

55from lsst.multiprofit.fitting.fit_source import ( 

56 CatalogExposureSourcesABC, 

57 CatalogSourceFitterABC, 

58 CatalogSourceFitterConfig, 

59 CatalogSourceFitterConfigData, 

60) 

61from lsst.multiprofit.modeller import Model 

62from lsst.multiprofit.utils import frozen_arbitrary_allowed_config, get_params_uniq, set_config_from_dict 

63import lsst.pex.config as pexConfig 

64from lsst.pex.config.configurableActions import ConfigurableAction, ConfigurableActionField 

65import lsst.pipe.base as pipeBase 

66import lsst.pipe.tasks.fit_coadd_multiband as fitMB 

67import lsst.utils.timer as utilsTimer 

68import numpy as np 

69import pydantic 

70 

71from .errors import IsParentError, NotPrimaryError 

72from .input_config import InputConfig 

73from .utils import get_spanned_image 

74 

75_LOG = logging.getLogger(__name__) 

76TWO_SQRT_PI = 2 * math.sqrt(np.pi) 

77 

78 

79class PsfFitSuccessActionBase(ConfigurableAction): 

80 """Base action to return whether a source had a succesful PSF fit.""" 

81 

82 def get_schema(self) -> list[str]: 

83 """Return the list of columns required to call this action.""" 

84 raise NotImplementedError("This method must be overloaded in subclasses") 

85 

86 def __call__(self, source: Mapping[str, Any], *args: Any, **kwargs: Any) -> bool: 

87 raise NotImplementedError("This method must be overloaded in subclasses") 

88 

89 

90class PsfComponentsActionBase(ConfigurableAction): 

91 """Base action to return a list of Gaussians from a source mapping. 

92 

93 This base class should be used as a sentinel when using a MultiProFit PSF 

94 fit table, and only needs to be specialized for external PSF fitters. 

95 """ 

96 

97 def get_schema(self) -> list[str]: 

98 """Return the list of columns required to call this action.""" 

99 raise NotImplementedError("This method must be overloaded in subclasses") 

100 

101 def __call__(self, source: Mapping[str, Any], *args: Any, **kwargs: Any) -> list[g2.Gaussian]: 

102 raise NotImplementedError("This method must be overloaded in subclasses") 

103 

104 

105class SourceTablePsfFitSuccessAction(PsfFitSuccessActionBase): 

106 """Action to return PSF fit status from a SourceTable row.""" 

107 

108 flag_format = pexConfig.Field[str]( 

109 doc="Format for the flag field; flag_prefix, flag_suffix and flag_sub are substituted", 

110 default="{flag_prefix}{flag_suffix}{flag_sub}", 

111 ) 

112 flag_prefix = pexConfig.Field[str]( 

113 doc="Prefix for the key for the summed flag field", 

114 default="modelfit_DoubleShapeletPsfApprox", 

115 ) 

116 flag_suffix = pexConfig.Field[str]( 

117 doc="Suffix for all flag fields", 

118 default="_flag", 

119 ) 

120 flags_sub = pexConfig.ListField[str]( 

121 doc="Suffixes for specific flag fields that must not be true", 

122 default=["_invalidPointForPsf", "_invalidMoments", "_maxIterations"], 

123 ) 

124 

125 def _format(self, flag_sub: str) -> str: 

126 return self.flag_format.format( 

127 flag_prefix=self.flag_prefix, 

128 flag_sub=flag_sub, 

129 flag_suffix=self.flag_suffix, 

130 ) 

131 

132 def get_schema(self) -> Iterable[str]: 

133 for flag_sub in self.flags_sub: 

134 yield self._format(flag_sub=flag_sub) 

135 

136 def __call__(self, source: Mapping[str, Any], *args: Any, **kwargs: Any) -> bool: 

137 good = True 

138 for flag_sub in self.flags_sub: 

139 good &= not source[self._format(flag_sub=flag_sub)] 

140 return good 

141 

142 

143class SourceTablePsfComponentsAction(PsfComponentsActionBase): 

144 """Action to return PSF components from a SourceTable. 

145 

146 This is anticipated to be a deepCoadd_meas with PSF fit parameters from a 

147 measurement plugin returning covariance matrix terms. 

148 """ 

149 

150 action_source = ConfigurableActionField[PsfFitSuccessActionBase]( 

151 doc="Action to return whether the PSF fit was successful for a single source row", 

152 default=SourceTablePsfFitSuccessAction, 

153 ) 

154 format = pexConfig.Field[str]( 

155 doc="Format for the field names, where {idx_comp} is the index of the component and {moment}" 

156 "is the name of the moment (xx, xy or yy, integral)", 

157 default="modelfit_DoubleShapeletPsfApprox_{idx_comp}_{moment}", 

158 ) 

159 name_moment_xx = pexConfig.Field[str](doc="Name of the xx (2nd x-axis) moment", default="xx") 

160 name_moment_xy = pexConfig.Field[str](doc="Name of the xy (covariance term) moment", default="xy") 

161 name_moment_yy = pexConfig.Field[str](doc="Name of the yy (2nd y-axis) moment", default="yy") 

162 name_moment_integral = pexConfig.Field[str](doc="Name of the integral (zeroth) moment", default="0") 

163 n_components = pexConfig.Field[int]( 

164 doc="Number of Gaussian components", 

165 default=2, 

166 check=lambda x: x >= 2, 

167 ) 

168 

169 @staticmethod 

170 def get_integral(moment_zero) -> float: 

171 """Get the total integrated flux from a zeroth moment value. 

172 

173 The zeroth moment is simply the integrated flux divided by a 

174 constant value of 2*sqrt(pi). 

175 

176 Parameters 

177 ---------- 

178 moment_zero 

179 The zeroth moment value. 

180 

181 Returns 

182 ------- 

183 integral 

184 The total integrated weight (flux). 

185 """ 

186 return moment_zero * TWO_SQRT_PI 

187 

188 def get_schema(self) -> list[str]: 

189 names_moments = ( 

190 self.name_moment_xx, 

191 self.name_moment_yy, 

192 self.name_moment_xy, 

193 self.name_moment_integral, 

194 ) 

195 columns = [ 

196 column 

197 for idx_comp in range(self.n_components) 

198 for column in ( 

199 self.format.format(name_moment=name_moment, idx_comp=idx_comp) 

200 for name_moment in names_moments 

201 ) 

202 ] + self.action_source.get_schema() 

203 return columns 

204 

205 def __call__(self, source: Mapping[str, Any], *args: Any, **kwargs: Any) -> list[g2.Gaussian]: 

206 if not self.action_source(source): 

207 raise PsfRebuildFitFlagError( 

208 f"PSF fit failed due to action based on schema: {self.action_source.get_schema()}" 

209 ) 

210 gaussians = [None] * self.n_components 

211 for idx_comp in range(self.n_components): 

212 gaussian = g2.Gaussian( 

213 ellipse=g2.Ellipse( 

214 g2.Covariance( 

215 sigma_x_sq=source[self.format.format(moment=self.name_moment_xx, idx_comp=idx_comp)], 

216 sigma_y_sq=source[self.format.format(moment=self.name_moment_yy, idx_comp=idx_comp)], 

217 cov_xy=source[self.format.format(moment=self.name_moment_xy, idx_comp=idx_comp)], 

218 ) 

219 ), 

220 integral=g2.GaussianIntegralValue( 

221 value=self.get_integral( 

222 source[self.format.format(moment=self.name_moment_integral, idx_comp=idx_comp)] 

223 ) 

224 ), 

225 ) 

226 gaussians[idx_comp] = gaussian 

227 return gaussians 

228 

229 

230class MagnitudeDependentSizePriorConfig(pexConfig.Config): 

231 """Configuration for a magnitude-dependent size prior. 

232 

233 Defaults are for ugrizy total mag and log10(r_eff/arcsec). 

234 """ 

235 

236 intercept_mag = pexConfig.Field[float]( 

237 doc="The magnitude at which no adjustment is applied", 

238 default=18.0, 

239 ) 

240 slope_median_per_mag = pexConfig.Field[float]( 

241 doc="The slope in the median size, in dex per mag", 

242 default=-0.15, 

243 ) 

244 slope_stddev_per_mag = pexConfig.Field[float]( 

245 doc="The slope in the standard deviation of the size, in dex per mag", 

246 default=0.0, 

247 ) 

248 

249 

250class ModelInitializer(ABC, pydantic.BaseModel): 

251 """An interface for a configurable model initializer based on priors 

252 and optional external data. 

253 """ 

254 

255 model_config: ClassVar[pydantic.ConfigDict] = frozen_arbitrary_allowed_config 

256 

257 inputs: dict[str, Any] = pydantic.Field( 

258 title="Additional external inputs used in initialization", 

259 default_factory=dict, 

260 ) 

261 priors_shape_mag: dict = pydantic.Field( 

262 title="Magnitude-dependent shape prior configurations", 

263 default_factory=dict, 

264 ) 

265 

266 @abstractmethod 

267 def initialize_model( 

268 self, 

269 model: Model, 

270 source: Mapping[str, Any], 

271 catexps: list[CatalogExposureSourcesABC], 

272 config_data: CatalogSourceFitterConfigData, 

273 values_init: Mapping[g2f.ParameterD, float] | None = None, 

274 **kwargs, 

275 ): 

276 """Initialize a MultiProFit model for a single object corresponding 

277 to a row in a catalog. 

278 

279 Parameters 

280 ---------- 

281 model 

282 The model to initialize parameter values for. 

283 source 

284 A mapping with fields expected to be populated in the 

285 corresponding source catalog for initialization. 

286 catexps 

287 Per-band catalog-exposure pairs. 

288 config_data 

289 Fitter configuration and data. 

290 values_init 

291 Default initial values for parameters. 

292 **kwargs 

293 Additional keyword arguments for any purpose. 

294 """ 

295 raise NotImplementedError(f"{self.__name__} must implement initialize_model") 

296 

297 

298class MakeInitializerActionBase(ConfigurableAction): 

299 """An interface for an action that creates an initializer.""" 

300 

301 def __call__( 

302 self, 

303 catalog_multi: Sequence, 

304 catexps: list[fitMB.CatalogExposureInputs], 

305 config_data: CatalogSourceFitterConfigData, 

306 **kwargs, 

307 ) -> ModelInitializer: 

308 """Make a ModelInitializer object that can initialize model 

309 parameter values for a given object in a catalog. 

310 

311 Parameters 

312 ---------- 

313 catalog_multi 

314 The multiband catalog with one row per object to fit. 

315 catexps 

316 Per-band catalog-exposure pairs. 

317 config_data 

318 Fitter configuration and data. 

319 **kwargs 

320 Additional arguments to pass to add to ModelInitializer.inputs. 

321 

322 Returns 

323 ------- 

324 initializer 

325 The configured ModelInitializer. 

326 """ 

327 raise NotImplementedError(f"{self.__name__} must implement __call__") 

328 

329 

330class BasicModelInitializer(ModelInitializer): 

331 """A generic model initializer that should work on most kinds of models 

332 with a single source. 

333 """ 

334 

335 def _get_params_init(self, model_sources: tuple[g2f.Source]) -> tuple[g2f.ParameterD]: 

336 """Return an ordered set of free parameters from a model's sources. 

337 

338 Parameters 

339 ---------- 

340 model_sources 

341 The sources in the model. 

342 

343 Returns 

344 ------- 

345 params_init 

346 The parameter objects for sources in the model. 

347 

348 Notes 

349 ----- 

350 Only free and/or centroid parameters are returned (centroids are 

351 always needed even if they are fixed). 

352 """ 

353 # TODO: There ought to be a better way to not get the PSF centroids 

354 # (those are part of model.data's fixed parameters) 

355 params_init = ( 

356 tuple( 

357 ( 

358 param 

359 for param in get_params_uniq(model_sources[0]) 

360 if param.free 

361 or ( 

362 isinstance(param, g2f.CentroidXParameterD) 

363 or isinstance(param, g2f.CentroidYParameterD) 

364 ) 

365 ) 

366 ) 

367 if (len(model_sources) == 1) 

368 else tuple( 

369 { 

370 param: None 

371 for source in model_sources 

372 for param in get_params_uniq(source) 

373 if param.free 

374 or ( 

375 isinstance(param, g2f.CentroidXParameterD) 

376 or isinstance(param, g2f.CentroidYParameterD) 

377 ) 

378 }.keys() 

379 ) 

380 ) 

381 return params_init 

382 

383 def _get_priors_type( 

384 self, 

385 priors: tuple[g2f.Prior], 

386 ) -> tuple[tuple[g2f.GaussianPrior], tuple[g2f.ShapePrior]]: 

387 """Return the list of priors of known type, by type. 

388 

389 Parameters 

390 ---------- 

391 priors 

392 A list of priors of any type, typically from a model. 

393 

394 Returns 

395 ------- 

396 priors_gauss 

397 A list of all of the Gaussian priors, in the order they occurred. 

398 priors_shape 

399 A list of all of the shape priors, in the order they occurred. 

400 """ 

401 priors_gauss: list[g2f.GaussianPrior] = [] 

402 priors_shape: list[g2f.ShapePrior] = [] 

403 for prior in priors: 

404 if isinstance(prior, g2f.GaussianPrior): 

405 priors_gauss.append(prior) 

406 elif isinstance(prior, g2f.ShapePrior): 

407 priors_shape.append(prior) 

408 return tuple(priors_gauss), tuple(priors_shape) 

409 

410 def get_centroid_and_shape( 

411 self, 

412 source: Mapping[str, Any], 

413 catexps: list[CatalogExposureSourcesABC], 

414 config_data: CatalogSourceFitterConfigData, 

415 values_init: Mapping[g2f.ParameterD, float] | None = None, 

416 ) -> tuple[tuple[float, float], tuple[float, float, float]]: 

417 """Get the centroid and shape for a source. 

418 

419 Parameters 

420 ---------- 

421 source 

422 A mapping with fields expected to be populated in the 

423 corresponding source catalog for initialization. 

424 catexps 

425 A list of (source and psf) catalog-exposure pairs. 

426 config_data 

427 Configuration settings and data for fitting and output. 

428 values_init 

429 Initial parameter values from the model configuration. 

430 

431 Returns 

432 ------- 

433 centroid 

434 The x- and y-axis centroid values. 

435 sig_x, sig_y, rho 

436 The x- and y-axis Gaussian sigma and rho values defining the 

437 estimated elliptical shape of the source. 

438 """ 

439 centroid = source["slot_Centroid_x"], source["slot_Centroid_y"] 

440 sig_x = math.sqrt(source["slot_Shape_xx"]) 

441 sig_y = math.sqrt(source["slot_Shape_yy"]) 

442 # TODO: Verify if there's a sign difference here 

443 rho = np.clip(source["slot_Shape_xy"] / (sig_x * sig_y), -0.5, 0.5) 

444 shape = sig_x, sig_y, rho 

445 return centroid, shape 

446 

447 def get_params_init(self, model: Model) -> tuple[g2f.ParameterD]: 

448 """Return the free and/or centroid parameters for a model. 

449 

450 Parameters 

451 ---------- 

452 model 

453 The model to return parameters for. 

454 

455 Returns 

456 ------- 

457 parameters 

458 The ordered list of parameters for the model. 

459 """ 

460 return self._get_params_init(model_sources=model.sources) 

461 

462 def get_priors_type(self, model: Model) -> tuple[tuple[g2f.GaussianPrior], tuple[g2f.ShapePrior]]: 

463 """Return the list of priors of known type, by type. 

464 

465 Parameters 

466 ---------- 

467 model 

468 The model to return priors for. 

469 

470 Returns 

471 ------- 

472 priors_gauss 

473 A list of all of the Gaussian priors, in the order they occurred. 

474 priors_shape 

475 A list of all of the shape priors, in the order they occurred. 

476 """ 

477 return self._get_priors_type(model.priors) 

478 

479 def initialize_model( 

480 self, 

481 model: Model, 

482 source: Mapping[str, Any], 

483 catexps: list[CatalogExposureSourcesABC], 

484 config_data: CatalogSourceFitterConfigData, 

485 values_init: Mapping[g2f.ParameterD, float] | None = None, 

486 **kwargs, 

487 ): 

488 if values_init is None: 

489 values_init = {} 

490 set_flux_limits = kwargs.pop("set_flux_limits", True) 

491 flux_init_min = kwargs.pop("value_init_min", 1e-10) 

492 flux_limit_min = kwargs.pop("flux_limit_min", 1e-12) 

493 if kwargs: 

494 raise ValueError(f"Unexpected {kwargs=}") 

495 centroid_pixel_offset = config_data.config.centroid_pixel_offset 

496 (cen_x, cen_y), (sig_x, sig_y, rho) = self.get_centroid_and_shape( 

497 source, 

498 catexps, 

499 config_data, 

500 values_init=values_init, 

501 ) 

502 # If we couldn't get a shape at all, make it small and roundish 

503 if not np.isfinite(rho): 

504 # Note rho=0 (circular) is generally disfavoured by shape priors 

505 # However, setting it to a non-zero value seems to make scipy 

506 # fail to move off initial conditions, as do sizes below 2 pixels 

507 sig_x, sig_y, rho = 2.0, 2.0, 0.0 

508 

509 # Make restrictive centroid limits (intersection, not union) 

510 x_min, y_min, x_max, y_max = -np.inf, -np.inf, np.inf, np.inf 

511 

512 fluxes_init = {} 

513 fluxes_limits = {} 

514 

515 # This is the maximum number of potential observations 

516 # They might not all have made it into the data 

517 n_catexps = len(catexps) 

518 n_components = len(model.sources[0].components) 

519 

520 # If not true, some bands must have no data to fit 

521 if len(catexps) != len(model.data): 

522 catexps_obs = [] 

523 for catexp in catexps: 

524 fluxes_init[catexp.channel] = flux_init_min 

525 fluxes_limits[catexp.channel] = (0, np.inf) 

526 # No associated catalog means we can't fit (and should be 

527 # because there's no exposure for this band in this patch) 

528 if len(catexp.get_catalog()) > 0: 

529 catexps_obs.append(catexp) 

530 else: 

531 catexps_obs = catexps 

532 

533 for idx_obs, observation in enumerate(model.data): 

534 coordsys = observation.image.coordsys 

535 catexp = catexps_obs[idx_obs] 

536 band = catexp.band 

537 

538 x_min = max(x_min, coordsys.x_min) 

539 y_min = max(y_min, coordsys.y_min) 

540 x_max = min(x_max, coordsys.x_min + float(observation.image.n_cols)) 

541 y_max = min(y_max, coordsys.y_min + float(observation.image.n_rows)) 

542 

543 flux_total = np.nansum(observation.image.data[observation.mask_inv.data]) 

544 

545 column_ref = f"merge_measurement_{band}" 

546 if column_ref in source.schema.getNames() and source[column_ref]: 

547 row = source 

548 else: 

549 row = catexp.catalog.find(source["id"]) 

550 

551 if not row["base_SdssShape_flag"]: 

552 flux_init = row["base_SdssShape_instFlux"] 

553 else: 

554 flux_init = row["slot_GaussianFlux_instFlux"] 

555 if not (flux_init > 0): 

556 flux_init = row["slot_PsfFlux_instFlux"] 

557 

558 calib = catexp.exposure.photoCalib 

559 flux_init = calib.instFluxToNanojansky(flux_init) if (flux_init > 0) else max(flux_total, 1.0) 

560 if set_flux_limits: 

561 flux_max = 10 * max((flux_init, flux_total)) 

562 flux_min = min(flux_limit_min, flux_max / 1000) 

563 else: 

564 flux_min, flux_max = 0, np.inf 

565 if not (flux_init > flux_min): 

566 flux_upper = flux_max if (flux_max < np.inf) else 10.*flux_min 

567 flux_init = flux_min + 0.01*(flux_upper - flux_min) 

568 fluxes_init[observation.channel] = flux_init / n_components 

569 fluxes_limits[observation.channel] = (flux_min, flux_max) 

570 

571 if not np.isfinite(cen_x): 

572 cen_x = observation.image.n_cols / 2.0 

573 else: 

574 cen_x -= centroid_pixel_offset 

575 if not np.isfinite(cen_y): 

576 # TODO: Add bbox coords or remove 

577 cen_y = observation.image.n_rows / 2.0 

578 else: 

579 cen_y -= centroid_pixel_offset 

580 

581 # An R_eff larger than the box size is problematic. This should also 

582 # stop unreasonable size proposals; a log10 transform isn't enough. 

583 # TODO: Try logit for r_eff? 

584 size_major = g2.EllipseMajor(g2.Ellipse(sigma_x=sig_x, sigma_y=sig_y, rho=rho)).r_major 

585 limits_size = max(5.0 * size_major, 2.0 * np.hypot(x_max - x_min, y_max - y_min)) 

586 limits_xy = (1e-5, limits_size) 

587 params_limits_init = { 

588 g2f.CentroidXParameterD: (cen_x, (x_min, x_max)), 

589 g2f.CentroidYParameterD: (cen_y, (y_min, y_max)), 

590 g2f.ReffXParameterD: (sig_x, limits_xy), 

591 g2f.ReffYParameterD: (sig_y, limits_xy), 

592 g2f.SigmaXParameterD: (sig_x, limits_xy), 

593 g2f.SigmaYParameterD: (sig_y, limits_xy), 

594 g2f.RhoParameterD: (rho, None), 

595 # TODO: get guess from configs? 

596 g2f.SersicMixComponentIndexParameterD: (1.0, None), 

597 } 

598 

599 fluxes_init_tuple = tuple(fluxes_init.values()) 

600 fluxes_limits_tuple = tuple(fluxes_limits.values()) 

601 idx_obs = 0 

602 for param in self.params_init: 

603 if param.linear: 

604 value_init = fluxes_init_tuple[idx_obs] 

605 limits_new = fluxes_limits_tuple[idx_obs] 

606 idx_obs += 1 

607 if idx_obs == n_catexps: 

608 idx_obs = 0 

609 else: 

610 type_param = type(param) 

611 value_init, limits_new = params_limits_init.get(type_param, (values_init.get(param), None)) 

612 if limits_new: 

613 param.limits = g2f.LimitsD(limits_new[0], limits_new[1]) 

614 if value_init is not None: 

615 param.value = value_init 

616 

617 priors_shape_mag = self.priors_shape_mag 

618 has_priors_mag = len(priors_shape_mag) > 0 

619 if has_priors_mag: 

620 mag_total = u.nJy.to(u.ABmag, np.nansum(fluxes_init_tuple)) 

621 

622 # TODO: Add centroid prior 

623 priors_gauss, priors_shape = self.get_priors_type(model) 

624 for prior in priors_shape: 

625 if has_priors_mag and ((prior_adjustments := priors_shape_mag.get(prior)) is not None): 

626 mag_dep_prior, prior_shape_new = prior_adjustments 

627 prior_size_new = prior_shape_new.prior_size 

628 # the size-apparent mag relation probably flattens 

629 # for very bright/faint objects - maybe not so 

630 # sharply, but clipping a broad mag range ought to be fine 

631 prior.prior_size.mean_parameter.value = prior_size_new.mean_parameter.value * 10 ** ( 

632 mag_dep_prior.slope_median_per_mag 

633 * np.clip( 

634 mag_total - mag_dep_prior.intercept_mag, 

635 -12.5, 

636 12.5, 

637 ) 

638 ) 

639 # it's uncertain how the intrinsic scatter behaves 

640 # educated guess is it doesn't change much, also 

641 # one runs out of bright galaxies to measure it anyway 

642 prior.prior_size.stddev_parameter.value = prior_size_new.stddev_parameter.value * 10 ** ( 

643 mag_dep_prior.slope_stddev_per_mag 

644 * np.clip( 

645 mag_total - mag_dep_prior.intercept_mag, 

646 -12.5, 

647 12.5, 

648 ) 

649 ) 

650 else: 

651 prior.prior_size.mean_parameter.value = size_major 

652 

653 

654class CachedBasicModelInitializer(BasicModelInitializer): 

655 """A basic initializer with a cached list of model sources and priors.""" 

656 

657 priors: tuple[g2f.Prior, ...] = pydantic.Field(title="The gauss2d_fit model priors") 

658 sources: tuple[g2f.Source, ...] = pydantic.Field(title="The gauss2d_fit model sources") 

659 

660 @cached_property 

661 def params_init(self) -> tuple[g2f.ParameterD]: 

662 """Return a cached reference to the result of _get_params_init.""" 

663 return self._get_params_init(model_sources=self.sources) 

664 

665 @cached_property 

666 def priors_type(self) -> tuple[tuple[g2f.GaussianPrior], tuple[g2f.ShapePrior]]: 

667 """Return a cached reference to the result of _get_priors_type.""" 

668 return self._get_priors_type(self.priors) 

669 

670 def get_params_init(self, model: Model) -> tuple[g2f.ParameterD]: 

671 assert tuple(model.sources) == self.sources 

672 return self.params_init 

673 

674 def get_priors_type(self, model: Model) -> tuple[tuple[g2f.GaussianPrior], tuple[g2f.ShapePrior]]: 

675 assert tuple(model.priors) == self.priors 

676 return self.priors_type 

677 

678 

679class InitialInputData(pydantic.BaseModel): 

680 """A configurable wrapper to retrieve formatted columns from a catalog. 

681 

682 This provides a common interface to typical MultiProFit table outputs. 

683 """ 

684 

685 model_config: ClassVar[pydantic.ConfigDict] = frozen_arbitrary_allowed_config 

686 

687 column_id: str | None = pydantic.Field( 

688 title="Override for id column specified in config_input", 

689 default=None, 

690 ) 

691 config_input: InputConfig = pydantic.Field(title="Configuration for the data table") 

692 data: Table = pydantic.Field(title="The data table") 

693 name_model: str = pydantic.Field(title="The name of the model in columns") 

694 prefix_column: str = pydantic.Field(title="The prefix for all fitted column names") 

695 size_column: str = pydantic.Field(title="The name of the size column", default="reff") 

696 

697 def get_column_id(self): 

698 """Return the name of the object ID column.""" 

699 return self.column_id or self.config_input.column_id 

700 

701 def get_column(self, name_column: str, data=None): 

702 """Get the values from a column. 

703 

704 Parameters 

705 ---------- 

706 name_column 

707 The name of the column to retrieve. 

708 data 

709 The catalog to retrieve the column from. Default is self.data. 

710 

711 Returns 

712 ------- 

713 values 

714 The column values. 

715 """ 

716 if data is None: 

717 data = self.data 

718 return data[f"{self.prefix_column}{name_column}"] 

719 

720 def model_post_init(self, __context: Any) -> None: 

721 # Initialize a mapping of the row number for a given object ID value 

722 # This is implemented in afw catalogs but not most other tabular types 

723 id_index = {idnum: idx for idx, idnum in enumerate(self.data[self.get_column_id()])} 

724 object.__setattr__(self, "id_index", id_index) 

725 

726 

727class MakeBasicInitializerAction(MakeInitializerActionBase): 

728 """An action to construct an initializer for a single-component, 

729 single-source model. 

730 """ 

731 

732 def _make_initializer( 

733 self, 

734 catalog_multi: Sequence, 

735 catexps: list[fitMB.CatalogExposureInputs], 

736 config_data: CatalogSourceFitterConfigData, 

737 ) -> ModelInitializer: 

738 return BasicModelInitializer() 

739 

740 def __call__( 

741 self, 

742 catalog_multi: Sequence, 

743 catexps: list[fitMB.CatalogExposureInputs], 

744 config_data: CatalogSourceFitterConfigData, 

745 **kwargs, 

746 ) -> ModelInitializer: 

747 initializer = self._make_initializer( 

748 catalog_multi=catalog_multi, 

749 catexps=catexps, 

750 config_data=config_data, 

751 ) 

752 for name, (config_input, data) in kwargs.items(): 

753 if not isinstance(data, Table) and hasattr(data, "meta"): 

754 _LOG.warning( 

755 f"Ignoring extra input {name=} because it is of type {type(data)} and is either not an" 

756 f" astropy.table.Table or missing a 'meta' attr" 

757 ) 

758 config_data = data.meta["config"] 

759 prefix_column = config_data["prefix_column"] 

760 config_source = next(iter(config_data["config_model"]["sources"].values())) 

761 config_group = next(iter(config_source["component_groups"].values())) 

762 is_sersic = len(config_group["components_sersic"]) > 0 

763 name_model = next( 

764 iter(config_group["components_sersic"] if is_sersic else config_group["components_gaussian"]) 

765 ) 

766 initializer.inputs[name] = InitialInputData( 

767 column_id=config_data.get("column_id"), 

768 config_input=config_input, 

769 data=data, 

770 name_model=name_model, 

771 prefix_column=prefix_column, 

772 size_column="reff" if is_sersic else "sig", 

773 ) 

774 return initializer 

775 

776 

777class MakeCachedBasicInitializerAction(MakeBasicInitializerAction): 

778 """A MakeBasicInitializerAction that caches references to the source 

779 and prior objects of the model. 

780 

781 This is solely a performance optimization and should be favored over 

782 MakeBasicInitializerAction unless the caching is shown to be slower. 

783 """ 

784 

785 def _make_initializer( 

786 self, 

787 catalog_multi: Sequence, 

788 catexps: list[fitMB.CatalogExposureInputs], 

789 config_data: CatalogSourceFitterConfigData, 

790 ) -> ModelInitializer: 

791 sources, priors = config_data.sources_priors 

792 return CachedBasicModelInitializer(priors=priors, sources=sources) 

793 

794 

795class MultiProFitSourceConfig(CatalogSourceFitterConfig, fitMB.CoaddMultibandFitSubConfig): 

796 """Configuration for the MultiProFit profile fitter.""" 

797 

798 action_initializer = ConfigurableActionField[MakeInitializerActionBase]( 

799 doc="The action to return an initializer", 

800 default=MakeCachedBasicInitializerAction, 

801 ) 

802 action_psf = ConfigurableActionField[PsfComponentsActionBase]( 

803 doc="The action to return PSF component values from catalogs, if implemented", 

804 default=None, 

805 ) 

806 columns_copy = pexConfig.DictField[str, str]( 

807 doc="Mapping of input/output column names to copy from the input" 

808 "multiband catalog to the output fit catalog.", 

809 default={}, 

810 dictCheck=lambda x: len(set(x.values())) == len(x.values()), 

811 ) 

812 mask_names_zero = pexConfig.ListField[str]( 

813 doc="Mask bits to mask out", 

814 default=["BAD", "EDGE", "SAT", "NO_DATA"], 

815 ) 

816 psf_sigma_subtract = pexConfig.Field[float]( 

817 doc="PSF x/y sigma value to subtract in quadrature from best-fit values", 

818 default=0.1, 

819 check=lambda x: np.isfinite(x) and (x >= 0), 

820 ) 

821 prefix_column = pexConfig.Field[str](default="mpf_", doc="Column name prefix") 

822 size_priors = pexConfig.ConfigDictField[str, MagnitudeDependentSizePriorConfig]( 

823 doc="Per-component magnitude-dependent size prior configurations." 

824 " Will be added to component with existing configs.", 

825 default={}, 

826 ) 

827 

828 def bands_read_only(self) -> set[str]: 

829 # TODO: Re-implement determination of prior-only bands once 

830 # data-driven priors are re-implemented (DM-4xxxx) 

831 return set() 

832 

833 def requires_psf(self): 

834 """Return whether the PSF action is not None.""" 

835 return type(self.action_psf) is PsfComponentsActionBase 

836 

837 def setDefaults(self): 

838 super().setDefaults() 

839 self.flag_errors = { 

840 IsParentError.column_name(): "IsParentError", 

841 NoDataError.column_name(): "NoDataError", 

842 NotPrimaryError.column_name(): "NotPrimaryError", 

843 PsfRebuildFitFlagError.column_name(): "PsfRebuildFitFlagError", 

844 } 

845 self.centroid_pixel_offset = -0.5 

846 self.naming_scheme = "lsst" 

847 self.prefix_column = "" 

848 self.suffix_error = "Err" 

849 

850 

851@pydantic.dataclasses.dataclass(frozen=True, kw_only=True, config=fitMB.CatalogExposureConfig) 

852class CatalogExposurePsfs(fitMB.CatalogExposureInputs, CatalogExposureSourcesABC): 

853 """Input data from lsst pipelines, parsed for MultiProFit.""" 

854 

855 channel: g2f.Channel = pydantic.Field(title="Channel for the image's band") 

856 config_fit: MultiProFitSourceConfig = pydantic.Field(title="Config for fitting options") 

857 

858 @cached_property 

859 def _psf_flux_params(self) -> tuple[list[g2f.ParameterD], bool]: 

860 psf_model = self.psf_model_data.psf_model 

861 n_comps = len(psf_model.components) 

862 params_flux = [None] * n_comps 

863 is_frac = [False] * n_comps 

864 for idx_comp, comp in enumerate(psf_model.components): 

865 # TODO: Change to comp.integralmodel when DM-44344 is fixed 

866 # integralmodels will still need to be handled differently 

867 params_all = get_params_uniq(comp) 

868 params_frac = [param for param in params_all if isinstance(param, g2f.ProperFractionParameterD)] 

869 if params_frac: 

870 is_last = idx_comp == (n_comps - 1) 

871 if len(params_frac) != (idx_comp + 1 - is_last): 

872 raise RuntimeError( 

873 f"Got unexpected {params_frac=} for" 

874 f" {self.psf_model_data.psf_model.components[idx_comp]=} ({idx_comp=});" 

875 f" len should be idx_comp+1" 

876 ) 

877 params_flux[idx_comp] = None if is_last else params_frac[idx_comp] 

878 is_frac[idx_comp] = True 

879 else: 

880 params_integral = [param for param in params_all if isinstance(param, g2f.IntegralParameterD)] 

881 if len(params_integral != 1): 

882 raise RuntimeError( 

883 f"Got unexpected {params_integral=} != 1 for" 

884 f" {self.psf_model_data.psf_model.components[idx_comp]=} ({idx_comp=})" 

885 ) 

886 params_flux[idx_comp] = params_integral[0] 

887 is_frac_any = any(is_frac) 

888 if is_frac_any and not all(is_frac): 

889 # TODO: This should work by iterating through componentgroups 

890 # But that's not trivial or supported now 

891 raise RuntimeError("Got PSF model with a mix of fractional and linear models; cannot initialize") 

892 

893 return params_flux, is_frac_any 

894 

895 def get_psf_model(self, source): 

896 psf_model = self.psf_model_data.psf_model 

897 # PsfComponentsActionBase is an abstract class, so check if the action 

898 # is a subclass that needs to be called 

899 if not self.config_fit.requires_psf(): 

900 try: 

901 gaussians = self.config_fit.action_psf(source) 

902 except PsfRebuildFitFlagError: 

903 return None 

904 n_comps = len(psf_model.components) 

905 fluxes = [0.0] * n_comps 

906 params_flux, is_frac = self._psf_flux_params 

907 for idx_comp, (comp, gaussian) in enumerate(zip(psf_model.components, gaussians)): 

908 ellipse_out = comp.ellipse 

909 ellipse_in = gaussian.ellipse 

910 ellipse_out.sigma_x = ellipse_in.sigma_x 

911 ellipse_out.sigma_y = ellipse_in.sigma_y 

912 ellipse_out.rho = ellipse_in.rho 

913 fluxes[idx_comp] = gaussian.integral.value 

914 # Apparently negative fluxes are possible. Not much can be done to 

915 # fix that but set them to a tiny value (zero might work) 

916 fluxes = np.clip(fluxes, 1e-3, np.inf) 

917 flux_total = sum(fluxes) 

918 if is_frac: 

919 flux_remaining = 1.0 

920 for flux, param_frac in zip(fluxes, params_flux[:-1]): 

921 flux_component = flux / flux_total 

922 param_frac.value = flux_component / flux_remaining 

923 flux_remaining -= flux_component 

924 else: 

925 for flux, param_flux in zip(fluxes, params_flux): 

926 param_flux.value = flux / flux_total 

927 else: 

928 # TODO: this should probably use .index or something 

929 match = np.argwhere( 

930 self.table_psf_fits[self.psf_model_data.config.column_id] == source[self.config_fit.column_id] 

931 )[0][0] 

932 psf_model = self.psf_model_data.psf_model 

933 self.psf_model_data.init_psf_model(self.table_psf_fits[match]) 

934 

935 sigma_subtract = self.config_fit.psf_sigma_subtract 

936 if sigma_subtract > 0: 

937 sigma_subtract_sq = sigma_subtract * sigma_subtract 

938 for param in self.psf_model_data.parameters.values(): 

939 if isinstance( 

940 param, 

941 g2f.SigmaXParameterD | g2f.SigmaYParameterD | g2f.ReffXParameterD | g2f.ReffYParameterD, 

942 ): 

943 param.value = math.sqrt(param.value**2 - sigma_subtract_sq) 

944 return psf_model 

945 

946 def get_source_observation(self, source, **kwargs) -> g2f.ObservationD: 

947 if not kwargs.get("skip_flags"): 

948 if (not source["detect_isPrimary"]) or source["merge_peak_sky"]: 

949 raise NotPrimaryError(f"source {source[self.config_fit.column_id]} has invalid flags for fit") 

950 footprint = source.getFootprint() 

951 bbox = footprint.getBBox() 

952 if not (bbox.getArea() > 0): 

953 return None 

954 bitmask = 0 

955 mask = self.exposure.mask[bbox] 

956 spans = footprint.spans.asArray() 

957 for bitname in self.config_fit.mask_names_zero: 

958 bitval = mask.getPlaneBitMask(bitname) 

959 bitmask |= bitval 

960 mask = ((mask.array & bitmask) != 0) & (spans != 0) 

961 mask = ~mask 

962 

963 is_deblended_child = source["parent"] != 0 

964 

965 img, _, sigma_inv = get_spanned_image( 

966 exposure=self.exposure, 

967 footprint=footprint if is_deblended_child else None, 

968 bbox=bbox, 

969 spans=spans, 

970 get_sig_inv=True, 

971 ) 

972 x_min_bbox, y_min_bbox = bbox.beginX, bbox.beginY 

973 # Crop to tighter box for deblended model if edges are unusable 

974 # ... this rarely ever seems to happen though 

975 if is_deblended_child: 

976 coords = np.argwhere(np.isfinite(img) & (sigma_inv > 0) & np.isfinite(sigma_inv)) 

977 x_min, y_min = coords.min(axis=0) 

978 x_max, y_max = coords.max(axis=0) 

979 x_max += 1 

980 y_max += 1 

981 

982 if (x_min > 0) or (y_min > 0) or (x_max < img.shape[0]) or (y_max < img.shape[1]): 

983 # Ensure the nominal centroid is still inside the box 

984 # ... although it's a bad sign if that row/column is all bad 

985 x_cen = source["slot_Centroid_x"] - x_min_bbox 

986 y_cen = source["slot_Centroid_y"] - y_min_bbox 

987 x_min = min(x_min, int(np.floor(x_cen))) 

988 x_max = max(x_max, int(np.ceil(x_cen))) 

989 y_min = min(y_min, int(np.floor(y_cen))) 

990 y_max = max(y_max, int(np.ceil(y_cen))) 

991 x_min_bbox += x_min 

992 y_min_bbox += y_min 

993 img = img[x_min:x_max, y_min:y_max] 

994 sigma_inv = sigma_inv[x_min:x_max, y_min:y_max] 

995 mask = mask[x_min:x_max, y_min:y_max] 

996 

997 sigma_inv[~mask] = 0 

998 

999 coordsys = g2.CoordinateSystem(1.0, 1.0, x_min_bbox, y_min_bbox) 

1000 

1001 obs = g2f.ObservationD( 

1002 image=g2.ImageD(img, coordsys), 

1003 sigma_inv=g2.ImageD(sigma_inv, coordsys), 

1004 mask_inv=g2.ImageB(mask, coordsys), 

1005 channel=self.channel, 

1006 ) 

1007 return obs 

1008 

1009 def __post_init__(self): 

1010 # TODO: Can/should this be the derived type (MultiProFitPsfConfig)? 

1011 config = CatalogPsfFitterConfig() 

1012 config_dict = self.table_psf_fits.meta.get("config") 

1013 if config_dict: 

1014 set_config_from_dict(config, config_dict) 

1015 else: 

1016 # TODO: How should this be set? 

1017 # If using external PSF fits, it needs to be configured normally 

1018 pass 

1019 config_data = CatalogPsfFitterConfigData(config=config) 

1020 object.__setattr__(self, "psf_model_data", config_data) 

1021 

1022 

1023class MultiProFitSourceFitter(CatalogSourceFitterABC): 

1024 """A MultiProFit source fitter. 

1025 

1026 Parameters 

1027 ---------- 

1028 wcs 

1029 A WCS solution that applies to all exposures. 

1030 errors_expected 

1031 A dictionary of exceptions that are expected to sometimes be raised 

1032 during processing (e.g. for missing data) keyed by the name of the 

1033 flag column used to record the failure. 

1034 add_missing_errors 

1035 Whether to add all of the standard MultiProFit errors with default 

1036 column names to errors_expected, if not already present. 

1037 **kwargs 

1038 Keyword arguments to pass to the superclass constructor. 

1039 """ 

1040 

1041 initializer: ModelInitializer = pydantic.Field( 

1042 title="The model parameter initializer", 

1043 default_factory=lambda: BasicModelInitializer(), 

1044 ) 

1045 wcs: lsst.afw.geom.SkyWcs = pydantic.Field( 

1046 title="The WCS object to use to convert pixel coordinates to RA/dec", 

1047 ) 

1048 

1049 def __init__( 

1050 self, 

1051 wcs: lsst.afw.geom.SkyWcs, 

1052 errors_expected: dict[str, Exception] | None = None, 

1053 add_missing_errors: bool = True, 

1054 **kwargs: Any, 

1055 ): 

1056 if errors_expected is None: 

1057 errors_expected = {} 

1058 if add_missing_errors: 

1059 for error_catalog in (IsParentError, NoDataError, NotPrimaryError, PsfRebuildFitFlagError): 

1060 if error_catalog not in errors_expected: 

1061 errors_expected[error_catalog] = error_catalog.column_name() 

1062 super().__init__(wcs=wcs, errors_expected=errors_expected, **kwargs) 

1063 

1064 def copy_centroid_errors( 

1065 self, 

1066 columns_cenx_err_copy: tuple[str], 

1067 columns_ceny_err_copy: tuple[str], 

1068 results: Table, 

1069 catalog_multi: Sequence, 

1070 catexps: list[CatalogExposureSourcesABC], 

1071 config_data: CatalogSourceFitterConfigData, 

1072 ): 

1073 for column in columns_cenx_err_copy: 

1074 results[column] = catalog_multi["slot_Centroid_xErr"] 

1075 for column in columns_ceny_err_copy: 

1076 results[column] = catalog_multi["slot_Centroid_yErr"] 

1077 

1078 def get_model_radec(self, source: Mapping[str, Any], cen_x: float, cen_y: float): 

1079 # no extra conversions are needed here - cen_x, cen_y are in catalog 

1080 # coordinates already 

1081 ra, dec = self.wcs.pixelToSky(cen_x, cen_y) 

1082 return ra.asDegrees(), dec.asDegrees() 

1083 

1084 def initialize_model( 

1085 self, 

1086 model: g2f.ModelD, 

1087 source: Mapping[str, Any], 

1088 catexps: list[CatalogExposureSourcesABC], 

1089 config_data: CatalogSourceFitterConfigData, 

1090 values_init: Mapping[g2f.ParameterD, float] | None = None, 

1091 **kwargs, 

1092 ): 

1093 self.initializer.initialize_model( 

1094 model=model, 

1095 source=source, 

1096 catexps=catexps, 

1097 config_data=config_data, 

1098 values_init=values_init, 

1099 **kwargs, 

1100 ) 

1101 

1102 def make_CatalogExposurePsfs( 

1103 self, 

1104 catexp: fitMB.CatalogExposureInputs, 

1105 config: MultiProFitSourceConfig, 

1106 ) -> CatalogExposurePsfs: 

1107 """Make a CatalogExposurePsfs from a list of inputs and a fit config. 

1108 

1109 Parameters 

1110 ---------- 

1111 catexp 

1112 The input catalog-exposure pairs. 

1113 config 

1114 The MultiProFit source fitting config. 

1115 

1116 Returns 

1117 ------- 

1118 catexp_psf 

1119 The resulting CatalogExposurePsfs. 

1120 """ 

1121 catexp_psf = CatalogExposurePsfs( 

1122 # dataclasses.asdict(catexp)_makes a recursive deep copy. 

1123 # That must be avoided. 

1124 **{key: getattr(catexp, key) for key in catexp.__dataclass_fields__.keys()}, 

1125 channel=g2f.Channel.get(catexp.band), 

1126 config_fit=config, 

1127 ) 

1128 return catexp_psf 

1129 

1130 def validate_fit_inputs( 

1131 self, 

1132 catalog_multi: Sequence, 

1133 catexps: list[CatalogExposurePsfs], 

1134 config_data: CatalogSourceFitterConfigData = None, 

1135 logger: logging.Logger = None, 

1136 **kwargs: Any, 

1137 ) -> None: 

1138 errors = [] 

1139 for idx, catexp in enumerate(catexps): 

1140 if not isinstance(catexp, CatalogExposurePsfs): 

1141 errors.append(f"catexps[{idx=} {type(catexp)=} !isinstance(CatalogExposurePsfs)") 

1142 # Pre-validate the model 

1143 config_sources = config_data.config.config_model.sources 

1144 model_sources, priors = config_data.sources_priors 

1145 priors_shape = [prior for prior in priors if isinstance(prior, g2f.ShapePrior)] 

1146 

1147 if len(config_sources.keys()) > 1: 

1148 errors.append(f"model config has multiple sources: {list(config_sources.keys())}") 

1149 elif len(priors_shape) > 0: 

1150 idx_prior_found = 0 

1151 name_source, config_source = next(iter(config_sources.items())) 

1152 source = model_sources[0] 

1153 config_groups = config_source.component_groups 

1154 if len(config_groups.keys()) > 1: 

1155 errors.append(f"model {name_source=} has multiple groups: {list(config_source.keys())}") 

1156 else: 

1157 name_group, config_group = next(iter(config_groups.items())) 

1158 for idx_comp, (name_comp, config_comp) in enumerate( 

1159 config_group.get_component_configs().items() 

1160 ): 

1161 ellipse = source.components[idx_comp].ellipse 

1162 # component.ellipse returns a const ref and must be copied 

1163 # The ellipse classes might need copy constructors 

1164 ellipse_copy = type(ellipse)( 

1165 size_x=ellipse.size_x, 

1166 size_y=ellipse.size_y, 

1167 rho=ellipse.rho, 

1168 ) 

1169 prior_shape_new = config_comp.make_shape_prior(ellipse_copy) 

1170 if prior_shape_new is not None: 

1171 if idx_prior_found == len(priors_shape): 

1172 errors.append( 

1173 f"Could not validate prior for {name_source=} {name_group=} {name_comp=}" 

1174 ) 

1175 break 

1176 prior_shape_old = priors_shape[idx_prior_found] 

1177 ll_new, ll_old = ( 

1178 prior.evaluate().loglike for prior in (prior_shape_new, prior_shape_old) 

1179 ) 

1180 # The necessary tolerance for this check is uncertain 

1181 if not np.isclose(ll_new, ll_old): 

1182 logger.warning( 

1183 f"shape prior for {name_comp=} got inconsistent {ll_new=} vs {ll_old}" 

1184 ) 

1185 if (prior_shape_mod := config_data.config.size_priors.get(name_comp)) is not None: 

1186 self.initializer.priors_shape_mag[prior_shape_old] = ( 

1187 prior_shape_mod, 

1188 prior_shape_new, 

1189 ) 

1190 

1191 if errors: 

1192 raise RuntimeError("\n".join(errors)) 

1193 

1194 

1195class MultiProFitSourceTask(fitMB.CoaddMultibandFitSubTask): 

1196 """Run MultiProFit on Exposure/SourceCatalog pairs in multiple bands. 

1197 

1198 This task uses MultiProFit to fit a single model to all sources in a coadd, 

1199 using a previously-fit PSF model for each exposure. The task may also use 

1200 prior measurements from single- or merged multiband catalogs for 

1201 initialization. 

1202 """ 

1203 

1204 ConfigClass: ClassVar = MultiProFitSourceConfig 

1205 _DefaultName: ClassVar = "multiProFitSource" 

1206 

1207 def make_default_fitter( 

1208 self, 

1209 catalog_multi: Sequence, 

1210 catexps: list[fitMB.CatalogExposureInputs], 

1211 config_data: CatalogSourceFitterConfigData, 

1212 **kwargs, 

1213 ) -> MultiProFitSourceFitter: 

1214 """Make a default MultiProFitSourceFitter. 

1215 

1216 Parameters 

1217 ---------- 

1218 catalog_multi 

1219 A multi-band, indexable source catalog. 

1220 catexps 

1221 Catalog-exposure-PSF model tuples to fit source models for. 

1222 config_data 

1223 Configuration and data for the initalizer. 

1224 **kwargs 

1225 Additional keyword arguments to pass to 

1226 self.config.action_initializer. 

1227 

1228 Returns 

1229 ------- 

1230 fitter 

1231 A MultiProFitSourceFitter using the first catexp's wcs. 

1232 """ 

1233 initializer = self.config.action_initializer( 

1234 catalog_multi=catalog_multi, catexps=catexps, config_data=config_data, **kwargs 

1235 ) 

1236 # Look for the first WCS - they ought to be identical 

1237 # If they are not, the patch coadd data model must have changed 

1238 wcs = None 

1239 for catexp in catexps: 

1240 if catexp.exposure is not None: 

1241 wcs = catexp.exposure.wcs 

1242 break 

1243 if wcs is None: 

1244 raise RuntimeError(f"Could not find valid wcs in any of {catexps=}") 

1245 fitter = MultiProFitSourceFitter(wcs=wcs, initializer=initializer) 

1246 return fitter 

1247 

1248 @utilsTimer.timeMethod 

1249 def run( 

1250 self, 

1251 catalog_multi: Sequence, 

1252 catexps: list[fitMB.CatalogExposureInputs], 

1253 fitter: MultiProFitSourceFitter | None = None, 

1254 **kwargs, 

1255 ) -> pipeBase.Struct: 

1256 """Run the MultiProFit source fit task on catalog-exposure pairs. 

1257 

1258 Parameters 

1259 ---------- 

1260 catalog_multi 

1261 A multi-band, indexable source catalog. 

1262 catexps 

1263 Catalog-exposure-PSF model tuples to fit source models for. 

1264 fitter 

1265 The fitter instance to use. Default-initialized if not provided. 

1266 **kwargs 

1267 Additional keyword arguments to pass to self.fit. 

1268 

1269 Returns 

1270 ------- 

1271 catalog : `astropy.Table` 

1272 A table with fit parameters for the PSF model at the location 

1273 of each source. 

1274 """ 

1275 n_catexps = len(catexps) 

1276 if n_catexps == 0: 

1277 raise ValueError("Must provide at least one catexp") 

1278 catexps_conv: list[CatalogExposurePsfs] = [None] * n_catexps 

1279 channels = [g2f.Channel.get(catexp.band) for catexp in catexps] 

1280 config_data = CatalogSourceFitterConfigData(channels=channels, config=self.config) 

1281 if fitter is None: 

1282 inputs_init = kwargs.get("inputs_init") 

1283 if inputs_init: 

1284 del kwargs["inputs_init"] 

1285 else: 

1286 inputs_init = {} 

1287 fitter = self.make_default_fitter( 

1288 catalog_multi=catalog_multi, catexps=catexps, config_data=config_data, **inputs_init 

1289 ) 

1290 for idx, catexp in enumerate(catexps): 

1291 if not isinstance(catexp, CatalogExposurePsfs): 

1292 catexp = fitter.make_CatalogExposurePsfs(catexp, config=self.config) 

1293 catexps_conv[idx] = catexp 

1294 catalog = fitter.fit( 

1295 catalog_multi=catalog_multi, catexps=catexps_conv, config_data=config_data, **kwargs 

1296 ) 

1297 for name_in, name_out in self.config.columns_copy.items(): 

1298 catalog[name_out] = catalog_multi[name_in] 

1299 catalog[name_out].description = catalog_multi.schema.find(name_in).field.getDoc() 

1300 return pipeBase.Struct(output=astropy_to_arrow(catalog))