Coverage for estats/catalog.py: 83%

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 

7import scipy 

8from scipy import stats 

9 

10from ekit import context 

11from estats import utils 

12from ekit import logger 

13 

14LOGGER = logger.init_logger(__name__) 

15 

16 

17class catalog: 

18 """ 

19 The catalog module handles a galaxy catalog consisting out of 

20 coordinates, ellipticities, weights and tomographic bins for each galaxy. 

21 Main functionalities: Can rotate the catalog on the sky, output 

22 survey mask, weight maps, shear and convergence maps. 

23 The convergence maps are calculated from the 

24 shear maps using spherical Kaiser-Squires. 

25 Also allows generation of shape noise catalogs 

26 or shape noise maps using random 

27 rotation of the galaxies in place. 

28 Can also pass additional scalar and obtain maps for these. 

29 

30 The most important functionalities are: 

31 

32 - rotate_catalog: 

33 

34 Allows to rotate the survey on the sky by given angles. 

35 The spin-2 ellipticity field can also be 

36 rotated by the appropriate angle. 

37 

38 - get_mask: 

39 

40 Returns the survey mask as a binary Healpix 

41 (`Gorski et al. 2005 

42 <https://iopscience.iop.org/article/10.1086/427976>`_) 

43 map. 

44 

45 - get_map: 

46 

47 Can return Healpix 

48 shear maps, convergence maps or a weight map 

49 (number of galaxies per pixel). The convergence map is calculated from 

50 the ellipticities of the galaxies using the spherical Kaiser-Squires 

51 routine (`Wallis et al. 2017 <https://arxiv.org/pdf/1703.09233.pdf>`_). 

52 The shear and convergence maps are weighted by the galaxy weights. 

53 

54 - generate_shape_noise_map: 

55 

56 Generates a shape noise Healpix shear map by random rotation of the 

57 galaxies ellipticities in place. The weights are considered in the 

58 generation of the noise. 

59 

60 The accepted keywords are: 

61 

62 - NSIDE: 

63 

64 default: 1024 

65 

66 choices: an integer being a power of 2 

67 

68 The Healpix resolution that is used to produce the map products. 

69 

70 - degree: 

71 

72 default: True 

73 

74 choices: True, False 

75 

76 If True the coordinates are assumed to be given in degrees otherwise 

77 in radians. 

78 

79 - colat: 

80 

81 default: False 

82 

83 choices: True, False 

84 

85 If True the second coordinate is assumed to be a co-latitude otherwise 

86 a normal latitude is assumed. 

87 

88 - prec: 

89 

90 default: 32 

91 

92 choices: 64, 32, 16 

93 

94 Size of the float values in the Healpix maps in bits. For less then 32 

95 hp.UNSEEN is mapped to -inf -> No RAM optimization anymore 

96 """ 

97 

98 def __init__(self, alphas=[], deltas=[], e1s=[], e2s=[], 

99 redshift_bins=[], weights=[], pixels=[], context_in={}, 

100 store_pix=False, 

101 verbosity=3, **kwargs): 

102 """ 

103 Initialization function for catalog class. 

104 At least the position of the objects on the 

105 sky is required (alphas, deltas). 

106 Also accepts galaxy ellipticity components (e1s, e2s). 

107 The ellipticities can be converted to a convergence signal. 

108 Further accepts galaxy weights and redshift bins. 

109 Instead of alphas, deltas can also pass HEALPIX pixel locations. 

110 Can pass additional fields by name. 

111 They will be treated as scalar galaxy properties. 

112 The fields are identified if the string 

113 'field' is contained in the name. 

114 

115 :param alphas: List of first coordinates, either in radians or degrees 

116 :param deltas: List of second coordinates, either in radians or degrees 

117 :param e1s: List of first ellipticity components 

118 :param e2s: List of second ellipticity components 

119 :param redshift_bins: List of integers indicating the redshift bin the 

120 objects belong to (should start from 1) 

121 :param weights: List of object weights. 

122 If not passed weights are set to 1 for each object. 

123 :param context_in: A dictionary containing parameters. 

124 :param store_pix: If true only pixel location of objects is stored 

125 :param verbosity: Verbosity level (0-4) 

126 """ 

127 

