Coverage for src/galsbi/ucat/galaxy_population_models/galaxy_light_profile.py: 96%

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 betaprime 

12 

13########################################################## 

14# 

15# Sampling sersic 

16# 

17########################################################## 

18 

19 

20def sample_sersic_berge(numgalaxies, int_mag, par): 

21 """ 

22 Sample the Sersic index-distribution parametrized as described in (Berge et al. 

23 2013) 

24 

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

26 :param int_mag: intrinsic, unmagnified magnitude of all galaxies 

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

28 :return: Sersic index 

29 """ 

30 

31 sersic_n = np.zeros(numgalaxies, dtype=par.catalog_precision) 

32 

33 while np.any(sersic_n <= 0) or np.any(sersic_n > 10): 

34 index = np.where((int_mag < 20) & ((sersic_n <= 0) | (sersic_n > 10)))[0] 

35 

36 if index.size > 0: 

37 sersic_n_temp = np.zeros_like(sersic_n[index]) 

38 midpoint = index.size // 2 

39 rest = index.size % 2 

40 sersic_n_temp[:midpoint] = ( 

41 np.exp( 

42 np.random.normal( 

43 par.sersic_n_mean_1_hi, par.sersic_n_sigma_1_hi, size=midpoint 

44 ) 

45 ) 

46 + np.ones_like(midpoint) * par.sersic_n_offset 

47 ) 

48 

49 sersic_n_temp[midpoint:] = ( 

50 np.exp( 

51 np.random.normal( 

52 par.sersic_n_mean_2_hi, 

53 par.sersic_n_sigma_2_hi, 

54 size=midpoint + rest, 

55 ) 

56 ) 

57 + np.ones_like(midpoint + rest) * par.sersic_n_offset 

58 ) 

59 

60 sersic_n[index] = np.random.permutation(sersic_n_temp).copy() 

61 

62 index = (int_mag >= 20) & ((sersic_n <= 0) | (sersic_n > 10)) 

63 

64 if np.any(index): 

65 sersic_n[index] = ( 

66 np.exp( 

67 np.random.normal( 

68 par.sersic_n_mean_low, 

69 par.sersic_n_sigma_low, 

70 size=np.sum(index), 

71 ) 

72 ) 

73 + np.ones_like(np.sum(index)) * par.sersic_n_offset 

74 ) 

75 

76 return sersic_n 

77 

78 

79def sample_sersic_betaprime(n_gal, mode, size, alpha=0, z=0.0, min_n=0.2, max_n=5.0): 

80 """ 

81 Sample the Sersic index-distribution parametrized as described in (Moser et al. 

82 2024). The parameter mode corresponds to the mode of the distribution, the size 

83 parameter is responsible for the scatter (with larger size the distribution becomes 

84 tighter). The alpha parameter is responsible for the redshift dependence of the 

85 mode. This was first introduced in Fischbacher et al. 2024. 

86 

87 :param n_gal: number of galaxies, i.e. number of samples 

88 :param mode: mode of the distribution 

89 :param size: size of the distribution (with larger size the distribution becomes 

90 tighter) 

91 :param alpha: redshift dependence of the mode 

92 :param z: redshift 

93 :param min_n: minimum Sersic index, raise exception if the sampled value is below 

94 this value 

95 :param max_n: maximum Sersic index, raise exception if the sampled value is above 

96 this value 

97 :return: Sersic index 

98 """ 

99 sersic_n = np.zeros(n_gal) 

100 mode_z = mode * (1 + z) ** alpha 

101 betaprime_a = mode_z * (size + 1) + 1 

102 betaprime_b = size 

103 sersic_n[:] = betaprime.rvs(a=betaprime_a, b=betaprime_b, size=n_gal) 

104 n_iter_max = 10 * n_gal 

105 for _ in range(n_iter_max): 

106 mask = (sersic_n > max_n) | (sersic_n < min_n) 

107 if np.any(mask): 

108 if isinstance(betaprime_a, np.ndarray): 

109 a = betaprime_a[mask] 

110 else: 

111 a = betaprime_a 

112 sersic_n[mask] = betaprime.rvs( 

113 a=a, b=betaprime_b, size=np.count_nonzero(mask) 

114 ) 

115 else: 

116 return sersic_n 

117 

118 raise RuntimeError( 

119 "Failed to get valid sersic indices after" 

120 f" {n_iter_max} iterations for mode {mode} and alpha {alpha}" 

121 ) 

122 

123 

124def sample_sersic_for_galaxy_type(n_gal, galaxy_type, app_mag, par, z=0.0): 

125 # Backward compability 

126 if par.sersic_sampling_method == "default": 

127 warn( 

128 DeprecationWarning( 

129 "sersic_sampling_method=default is deprecated, using berge instead" 

130 ), 

131 stacklevel=1, 

132 ) 

133 par.sersic_sampling_method = "berge" 

134 

135 if par.sersic_sampling_method == "berge": 

136 sersic_n = sample_sersic_berge(n_gal, app_mag, par).astype( 

137 par.catalog_precision 

138 ) 

139 

140 elif par.sersic_sampling_method == "blue_red_fixed": 

141 assert galaxy_type in [ 

142 "blue", 

143 "red", 

144 ], f"sersic index for galaxy_type {galaxy_type} not implemented" 

145 

146 sersic_n = np.full( 

147 n_gal, 

148 fill_value=getattr(par, f"sersic_index_{galaxy_type}"), 

149 dtype=par.catalog_precision, 

150 ) 

151 

152 elif par.sersic_sampling_method == "single": 

153 sersic_n = np.full( 

154 n_gal, fill_value=par.sersic_single_value, dtype=par.catalog_precision 

155 ) 

156 

157 elif par.sersic_sampling_method == "blue_red_betaprime": 

158 assert galaxy_type in [ 

159 "blue", 

160 "red", 

161 ], f"sersic index for galaxy_type {galaxy_type} not implemented" 

162 

163 sersic_n = sample_sersic_betaprime( 

164 n_gal=n_gal, 

165 mode=getattr(par, f"sersic_betaprime_{galaxy_type}_mode"), 

166 size=getattr(par, f"sersic_betaprime_{galaxy_type}_size"), 

167 alpha=getattr(par, f"sersic_betaprime_{galaxy_type}_mode_alpha"), 

168 z=z, 

169 min_n=par.sersic_n_min, 

170 max_n=par.sersic_n_max, 

171 ) 

172 

173 else: 

174 raise ValueError("unknown sersic_sampling_method {par.sersic_sampling_method}") 

175 

176 return sersic_n.astype(par.catalog_precision)