Coverage for python / lsst / multiprofit / fitting / fit_catalog.py: 50%

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/>. 

21 

22__all__ = ["CatalogExposureABC", "ColumnInfo", "CatalogFitterConfig"] 

23 

24from abc import ABC, abstractmethod 

25from collections.abc import Iterable 

26from typing import ClassVar 

27 

28from astropy.table import Table 

29import astropy.units as u 

30import lsst.pex.config as pexConfig 

31import numpy as np 

32import pydantic 

33 

34from ..componentconfig import GaussianComponentConfig, SersicComponentConfig 

35from ..modeller import ModelFitConfig 

36from ..utils import frozen_arbitrary_allowed_config 

37 

38 

39class CatalogExposureABC(ABC): 

40 """Interface for catalog-exposure pairs.""" 

41 

42 # TODO: add get_exposure (with Any return type?) 

43 

44 @abstractmethod 

45 def get_catalog(self) -> Iterable: 

46 """Return a row-iterable catalog covering an exposure.""" 

47 

48 

49class ColumnInfo(pydantic.BaseModel): 

50 """Metadata for a column in a catalog.""" 

51 

52 model_config: ClassVar[pydantic.ConfigDict] = frozen_arbitrary_allowed_config 

53 

54 dtype: str = pydantic.Field(title="Column data type name (numpy or otherwise)") 

55 key: str = pydantic.Field(title="Column key (name)") 

56 description: str = pydantic.Field("", title="Column description") 

57 unit: u.UnitBase | None = pydantic.Field(None, title="Column unit (astropy)") 

58 

59 

60class CatalogFitterConfig(pexConfig.Config): 

61 """Configuration for generic MultiProFit fitting tasks.""" 

62 

63 column_id = pexConfig.Field[str](default="id", doc="Catalog index column key") 

64 compute_errors = pexConfig.ChoiceField[str]( 

65 default="INV_HESSIAN_BESTFIT", 

66 doc="Whether/how to compute sqrt(variances) of each free parameter", 

67 allowed={ 

68 "NONE": "no errors computed", 

69 "INV_HESSIAN": "inverse hessian using noisy image as data", 

70 "INV_HESSIAN_BESTFIT": "inverse hessian using best-fit model as data", 

71 }, 

72 ) 

73 compute_errors_from_jacobian = pexConfig.Field[bool]( 

74 default=True, 

75 doc="Whether to estimate the Hessian from the Jacobian first, with finite differencing as a backup", 

76 ) 

77 compute_errors_no_covar = pexConfig.Field[bool]( 

78 default=True, 

79 doc="Whether to compute parameter errors independently, ignoring covariances", 

80 ) 

81 config_fit = pexConfig.ConfigField[ModelFitConfig](default=ModelFitConfig, doc="Fitter configuration") 

82 fit_centroid = pexConfig.Field[bool](default=True, doc="Fit centroid parameters") 

83 fit_linear_init = pexConfig.Field[bool](default=True, doc="Fit linear parameters after initialization") 

84 fit_linear_final = pexConfig.Field[bool](default=True, doc="Fit linear parameters after optimization") 

85 float_fill_value = pexConfig.Field[float]( 

86 default=np.nan, doc="Fill value for float fields when creating the output table." 

87 ) 

88 flag_errors = pexConfig.DictField( 

89 default={}, 

90 keytype=str, 

91 itemtype=str, 

92 doc="Flag column names to set, keyed by name of exception to catch", 

93 ) 

94 integer_fill_value = pexConfig.Field[int]( 

95 default=-1, doc="Fill value for integer fields when creating the output table." 

96 ) 

97 naming_scheme = pexConfig.ChoiceField[str]( 

98 doc="Naming scheme for column names", 

99 allowed={ 

100 "default": "snake_case with {component_name}[_{band}]_{parameter}[_err]", 

101 "camel": "CamelCase with {component_name}[_{band}]_{parameter}[Err]", 

102 "lsst": "snake_case with [{band}_]{component_name}_{parameter}[Err]", 

103 }, 

104 default="default", 

105 ) 

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