128 logger.set_logger_level(LOGGER, verbosity) 

129 

130 LOGGER.debug("Initializing catalog object") 

131 

132 # setup context 

133 allowed = ['NSIDE', 'degree', 'colat', 'prec'] 

134 types = ['int', 'bool', 'bool', 'int'] 

135 defaults = [1024, True, False, 32] 

136 

137 # separate kwargs into additional fields and context variables 

138 extra_variables = {} 

139 self.extra_fields = [] 

140 for item in kwargs: 

141 if 'field' in item: 

142 LOGGER.debug(f"Found additional field {item}") 

143 setattr(self, item, kwargs[item][:]) 

144 self.extra_fields.append(item) 

145 else: 

146 LOGGER.debug(f"Found additional context variable {item}") 

147 extra_variables[item] = kwargs[item] 

148 

149 self.ctx = context.setup_context( 

150 {**context_in, **extra_variables}, allowed, types, defaults, 

151 verbosity=verbosity) 

152 

153 self.ctx['prec'] = 'float{}'.format(self.ctx['prec']) 

154 

155 # assign objects 

156 self.pixels = {} 

157 self.set_coordinates(alphas, deltas, pixels) 

158 self.set_ellipticities(e1s, e2s) 

159 self.set_weights(weights) 

160 self.set_redhift_bins(redshift_bins) 

161 if store_pix: 

162 LOGGER.debug("Only storing pixels instead of ra,dec") 

163 self.pixels[0] = self._pixelize(alphas, deltas) 

164 self.alpha = [] 

165 self.delta = [] 

166 self.store_pix = True 

167 else: 

168 self.store_pix = False 

169 

170 self._assert_lengths() 

171 

172 ################## 

173 # ACCESING OBJECTS 

174 ################## 

175 

176 def get_coordinates(self, bin=0): 

177 """ 

178 Returns the coordinates of the catalog objects 

179 

180 :param bin: Indicate the redshift bin of the objects which should be 

181 returned. If 0 all bins are used. 

182 """ 

183 idx = self._get_redshift_bool(bin) 

184 if len(self.alpha) > 0: 

185 return self.alpha[idx], self.delta[idx] 

186 else: 

187 LOGGER.info( 

188 "Found only pixel coordinates. " 

189 "Converting to angular coordinates.") 

190 delta, alpha = hp.pixelfunc.pix2ang(self.ctx['NSIDE'], 

191 self.pixels[0][idx]) 

192 if self.ctx['degree']: 

193 LOGGER.debug("Returning coordinates in degrees") 

194 alpha = np.degrees(alpha) 

195 delta = np.degrees(delta) 

196 if not self.ctx['colat']: 

197 delta = 90. - delta 

198 else: 

199 LOGGER.debug("Returning coordinates in co-latitude") 

200 else: 

201 if not self.ctx['colat']: 

202 delta = np.pi / 2. - delta 

203 else: 

204 LOGGER.debug("Returning coordinates in co-latitude") 

205 return alpha, delta 

206 

207 def get_ellipticities(self, bin=0): 

208 """ 

209 Returns the ellipticities of the catalog objects 

210 

211 :param bin: Indicate the redshift bin of the objects which should be 

212 returned. If 0 all bins are used. 

213 """ 

214 idx = self._get_redshift_bool(bin) 

215 return self.e1s[idx], self.e2s[idx] 

216 

217 def get_redshift_bins(self, bin=0): 

218 """ 

219 Returns the redshift bins of the catalog objects 

220 

221 :param bin: Indicate the redshift bin of the objects which should be 

222 returned. If 0 all bins are used. 

223 """ 

224 bin = int(bin) 

225 idx = self._get_redshift_bool(bin) 

226 return self.redshift_bin[idx] 

227 

228 def get_weights(self, bin=0): 

229 """ 

230 Returns the weights of the catalog objects 

231 

232 :param bin: Indicate the redshift bin of the objects which should be 

233 returned. If 0 all bins are used. 

234 """ 

235 bin = int(bin) 

236 idx = self._get_redshift_bool(bin) 

237 return self.weights[idx] 

238 

239 def get_pixels(self, pix_key=0, bin=0): 

