getting_started_3.py 4.82 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
14
# Copyright(C) 2013-2019 Max-Planck-Society
15
#
16
17
18
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.

############################################################
Philipp Arras's avatar
Philipp Arras committed
19
# Non-linear tomography
Torsten Ensslin's avatar
Torsten Ensslin committed
20
21
22
23
24
25
#
# The signal is a sigmoid-normal distributed field.
# The data is the field integrated along lines of sight that are
# randomly (set mode=0) or radially (mode=1) distributed
#
# Demo takes a while to compute
26
#############################################################
27

Philipp Arras's avatar
Philipp Arras committed
28
29
import sys

Jakob Knollmueller's avatar
Jakob Knollmueller committed
30
31
import numpy as np

Martin Reinecke's avatar
5->6    
Martin Reinecke committed
32
import nifty6 as ift
Philipp Arras's avatar
Philipp Arras committed
33

Jakob Knollmueller's avatar
Jakob Knollmueller committed
34

Philipp Arras's avatar
Philipp Arras committed
35
def random_los(n_los):
36
    starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
Torsten Ensslin's avatar
Torsten Ensslin committed
37
    ends = list(np.random.uniform(0, 1, (n_los, 2)).T)
Philipp Arras's avatar
Philipp Arras committed
38
39
40
41
42
    return starts, ends


def radial_los(n_los):
    starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
Torsten Ensslin's avatar
Torsten Ensslin committed
43
    ends = list(0.5 + 0*np.random.uniform(0, 1, (n_los, 2)).T)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
44
45
    return starts, ends

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
46

Jakob Knollmueller's avatar
Jakob Knollmueller committed
47
if __name__ == '__main__':
Philipp Arras's avatar
Philipp Arras committed
48
    np.random.seed(420)
Philipp Arras's avatar
Philipp Arras committed
49

Torsten Ensslin's avatar
Torsten Ensslin committed
50
51
    # Choose between random line-of-sight response (mode=0) and radial lines
    # of sight (mode=1)
Lukas Platz's avatar
Lukas Platz committed
52
53
54
55
    if len(sys.argv) == 2:
        mode = int(sys.argv[1])
    else:
        mode = 0
Philipp Arras's avatar
Philipp Arras committed
56
    filename = "getting_started_3_mode_{}_".format(mode) + "{}.png"
Philipp Arras's avatar
Philipp Arras committed
57

58
    position_space = ift.RGSpace([128, 128])
59

60
    cfmaker = ift.CorrelatedFieldMaker.make(1e-3, 1e-6, '')
61
    cfmaker.add_fluctuations(position_space,
62
63
                             1., 1e-2, 1, .5, .1, .5, -3, 0.5, '')
    correlated_field = cfmaker.finalize()
64
    A = cfmaker.amplitude
Jakob Knollmueller's avatar
Jakob Knollmueller committed
65

Philipp Arras's avatar
Philipp Arras committed
66
    # Apply a nonlinearity
Jakob Knollmueller's avatar
Jakob Knollmueller committed
67
    signal = ift.sigmoid(correlated_field)
Martin Reinecke's avatar
Martin Reinecke committed
68

Philipp Arras's avatar
Philipp Arras committed
69
    # Build the line-of-sight response and define signal response
Torsten Ensslin's avatar
Torsten Ensslin committed
70
    LOS_starts, LOS_ends = random_los(100) if mode == 0 else radial_los(100)
Philipp Arras's avatar
Philipp Arras committed
71
    R = ift.LOSResponse(position_space, starts=LOS_starts, ends=LOS_ends)
Martin Reinecke's avatar
Martin Reinecke committed
72
    signal_response = R(signal)
Philipp Arras's avatar
Philipp Arras committed
73
74

    # Specify noise
Jakob Knollmueller's avatar
Jakob Knollmueller committed
75
    data_space = R.target
Jakob Knollmueller's avatar
Jakob Knollmueller committed
76
    noise = .001
77
    N = ift.ScalingOperator(data_space, noise)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
78

Philipp Arras's avatar
Philipp Arras committed
79
80
81
    # Generate mock signal and data
    mock_position = ift.from_random('normal', signal_response.domain)
    data = signal_response(mock_position) + N.draw_sample()
Jakob Knollmueller's avatar
Jakob Knollmueller committed
82

Philipp Arras's avatar
Philipp Arras committed
83
    # Minimization parameters
