Coverage for estats/stats/ShearCLs.py: 97%
Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# Copyright (C) 2019 ETH Zurich
2# Institute for Particle Physics and Astrophysics
3# Author: Dominik Zuercher
5import numpy as np
6import healpy as hp
9def context():
10 """
11 Defines the paramters used by the plugin
12 """
13 stat_type = 'shear'
15 required = ['lmax', 'ShearCLs_sliced_bins', 'ShearCLs_min',
16 'ShearCLs_max', 'cl_engine']
17 defaults = [1000, 20, 100, 1000, 'anafast']
18 types = ['int', 'int', 'int', 'int', 'str']
19 return required, defaults, types, stat_type
22def ShearCLs(map_1, map_2, weight_map, ctx):
23 """
24 Calculates cross angular power spectrum of two sets of shear maps.
25 :param map_1: A Healpix shear map, first shear component.
26 :param map_2: A Healpix shear map, second shear component.
27 :param weight_map: A Healpix map with pixel weights
28 for the first shear maps
29 """
30 # anafast uses HealPix polarisation -> negate map_1!
31 # not sure if this is also true for polspice
32 shear_map_1 = [-map_1, map_2]
34 lmax = ctx['lmax'] - 1
36 # check if weights are set
37 if len(weight_map) == 0:
38 weight_map = np.ones_like(map_1, dtype=ctx['prec'])
40 if ctx['cl_engine'] == 'anafast':
41 # proper weight masking
42 select_seen = weight_map > 0
43 select_unseen = ~select_seen
44 select_unseen |= shear_map_1[0] == hp.UNSEEN
45 select_unseen |= shear_map_1[1] == hp.UNSEEN
46 maps_0_w = (shear_map_1[0] * weight_map).astype(np.float64)
47 maps_1_w = (shear_map_1[1] * weight_map).astype(np.float64)
48 maps_0_w[select_unseen] = hp.UNSEEN
49 maps_1_w[select_unseen] = hp.UNSEEN
50 shear_map_1_w = [maps_0_w, maps_1_w]
52 dummie_map = np.zeros_like(shear_map_1[0])
54 _, cl_e1e2, cl_b1b2, _, _, _ = np.array(
55 hp.sphtfunc.anafast(
56 (dummie_map, shear_map_1_w[0], shear_map_1_w[1]),
57 (dummie_map, shear_map_1_w[0], shear_map_1_w[1]),
58 lmax=lmax))
59 Cl = {'cl_EE': cl_e1e2, 'cl_BB': cl_b1b2}
61 else:
62 raise Exception("Unknown cl_engine {}".format(ctx['cl_engine']))
64 Cl_EE = Cl['cl_EE']
65 Cl_EE = Cl_EE.reshape(1, Cl_EE.size)
66 Cl_BB = Cl['cl_BB']
67 Cl_BB = Cl_BB.reshape(1, Cl_BB.size)
68 return (Cl_EE, Cl_BB)
71def process(data, ctx, scale_to_unity=False):
72 if scale_to_unity:
73 data *= 1e4
74 return data
77def slice(ctx):
78 # number of datavectors for each scale
79 mult = 1
80 # number of scales
81 num_of_scales = 1
82 # either mean or sum, for how to assemble the data into the bins
83 operation = 'mean'
85 n_bins_sliced = ctx['ShearCLs_sliced_bins']
87 return num_of_scales, n_bins_sliced, operation, mult
90def decide_binning_scheme(data, meta, bin, ctx):
91 # Perform simple equal width splitting for Cls.
93 range_edges = [ctx['ShearCLs_min'], ctx['ShearCLs_max']]
94 num_of_scales = 1
95 n_bins_sliced = ctx['ShearCLs_sliced_bins']
96 bin_centers = np.zeros((num_of_scales, n_bins_sliced))
97 bin_edge_indices = np.zeros((num_of_scales, n_bins_sliced + 1))
99 # Cut out desired l range
100 minimum = range_edges[0]
101 maximum = range_edges[1]
102 diff = maximum - minimum
104 per_bin = diff // n_bins_sliced
105 remain = diff % n_bins_sliced
106 remain_front = remain // 2
107 remain_back = remain_front + remain % 2
109 # Decide on edge indices
110 bin_edge_indices_temp = np.arange(
111 remain_front + minimum, maximum - remain_back, per_bin)
112 bin_edge_indices_temp[0] -= remain_front
113 bin_edge_indices_temp = np.append(bin_edge_indices_temp, maximum)
115 # Decide on central bin values
116 bin_centers_temp = np.zeros(0)
117 for jj in range(len(bin_edge_indices_temp) - 1):
118 bin_centers_temp = np.append(
119 bin_centers_temp,
120 np.nanmean(bin_edge_indices_temp[jj:jj + 2]))
122 # For consistency with other binning scheme
123 # assign same binning to all scales
124 for scale in range(num_of_scales):
125 bin_centers[scale, :] = bin_centers_temp
126 bin_edge_indices[scale, :] = bin_edge_indices_temp
128 return bin_edge_indices, bin_centers
131def filter(ctx):
132 return [True] * ctx['ShearCLs_sliced_bins']