Coverage for metadetect / lsst / tests / test_lsst_metacal_exposures.py: 12%

1import numpy as np 

2from metadetect.lsst.util import exp2obs 

3from metadetect.lsst.util import get_integer_center, get_jacobian, get_mbexp 

4from metadetect.lsst.metacal_exposures import ( 

5 get_metacal_exps_fixnoise, get_metacal_mbexps_fixnoise, 

6) 

7from copy import deepcopy 

8import ngmix 

9from ngmix.metacal import get_all_metacal 

10import pytest 

11from lsst.utils import getPackageDir 

12 

13try: 

14 getPackageDir('descwl_shear_sims') 

15 skip_tests_on_simulations = False 

16except LookupError: 

17 skip_tests_on_simulations = True 

18 

19 

20@pytest.mark.skipif( 

21 skip_tests_on_simulations, 

22 reason='descwl_shear_sims not available' 

23) 

24def test_metacal_exps(ntrial=10, show=False): 

25 from descwl_shear_sims.galaxies import make_galaxy_catalog 

26 from descwl_shear_sims.sim import make_sim 

27 from descwl_shear_sims.psfs import make_fixed_psf 

28 

29 seed = 20251117 

30 dim = 250 

31 buff = 50 

32 

33 rng = np.random.RandomState(seed) 

34 

35 for itrial in range(ntrial): 

36 galaxy_catalog = make_galaxy_catalog( 

37 rng=rng, 

38 gal_type='fixed', 

39 layout='grid', 

40 # gal_config={'mag': 22}, 

41 coadd_dim=dim, 

42 buff=buff, 

43 ) 

44 psf = make_fixed_psf(psf_type='gauss') 

45 

46 sim_data = make_sim( 

47 rng=rng, 

48 galaxy_catalog=galaxy_catalog, 

49 coadd_dim=dim, 

50 se_dim=dim, 

51 g1=0.02, 

52 g2=0.00, 

53 psf=psf, 

54 ) 

55 exp = sim_data['band_data']['i'][0] 

56 nexp = deepcopy(exp) 

57 # copy doesn't copy the psf 

58 nexp.setPsf(exp.getPsf()) 

59 assert nexp.getWcs() == exp.getWcs() 

60 

61 obs = exp2obs(exp) 

62 assert np.all(obs.image == exp.image.array) 

63 assert np.all(obs.weight == 1/exp.variance.array) 

64 cen, _ = get_integer_center(exp.getWcs(), exp.getBBox(), as_double=True) 

65 

66 psf_obj = exp.getPsf() 

67 psf_image_array = psf_obj.computeKernelImage(cen).array 

68 

69 assert np.all(obs.psf.image == psf_image_array) 

70 

71 # gwcs = get_galsim_jacobian_wcs(exp, cen=cen) 

72 jdata = obs.jacobian._data[0] 

73 ejac = get_jacobian(exp=exp, cen=cen) 

74 edata = ejac._data[0] 

75 for n in jdata.dtype.names: 

76 assert jdata[n] == edata[n] 

77 

78 if show: 

79 compare_images( 

80 obs.image, exp.image.array, label1='obs', label2='exp', 

81 ) 

82 compare_images( 

83 obs.psf.image, psf_image_array, 

84 label1='obs psf', label2='exp psf', 

85 ) 

86 

87 noise = np.sqrt(exp.variance.array) 

88 nexp.image.array[:, :] = noise * rng.normal(size=exp.image.array.shape) 

89 assert np.all(exp.image.array != nexp.image.array) 

90 

91 obs.noise = nexp.image.array 

92 

93 types = ('noshear', '1p', '1m', '2p', '2m') 

94 mdict_obs = get_all_metacal( 

95 obs, use_noise_image=True, types=types, 

96 ) 

97 mdict_exp, noise_mdict = get_metacal_exps_fixnoise( 

98 exp, nexp, types=types, 

99 ) 

100 

101 for key in mdict_obs: 

102 print('checking:', key) 

103 assert key in mdict_exp, f'checking for {key}' 

104 

105 timage = mdict_obs[key].image 

106 eimage = mdict_exp[key].image.array 

107 

108 tweight = mdict_obs[key].weight 