84
85
86
87
    ic_sampling = ift.AbsDeltaEnergyController(
        name='Sampling', deltaE=0.05, iteration_limit=100)
    ic_newton = ift.AbsDeltaEnergyController(
        name='Newton', deltaE=0.5, iteration_limit=35)
88
    minimizer = ift.NewtonCG(ic_newton)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
89

Philipp Arras's avatar
Philipp Arras committed
90
    # Set up likelihood and information Hamiltonian
91
92
    likelihood = (ift.GaussianEnergy(mean=data, inverse_covariance=N.inverse) @
                  signal_response)
93
    H = ift.StandardHamiltonian(likelihood, ic_sampling)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
94

95
96
    initial_mean = ift.MultiField.full(H.domain, 0.)
    mean = initial_mean
Jakob Knollmueller's avatar
Jakob Knollmueller committed
97

98
    plot = ift.Plot()
Philipp Arras's avatar
Philipp Arras committed
99
    plot.add(signal(mock_position), title='Ground Truth')
Martin Reinecke's avatar
merge    
Martin Reinecke committed
100
    plot.add(R.adjoint_times(data), title='Data')
Philipp Arras's avatar
Philipp Arras committed
101
    plot.add([A.force(mock_position)], title='Power Spectrum')
Lukas Platz's avatar
Lukas Platz committed
102
    plot.output(ny=1, nx=3, xsize=24, ysize=6, name=filename.format("setup"))
Jakob Knollmueller's avatar
Jakob Knollmueller committed
103

Jakob Knollmueller's avatar
Jakob Knollmueller committed
104
    # number of samples used to estimate the KL
105
    N_samples = 20
Philipp Arras's avatar
Philipp Arras committed
106
107

    # Draw new samples to approximate the KL five times
108
    for i in range(5):
Philipp Arras's avatar
Philipp Arras committed
109
        # Draw new samples and minimize KL
110
        KL = ift.MetricGaussianKL(mean, H, N_samples)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
111
        KL, convergence = minimizer(KL)
112
        mean = KL.position
Philipp Arras's avatar
Philipp Arras committed
113
114

        # Plot current reconstruction
115
        plot = ift.Plot()
Martin Reinecke's avatar
merge    
Martin Reinecke committed
116
        plot.add(signal(KL.position), title="reconstruction")
Philipp Arras's avatar
Philipp Arras committed
117
        plot.add([A.force(KL.position), A.force(mock_position)], title="power")
Lukas Platz's avatar
Lukas Platz committed
118
        plot.output(ny=1, ysize=6, xsize=16,
Philipp Arras's avatar
Philipp Arras committed
119
                    name=filename.format("loop_{:02d}".format(i)))
Jakob Knollmueller's avatar
Jakob Knollmueller committed
120

Philipp Arras's avatar
Philipp Arras committed
121
    # Draw posterior samples
122
    KL = ift.MetricGaussianKL(mean, H, N_samples)
Martin Reinecke's avatar
Martin Reinecke committed
123
    sc = ift.StatCalculator()
124
    for sample in KL.samples:
Philipp Arras's avatar
Philipp Arras committed
125
        sc.add(signal(sample + KL.position))
Philipp Arras's avatar
Philipp Arras committed
126
127

    # Plotting
Lukas Platz's avatar
Lukas Platz committed
128
    filename_res = filename.format("results")
Philipp Arras's avatar
Philipp Arras committed
129
    plot = ift.Plot()
Martin Reinecke's avatar
merge    
Martin Reinecke committed
130
131
    plot.add(sc.mean, title="Posterior Mean")
    plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation")
Martin Reinecke's avatar
Martin Reinecke committed
132

Philipp Arras's avatar
Philipp Arras committed
133
    powers = [A.force(s + KL.position) for s in KL.samples]
Martin Reinecke's avatar
merge    
Martin Reinecke committed
134
    plot.add(
Philipp Arras's avatar
Philipp Arras committed
135
136
        powers + [A.force(KL.position),
                  A.force(mock_position)],
Lukas Platz's avatar
Lukas Platz committed
137
138
        title="Sampled Posterior Power Spectrum",
        linewidth=[1.]*len(powers) + [3., 3.])
Lukas Platz's avatar
Lukas Platz committed
139
    plot.output(ny=1, nx=3, xsize=24, ysize=6, name=filename_res)
Philipp Arras's avatar
Philipp Arras committed
140
    print("Saved results as '{}'.".format(filename_res))