Coverage for estats/stats/CrossStarletL1NormDi.py: 16%

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 

7from astropy.stats import mad_std 

8 

9 

10def context(): 

11 """ 

12 Defines the paramters used by the plugin 

13 """ 

14 stat_type = 'convergence-cross' 

15 

16 required = ['Starlet_L1_steps', 'Starlet_scalesDi', 

17 'Starlet_L1_selected_scalesDi', 

18 'Starlet_L1_sliced_bins', 'NSIDE', 

19 'min_SL1_SNR', 'max_SL1_SNR'] 

20 defaults = [1000, [8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096], 

21 [8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096], 

22 15, 1024, -4., 4.] 

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

24 'float', 'float'] 

25 return required, defaults, types, stat_type 

26 

27 

28def CrossStarletL1NormDi(map_w, weights, ctx): 

29 """ 

30 Performs Starlet decompostion of map and calculates the L1 norm of 

31 each filter band. Uses the dyadic scheme. 

32 :param map: A Healpix convergence map 

33 :param weights: A Healpix map with pixel weights 

34 :param ctx: Context instance 

35 :return: Starlet L1 norm (num filter bands * Starlet_steps) 

36 """ 

37 

38 try: 

39 from esd import esd 

40 except ImportError: 

41 raise ImportError( 

42 "Did not find esd package. " 

43 "It is required for this module to work properly. " 

44 "Download from: " 

45 "https://cosmo-gitlab.phys.ethz.ch/cosmo_public/esd") 

46 

47 l1_coll = np.zeros((len(ctx['Starlet_scalesDi']), 

48 ctx['Starlet_L1_steps'] + 1)) 

49 

50 # noise map to estimate std of wavelet coeffs 

51 std = mad_std(map_w[map_w > hp.UNSEEN]) 

52 rands = np.random.randn(np.sum(map_w > hp.UNSEEN)) 

53 noise_map = np.full(map_w.size, hp.UNSEEN) 

54 noise_map[map_w > hp.UNSEEN] = rands 

55 

56 # generators for wavelet decompositions 

57 wave_iter = esd.calc_wavelet_decomp_iter( 

58 map_w, l_bins=ctx['Starlet_scalesDi']) 

59 noise_iter = esd.calc_wavelet_decomp_iter( 

60 noise_map, l_bins=ctx['Starlet_scalesDi']) 

61 

62 counter = 0 

63 for ii, maps in enumerate(zip(wave_iter, noise_iter)): 

64 if ii == 0: 

65 continue 

66 

67 wmap = maps[0] 

68 nmap = maps[1] 

69 

70 # redo masking 

71 wmap = wmap[weights > 0.0] 

72 nmap = nmap[weights > 0.0] 

73 

74 noise_est = np.std(nmap) * std 

75 snr = wmap / noise_est 

76 minimum = np.min(snr) 

77 maximum = np.max(snr) 

78 

79 thresholds_snr = np.linspace(minimum, maximum, 

80 ctx['Starlet_L1_steps'] - 1) 

81 digitized = np.digitize(snr, thresholds_snr) 

82 snr_abs = np.abs(snr) 

83 bin_l1_norm = [np.sum(snr_abs[digitized == i]) 

84 for i in range(1, len(thresholds_snr))] 

85 

86 # append min, max and std of the map 

87 bin_l1_norm = np.hstack((np.asarray([minimum]), 

88 bin_l1_norm, np.asarray([maximum]), 

89 np.asarray([std]))) 

90 bin_l1_norm = bin_l1_norm.reshape(1, bin_l1_norm.size) 

91 l1_coll[counter] = bin_l1_norm 

92 counter += 1 

93 return l1_coll 

94 

95 

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

97 # backwards compatibility for data without map std 

98 if data.shape[1] > ctx['Starlet_L1_steps']: 

99 data = data[:, :-1] 

100 

101 num_of_scales = len(ctx['Starlet_scalesDi']) 

102 

103 new_data = np.zeros( 

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

105 * num_of_scales)) 

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

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

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

109 return new_data 

110 

111 

112def slice(ctx): 

113 # number of datavectors for each scale 

114 mult = 1 

115 # number of scales 

116 num_of_scales = len(ctx['Starlet_scalesDi']) 

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

118 operation = 'sum' 

119 

120 n_bins_sliced = ctx['Starlet_L1_sliced_bins'] 

121 

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

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

124 # bin edges rather than their indices 

125 range_mode = True 

126 

127 return num_of_scales, n_bins_sliced, operation, mult, range_mode 

128 

129 

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

131 num_of_scales = len(ctx['Starlet_scalesDi']) 

132 n_bins_original = ctx['Starlet_L1_steps'] 

133 n_bins_sliced = ctx['Starlet_L1_sliced_bins'] 

134 

135 # get the correct tomographic bins 

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

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

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

139 data = data[bin_idx, :] 

140 

141 # Get bins for each smooting scale 

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

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

144 for scale in range(num_of_scales): 

145 # cut correct scale and minimum and maximum kappa values 

146 data_act = data[:, 

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

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

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

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

151 bin_edges[scale, :] = new_kappa_bins 

152 

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

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

155 bin_centers[scale, :] = bin_centers_act 

156 return bin_edges, bin_centers 

157 

158 

159def filter(ctx): 

160 filter = np.zeros(0) 

161 for scale in reversed(ctx['Starlet_scalesDi']): 

162 if scale in ctx['Starlet_L1_selected_scalesDi']: 

163 f = [True] * \ 

164 ctx['Starlet_L1_sliced_bins'] 

165 f = np.asarray(f) 

166 else: 

167 f = [False] * \ 

168 ctx['Starlet_L1_sliced_bins'] 

169 f = np.asarray(f) 

170 filter = np.append(filter, f) 

171 return filter