getting_started_2.py 3.55 KB
Newer Older
Natalia Porqueres's avatar
Natalia Porqueres committed
1
2
import nifty5 as ift
import numpy as np
3
# from nifty5.library.nonlinearities import Exponential
Natalia Porqueres's avatar
Natalia Porqueres committed
4

Natalia Porqueres's avatar
Natalia Porqueres committed
5
#DEFINE THE SIGNAL:
Natalia Porqueres's avatar
Natalia Porqueres committed
6

Natalia Porqueres's avatar
Natalia Porqueres committed
7
#Define signal space as a regular grid
8
9
10
#s_space = ift.RGSpace(1024)
# s_space = ift.RGSpace([128,128])
s_space = ift.HPSpace(128)
Natalia Porqueres's avatar
Natalia Porqueres committed
11
#Define the harmonic space
Natalia Porqueres's avatar
Natalia Porqueres committed
12
13
h_space = s_space.get_default_codomain()

Natalia Porqueres's avatar
Natalia Porqueres committed
14
15
16
17
#Prepare Harmonic transformation between the two spaces
HT = ift.HarmonicTransformOperator(h_space, s_space)

#Define domain
Natalia Porqueres's avatar
Natalia Porqueres committed
18
domain = ift.MultiDomain.make({'xi': h_space})
Natalia Porqueres's avatar
Natalia Porqueres committed
19
20

#Define positions from a Gaussian distribution
Natalia Porqueres's avatar
Natalia Porqueres committed
21
position = ift.from_random('normal', domain)
22
Nsamples = 5
Natalia Porqueres's avatar
Natalia Porqueres committed
23
#Define a power spectrum
Natalia Porqueres's avatar
Natalia Porqueres committed
24
def sqrtpspec(k):
25
    return 10. / (20.+k**2)
Natalia Porqueres's avatar
Natalia Porqueres committed
26

Natalia Porqueres's avatar
Natalia Porqueres committed
27
#Define a power space
Natalia Porqueres's avatar
Natalia Porqueres committed
28
p_space = ift.PowerSpace(h_space)
Natalia Porqueres's avatar
Natalia Porqueres committed
29
30

#Define the power distribution between the harmonic and the power spaces
Natalia Porqueres's avatar
Natalia Porqueres committed
31
pd = ift.PowerDistributor(h_space, p_space)
Natalia Porqueres's avatar
Natalia Porqueres committed
32
33

#Create a field with the defined power spectrum
Natalia Porqueres's avatar
Natalia Porqueres committed
34
a = ift.PS_field(p_space, sqrtpspec)
Natalia Porqueres's avatar
Natalia Porqueres committed
35
36

#Define the amplitudes
Natalia Porqueres's avatar
Natalia Porqueres committed
37
38
A = pd(a)

Natalia Porqueres's avatar
Natalia Porqueres committed
39
#Unpack the positions xi from the Multifield
Natalia Porqueres's avatar
Natalia Porqueres committed
40
xi = ift.Variable(position)['xi']
Natalia Porqueres's avatar
Natalia Porqueres committed
41
42

#Multiply the positions by the amplitudes in the harmonic domain
Natalia Porqueres's avatar
Natalia Porqueres committed
43
logsky_h = xi * A
Natalia Porqueres's avatar
Natalia Porqueres committed
44
45

#Transform to the real domain
Natalia Porqueres's avatar
Natalia Porqueres committed
46
logsky = HT(logsky_h)
Natalia Porqueres's avatar
Natalia Porqueres committed
47
48

#Create a sky model by applying the exponential (Poisson)
Natalia Porqueres's avatar
Natalia Porqueres committed
49
sky = ift.PointwiseExponential(logsky)
Natalia Porqueres's avatar
Natalia Porqueres committed
50
51
52
53

#DEFINE THE RESPONSE OPERATOR:

#Define a mask to cover a patch of the real space
54
55
exposure = 1*np.ones(s_space.shape)
# exposure[int(s_space.shape[0]/3):int(s_space.shape[0]/3+10)] = 10.
Natalia Porqueres's avatar
Natalia Porqueres committed
56
57

#Convert the mask into a field
58
exposure = ift.Field(s_space,val=exposure)
Natalia Porqueres's avatar
Natalia Porqueres committed
59
60