240 """ 

241 Returns the shears of the catalog objects 

242 

243 :param pix_key: If multiple rotated versions of the original catalog 

244 are stored can access them with the pix_key index. 

245 Only works if store_pix was set to True. 

246 :param bin: Indicate the redshift bin of the objects which should be 

247 returned. If 0 all bins are used. 

248 """ 

249 idx = self._get_redshift_bool(bin) 

250 # if pixels not set 

251 if (len(self.pixels.keys()) == 0) & (pix_key == 0): 

252 LOGGER.debug("Converting coordinates to pixels") 

253 pixs = self._pixelize(self.alpha, self.delta, bin=bin) 

254 return pixs 

255 else: 

256 if pix_key in self.pixels.keys(): 

257 return self.pixels[pix_key][idx] 

258 else: 

259 raise KeyError( 

260 f"No pixels with key number {pix_key} found. " 

261 "Rotate catalog first") 

262 

263 # SETTING OBJECTS 

264 def set_coordinates(self, alphas, deltas, pixels=[]): 

265 """ 

266 Sets coordinate arrays 

267 :param alphas: First coordinates 

268 :param deltas: Second coordinates 

269 :param pixels: Alternatively can pass list of 

270 pixels instead of angular coordinates 

271 """ 

272 if len(pixels) > 0: 

273 LOGGER.debug("Go list of pixels. Ignoring passed ra,decs.") 

274 delta, alpha = hp.pixelfunc.pix2ang(self.ctx['NSIDE'], pixels) 

275 

276 # converting coordinates 

277 if self.ctx['degree']: 

278 LOGGER.debug("Setting coordinates in degrees") 

279 alpha = np.degrees(alpha) 

280 delta = np.degrees(delta) 

281 if not self.ctx['colat']: 

282 delta = 90. - delta 

283 else: 

284 LOGGER.debug("Setting coordinates in co-latitude") 

285 else: 

286 if not self.ctx['colat']: 

287 delta = np.pi / 2. - delta 

288 else: 

289 LOGGER.debug("Setting coordinates in co-latitude") 

290 

291 self.alpha = alpha 

292 self.delta = delta 

293 else: 

294 self.alpha = alphas[:] 

295 self.delta = deltas[:] 

296 

297 def set_ellipticities(self, e1s, e2s): 

298 """ 

299 Sets elliticity arrays 

300 

301 :param e1s: First ellipticity component 

302 :param e2s: Second ellipticity component 

303 """ 

304 if (len(e1s) == 0): 

305 LOGGER.warning( 

306 "Passed empty list for first ellipicities. Setting all to 0") 

307 self.e1s = np.zeros_like(self.alpha) 

308 if (len(e2s) == 0): 

309 LOGGER.warning( 

310 "Passed empty list for second ellipicities. Setting all to 0") 

311 self.e2s = np.zeros_like(self.alpha) 

312 else: 

313 self.e1s = e1s[:] 

314 self.e2s = e2s[:] 

315 

316 def set_redhift_bins(self, redshift_bins): 

317 """ 

318 Sets redshift bin array 

319 

320 :param redshift_bins: Redshift bins for the stored objects. 

321 If empty set to 0 

322 """ 

323 if len(redshift_bins) == 0: 

324 LOGGER.warning("Passed empty list for redshifts. Setting all to 0") 

325 self.redshift_bin = np.zeros_like(self.alpha, dtype=int) 

326 else: 

327 self.redshift_bin = redshift_bins[:] 

328 

329 def set_weights(self, weights): 

330 """ 

331 Sets weight array 

332 

333 :param weights: weights for the stored objects. If empty set to 1. 

334 """ 

335 if len(weights) == 0: 

336 LOGGER.warning("Passed empty list for weights. Setting all to 1") 

337 self.weights = np.ones_like(self.alpha, dtype=int) 

338 else: 

339 self.weights = weights[:] 

340 

341 # METHODS 

342 def clip_ellipticities(self, clip_value): 

343 """ 

344 Clips ellipticities at a provided value. 

345 All objects exceeding sqrt(e1^2+e2^2) > clip_value are removed. 

346 

347 :param clip_value: The maximum absolute values 

348 of the ellipticities that is allowed 

349 """ 

350 

351 if len(self.e1s) == 0: 

352 raise Exception("Cannot clip ellipticities because none are set.") 

