getting_started_2.py 4.24 KB
Newer Older
Natalia Porqueres's avatar
Natalia Porqueres committed
1
2
3
4
5
6
import nifty5 as ift
import sys
import numpy as np
import global_newton as gn
from nifty5.library.nonlinearities import Exponential

Natalia Porqueres's avatar
Natalia Porqueres committed
7
#DEFINE THE SIGNAL:
Natalia Porqueres's avatar
Natalia Porqueres committed
8

Natalia Porqueres's avatar
Natalia Porqueres committed
9
10
11
#Define signal space as a regular grid
#s_space = ift.RGSpace(10024)
s_space = ift.RGSpace([128,128])
Natalia Porqueres's avatar
Natalia Porqueres committed
12

Natalia Porqueres's avatar
Natalia Porqueres committed
13
#Define the harmonic space
Natalia Porqueres's avatar
Natalia Porqueres committed
14
15
h_space = s_space.get_default_codomain()

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

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

#Define positions from a Gaussian distribution
Natalia Porqueres's avatar
Natalia Porqueres committed
23
24
position = ift.from_random('normal', domain)

Natalia Porqueres's avatar
Natalia Porqueres committed
25
#Define a power spectrum
Natalia Porqueres's avatar
Natalia Porqueres committed
26
27
28
def sqrtpspec(k):
    return 16. / (20.+k**2)

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

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

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

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

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

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

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

#Create a sky model by applying the exponential (Poisson)
Natalia Porqueres's avatar
Natalia Porqueres committed
51
sky = ift.PointwiseExponential(logsky)
Natalia Porqueres's avatar
Natalia Porqueres committed
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

#DEFINE THE RESPONSE OPERATOR:

#Define a mask to cover a patch of the real space
mask = np.ones(s_space.shape)
mask[int(s_space.shape[0]/3):int(s_space.shape[0]/3+10)] = 0.

#Convert the mask into a field
mask = ift.Field(s_space,val=mask)

#Create a diagonal matrix corresponding to the mask
M = ift.DiagonalOperator(mask)

#Create the response operator and apply the mask on it
R = ift.ScalingOperator(1., s_space) * M

#CREATE THE MOCK DATA:

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

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

Natalia Porqueres's avatar
Natalia Porqueres committed
77
78
79
80
81
82
83
#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
84
85
data = ift.Field.from_local_data(d_space, data)

Natalia Porqueres's avatar
Natalia Porqueres committed
86
87
88
#RECONSTRUCT THE SIGNAL:

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

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

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

#Define a iteration controller with convergence criteria
Natalia Porqueres's avatar
Natalia Porqueres committed
98
99
ic_samps = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-4)

Natalia Porqueres's avatar
Natalia Porqueres committed
100
#Define a iteration controller for the minimizer
Natalia Porqueres's avatar
Natalia Porqueres committed
101
102
103
ic_newton = ift.GradientNormController(name='Newton', tol_abs_gradnorm=1e-3)
minimizer = ift.RelaxedNewton(ic_newton)

Natalia Porqueres's avatar
Natalia Porqueres committed
104
#Build the Hamiltonian
Natalia Porqueres's avatar
Natalia Porqueres committed
105
106
H = gn.Hamiltonian(likelihood, ic_cg, ic_samps)

Natalia Porqueres's avatar
Natalia Porqueres committed
107
#Iterate the reconstruction
Natalia Porqueres's avatar
Natalia Porqueres committed
108
for _ in range(N):
Natalia Porqueres's avatar
Natalia Porqueres committed
109
    #Draw samples from the curvature
Natalia Porqueres's avatar
Natalia Porqueres committed
110
111
    samples = [H.curvature.draw_sample(from_inverse=True)
               for _ in range(Nsamples)]
Natalia Porqueres's avatar
Natalia Porqueres committed
112
                          
Natalia Porqueres's avatar
Natalia Porqueres committed
113
114
115
116
117
    sc_samplesky = ift.StatCalculator()
    for s in samples:
        sc_samplesky.add(sky.at(s+position).value)
    ift.plot(sc_samplesky.mean, name='sample_mean.png')

Natalia Porqueres's avatar
Natalia Porqueres committed
118
    #Compute the Kullbachh Leibler divergence
Natalia Porqueres's avatar
Natalia Porqueres committed
119
    KL = gn.SampledKullbachLeiblerDivergence(H, samples, ic_cg)
Natalia Porqueres's avatar
Natalia Porqueres committed
120
121
    
    #Update the position
Natalia Porqueres's avatar
Natalia Porqueres committed
122
    position = KL.position
Natalia Porqueres's avatar
Natalia Porqueres committed
123
124
125
126
127
128
129
130
131
132
133
    
    #Obtain the value of the KL and the convergence at the new position
    KL, convergence = minimizer(KL)

#PLOT RESULTS:
   
#Evaluate lambda at final position
lamb_recontructed = lamb.at(KL.position).value

#Evaluate lambda at the data position for comparison
lamb_mock = lamb.at(mock_position).value
Natalia Porqueres's avatar
Natalia Porqueres committed
134

Natalia Porqueres's avatar
Natalia Porqueres committed
135
136
#Plot the data, reconstruction and underlying signal
ift.plot([data, lamb_mock, lamb_recontructed], name="poisson.png",
Natalia Porqueres's avatar
Natalia Porqueres committed
137
138
         label=['Data', 'Mock signal', 'Reconstruction'],
         alpha=[.5, 1, 1])
Natalia Porqueres's avatar
Natalia Porqueres committed
139
140
141
142
143
         
#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
144
145
         name='power_spectrum.png', label=['Mock', 'Reconstruction', 'power_analyze'])