107 suffix_error = pexConfig.Field[str]( 

108 default="_err", 

109 doc="Default suffix for error columns. Can be overridden by naming_scheme.", 

110 ) 

111 

112 _format_flux = { 

113 "default": "{label}{band}_flux", 

114 "lsst": "{band}_{label}Flux", 

115 "camel": "{label}{band}Flux", 

116 } 

117 _key_cen = {"default": "_cen", "lsst": "_cen", "camel": "Cen"} 

118 _key_reff = {"default": f"_{SersicComponentConfig._size_label}", "lsst": "_reff", "camel": "Reff"} 

119 _key_rho = {"default": "_rho", "lsst": "_rho", "camel": "Rho"} 

120 _key_sigma = {"default": f"_{GaussianComponentConfig._size_label}", "lsst": "_sigma", "camel": "Sigma"} 

121 _key_sersicindex = {"default": "_sersic_index", "lsst": "sersic_index", "camel": "SersicIndex"} 

122 _suffix_dec = {"default": "_dec", "lsst": "_dec", "camel": "Dec"} 

123 _suffix_ra = {"default": "_ra", "lsst": "_ra", "camel": "Ra"} 

124 _suffix_x = {"default": "_x", "lsst": "_x", "camel": "X"} 

125 _suffix_y = {"default": "_y", "lsst": "_y", "camel": "Y"} 

126 

127 def _get_label(self, format_name: str, values: dict[str, str]) -> str: 

128 """Get the label for part of a column name for a given format. 

129 

130 Parameters 

131 ---------- 

132 format_name 

133 The name of the format to get the label for. 

134 values 

135 The values of the name by format. 

136 

137 Returns 

138 ------- 

139 label 

140 The formatted label, if specified for that format, else the 

141 value for the default format. 

142 """ 

143 return values.get(format_name, values["default"]) 

144 

145 def get_key_cen(self) -> str: 

146 """Get the key for centroid columns.""" 

147 return self._get_label(self.naming_scheme, self._key_cen) 

148 

149 def get_key_flux(self, band: str, label: str = "") -> str: 

150 """Get the key for a flux column. 

151 

152 Parameters 

153 ---------- 

154 band 

155 The band of the flux column. 

156 label 

157 A label for this flux, e.g. a component name. 

158 

159 Returns 

160 ------- 

161 key_flux 

162 The flux column key. 

163 """ 

164 return self._get_label(self.naming_scheme, self._format_flux).format(band=band, label=label) 

165 

166 def get_key_reff(self) -> str: 

167 """Get the key for Sersic effective radius columns.""" 

168 return self._get_label(self.naming_scheme, self._key_reff) 

169 

170 def get_key_rho(self) -> str: 

171 """Get the key for ellipse rho columns.""" 

172 return self._get_label(self.naming_scheme, self._key_rho) 

173 

174 def get_key_sersicindex(self) -> str: 

175 """Get the key for Sersic index columns.""" 

176 return self._get_label(self.naming_scheme, self._key_sersicindex) 

177 

178 def get_key_sigma(self) -> str: 

179 """Get the key for Gaussian sigma columns.""" 

180 return self._get_label(self.naming_scheme, self._key_sigma) 

181 

182 def get_key_size(self, label_size: str) -> str: 

183 """Get the key for a size column by its label. 

184 

185 Parameters 

186 ---------- 

187 label_size 

188 The label of the size, usually specified in a ComponentConfig. 

189 

190 Returns 

191 ------- 

192 key_size 

193 The size column key. 

194 """ 

195 if label_size == GaussianComponentConfig._size_label: 

196 return self._get_label(self.naming_scheme, self._key_sigma) 

197 elif label_size == SersicComponentConfig._size_label: 

198 return self._get_label(self.naming_scheme, self._key_reff) 

199 return label_size 

200 

201 def get_prefixed_label(self, label: str, prefix: str) -> str: 

202 """Get a prefixed label with redundant underscores removed. 

203 

204 Parameters 

205 ---------- 

206 label 

207 The label to format. 

208 prefix 

209 The prefix to prepend. 

210 

211 Returns 

212 ------- 

213 label_prefixed 

214 The prefixed label, with redundant underscores removed. 

215 """ 

