Coverage for estats/stats/CrossVoids.py: 93%

1# Copyright (C) 2019 ETH Zurich 

2# Institute for Particle Physics and Astrophysics 

3# Author: Dominik Zuercher 

4 

5import numpy as np 

6import healpy as hp 

7 

8 

9def context(): 

10 """ 

11 Defines the paramters used by the plugin 

12 """ 

13 stat_type = 'multi' 

14 

15 required = ['CrossVoids_steps', 

16 'void_upper_threshold', 'CrossVoids_sliced_bins', 'NSIDE', 

17 'scales', 'selected_scales', 'min_count', 'SNR_voids', 

18 'min_SNR'] 

19 defaults = [1000, -2.5, 15, 1024, 

20 [31.6, 29.0, 26.4, 23.7, 21.1, 18.5, 15.8, 13.2, 

21 10.5, 7.9, 5.3, 2.6], [31.6, 29.0, 26.4, 23.7, 21.1, 18.5, 

22 15.8, 13.2, 10.5], 30, False, -100.] 

23 types = ['int', 'float', 'int', 'int', 'list', 'list', 

24 'int', 'bool', 'float'] 

25 return required, defaults, types, stat_type 

26 

27 

28def CrossVoids(map_w, weights, ctx): 

29 """ 

30 Calculates Voids on a convergence map. 

31 :param map_w: A Healpix convergence map 

32 :param weights: A Healpix map with pixel weights 

33 :param ctx: Context instance 

34 :return: Void abundance function 

35 """ 

36 

37 # standard devaition of map 

38 sig = np.std(map_w[map_w > hp.pixelfunc.UNSEEN]) 

39 

40 ell = map_w.size 

41 

42 # Exclude last element 

43 nside = hp.get_nside(map_w) 

44 ran = np.arange(ell - 1) 

45 

46 # restrict to mask 

47 ran = ran[weights[:-1] > 0.0] 

48 

49 temp = map_w[-1] 

50 map_w[-1] = hp.pixelfunc.UNSEEN 

51 

52 voids = np.zeros(0, dtype=np.int64) 

53 

54 # calculate all the neighbours, chunked (slow but slim memory) 

55 num_chunks = 1 

56 for r in np.array_split(ran, num_chunks): 

57 neighbours = hp.get_all_neighbours(nside, r) 

58 

59 edges = np.any(map_w[neighbours] == hp.pixelfunc.UNSEEN, axis=0) 

60 

61 # Get Void positions (minima) 

62 ps = np.all(map_w[neighbours] > map_w[r], axis=0) 

63 

64 # Remove pixels which are next to an UNSEEN pixel 

65 ps = np.logical_and(ps, np.logical_not(edges)) 

66 voids = np.append(voids, r[ps]) 

67 

68 # Last Value treatment 

69 map_w[-1] = temp 

70 n0 = hp.get_all_neighbours(nside, ell - 1) 

71 if (np.all(map_w[n0] < map_w[ell - 1])): 

72 voids = np.append(voids, ell - 1) 

73 n1 = np.reshape(hp.get_all_neighbours(nside, n0).T, (-1)) 

74 voids2 = np.all(np.reshape(map_w[n1], (-1, 8)) 

75 > map_w[n0].reshape((-1, 1)), axis=1) 

76 voids2 = n0[voids2] 

77 voids = _select(n1, voids, voids2) 

78 

79 # Remove UNSEENS labeled as Voids 

80 valids = map_w[voids] > -1e+20 

81 voids = voids[valids] 

82 

83 # get values 

84 # and cutting off below threshold coordinates 

85 void_vals = map_w[voids] 

86 

87 if ctx['SNR_voids']: 

88 void_vals /= sig 

89 

90 # restrict to peaks with SNR < 4.0 

91 void_vals = void_vals[void_vals / sig >= ctx['min_SNR']] 

92 

93 # edge case (nothing found) 

94 if len(void_vals) == 0: 

95 return np.full(ctx['Voids_steps'] + 1, np.nan) 

96 

97 # Binning for values 

98 if len(void_vals) > (ctx['min_count'] * 200.): 

99 minimum = np.max( 

100 np.partition(void_vals, ctx['min_count'])[:ctx['min_count']]) 

101 maximum = np.min( 

102 np.partition(void_vals, -ctx['min_count'])[-ctx['min_count']:]) 