353 

354 # clipping 

355 clip_value = float(clip_value) 

356 gamma_mag = np.sqrt(self.e1s**2. + self.e2s**2.) 

357 flagged = np.where(gamma_mag > clip_value)[0] 

358 

359 LOGGER.info( 

360 f"Clipped ellipticity to {clip_value} -> " 

361 f"Removed {len(flagged)/len(self.e1s) * 100.}% of objects.") 

362 

363 # remove flagged elements 

364 self.e1s = np.delete(self.e1s, flagged) 

365 self.e2s = np.delete(self.e2s, flagged) 

366 if len(self.alpha) > 0: 

367 self.alpha = np.delete(self.alpha, flagged) 

368 self.delta = np.delete(self.delta, flagged) 

369 for k in self.pixels.keys(): 

370 self.pixels[k] = np.delete(self.pixels[k], flagged) 

371 self.redshift_bin = np.delete(self.redshift_bin, flagged) 

372 self.weights = np.delete(self.weights, flagged) 

373 for field in self.extra_fields: 

374 foo = np.delete(getattr(self, field), flagged) 

375 setattr(self, field, foo) 

376 

377 def rotate_catalog(self, alpha_rot, delta_rot, store_pixels=False, 

378 mirror=False, rotate_ellipticities=False): 

379 """ 

380 Rotates the galaxy coordinates on the sky. Can overwrite coordinates 

381 or just store the pixel coordinates of the rotated survey. 

382 

383 :param alpha_rot: The angle along the first coordinates 

384 to rotate (in radians) 

385 :param delta_rot: The angle along the second coordinates 

386 to rotate (in radians) 

387 :param store_pixels: If False the old coordinates are overwritten. 

388 If True old coordinates are conserved and only the 

389 pixel numbers of the rotated objects are stored in 

390 self.pixels. Allows holding many rotated catalogs 

391 without requiring too much memory. 

392 :param mirror: If True mirrors all coordinates on the equator 

393 :param rotate_ellipticities: If True rotates the ellipticities as well. 

394 Only supported for store_pixels=False 

395 and mirror=False. 

396 """ 

397 

398 alpha, delta = self._rotate_coordinates( 

399 alpha_rot, delta_rot, mirror) 

400 

401 if store_pixels: 

402 pixs = self._pixelize(alpha, delta, bin=0, full_output=False) 

403 key = len(self.pixels.keys()) 

404 self.pixels[key] = pixs 

405 LOGGER.info( 

406 "The pixel locations of the rotated catalog have been " 

407 "stored but the coordinates were rotated back") 

408 else: 

409 self.alpha, self.delta = alpha, delta 

410 if rotate_ellipticities: 

411 if store_pixels: 

412 raise Exception( 

413 "store_pixel is not supported " 

414 "if rotate_ellipticities is used.") 

415 if len(self.e1s) > 0: 

416 self.e1s, self.e2s = self._rotate_ellipticities( 

417 alpha_rot, delta_rot, mirror) 

418 

419 def get_object_pixels(self, bin=0, alpha_rot=0.0, delta_rot=0.0, pix=-1, 

420 mirror=False): 

421 """ 

422 Returns the pixel of each object, all unique pixels where the 

423 objects are located and the object indices in each 

424 pixel. Can optionally rotate the catalog before returning products. 

425 

426 :param bin: Redshift bin to use (0=full sample) 

427 :param alpha_rot: The angle along the first coordinates to rotate 

428 :param delta_rot: The angle along the second coordinates to rotate 

429 :param pix: Alternative operation mode. Uses the stored pixel catalogs 

430 instead of the actual coordinates (-1 not used, greater or 

431 equal than 0 indicates the map to use) 

432 :param mirror: If True mirrors coordinates on equator 

433 :return: Three lists containing the different pixelization information 

434 """ 

435 

436 # pixelize the coordinates 

437 if pix < -0.5: 

438 if len(self.alpha) == 0: 

439 raise Exception( 

440 "No coordinates are set. " 

441 "Cannot convert to pixel locations.") 

442 

443 # rotate catalog if appicable 

444 if (not np.isclose(alpha_rot, 0.0)) | \ 

445 (not np.isclose(delta_rot, 0.0)) | \ 