109 eweight = 1/mdict_exp[key].variance.array 

110 if show: 

111 compare_images(timage, eimage, label1='obs', label2='exp') 

112 assert np.all(timage == eimage) 

113 assert np.all(tweight == eweight) 

114 

115 

116@pytest.mark.skipif( 

117 skip_tests_on_simulations, 

118 reason='descwl_shear_sims not available' 

119) 

120def test_metacal_mbexp(ntrial=10, show=False): 

121 from descwl_shear_sims.galaxies import make_galaxy_catalog 

122 from descwl_shear_sims.sim import make_sim 

123 from descwl_shear_sims.psfs import make_fixed_psf 

124 

125 seed = 20251117 

126 dim = 250 

127 buff = 50 

128 bands = ['r', 'i'] 

129 

130 rng = np.random.RandomState(seed) 

131 

132 for itrial in range(ntrial): 

133 galaxy_catalog = make_galaxy_catalog( 

134 rng=rng, 

135 gal_type='fixed', 

136 layout='grid', 

137 # gal_config={'mag': 22}, 

138 coadd_dim=dim, 

139 buff=buff, 

140 ) 

141 psf = make_fixed_psf(psf_type='gauss') 

142 

143 sim_data = make_sim( 

144 bands=bands, 

145 rng=rng, 

146 galaxy_catalog=galaxy_catalog, 

147 coadd_dim=dim, 

148 se_dim=dim, 

149 g1=0.02, 

150 g2=0.00, 

151 psf=psf, 

152 ) 

153 

154 exps = [sim_data['band_data'][band][0] for band in bands] 

155 noise_exps = [] 

156 mbobs = ngmix.MultiBandObsList() 

157 for exp in exps: 

158 nexp = deepcopy(exp) 

159 nexp.setPsf(exp.getPsf()) 

160 

161 noise = np.sqrt(exp.variance.array) 

162 nexp.image.array[:, :] = noise * rng.normal(size=exp.image.array.shape) 

163 

164 noise_exps.append(nexp) 

165 

166 obs = exp2obs(exp) 

167 obs.noise = nexp.image.array 

168 obslist = ngmix.ObsList() 

169 obslist.append(obs) 

170 mbobs.append(obslist) 

171 

172 mbexp = get_mbexp(exps) 

173 noise_mbexp = get_mbexp(noise_exps) 

174 

175 types = ('noshear', '1p', '1m', '2p', '2m') 

176 mdict_obs = get_all_metacal( 

177 mbobs, use_noise_image=True, types=types, 

178 ) 

179 mdict_mbexp, noise_mdict = get_metacal_mbexps_fixnoise( 

180 mbexp, noise_mbexp, types=types, 

181 ) 

182 

183 for key in mdict_obs: 

184 print('checking:', key) 

185 assert key in mdict_mbexp, f'checking for {key}' 

186 

187 mbobs = mdict_obs[key] 

188 mbexp = mdict_mbexp[key] 

189 for iband, band in enumerate(bands): 

190 print(f'checking band {band}') 

191 obs = mbobs[iband][0] 

192 exp = mbexp[band] 

193 

194 timage = obs.image 

195 eimage = exp.image.array 

196 

197 tweight = obs.weight 

198 eweight = 1/exp.variance.array 

199 if show: 

200 compare_images(timage, eimage, label1='obs', label2='exp') 

201 assert np.all(timage == eimage) 

202 assert np.all(tweight == eweight) 

203 

204 

205def compare_images(im1, im2, label1='im1', label2='im2'): 

206 import matplotlib.pyplot as mplt 

207 fig, axs = mplt.subplots(nrows=2, ncols=2) 

208 axs[1, 1].axis('off') 

209 

210 axs[0, 0].imshow(im1) 

211 axs[0, 0].set_title(label1) 

212 

213 axs[0, 1].imshow(im1) 

214 axs[0, 1].set_title(label2) 

215 

216 axs[1, 0].imshow(im1 - im2) 

217 axs[1, 0].set_title(f'{label1} - {label2}') 

218 mplt.show() 

219 

220 

221if __name__ == '__main__': 

222 # test_metacal_exps(show=False) 

223 test_metacal_mbexp(show=True)