Coverage for estats/stats/CrossStarletPeaksDi.py: 20%

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 estats.stats import CrossPeaks 

8 

9 

10def context(): 

11 """ 

12 Defines the paramters used by the plugin 

13 """ 

14 stat_type = 'convergence-cross' 

15 

16 required = ['Starlet_steps', 'Starlet_scalesDi', 

17 'Starlet_selected_scalesDi', 

18 'peak_lower_threshold', 'Starlet_sliced_bins', 

19 'NSIDE', 'min_count', 'SNR_peaks', 'max_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 2.5, 15, 1024, 30, False, 100.] 

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

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

25 return required, defaults, types, stat_type 

26 

27 

28def CrossStarletPeaksDi(map_w, weights, ctx): 

29 """ 

30 Performs the starlet-wavelet decomposition of map and counts the local 

31 maxima in each filter band. 

32 :param map: A Healpix convergence map 

33 :param weights: A Healpix map with pixel weights (integer >=0) 

34 :param ctx: Context instance 

35 :return: Starlet counts (num filter bands, Starlet_steps + 1) 

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 # build decomposition 

48 # (remove first map that contains remaining small scales) 

49 wavelet_counts = np.zeros((len(ctx['Starlet_scalesDi']), 

50 ctx['Starlet_steps'] + 1)) 

51 

52 # count peaks in each filter band 

53 wave_iter = esd.calc_wavelet_decomp_iter( 

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

55 counter = 0 

56 for ii, wmap in enumerate(wave_iter): 

57 if ii == 0: 

58 continue 

59 # reapply mask 

60 wmap[np.isclose(weights, 0)] = hp.UNSEEN 

61 

62 peak_vals = CrossPeaks.CrossPeaks(wmap, weights, ctx) 

63 wavelet_counts[counter] = peak_vals 

64 counter += 1 

65 

66 return wavelet_counts 

67 

68 

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

70 # backwards compatibility for data without map std 

71 if data.shape[1] > ctx['CrossPeaks_steps']: 

72 data = data[:, :-1] 

73 

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

75 

76 new_data = np.zeros( 

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

78 * num_of_scales)) 

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

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

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

82 return new_data 

83 

84 

85def slice(ctx): 

86 # number of datavectors for each scale 

87 mult = 1 

88 # number of scales 

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

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

91 operation = 'sum' 

92 

93 n_bins_sliced = ctx['Starlet_sliced_bins'] 

94 

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

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

97 # bin edges rather than their indices 

98 range_mode = True 

99 

100 return num_of_scales, n_bins_sliced, operation, mult, range_mode 

101 

102 

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

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

105 n_bins_original = ctx['Starlet_steps'] 

106 n_bins_sliced = ctx['Starlet_sliced_bins'] 

107 

108 # get the correct tomographic bins 

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

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

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

112 data = data[bin_idx, :] 

113 

114 # Get bins for each smooting scale 

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

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

117 for scale in range(num_of_scales): 

118 # cut correct scale and minimum and maximum kappa values 

119 data_act = data[:, 

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

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

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

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

124 bin_edges[scale, :] = new_kappa_bins 

125 

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

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

128 bin_centers[scale, :] = bin_centers_act 

129 return bin_edges, bin_centers 

130 

131 

132def filter(ctx): 

133 filter = np.zeros(0) 

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

135 if scale in ctx['Starlet_selected_scalesDi']: 

136 f = [True] * \ 

137 ctx['Starlet_sliced_bins'] 

138 f = np.asarray(f) 

139 else: 

140 f = [False] * \ 

141 ctx['Starlet_sliced_bins'] 

142 f = np.asarray(f) 

143 filter = np.append(filter, f) 

144 return filter