446 (mirror): 

447 alpha, delta = self._rotate_coordinates( 

448 alpha_rot, delta_rot, mirror) 

449 else: 

450 alpha = self.alpha 

451 delta = self.delta 

452 unis, indices, pixs = self._pixelize( 

453 alpha, delta, bin, full_output=True) 

454 else: 

455 pixs = self.pixels[pix] 

456 unis, indices = np.unique( 

457 pixs, return_inverse=True) 

458 

459 return (pixs, unis, indices) 

460 

461 def get_map(self, type, bin=0, alpha_rot=0.0, delta_rot=0.0, 

462 pix=-1, trimming=True, mirror=False, normalize=True, 

463 method='KS'): 

464 """ 

465 Creates number count, weight, shear and convergence maps 

466 or custom field map from the catalog. 

467 Can optionally rotate the catalog before map creation. 

468 :param type: Can be counts, weights, ellipticities, convergence, 

469 a custom field or a list of them. 

470 Indicates types of the maps that are returned. 

471 The object weights are applied in the creation 

472 of all the maps except for the counts and weights map. 

473 Note that ellipticities and convergence return 

474 two maps each (e1, e2) and (E, B) respectively. 

475 :param bin: Redshift bin to use (0=full sample) 

476 :param alpha_rot: The angle along the first coordinates to rotate 

477 :param delta_rot: The angle along the second coordinates to rotate 

478 :param pix: Alternative operation mode. Uses the stored pixel catalogs 

479 instead of the actual coordinates (-1 not used, greater or 

480 equal than 0 indicates the instance to use) 

481 :param trimming: If True applies stricter masking to convergence map 

482 (some values on edges of survey mask can be 

483 very off due to spherical harmonics transformation). 

484 :param mirror: If True mirrors coordinates on equator 

485 :param normalize: If True all maps except for the counts and weights 

486 maps are normalized by the sum of all object weights per pixel. 

487 :param method: Mass mapping method for convergence maps. 

488 At the moment can only be KS (Kaiser-Squires). 

489 :return: The desired maps in a list (map1, map2, ...) 

490 """ 

491 

492 pixs, unis, indices = self.get_object_pixels( 

493 bin, alpha_rot, delta_rot, pix, mirror) 

494 

495 # Do not touch! 

496 sign_flip = True 

497 

498 if not isinstance(type, list): 

499 type = [type] 

500 

501 # check if all desired fields are present 

502 output_type = "Returning maps in order: " 

503 for field in type: 

504 if (field != 'convergence') & (field != 'ellipticities') \ 

505 & (field != 'weights') \ 

506 & (field != 'counts'): 

507 try: 

508 getattr(self, field) 

509 output_type += f"{field}, " 

510 except AttributeError: 

511 type.remove(field) 

512 LOGGER.warning( 

513 "Did not find field {}. Skipping...".format(field)) 

514 else: 

515 output_type += f"{field}, " 

516 output_type = output_type[:-2] 

517 

518 output = [] 

519 

520 idx = self._get_redshift_bool(bin) 

521 

522 # calculate the weights map 

523 weights = np.zeros(hp.pixelfunc.nside2npix(self.ctx['NSIDE']), 

524 dtype=float) 

525 weights[unis] += np.bincount(indices, weights=self.weights[idx]) 

526 

527 if 'counts' in type: 

528 # calculate the count map 

529 LOGGER.debug("Calculating object count map") 

530 count_weights = np.zeros(pixs.size) 

531 count_weights[idx] = 1.0 

532 counts = np.zeros(hp.pixelfunc.nside2npix(self.ctx['NSIDE']), 

533 dtype=float) 

534 counts[unis] += np.bincount(indices, weights=count_weights) 

535 output.append(counts) 

536 

537 if 'weights' in type: 

538 LOGGER.debug("Calculating weight map") 

539 output.append(weights) 

540 

541 if 'ellipticities' in type: 

542 LOGGER.debug("Calculating ellipticity maps") 

543 # calculate the shear maps 

544 ellipticity_map_1 = np.zeros( 

545 hp.pixelfunc.nside2npix(self.ctx['NSIDE']), 

546 dtype=self.ctx['prec']) 

