Coverage for src/galsbi/ucat/magnitude_calculator.py: 95%

1# Copyright (C) 2018 ETH Zurich, Institute for Particle Physics and Astrophysics 

2 

3""" 

4Created on Sept 2021 

5author: Tomasz Kacprzak 

6Code author: Joerg Herbel 

7""" 

8 

9from collections import OrderedDict 

10 

11import numpy as np 

12from cosmic_toolbox import logger 

13 

14from galsbi.ucat import sed_templates_util, spectrum_util 

15from galsbi.ucat.galaxy_population_models.galaxy_luminosity_function import ( 

16 find_closest_ind, 

17) 

18 

19LOGGER = logger.get_logger(__file__) 

20 

21speed_of_light_c = 3 * 1e8 * 1e6 # speed of light, units: micrometre/s 

22 

23TEMPL_INTEGRALS_CACHE = OrderedDict() 

24 

25 

26def flux_to_AB_mag(flux): 

27 mag = -2.5 * np.log10(flux) - 48.6 

28 mag[np.isnan(mag)] = np.inf 

29 return mag 

30 

31 

32def AB_mag_to_flux(mag): 

33 return 10 ** ((mag + 48.6) / (-2.5)) 

34 

35 

36class MagCalculatorDirect: 

37 """ 

38 Computes galaxy magnitudes by integrating redshifted and reddened galaxy spectra 

39 which are given by a sum of five template spectra. See also 

40 docs/jupyter_notebooks/spectra_to_magnitudes.ipynb and 

41 docs/jupyter_notebooks/extinction.ipynb. 

42 """ 

43 

44 def __init__(self, filters, sed_templates): 

45 self.filters = filters 

46 self.sed_templates = sed_templates 

47 self.extinction_spline = spectrum_util.spline_ext_coeff() 

48 self.n_templates = sed_templates["amp"].shape[0] 

49 

50 def __call__( 

51 self, redshifts, excess_b_v, coeffs, filter_names, return_fluxes=False 

52 ): 

53 # check if iterable 

54 from collections.abc import Iterable 

55 

56 if not isinstance(filter_names, Iterable): 

57 filter_names = [filter_names] 

58 

59 n_obj = redshifts.size 

60 fluxes = {f: np.empty(n_obj) for f in filter_names} 

61 

62 for i in LOGGER.progressbar( 

63 range(n_obj), desc="getting magnitudes", at_level="debug" 

64 ): 

65 lam_obs = self.sed_templates["lam"] * (1 + redshifts[i]) 

66 

67 spec = spectrum_util.construct_reddened_spectrum( 

68 lam_obs=lam_obs, 

69 templates_amp=self.sed_templates["amp"], 

70 coeff=coeffs[i], 

71 excess_b_v=excess_b_v[i], 

72 extinction_spline=self.extinction_spline, 

73 ) 

74 

75 for f in filter_names: 

76 filt_amp = np.interp( 

77 lam_obs, self.filters[f]["lam"], self.filters[f]["amp"] 

78 ) 

79 fluxes[f][i] = np.trapz(lam_obs * spec * filt_amp, x=lam_obs).item() / ( 

80 speed_of_light_c * self.filters[f]["integ"] 

81 ) 

82 

83 if return_fluxes: 

84 return fluxes 

85 

86 else: 

87 mags = {} 

88 for f in filter_names: 

89 mags[f] = flux_to_AB_mag(fluxes[f]) 

90 mags[f][np.isnan(mags[f])] = np.inf 

91 del fluxes[f] 

92 

93 return mags 

94 

95 

96class MagCalculatorTable: 

97 """ 

98 Computes galaxy magnitudes by looking up pre-computed values of the integrals of our 

99 template spectra as a function of redshift and extinction E(B-V). The integrals need 

100 to be pre-computed for every filter band separately, such that for every filter 

101 band, we have five (number of template spectra) 2-dim. tables of integrals with 

102 redshift on the one and E(B-V) on the other axis. See also 

103 docs/jupyter_notebooks/tabulate_template_integrals.ipynb and 

104 docs/jupyter_notebooks/extinction.ipynb. 

105 """ 

106 

107 def __init__( 

108 self, 

109 filter_names, 

110 filepath_sed_integrals, 

111 reload_cache=False, 

112 copy_to_cwd=False, 

113 ): 

114 # TODO: check if this cache helped at all and if it did, fix the memory leak 

115 """ 

116 if reload_cache: 

117 for key in TEMPL_INTEGRALS_CACHE.keys(): 

118 del TEMPL_INTEGRALS_CACHE[key] 

119 sed_templates_util.load_sed_integrals( 

120 filepath_sed_integrals, 

121 filter_names, 

122 crop_negative=True, 

123 sed_templates=TEMPL_INTEGRALS_CACHE, 

124 ) 

125 self.templates_int_table_dict = TEMPL_INTEGRALS_CACHE 

126 self.z_grid = TEMPL_INTEGRALS_CACHE.z_grid 

127 self.excess_b_v_grid = TEMPL_INTEGRALS_CACHE.excess_b_v_grid 

128 

129 """ 

130 self.templates_int_table_dict = sed_templates_util.load_sed_integrals( 

131 filepath_sed_integrals, 

132 filter_names, 

133 crop_negative=True, 

134 sed_templates=None, 

135 copy_to_cwd=copy_to_cwd, 

136 ) 

137 self.z_grid = self.templates_int_table_dict.z_grid 

138 self.excess_b_v_grid = self.templates_int_table_dict.excess_b_v_grid 

139 

140 self.n_templates = self.templates_int_table_dict.n_templates 

141 

142 def __call__( 

143 self, redshifts, excess_b_v, coeffs, filter_names, return_fluxes=False 

144 ): 

145 # check if iterable 

146 from collections.abc import Iterable 

147 

148 if not isinstance(filter_names, Iterable): 

149 filter_names = [filter_names] 

150 

151 # find the values of redshift and extinction in the grid for all objects 

152 z_ind = find_closest_ind(self.z_grid, redshifts) 

153 excess_b_v_ind = find_closest_ind(self.excess_b_v_grid, excess_b_v) 

154 

155 fluxes = {} 

156 

157 for f in filter_names: 

158 templates_int_tables = self.templates_int_table_dict[f] 

159 fluxes[f] = np.zeros(redshifts.size, dtype=np.float64) 

160 

161 # create fluxes from templates 

162 for i in range(coeffs.shape[1]): 

163 fluxes[f] += ( 

164 coeffs[:, i] * templates_int_tables[i][z_ind, excess_b_v_ind] 

165 ) 

166 

167 if return_fluxes: 

168 return fluxes 

169 else: 

170 # convert flux to AB mag 

171 magnitudes = {} 

172 for f in filter_names: 

173 magnitudes[f] = flux_to_AB_mag(fluxes[f]) 

174 del fluxes[f] 

175 return magnitudes