#Create a diagonal matrix corresponding to the mask
61
E = ift.DiagonalOperator(exposure)
Natalia Porqueres's avatar
Natalia Porqueres committed
62
63

#Create the response operator and apply the mask on it
64
R = ift.GeometryRemover(s_space) * E
Natalia Porqueres's avatar
Natalia Porqueres committed
65
66
67
68

#CREATE THE MOCK DATA:

#Define the data space
Natalia Porqueres's avatar
Natalia Porqueres committed
69
d_space = R.target[0]
Natalia Porqueres's avatar
Natalia Porqueres committed
70
71
72

#Apply the response operator to the signal
#lamb corresponds to the mean in the Poisson distribution
Natalia Porqueres's avatar
Natalia Porqueres committed
73
74
lamb = R(sky)

Natalia Porqueres's avatar
Natalia Porqueres committed
75
76
77
78
79
80
81
#Draw coordinates of the mock data from a Gaussian distribution
mock_position = ift.from_random('normal', lamb.position.domain)

#Generate mock data from a Poisson distribution using lamb as a mean
data = np.random.poisson(lamb.at(mock_position).value.val.astype(np.float64))

#Store the data as a field
Natalia Porqueres's avatar
Natalia Porqueres committed
82
83
data = ift.Field.from_local_data(d_space, data)

Natalia Porqueres's avatar
Natalia Porqueres committed
84
85
86
#RECONSTRUCT THE SIGNAL:

#Define the positions where we perform the analysis from a Gaussian distribution
Natalia Porqueres's avatar
Natalia Porqueres committed
87
position = ift.from_random('normal', lamb.position.domain)
Natalia Porqueres's avatar
Natalia Porqueres committed
88
89

#Define the Poisson likelihood knowing the mean and the data
Natalia Porqueres's avatar
Natalia Porqueres committed
90
91
likelihood = ift.library.PoissonLogLikelihood(lamb, data)

Natalia Porqueres's avatar
Natalia Porqueres committed
92
#Define a iteration controller with a maximum number of iterations
Natalia Porqueres's avatar
Natalia Porqueres committed
93
ic_cg = ift.GradientNormController(iteration_limit=50)
Natalia Porqueres's avatar
Natalia Porqueres committed
94

Natalia Porqueres's avatar
Natalia Porqueres committed
95

Natalia Porqueres's avatar
Natalia Porqueres committed
96
#Define a iteration controller for the minimizer
Natalia Porqueres's avatar
Natalia Porqueres committed
97
98
99
ic_newton = ift.GradientNormController(name='Newton', tol_abs_gradnorm=1e-3)
minimizer = ift.RelaxedNewton(ic_newton)

Natalia Porqueres's avatar
Natalia Porqueres committed
100
#Build the Hamiltonian
101
102
H = ift.Hamiltonian(likelihood, ic_cg)
H, convergence = minimizer(H)
Natalia Porqueres's avatar
Natalia Porqueres committed
103
104
105
106

#PLOT RESULTS:
   
#Evaluate lambda at final position
107
108
lamb_recontructed = lamb.at(H.position).value
sky_reconstructed = sky.at(H.position).value
Natalia Porqueres's avatar
Natalia Porqueres committed
109
110
#Evaluate lambda at the data position for comparison
lamb_mock = lamb.at(mock_position).value
Natalia Porqueres's avatar
Natalia Porqueres committed
111

Natalia Porqueres's avatar
Natalia Porqueres committed
112
113
#Plot the data, reconstruction and underlying signal
ift.plot([data, lamb_mock, lamb_recontructed], name="poisson.png",
Natalia Porqueres's avatar
Natalia Porqueres committed
114
115
         label=['Data', 'Mock signal', 'Reconstruction'],
         alpha=[.5, 1, 1])
Natalia Porqueres's avatar
Natalia Porqueres committed
116
117
118
119
120
         
#Plot power spectrum for posterior test         
a_mock = a
a_recon = a
ift.plot([a_mock**2, a_recon**2, ift.power_analyze(logsky_h.at(KL.position).value)],
Natalia Porqueres's avatar
Natalia Porqueres committed
121
122
         name='power_spectrum.png', label=['Mock', 'Reconstruction', 'power_analyze'])