103 else: 

104 maximum = -1. 

105 minimum = 1. 

106 if maximum < minimum: 

107 minimum = np.min(void_vals) 

108 maximum = np.max(void_vals) 

109 void_bins = np.linspace( 

110 minimum, maximum, ctx['Voids_steps'] - 1) 

111 

112 void_vals = np.histogram(void_vals, bins=void_bins)[0] 

113 # first and last bins indicate maximum and minmum of the kappa range 

114 void_vals = np.hstack((minimum, void_vals, maximum, np.asarray([sig]))) 

115 void_vals = void_vals.reshape(1, void_vals.size) 

116 return void_vals 

117 

118 

119def _select(n, p1, p2): 

120 for i in n: 

121 if ((i in p1) and (i not in p2)): 

122 p1 = p1[p1 != i] 

123 if ((i in p2) and (i not in p1)): 

124 p1 = np.append(p1, i) 

125 return p1 

126 

127 

128def process(data, ctx, scale_to_unity=False): 

129 # backwards compatibility for data without map std 

130 if data.shape[1] > ctx['CrossVoids_steps']: 

131 data = data[:, :-1] 

132 

133 num_of_scales = len(ctx['scales']) 

134 

135 new_data = np.zeros( 

136 (int(data.shape[0] / num_of_scales), data.shape[1] 

137 * num_of_scales)) 

138 for jj in range(int(data.shape[0] / num_of_scales)): 

139 new_data[jj, :] = data[jj * num_of_scales: 

140 (jj + 1) * num_of_scales, :].ravel() 

141 return new_data 

142 

143 

144def slice(ctx): 

145 # number of datavectors for each scale 

146 mult = 1 

147 # number of scales 

148 num_of_scales = len(ctx['scales']) 

149 # either mean or sum, for how to assemble the data into the bins 

150 operation = 'sum' 

151 

152 n_bins_sliced = ctx['CrossVoids_sliced_bins'] 

153 

154 # if True assumes that first and last entries of the data vector indicate 

155 # the upper and lower boundaries and that binning scheme indicates 

156 # bin edges rather than their indices 

157 range_mode = True 

158 

159 return num_of_scales, n_bins_sliced, operation, mult, range_mode 

160 

161 

162def decide_binning_scheme(data, meta, bin, ctx): 

163 num_of_scales = len(ctx['scales']) 

164 n_bins_original = ctx['CrossVoids_steps'] 

165 n_bins_sliced = ctx['CrossVoids_sliced_bins'] 

166 

167 # get the correct tomographic bins 

168 bin_idx = np.zeros(meta.shape[0], dtype=bool) 

169 bin_idx[np.where(meta[:, 0] == bin)[0]] = True 

170 bin_idx = np.repeat(bin_idx, meta[:, 1].astype(int)) 

171 data = data[bin_idx, :] 

172 

173 # Get bins for each smooting scale 

174 bin_centers = np.zeros((num_of_scales, n_bins_sliced)) 

175 bin_edges = np.zeros((num_of_scales, n_bins_sliced + 1)) 

176 for scale in range(num_of_scales): 

177 # cut correct scale and minimum and maximum kappa values 

178 data_act = data[:, 

179 n_bins_original * scale:n_bins_original * (scale + 1)] 

180 minimum = np.max(data_act[:, 0]) 

181 maximum = np.min(data_act[:, -1]) 

182 new_kappa_bins = np.linspace(minimum, maximum, n_bins_sliced + 1) 

183 bin_edges[scale, :] = new_kappa_bins 

184 

185 bin_centers_act = new_kappa_bins[:-1] + 0.5 * \ 

186 (new_kappa_bins[1:] - new_kappa_bins[:-1]) 

187 bin_centers[scale, :] = bin_centers_act 

188 return bin_edges, bin_centers 

189 

190 

191def filter(ctx): 

192 filter = np.zeros(0) 

193 for scale in ctx['scales']: 

194 if scale in ctx['selected_scales']: 

195 f = [True] * \ 

196 ctx['CrossVoids_sliced_bins'] 

197 f = np.asarray(f) 

198 else: 

199 f = [False] * \ 

200 ctx['CrossVoids_sliced_bins'] 

201 f = np.asarray(f) 

202 filter = np.append(filter, f) 

203 return filter