216 if label.startswith("_") and ((prefix == "") or (prefix[-1] == "_")): 

217 return f"{prefix}{label[1:]}" 

218 return f"{prefix}{label}" 

219 

220 def get_suffix_dec(self) -> str: 

221 """Get the suffix for declination columns.""" 

222 return self._get_label(self.naming_scheme, self._suffix_dec) 

223 

224 def get_suffix_ra(self) -> str: 

225 """Get the suffix for right ascension columns.""" 

226 return self._get_label(self.naming_scheme, self._suffix_ra) 

227 

228 def get_suffix_x(self) -> str: 

229 """Get the suffix for x-axis columns.""" 

230 return self._get_label(self.naming_scheme, self._suffix_x) 

231 

232 def get_suffix_y(self) -> str: 

233 """Get the suffix for y-axis columns.""" 

234 return self._get_label(self.naming_scheme, self._suffix_y) 

235 

236 def make_catalog(self, n_rows: int, **kwargs): 

237 """Make a catalog with default-initialized column values. 

238 

239 Parameters 

240 ---------- 

241 n_rows 

242 The number of rows to create. 

243 **kwargs 

244 Keyword arguments to pass to self.schema. 

245 

246 Returns 

247 ------- 

248 catalog 

249 The initialized catalog. 

250 columns 

251 The columns as returned by self.schema. 

252 """ 

253 columns = self.schema(**kwargs) 

254 keys = [column.key for column in columns] 

255 prefix = self.prefix_column 

256 

257 idx_flag_first = keys.index("unknown_flag") 

258 idx_flag_last = idx_flag_first + len(self.flag_errors) 

259 dtypes = [(f'{prefix if col.key != self.column_id else ""}{col.key}', col.dtype) for col in columns] 

260 

261 results = Table(np.empty(n_rows, dtype=dtypes)) 

262 for colname in results.colnames: 

263 column = results[colname] 

264 dtype = column.dtype 

265 if ( 

266 value := ( 

267 self.float_fill_value 

268 if np.issubdtype(dtype, np.floating) 

269 else (self.integer_fill_value if np.issubdtype(dtype, np.integer) else None) 

270 ) 

271 ) is not None: 

272 column[:] = value 

273 

274 # Set nan-default flags to False instead 

275 errors = [] 

276 for flag in columns[idx_flag_first : (idx_flag_last + 1)]: 

277 column = results[f"{prefix}{flag.key}"] 

278 column[:] = False 

279 if not ((column.dtype == bool) and column.name.endswith("_flag")): 

280 errors.append(f"{column.name=} should end with _flag and {column.dtype=} must be bool") 

281 if errors: 

282 errors.append(f"These may be logic errors in {self=}") 

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

284 results.meta["config"] = self.toDict() 

285 

286 return results, columns 

287 

288 def schema( 

289 self, 

290 bands: list[str] | None = None, 

291 ) -> list[ColumnInfo]: 

292 """Return the schema as an ordered list of columns. 

293 

294 Parameters 

295 ---------- 

296 bands 

297 A list of band names to prefix band-dependent columns with. 

298 Band prefixes should not be used if None. 

299 

300 Returns 

301 ------- 

302 schema 

303 An ordered list of ColumnInfo instances. 

304 """ 

305 schema = [ 

306 ColumnInfo(key=self.column_id, dtype="i8"), 

307 ColumnInfo(key="n_iter", dtype="i4"), 

308 ColumnInfo(key="time_eval", dtype="f8", unit=u.s), 

309 ColumnInfo(key="time_fit", dtype="f8", unit=u.s), 

310 ColumnInfo(key="time_full", dtype="f8", unit=u.s), 

311 ColumnInfo(key="chisq_reduced", dtype="f8"), 

312 ColumnInfo(key="unknown_flag", dtype="bool"), 

313 ] 

314 schema.extend([ColumnInfo(key=key, dtype="bool") for key in self.flag_errors.keys()]) 

315 # Subclasses should always write out centroids even if not fitting 

316 # They are helpful for reconstructing models 

317 return schema