547 ellipticity_map_1[unis] += np.bincount( 

548 indices, weights=self.weights[idx] * self.e1s[idx]) 

549 ellipticity_map_2 = np.zeros( 

550 hp.pixelfunc.nside2npix(self.ctx['NSIDE']), 

551 dtype=self.ctx['prec']) 

552 ellipticity_map_2[unis] += np.bincount( 

553 indices, weights=self.weights[idx] * self.e2s[idx]) 

554 

555 mask = np.logical_not(weights > 0.0) 

556 

557 if normalize: 

558 widx = np.logical_not(mask) 

559 ellipticity_map_1[widx] = \ 

560 ellipticity_map_1[widx] / weights[widx] 

561 ellipticity_map_2[widx] = \ 

562 ellipticity_map_2[widx] / weights[widx] 

563 

564 # masking 

565 ellipticity_map_1[mask] = hp.pixelfunc.UNSEEN 

566 ellipticity_map_2[mask] = hp.pixelfunc.UNSEEN 

567 

568 output.append(ellipticity_map_1) 

569 output.append(ellipticity_map_2) 

570 

571 if 'convergence' in type: 

572 LOGGER.debug("Calculating convergence maps") 

573 if 'ellipticities' not in type: 

574 # calculate the shear maps 

575 ellipticity_map_1 = np.zeros( 

576 hp.pixelfunc.nside2npix(self.ctx['NSIDE']), 

577 dtype=self.ctx['prec']) 

578 ellipticity_map_1[unis] += np.bincount( 

579 indices, weights=self.weights[idx] * self.e1s[idx]) 

580 ellipticity_map_2 = np.zeros( 

581 hp.pixelfunc.nside2npix(self.ctx['NSIDE']), 

582 dtype=self.ctx['prec']) 

583 ellipticity_map_2[unis] += np.bincount( 

584 indices, weights=self.weights[idx] * self.e2s[idx]) 

585 

586 mask = np.logical_not(weights > 0.0) 

587 

588 if normalize: 

589 widx = np.logical_not(mask) 

590 ellipticity_map_1[widx] = \ 

591 ellipticity_map_1[widx] / weights[widx] 

592 ellipticity_map_2[widx] = \ 

593 ellipticity_map_2[widx] / weights[widx] 

594 

595 # masking 

596 ellipticity_map_1[mask] = hp.pixelfunc.UNSEEN 

597 ellipticity_map_2[mask] = hp.pixelfunc.UNSEEN 

598 # calculate convergence maps 

599 m_kappa_E, m_kappa_B = self._calc_convergence_map( 

600 ellipticity_map_1, ellipticity_map_2, weights, bin, 

601 trimming, sign_flip, method, 

602 unis, indices, idx, normalize) 

603 

604 output.append(m_kappa_E) 

605 output.append(m_kappa_B) 

606 

607 for field in type: 

608 if (field != 'convergence') & (field != 'ellipticities') \ 

609 & (field != 'weights') \ 

610 & (field != 'counts'): 

611 LOGGER.debug(f"Calculating {field} map") 

612 map_1 = np.zeros(hp.pixelfunc.nside2npix(self.ctx['NSIDE']), 

613 dtype=self.ctx['prec']) 

614 map_1[unis] += np.bincount( 

615 indices, 

616 weights=self.weights[idx] * getattr(self, field)[idx]) 

617 

618 mask = np.logical_not(weights > 0.0) 

619 if normalize: 

620 widx = np.logical_not(mask) 

621 map_1[widx] = map_1[widx] / weights[widx] 

622 

623 # masking 

624 map_1[mask] = hp.pixelfunc.UNSEEN 

625 output.append(map_1) 

626 

627 LOGGER.debug(output_type) 

628 return output 

629 

630 def get_mask(self, bin=0, alpha_rot=0.0, delta_rot=0.0, mirror=False): 

631 """ 

632 Returns binary survey mask. Can optionally rotate catalog on sky. 

633 

634 :param bin: Redshift bin to use (0=full sample) 

635 :param alpha_rot: The angle along the first coordinates to rotate 

636 :param delta_rot: The angle along the second coordinates to rotate 

637 :param mirror: If True mirrors coordinates on equator 

638 :return: A Healpix mask (0 outside and 1 inside of mask) 

639 """ 

