Coverage for src/galsbi/galsbi.py: 99%
139 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-03 17:18 +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.catalog_name = file_name
100 self._run(config, model_parameters[index], **kwargs)
102 __call__ = generate_catalog
104 def _run(self, config, model_parameters, **kwargs):
105 """
106 Runs the workflow with the given configuration and model parameters
108 :param config: configuration to use for generating the catalog
109 :param model_parameters: model parameters to use for generating the catalog
110 :param kwargs: additional keyword arguments to pass to the workflow (overwrites
111 the values from the model parameters and config file)
112 """
113 kargs = {}
114 for col in model_parameters.dtype.names:
115 kargs[col] = model_parameters[col]
116 if ("moffat_beta1" in model_parameters.dtype.names) and (
117 "moffat_beta2" in model_parameters.dtype.names
118 ):
119 kargs["psf_beta"] = [
120 model_parameters["moffat_beta1"][0],
121 model_parameters["moffat_beta2"][0],
122 ]
123 kargs.update(kwargs)
124 if "filters" in kargs:
125 self.filters = kargs["filters"]
126 else:
127 config_module = importlib.import_module(config)
128 self.filters = config_module.filters
130 run_util.run_ufig_from_config(config, **kargs)
132 def cite(self):
133 """
134 Prints all the papers that should be cited when using the configuration
135 specified
136 """
137 print("\033[1mPlease cite the following papers\033[0m")
138 print("=================================")
139 print("\033[1mFor using the GalSBI model:\033[0m")
140 citations.cite_galsbi_release()
141 print("\033[1mFor using the galsbi python package:\033[0m")
142 citations.cite_code_release(self.mode)
143 print("")
145 print(
146 "\033[1mFor the galaxy population model and redshift distribution:\033[0m"
147 )
148 citations.cite_abc_posterior(self.name)
149 print("")
150 print("Example:")
151 print("--------")
152 print(
153 "We use the GalSBI framework (PAPERS GalSBI release) to generate mock"
154 " galaxy catalogs. The galaxy population model corresponds to the"
155 " posterior from (PAPER model). (...) "
156 "Acknowledgements: We acknowledge the use of the following software:"
157 "(numpy), (scipy), (PAPERS code release), (...)"
158 )
160 def load_catalogs(self, output_format="rec", model_index=0, combine=False):
161 """
162 Loads the catalogs generated by the model.
164 :param output_format: format of the output, either "rec", "df" or "fits"
165 :param model_index: index of the model parameters to use for loading the
166 catalogs
167 :param combine: if True, combines the catalogs from all bands into a single
168 catalog
169 :return: catalogs in the specified format
170 """
171 if self.filters is None:
172 raise RuntimeError("please generate catalogs first")
174 if output_format == "rec":
175 convert = lambda x: x # noqa: E731
176 elif output_format == "df":
177 convert = at.rec2pd
178 elif output_format == "fits":
179 convert = Table
180 else:
181 raise ValueError(f"Unknown output format {output_format}")
183 output = {}
184 for band in self.filters:
185 catalog_name = f"{self.catalog_name}_{model_index}_{band}_ucat.gal.cat"
186 with contextlib.suppress(FileNotFoundError):
187 output[f"ucat galaxies {band}"] = at.load_hdf(catalog_name)
188 catalog_name = f"{self.catalog_name}_{model_index}_{band}_ucat.star.cat"
189 with contextlib.suppress(FileNotFoundError):
190 output[f"ucat stars {band}"] = at.load_hdf(catalog_name)
191 catalog_name = f"{self.catalog_name}_{model_index}_{band}_se.cat"
192 with contextlib.suppress(FileNotFoundError):
193 output[f"sextractor {band}"] = at.load_hdf(catalog_name)
194 if len(output) == 0:
195 LOGGER.warning(
196 "No catalogs found. Did you already generate catalogs? Does the "
197 "model_index match the one used for generating the catalogs?"
198 )
199 if not combine:
200 catalogs = {key: convert(value) for key, value in output.items()}
201 return catalogs
203 combined_catalogs = self._build_combined_catalogs(output)
204 return {key: convert(value) for key, value in combined_catalogs.items()}
206 def load_images(self, model_index=0):
207 """
208 Loads the images generated by the model. This include the actual image,
209 the segmentation map and the background map.
211 param model_index: index of the model parameters to use for loading the images
212 return: images as numpy arrays
213 """
214 output = {}
215 for band in self.filters:
216 image_name = f"{self.catalog_name}_{model_index}_{band}_image.fits"
217 try:
218 hdul = fits.open(image_name)
219 image = hdul[0].data
220 hdul.close()
221 output[f"image {band}"] = image
222 except FileNotFoundError:
223 pass
224 segmap_name = f"{self.catalog_name}_{model_index}_{band}_se_seg.h5"
225 with contextlib.suppress(FileNotFoundError):
226 output[f"segmentation {band}"] = at.load_hdf_cols(segmap_name)[
227 "SEGMENTATION"
228 ]
229 bkgmap_name = f"{self.catalog_name}_{model_index}_{band}_se_bkg.h5"
230 with contextlib.suppress(FileNotFoundError):
231 output[f"background {band}"] = at.load_hdf_cols(bkgmap_name)[
232 "BACKGROUND"
233 ]
234 return output
236 def _build_combined_catalogs(self, catalogs):
237 band_dep_params = ["int_mag", "mag", "abs_mag", "bkg_noise_amp"]
238 combined_catalogs = {}
239 filter = self.filters[0]
240 if f"ucat galaxies {filter}" in catalogs:
241 new_cat = {}
242 for f in self.filters:
243 cat = catalogs[f"ucat galaxies {f}"]
244 for par in cat.dtype.names:
245 if par not in band_dep_params:
246 new_cat[par] = cat[par]
247 else:
248 new_cat[f"{par} {f}"] = cat[par]
249 combined_catalogs["ucat galaxies"] = at.dict2rec(new_cat)
250 if f"ucat stars {filter}" in catalogs:
251 new_cat = {}
252 for f in self.filters:
253 cat = catalogs[f"ucat stars {f}"]
254 for par in cat.dtype.names:
255 if par not in band_dep_params:
256 new_cat[par] = cat[par]
257 else:
258 new_cat[f"{par} {f}"] = cat[par]
259 combined_catalogs["ucat stars"] = at.dict2rec(new_cat)
260 if f"sextractor {filter}" in catalogs:
261 band_ind_params = [
262 "dec",
263 "ra",
264 "z",
265 "e1",
266 "e2",
267 "r50",
268 "r50_arcsec",
269 "r50_phys",
270 "sersic_n",
271 "galaxy_type",
272 "id",
273 "x",
274 "y",
275 ]
276 new_cat = {}
277 for f in self.filters:
278 cat = catalogs[f"sextractor {f}"]
279 for par in cat.dtype.names:
280 if par in band_ind_params:
281 new_cat[par] = cat[par]
282 else:
283 new_cat[f"{par} {f}"] = cat[par]
284 combined_catalogs["sextractor"] = at.dict2rec(new_cat)
285 return combined_catalogs