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'])