Coverage for src/galsbi/galsbi.py: 99%
140 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-22 16:58 +0000
1# Copyright (C) 2024 ETH Zurich
2# Institute for Particle Physics and Astrophysics
3# Author: Silvan Fischbacher
4# created: Thu Aug 08 2024
6import contextlib
7import importlib
9from astropy.io import fits
10from astropy.table import Table
11from cosmic_toolbox import arraytools as at
12from cosmic_toolbox import logger
13from ufig import run_util
15from . import citations, load
17LOGGER = logger.get_logger(__name__)
20class GalSBI:
21 """
22 This class is the main interface to the model. It provides methods to generate
23 mock galaxy catalogs and to cite the model.
24 """
26 def __init__(self, name, verbosity="info"):
27 """
28 :param name: name of the model to use
29 :param verbosity: verbosity level of the logger, either "debug", "info",
30 "warning", "error" or "critical"
31 """
32 self.name = name
33 self.mode = None
34 self.filters = None
35 self.verbosity = verbosity
37 def generate_catalog(
38 self,
39 mode="intrinsic",
40 config_file=None,
41 model_index=0,
42 file_name="GalSBI_sim",
43 verbosity=None,
44 **kwargs,
45 ):
46 """
47 Generates a mock galaxy catalog using the model and configuration specified. The
48 parameter model_index is used to select a specific set of model parameters from
49 the ABC posterior. If a list of model parameters is provided, catalogs are
50 generated for each set of parameters. The saved catalogs and images are named
51 according to the file_name and model_index.
53 Names of the files
54 ------------------
55 Intrinsic ucat galaxy catalog: f"{file_name}_{index}_{band}_ucat.gal.cat"
56 Intrinsic ucat star catalog: f"{file_name}_{index}_{band}_ucat.star.cat"
57 Output catalog: f"{file_name}_{index}_{band}_se.cat"
58 Output image: f"{file_name}_{index}_{band}_image.fits"
59 Segmentation map: f"{file_name}_{index}_{band}_se_seg.h5"
60 Background map: f"{file_name}_{index}_{band}_se_bkg.h5"
62 :param mode: mode to use for generating the catalog, either "intrinsic", "emu",
63 "image", "image+SE", "config_file"
64 :param config_file: dictionary or path to a configuration file to use for
65 generating the catalog (only used if mode="config_file")
66 :param model_index: index of the model parameters to use for generating the
67 catalog
68 :param file_name: filename of the catalog and images to generate
69 :param verbosity: verbosity level of the logger, either "debug", "info",
70 "warning", "error" or "critical"
71 :param kwargs: additional keyword arguments to pass to the workflow (overwrites
72 the values from the model parameters and config file)
73 """
74 if verbosity is not None:
75 self.verbosity = verbosity
76 logger.set_all_loggers_level(self.verbosity)
77 self.mode = mode
78 model_parameters = load.load_abc_posterior(self.name)
79 config = load.load_config(self.name, mode, config_file)
81 if isinstance(model_index, int):
82 model_index = [model_index]
83 for index in model_index:
84 LOGGER.info(
85 "Generating catalog for model"
86 f" {self.name} and mode {mode} with index {index}"
87 )
88 kwargs["galaxy_catalog_name_format"] = (
89 f"{file_name}_{index}_{ } { } _ucat.gal.cat"
90 )
91 kwargs["star_catalog_name_format"] = (
92 f"{file_name}_{index}_{ } { } _ucat.star.cat"
93 )
94 kwargs["sextractor_forced_photo_catalog_name_format"] = (
95 f"{file_name}_{index}_{ } { } _se.cat"
96 )
97 kwargs["image_name_format"] = f"{file_name}_{index}_{ } { } _image.fits"
98 kwargs["tile_name"] = ""
99 self.file_name = file_name
100 self.catalog_name = file_name # for backward compatibility
101 self._run(config, model_parameters[index], **kwargs)
103 __call__ = generate_catalog
105 def _run(self, config, model_parameters, **kwargs):
106 """
107 Runs the workflow with the given configuration and model parameters
109 :param config: configuration to use for generating the catalog
110 :param model_parameters: model parameters to use for generating the catalog
111 :param kwargs: additional keyword arguments to pass to the workflow (overwrites
112 the values from the model parameters and config file)
113 """
114 kargs = {}
115 for col in model_parameters.dtype.names:
116 kargs[col] = model_parameters[col]
117 if ("moffat_beta1" in model_parameters.dtype.names) and (
118 "moffat_beta2" in model_parameters.dtype.names
119 ):
120 kargs["psf_beta"] = [
121 model_parameters["moffat_beta1"][0],
122 model_parameters["moffat_beta2"][0],
123 ]
124 kargs.update(kwargs)
125 if "filters" in kargs:
126 self.filters = kargs["filters"]
127 else:
128 config_module = importlib.import_module(config)
129 self.filters = config_module.filters
131 run_util.run_ufig_from_config(config, **kargs)
133 def cite(self):
134 """
135 Prints all the papers that should be cited when using the configuration
136 specified
137 """
138 print("\033[1mPlease cite the following papers\033[0m")
139 print("=================================")
140 print("\033[1mFor using the GalSBI model:\033[0m")
141 citations.cite_galsbi_release()
142 print("\033[1mFor using the galsbi python package:\033[0m")
143 citations.cite_code_release(self.mode)
144 print("")
146 print(
147 "\033[1mFor the galaxy population model and redshift distribution:\033[0m"
148 )
149 citations.cite_abc_posterior(self.name)
150 print("")
151 print("Example:")
152 print("--------")
153 print(
154 "We use the GalSBI framework (PAPERS GalSBI release) to generate mock"
155 " galaxy catalogs. The galaxy population model corresponds to the"
156 " posterior from (PAPER model). (...) "
157 "Acknowledgements: We acknowledge the use of the following software:"
158 "(numpy), (scipy), (PAPERS code release), (...)"
159 )
161 def load_catalogs(self, output_format="rec", model_index=0, combine=False):
162 """
163 Loads the catalogs generated by the model.
165 :param output_format: format of the output, either "rec", "df" or "fits"
166 :param model_index: index of the model parameters to use for loading the
167 catalogs
168 :param combine: if True, combines the catalogs from all bands into a single
169 catalog
170 :return: catalogs in the specified format
171 """
172 if self.filters is None:
173 raise RuntimeError("please generate catalogs first")
175 if output_format == "rec":
176 convert = lambda x: x # noqa: E731
177 elif output_format == "df":
178 convert = at.rec2pd
179 elif output_format == "fits":
180 convert = Table
181 else:
182 raise ValueError(f"Unknown output format {output_format}")
184 output = {}
185 for band in self.filters:
186 catalog_name = f"{self.file_name}_{model_index}_{band}_ucat.gal.cat"
187 with contextlib.suppress(FileNotFoundError):
188 output[f"ucat galaxies {band}"] = at.load_hdf(catalog_name)
189 catalog_name = f"{self.file_name}_{model_index}_{band}_ucat.star.cat"
190 with contextlib.suppress(FileNotFoundError):
191 output[f"ucat stars {band}"] = at.load_hdf(catalog_name)
192 catalog_name = f"{self.file_name}_{model_index}_{band}_se.cat"
193 with contextlib.suppress(FileNotFoundError):
194 output[f"sextractor {band}"] = at.load_hdf(catalog_name)
195 if len(output) == 0:
196 LOGGER.warning(
197 "No catalogs found. Did you already generate catalogs? Does the "
198 "model_index match the one used for generating the catalogs?"
199 )
200 if not combine:
201 catalogs = {key: convert(value) for key, value in output.items()}
202 return catalogs
204 combined_catalogs = self._build_combined_catalogs(output)
205 return {key: convert(value) for key, value in combined_catalogs.items()}
207 def load_images(self, model_index=0):
208 """
209 Loads the images generated by the model. This include the actual image,
210 the segmentation map and the background map.
212 param model_index: index of the model parameters to use for loading the images
213 return: images as numpy arrays
214 """
215 output = {}
216 for band in self.filters:
217 image_name = f"{self.file_name}_{model_index}_{band}_image.fits"
218 try:
219 hdul = fits.open(image_name)
220 image = hdul[0].data
221 hdul.close()
222 output[f"image {band}"] = image
223 except FileNotFoundError:
224 pass
225 segmap_name = f"{self.file_name}_{model_index}_{band}_se_seg.h5"
226 with contextlib.suppress(FileNotFoundError):
227 output[f"segmentation {band}"] = at.load_hdf_cols(segmap_name)[
228 "SEGMENTATION"
229 ]
230 bkgmap_name = f"{self.file_name}_{model_index}_{band}_se_bkg.h5"
231 with contextlib.suppress(FileNotFoundError):
232 output[f"background {band}"] = at.load_hdf_cols(bkgmap_name)[
233 "BACKGROUND"
234 ]
235 return output
237 def _build_combined_catalogs(self, catalogs):
238 band_dep_params = ["int_mag", "mag", "abs_mag", "bkg_noise_amp"]
239 combined_catalogs = {}
240 filter = self.filters[0]
241 if f"ucat galaxies {filter}" in catalogs:
242 new_cat = {}
243 for f in self.filters:
244 cat = catalogs[f"ucat galaxies {f}"]
245 for par in cat.dtype.names:
246 if par not in band_dep_params:
247 new_cat[par] = cat[par]
248 else:
249 new_cat[f"{par} {f}"] = cat[par]
250 combined_catalogs["ucat galaxies"] = at.dict2rec(new_cat)
251 if f"ucat stars {filter}" in catalogs:
252 new_cat = {}
253 for f in self.filters:
254 cat = catalogs[f"ucat stars {f}"]
255 for par in cat.dtype.names:
256 if par not in band_dep_params:
257 new_cat[par] = cat[par]
258 else:
259 new_cat[f"{par} {f}"] = cat[par]
260 combined_catalogs["ucat stars"] = at.dict2rec(new_cat)
261 if f"sextractor {filter}" in catalogs:
262 band_ind_params = [
263 "dec",
264 "ra",
265 "z",
266 "e1",
267 "e2",
268 "r50",
269 "r50_arcsec",
270 "r50_phys",
271 "sersic_n",
272 "galaxy_type",
273 "id",
274 "x",
275 "y",
276 ]
277 new_cat = {}
278 for f in self.filters:
279 cat = catalogs[f"sextractor {f}"]
280 for par in cat.dtype.names:
281 if par in band_ind_params:
282 new_cat[par] = cat[par]
283 else:
284 new_cat[f"{par} {f}"] = cat[par]
285 combined_catalogs["sextractor"] = at.dict2rec(new_cat)
286 return combined_catalogs