/
gaussianity.py
77 lines (47 loc) · 2.03 KB
/
gaussianity.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
"""
Tests the gaussianity of simulated noise maps by measuring their cubic and quartic moments
"""
import sys
from lenstools import ConvergenceMap,Ensemble,GaussianNoiseGenerator
from lenstools.utils.defaults import load_fits_default_convergence
from lenstools.simulations import IGS1
import numpy as np
from astropy.units import deg,arcmin
import logging
from emcee.utils import MPIPool
def generate_and_measure(map_id,generator,power_func):
#Generate the noise map
logging.debug("Processing mock map {0}".format(map_id))
conv_map = generator.fromConvPower(power_func=power_func,seed=map_id,bounds_error=False,fill_value=0.0)
#Measure its moments
return conv_map.moments(connected=True,dimensionless=True)
def measure_from_IGS1(filename):
#Read in the map
logging.debug("Processing IGS1 map {0}".format(filename))
conv_map = ConvergenceMap.load(filename)
#Smooth 1 arcmin
conv_map.smooth(1.0*arcmin,inplace=True)
#Measure the moments
return conv_map.moments(connected=True,dimensionless=True)
logging.basicConfig(level=logging.DEBUG)
try:
pool = MPIPool()
except ValueError:
pool = None
if (pool is not None) and not(pool.is_master()):
pool.wait()
sys.exit(0)
map_mock_ids = range(int(sys.argv[1]))
igs1_set = IGS1(root_path="/Users/andreapetri/Documents/Columbia/spurious_shear/convergence_maps")
map_igs1_ids = igs1_set.getNames(z=1.0,realizations=range(1,int(sys.argv[1])+1))
gen = GaussianNoiseGenerator(shape=(2048,2048),side_angle=3.41*deg,label="convergence")
power_func = np.loadtxt("Data/ee4e-7.txt",unpack=True)
ens_mock = Ensemble.fromfilelist(map_mock_ids)
ens_igs1 = Ensemble.fromfilelist(map_igs1_ids)
ens_mock.load(callback_loader=generate_and_measure,pool=pool,generator=gen,power_func=power_func)
ens_igs1.load(callback_loader=measure_from_IGS1,pool=pool)
if pool is not None:
pool.close()
np.savetxt("moments_mock.txt",np.array([ens_mock.mean(),np.sqrt(ens_mock.covariance().diagonal())]))
np.savetxt("moments_igs1.txt",np.array([ens_igs1.mean(),np.sqrt(ens_igs1.covariance().diagonal())]))
logging.info("Done!")