Coverage for src/galsbi/ucat/galaxy_population_models/galaxy_shape.py: 87%

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

2 

3""" 

4Created 2021 

5author: Tomasz Kacprzak 

6""" 

7 

8from warnings import warn 

9 

10import numpy as np 

11from scipy.stats import beta 

12 

13LOG_0_2539 = np.log(0.2539) 

14 

15 

16def distortion_to_shear(distortion_x): 

17 q_sq = (1 - distortion_x) / (1 + distortion_x) 

18 q = np.sqrt(q_sq) 

19 shear_x = (1 - q) / (1 + q) 

20 return shear_x 

21 

22 

23########################################################## 

24# 

25# Sample galaxy ellipticities 

26# 

27########################################################## 

28 

29 

30def sample_ellipticities_gaussian(numgalaxies, e1_mean, e2_mean, e1_sigma, e2_sigma): 

31 """ 

32 Sample Gaussian distributions for the intrinsic e1 and e2 while enforcing that 

33 e1**2 + e2**2 <= 1 

34 

35 :param numgalaxies: number of galaxies, i.e. number of samples 

36 :param par: ctx.parameters; part of ctx containing parameters 

37 :return: e1 values 

38 :return: e2 values 

39 """ 

40 

41 gal_e1 = np.ones(numgalaxies) 

42 gal_e2 = np.ones(numgalaxies) 

43 

44 while np.any(gal_e1**2 + gal_e2**2 >= 1.0): 

45 index = gal_e1**2 + gal_e2**2 >= 1.0 

46 

47 if np.any(index): 

48 gal_e1[index] = np.random.normal(e1_mean, e1_sigma, size=np.sum(index)) 

49 gal_e2[index] = np.random.normal(e2_mean, e2_sigma, size=np.sum(index)) 

50 

51 return gal_e1, gal_e2 

52 

53 

54def pe_disc(ngal, log_a=LOG_0_2539, emax=0.8, emin=0.0256, pow_alpha=1): 

55 """ 

56 From miller2013. 

57 """ 

58 

59 e_distortion = np.linspace(0, 1, 100000) 

60 e_shear = distortion_to_shear(e_distortion) 

61 pe = ( 

62 (e_shear**pow_alpha) 

63 * (1 - np.exp((e_shear - emax) / np.exp(log_a))) 

64 / (1.0 + e_shear) 

65 / np.sqrt(e_shear**2 + emin**2) 

66 ) 

67 pe[e_shear > 1] = 0 

68 pe[pe < 0] = 0 

69 pe = pe / np.sum(pe) 

70 pep = np.random.choice(a=e_distortion, size=ngal, p=pe) 

71 ea = pep * np.exp(1j * np.random.uniform(0.0, 2 * np.pi, len(pep))) 

72 e1, e2 = ea.real, ea.imag 

73 return e1, e2 

74 

75 

76def pe_bulge(ngal, b=2.368, c=6.691): 

77 """ 

78 From miller2013. 

79 """ 

80 

81 e_distortion = np.linspace(0, 1, 100000) 

82 e_shear = distortion_to_shear(e_distortion) 

83 pe = e_shear * np.exp(-b * e_shear - c * e_shear**2) 

84 pe[e_shear > 1] = 0 

85 pe[pe < 0] = 0 

86 pe = pe / np.sum(pe) 

87 pep = np.random.choice(a=e_distortion, size=ngal, p=pe) 

88 ea = pep * np.exp(1j * np.random.uniform(0.0, 2 * np.pi, len(pep))) 

89 e1, e2 = ea.real, ea.imag 

90 return e1, e2 

91 

92 

93def sample_ellipticities_beta(n_gal, par): 

94 grid_e = np.linspace(0, 1, 100000) 

95 beta_a = par.ell_beta_ab_sum * (par.ell_beta_ab_ratio) 

96 beta_b = par.ell_beta_ab_sum * (1 - par.ell_beta_ab_ratio) 

97 rv = beta(beta_a, beta_b) 

98 pe = (grid_e / par.ell_beta_emax) ** 0.2 * rv.pdf(grid_e / par.ell_beta_emax) 

99 pe[~np.isfinite(pe)] = 0 

100 pe = pe / np.sum(pe) 

101 pep = np.random.choice(a=grid_e, size=n_gal, p=pe) 

102 ea = pep * np.exp(1j * np.random.uniform(0.0, 2 * np.pi, len(pep))) 

103 e1, e2 = ea.real, ea.imag 

104 return e1, e2 

105 

106 

107def sample_ellipticities_beta_mode( 

108 n_gal, ell_beta_ab_sum, ell_beta_mode, ell_beta_emax 

109): 

110 beta_a = ell_beta_ab_sum * ell_beta_mode - 2 * ell_beta_mode + 1 

111 beta_b = ell_beta_ab_sum - beta_a 

112 pep = beta.rvs(beta_a, beta_b, scale=ell_beta_emax, size=n_gal) 

113 ea = pep * np.exp(1j * np.random.uniform(0.0, 2 * np.pi, len(pep))) 

114 e1, e2 = ea.real, ea.imag 

115 return e1, e2 

116 

117 

118def sample_ellipticities_for_galaxy_type(n_gal, galaxy_type, par): 

119 backwards_compatibility(par) 

120 

121 if par.ellipticity_sampling_method == "gaussian": 

122 int_e1, int_e2 = sample_ellipticities_gaussian( 

123 numgalaxies=n_gal, 

124 e1_mean=par.e1_mean, 

125 e2_mean=par.e2_mean, 

126 e1_sigma=par.e1_sigma, 

127 e2_sigma=par.e2_sigma, 

128 ) 

129 

130 elif par.ellipticity_sampling_method == "gaussian_blue_red": 

131 int_e1, int_e2 = sample_ellipticities_gaussian( 

132 numgalaxies=n_gal, 

133 e1_mean=getattr(par, f"e1_mean_{galaxy_type}"), 

134 e2_mean=getattr(par, f"e2_mean_{galaxy_type}"), 

135 e1_sigma=getattr(par, f"ell_sigma_{galaxy_type}"), 

136 e2_sigma=getattr(par, f"ell_sigma_{galaxy_type}"), 

137 ) 

138 

139 elif par.ellipticity_sampling_method == "blue_red_miller2013": 

140 if galaxy_type == "blue": 

141 int_e1, int_e2 = pe_disc( 

142 ngal=n_gal, 

143 log_a=par.ell_disc_log_a, 

144 emin=par.ell_disc_min_e, 

145 pow_alpha=par.ell_disc_pow_alpha, 

146 ) 

147 

148 elif galaxy_type == "red": 

149 int_e1, int_e2 = pe_bulge(ngal=n_gal, b=par.ell_bulge_b) 

150 

151 else: 

152 raise ValueError( 

153 f"galaxy type {galaxy_type} not supported for" 

154 " par.ellipticity_sampling_method=blue_red_miller2013" 

155 ) 

156 

157 elif par.ellipticity_sampling_method == "beta_ratio": 

158 int_e1, int_e2 = sample_ellipticities_beta(n_gal=n_gal, par=par) 

159 

160 elif par.ellipticity_sampling_method == "beta_mode": 

161 int_e1, int_e2 = sample_ellipticities_beta_mode( 

162 n_gal=n_gal, 

163 ell_beta_ab_sum=par.ell_beta_ab_sum, 

164 ell_beta_mode=par.ell_beta_mode, 

165 ell_beta_emax=par.ell_beta_emax, 

166 ) 

167 

168 elif par.ellipticity_sampling_method == "beta_mode_red_blue": 

169 int_e1, int_e2 = sample_ellipticities_beta_mode( 

170 n_gal=n_gal, 

171 ell_beta_ab_sum=getattr(par, f"ell_beta_ab_sum_{galaxy_type}"), 

172 ell_beta_mode=getattr(par, f"ell_beta_mode_{galaxy_type}"), 

173 ell_beta_emax=par.ell_beta_emax, 

174 ) 

175 

176 else: 

177 raise ValueError( 

178 f"unknown ellipticity_sampling_method {par.ellipticity_sampling_method}" 

179 ) 

180 

181 return int_e1, int_e2 

182 

183 

184def backwards_compatibility(par): 

185 if par.ellipticity_sampling_method == "default": 

186 warn( 

187 DeprecationWarning( 

188 "ellipticity_sampling_method=default is deprecated, " 

189 "using gaussian instead" 

190 ), 

191 stacklevel=1, 

192 ) 

193 par.ellipticity_sampling_method = "gaussian" 

194 

195 if par.ellipticity_sampling_method == "blue_red": 

196 warn( 

197 DeprecationWarning( 

198 "ellipticity_sampling_method=blue_red is deprecated, " 

199 "using gaussian_blue_red instead" 

200 ), 

201 stacklevel=1, 

202 ) 

203 par.ellipticity_sampling_method = "gaussian_blue_red" 

204 

205 if par.ellipticity_sampling_method == "beta_function": 

206 warn( 

207 DeprecationWarning( 

208 "ellipticity_sampling_method=beta_function is deprecated, " 

209 "using beta_ratio instead" 

210 ), 

211 stacklevel=1, 

212 ) 

213 par.ellipticity_sampling_method = "beta_ratio" 

214 

215 if par.ellipticity_sampling_method == "beta_function_mode": 

216 warn( 

217 DeprecationWarning( 

218 "ellipticity_sampling_method=beta_function_mode is deprecated, " 

219 "using beta_mode instead" 

220 ), 

221 stacklevel=1, 

222 ) 

223 par.ellipticity_sampling_method = "beta_mode" 

224 

225 if par.ellipticity_sampling_method == "beta_function_mode_red_blue": 

226 warn( 

227 DeprecationWarning( 

228 "ellipticity_sampling_method=beta_function_mode_red_blue is " 

229 "deprecated, using beta_mode_red_blue instead" 

230 ), 

231 stacklevel=1, 

232 ) 

233 par.ellipticity_sampling_method = "beta_mode_red_blue"