Coverage for src/galsbi/ucat/galaxy_sampling_util.py: 98%

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

2 

3""" 

4Created on Mar 6, 2018 

5author: Joerg Herbel 

6""" 

7 

8import warnings 

9 

10import numpy as np 

11import scipy.interpolate 

12import scipy.optimize 

13import scipy.special 

14import scipy.stats 

15 

16warnings.filterwarnings("once") 

17 

18 

19class UCatNumGalError(ValueError): 

20 """ 

21 Raised when more galaxies than allowed by the input parameters are sampled 

22 """ 

23 

24 

25class UCatMZInterpError(ValueError): 

26 """ 

27 Raised when lum_fct_m_max is too low for the given redshift range 

28 """ 

29 

30 

31class Catalog: 

32 pass 

33 

34 

35def apply_pycosmo_distfun(dist_fun, z): 

36 z_unique, inv_ind = np.unique(z, return_inverse=True) 

37 dist = dist_fun(a=1 / (1 + z_unique))[inv_ind] # Unit: Mpc 

38 return dist 

39 

40 

41def intp_z_m_cut(cosmo, mag_calc, par): 

42 """ 

43 This function returns an interpolator which predicts the maximum absolute magnitude 

44 for a given redshift such that a galaxy will still fall below the limiting apparent 

45 magnitude in the reference band (par.gals_mag_max). To do this, we do a check for a 

46 number of grid points in redshift. The check consists of finding the template with 

47 the smallest ratio of flux in the ref-band at that redshift (apparent flux) and flux 

48 in the luminosity function band at redshift zero (absolute flux). We then compute 

49 the absolute magnitude a galaxy would need to have to still pass the cut assuming 

50 that its spectrum is given by only that one template which we found. This works 

51 because the mixing-in of other templates will only make the object brighter in the 

52 r-band at this redshift. See also docs/jupyter_notebooks/z_m_cut.ipynb. 

53 """ 

54 

55 def find_max_template_ind(z, n_templates): 

56 coeffs = np.eye(n_templates) 

57 

58 m_lf = mag_calc( 

59 redshifts=np.zeros(n_templates), 

60 excess_b_v=np.zeros(n_templates), 

61 coeffs=coeffs, 

62 filter_names=[par.lum_fct_filter_band], 

63 )[par.lum_fct_filter_band] 

64 

65 m_ref = mag_calc( 

66 redshifts=np.full(n_templates, z), 

67 excess_b_v=np.zeros(n_templates), 

68 coeffs=coeffs, 

69 filter_names=[par.reference_band], 

70 )[par.reference_band] 

71 

72 ind = np.argmin(m_ref - m_lf) 

73 

74 return ind 

75 

76 def app_mag_ref(z, temp_i, m_abs, n_templates): 

77 coeff = np.zeros((1, n_templates)) 

78 coeff[0, temp_i] = 1 

79 coeff[0, temp_i] = 10 ** ( 

80 0.4 

81 * ( 

82 mag_calc( 

83 redshifts=np.zeros(1), 

84 excess_b_v=np.zeros(1), 

85 coeffs=coeff, 

86 filter_names=[par.lum_fct_filter_band], 

87 )[par.lum_fct_filter_band][0] 

88 - m_abs 

89 ) 

90 ) 

91 coeff[0, temp_i] *= (1e-5 / cosmo.background.dist_lum_a(a=1 / (1 + z))) ** 2 / ( 

92 1 + z 

93 ) 

94 m = mag_calc( 

95 redshifts=np.atleast_1d(z), 

96 excess_b_v=np.zeros(1), 

97 coeffs=coeff, 

98 filter_names=[par.reference_band], 

99 )[par.reference_band] 

100 

101 return m[0] 

102 

103 n_templates = mag_calc.n_templates 

104 z_intp = [par.lum_fct_z_res] 

105 temp_ind = find_max_template_ind(z_intp[0], n_templates) 

106 abs_mag = [ 

107 scipy.optimize.brentq( 

108 f=lambda m: app_mag_ref(z_intp[0], temp_ind, m, n_templates) 

109 - par.gals_mag_max, 

110 a=-100, 

111 b=par.gals_mag_max, 

112 ) 

113 ] 

114 

115 i_loop = 0 

116 while True: 

117 i_loop += 1 

118 try: 

119 z_ = z_intp[-1] + 0.02 

120 

121 if z_ > par.lum_fct_z_max: 

122 break 

123 

124 temp_ind = find_max_template_ind(z_, n_templates) 

125 abs_mag_ = scipy.optimize.brentq( 

126 f=lambda m, 

127 temp_ind=temp_ind, 

128 z_=z_, 

129 n_templates=n_templates: app_mag_ref(z_, temp_ind, m, n_templates) 

130 - par.gals_mag_max, 

131 a=-100, 

132 b=abs_mag[-1], 

133 ) 

134 z_intp.append(z_) 

135 abs_mag.append(abs_mag_) 

136 except ValueError: 

137 break 

138 

139 intp = scipy.interpolate.interp1d( 

140 x=z_intp, 

141 y=abs_mag, 

142 kind="cubic", 

143 bounds_error=False, 

144 fill_value=(abs_mag[0], abs_mag[-1]), 

145 ) 

146 

147 if np.any(intp(intp.x) > par.lum_fct_m_max): 

148 msg = ( 

149 "par.lum_fct_m_max is too low according to z-m-interpolation," 

150 " some galaxies may be missing\n" 

151 f"gals_mag_max={par.gals_mag_max:2.2f}" 

152 f" lum_fct_z_max={par.lum_fct_z_max:2.2f}" 

153 f" m_max_sample={np.max(intp(intp.x)):2.2f}" 

154 f" lum_fct_m_max={par.lum_fct_m_max:2.2f}" 

155 ) 

156 if par.raise_z_m_interp_error: 

157 raise UCatMZInterpError(msg) 

158 else: 

159 warnings.warn(msg, stacklevel=1) 

160 

161 return intp