640 

641 weights = self.get_map('weights', bin, alpha_rot, 

642 delta_rot, mirror=mirror)[0] 

643 mask = weights > 0.0 

644 

645 return mask 

646 

647 def get_trimmed_mask(self, bin=0, alpha_rot=0.0, delta_rot=0.0, 

648 mirror=False, accepted_error=0.05, smoothing=21.1): 

649 """ 

650 Given a Healpix mask creates a trimmed mask, meaning it removes some 

651 pixels which are distored by edge 

652 effects in shear->convergence procedure. 

653 Optionally allows rotation of catalog. 

654 

655 :param bin: Redshift bin to use (0=full sample) 

656 :param alpha_rot: The angle along the first coordinates to rotate 

657 :param delta_rot: The angle along the second coordinates to rotate 

658 :param mirror: If True mirrors coordinates on equator 

659 :param accepted_error: Maximum fractional error of pixels after 

660 smoothing for trimming 

661 :param smoothing: FWHM of Gaussian kernel for smoothing in arcmin. 

662 :return: The trimmed Healpix mask 

663 """ 

664 mask = self.get_mask(bin=bin, alpha_rot=alpha_rot, 

665 delta_rot=delta_rot, mirror=mirror) 

666 mask = self._trimming(mask, accepted_error=accepted_error, 

667 smoothing=smoothing) 

668 

669 return mask 

670 

671 def generate_shape_noise(self, seed=None, bin=0, fields=None, 

672 input_fields=None): 

673 """ 

674 Generates shape noise which resembles the noise in the original 

675 catalog by randomly rotating the ellipticities in place. 

676 

677 :param seed: Seeding for the random generation of the rotation. 

678 :param bin: Redshift bin to use (0=full sample) 

679 :param fields: Can also provide a list of fields to randomize. 

680 If two field names are given as a list they are 

681 intrepreted as the two components of a spin-2 

682 field and rotated like the ellipticities, 

683 otherwise a scalar field is assumed and 

684 the noise is drawn from its global distribution. 

685 :param input_fields: Provide list of fields directly instead of using 

686 internal fields. 

687 :return: The new randomized fields 

688 """ 

689 

690 # Seeding 

691 if seed is not None: 

692 np.random.seed(int(seed)) 

693 LOGGER.debug(f"Seeded random generater with seed {seed}") 

694 

695 if input_fields is not None: 

696 LOGGER.debug("Not using internal fields but input_fields!") 

697 field_1 = input_fields[0] 

698 idx = np.ones(len(field_1), dtype=bool) 

699 try: 

700 field_2 = input_fields[1] 

701 except KeyError: 

702 pass 

703 else: 

704 idx = self._get_redshift_bool(bin) 

705 if (not isinstance(fields, list)) & (fields is not None): 

706 fields = [fields] 

707 

708 if fields is None: 

709 LOGGER.debug("Using ellipticities") 

710 field_1 = self.e1s 

711 field_2 = self.e2s 

712 elif isinstance(fields, list): 

713 if isinstance(fields[0], str): 

714 field_1 = getattr(self, fields[0]) 

715 else: 

716 raise Exception( 

717 "First field must be of type string, " 

718 "but received {}".format(type(fields[0]))) 

719 try: 

720 if isinstance(fields[1], str): 

721 field_2 = getattr(self, fields[1]) 

722 else: 

723 raise Exception( 

724 "Second field must be of type string, " 

725 "but received {}".format(type(fields[1]))) 

726 except IndexError: 

727 field_2 = None 

728 else: 

729 raise Exception( 

730 "Fields argument must be of type string or " 

731 "list of strings, but received {}".format(type(fields))) 

732 

733 if field_2 is not None: 

734 LOGGER.debug("Assuming spin2 field") 

735 # Draw random phase 

736 ell = len(field_1[idx]) 

737 

738 # Generate random phase 

739 rad, ang = self._from_rec_to_polar( 

740 field_1[idx] + field_2[idx] * 1j) 

741 ang += np.pi 

742 ang = (ang + np.random.uniform(0.0, 2.0 * np.pi, ell) 

743 ) % (2.0 * np.pi) - np.pi 

744 

745 noise_shear = self._from_polar_to_rec(rad, ang) 

746 

747 noise_shear_1 = np.real(noise_shear) 

748 noise_shear_2 = np.imag(noise_shear) 

749 return noise_shear_1, noise_shear_2 

750 else: 

751 LOGGER.debug("Assuming scalar field") 

752 # If only one field passed, draw from global distribution of pixels 

753 interp_steps = 10000 

754 min = np.min(field_1) 

755 max = np.min(field_1) 

756 

757 pdf = stats.binned_statistic( 

758 field_1[idx], 

759 np.ones_like(field_1[idx]), 

760 statistic=sum, bins=interp_steps, range=(min, max))[0] 

761 

762 # normalize 

763 pdf /= np.sum(pdf) 

764 # get cumulative function 

765 pdf = np.cumsum(pdf) 

766 

767 # create interplator for inverted function 

768 steps = np.linspace(min, max, interp_steps) 

769 

770 pdf_inv = scipy.interpolate.interp1d(pdf, steps) 

771 

772 noise = pdf_inv(np.random.random(len(field_1[idx]))) 

773 return [noise] 

774 

775 def generate_shape_noise_map(self, seed=None, pix=-1, bin=0, 

776 normalize=True, fields=None): 

777 """ 

778 Generates shape noise maps. 

779 :param seed: Seeding for the random generation of the rotation 

780 :param pix: Alternative operation mode. Uses the stored pixel catalogs 

781 instead of the actual coordinates(-1 not used, greater or 

782 equal than 0 indicates the map to use) 

783 :param bin: Redshift bin to use (0=full sample) 

784 :param normalize: If True the maps 

785 are normalized by the sum of all object weights per pixel. 

786 :param fields: Can also provide a list of fields to randomize. 

787 If two field names are given as a list they are 

788 intrepreted as the two components of a spin-2 

789 field and rotated like the ellipticities, 

790 otherwise a scalar field is assumed and 

791 the noise is drawn from its global distribution. 

792 :return: Two shape noise maps. 

793 """ 

794 

795 out = self.generate_shape_noise(seed, bin, fields=fields) 

796 

797 pixs, unis, indices = self.get_object_pixels( 

798 bin, alpha_rot=0.0, delta_rot=0.0, pix=pix, mirror=False) 

799 

800 idx = self._get_redshift_bool(bin) 

801 

802 # calulate weights maps 

803 weights = np.zeros(hp.pixelfunc.nside2npix(self.ctx['NSIDE']), 

804 dtype=float) 

805 weights[unis] += np.bincount(indices, weights=self.weights[idx]) 

806 mask = np.logical_not(weights > 0.0) 

807 

808 output = [] 

809 for field in out: 

810 noise_field_map = np.zeros( 

811 hp.pixelfunc.nside2npix(self.ctx['NSIDE']), 

812 dtype=self.ctx['prec']) 

813 

814 noise_field_map[unis] += np.bincount( 

815 indices, weights=self.weights[idx] * field) 

816 

817 if normalize: 

818 widx = np.logical_not(mask) 

819 noise_field_map[widx] = noise_field_map[widx] / weights[widx] 

820 

821 # masking 

822 noise_field_map[mask] = hp.pixelfunc.UNSEEN 

823 

824 output.append(noise_field_map) 

825 

826 return output 

827 

828 ################################## 

829 # HELPER FUNCTIONS 

830 ################################## 

831 def _pixelize(self, alpha, delta, bin=0, full_output=False): 

832 """ 

833 Converts angular coordinates into HEalpix map 

834 """ 

835 

836 # Handling for 1 sized arrays 

837 if isinstance(alpha, float): 

838 alpha = [alpha] 

839 if isinstance(delta, float): 

840 delta = [delta] 

841 

842 idx = self._get_redshift_bool(bin) 

843 

844 # converting coordinates to HealPix convention 

845 if self.ctx['degree']: 

846 alpha = np.radians(alpha) 

847 delta = np.radians(delta) 

848 if not self.ctx['colat']: 

849 delta = np.pi / 2. - delta 

850 pix = hp.pixelfunc.ang2pix(self.ctx['NSIDE'], delta[idx], alpha[idx]) 

851 pix = pix.astype(int) 

852 

853 # converting coordinates back 

854 if self.ctx['degree']: