getting_started_3.py 5.53 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):
Martin Reinecke's avatar
Martin Reinecke committed
36
37
    starts = list(ift.random.current_rng().random((n_los, 2)).T)
    ends = list(ift.random.current_rng().random((n_los, 2)).T)
Philipp Arras's avatar
Philipp Arras committed
38
39
40
41
    return starts, ends


def radial_los(n_los):
Martin Reinecke's avatar
Martin Reinecke committed
42
43
    starts = list(ift.random.current_rng().random((n_los, 2)).T)
    ends = list(0.5 + 0*ift.random.current_rng().random((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__':
Torsten Ensslin's avatar
Torsten Ensslin committed
48
49
    # Choose between random line-of-sight response (mode=0) and radial lines
    # of sight (mode=1)
Lukas Platz's avatar
Lukas Platz committed
50
51
52
53
    if len(sys.argv) == 2:
        mode = int(sys.argv[1])
    else:
        mode = 0
Philipp Arras's avatar
Philipp Arras committed
54
    filename = "getting_started_3_mode_{}_".format(mode) + "{}.png"
Philipp Arras's avatar
Philipp Arras committed
55

56
    position_space = ift.RGSpace([128, 128])
57

Lukas Platz's avatar
Lukas Platz committed
58
    cfmaker = ift.CorrelatedFieldMaker.make(
Lukas Platz's avatar
Lukas Platz committed
59
60
61
        offset_mean =      0.0,  # 0.
        offset_std_mean = 1e-3,  # 1e-3
        offset_std_std =  1e-6,  # 1e-6
Lukas Platz's avatar
Lukas Platz committed
62
63
64
65
        prefix = '')

    fluctuations_dict = {
        # Amplitude of the fluctuations
Lukas Platz's avatar
Lukas Platz committed
66
        'fluctuations_mean':   2.0,  # 1.0
Lukas Platz's avatar
Lukas Platz committed
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
        'fluctuations_stddev': 1.0,  # 1e-2

        # Smooth variation speed
        'flexibility_mean':   2.5,  # 1.0
        'flexibility_stddev': 1.0,  # 0.5

        # How strong the ragged component of the spectrum is
        # (Ratio of Wiener process and integrated Wiener process ?)
        'asperity_mean':   0.5,  # 0.1
        'asperity_stddev': 0.5,  # 0.5

        # Slope of linear spectrum component
        'loglogavgslope_mean': -2.0,  # -3.0
        'loglogavgslope_stddev': 0.5  #  0.5
    }
    cfmaker.add_fluctuations(position_space, **fluctuations_dict)

84
    correlated_field = cfmaker.finalize()
85
    A = cfmaker.amplitude
Jakob Knollmueller's avatar
Jakob Knollmueller committed
86

Philipp Arras's avatar
Philipp Arras committed
87
    # Apply a nonlinearity
Jakob Knollmueller's avatar
Jakob Knollmueller committed
88
    signal = ift.sigmoid(correlated_field)
Martin Reinecke's avatar
Martin Reinecke committed
89

Philipp Arras's avatar
Philipp Arras committed
90
    # Build the line-of-sight response and define signal response
Torsten Ensslin's avatar
Torsten Ensslin committed
91
    LOS_starts, LOS_ends = random_los(100) if mode == 0 else radial_los(100)
Philipp Arras's avatar
Philipp Arras committed
92
    R = ift.LOSResponse(position_space, starts=LOS_starts, ends=LOS_ends)
Martin Reinecke's avatar
Martin Reinecke committed
93
    signal_response = R(signal)
Philipp Arras's avatar
Philipp Arras committed
94
95

    # Specify noise
Jakob Knollmueller's avatar
Jakob Knollmueller committed
96
    data_space = R.target
Jakob Knollmueller's avatar
Jakob Knollmueller committed
97
    noise = .001
98
    N = ift.ScalingOperator(data_space, noise)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
99

Philipp Arras's avatar
Philipp Arras committed
100
101
    # Generate mock signal and data
    mock_position = ift.from_random('normal', signal_response.domain)
102
    data = signal_response(mock_position) + N.draw_sample(dtype=np.float64)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
103

Philipp Arras's avatar
Philipp Arras committed
104
    # Minimization parameters
105
106
107
108
    ic_sampling = ift.AbsDeltaEnergyController(
        name='Sampling', deltaE=0.05, iteration_limit=100)
    ic_newton = ift.AbsDeltaEnergyController(
        name='Newton', deltaE=0.5, iteration_limit=35)
109
    minimizer = ift.NewtonCG(ic_newton)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
110

Philipp Arras's avatar
Philipp Arras committed
111
    # Set up likelihood and information Hamiltonian
112
113
    likelihood = (ift.GaussianEnergy(mean=data, inverse_covariance=N.inverse) @
                  signal_response)
114
    H = ift.StandardHamiltonian(likelihood, ic_sampling)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
115

116
117
    initial_mean = ift.MultiField.full(H.domain, 0.)
    mean = initial_mean
Jakob Knollmueller's avatar
Jakob Knollmueller committed
118

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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
125
    # number of samples used to estimate the KL
126
    N_samples = 20
Philipp Arras's avatar
Philipp Arras committed
127
128

    # Draw new samples to approximate the KL five times
129
    for i in range(5):
Philipp Arras's avatar
Philipp Arras committed
130
        # Draw new samples and minimize KL
131
        KL = ift.MetricGaussianKL(mean, H, N_samples)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
132
        KL, convergence = minimizer(KL)
133
        mean = KL.position
Philipp Arras's avatar
Philipp Arras committed
134
135

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

Philipp Arras's avatar
Philipp Arras committed
142
    # Draw posterior samples
143
    KL = ift.MetricGaussianKL(mean, H, N_samples)
Martin Reinecke's avatar
Martin Reinecke committed
144
    sc = ift.StatCalculator()
145
    for sample in KL.samples:
Philipp Arras's avatar
Philipp Arras committed
146
        sc.add(signal(sample + KL.position))
Philipp Arras's avatar
Philipp Arras committed
147
148

    # Plotting
Lukas Platz's avatar
Lukas Platz committed
149
    filename_res = filename.format("results")
Philipp Arras's avatar
Philipp Arras committed
150
    plot = ift.Plot()
Martin Reinecke's avatar
merge    
Martin Reinecke committed
151
152
    plot.add(sc.mean, title="Posterior Mean")
    plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation")
Martin Reinecke's avatar
Martin Reinecke committed
153

Philipp Arras's avatar
Philipp Arras committed
154
    powers = [A.force(s + KL.position) for s in KL.samples]
Martin Reinecke's avatar
merge    
Martin Reinecke committed
155
    plot.add(
Philipp Arras's avatar
Philipp Arras committed
156
157
        powers + [A.force(KL.position),
                  A.force(mock_position)],
Lukas Platz's avatar
Lukas Platz committed
158
159
        title="Sampled Posterior Power Spectrum",
        linewidth=[1.]*len(powers) + [3., 3.])
Lukas Platz's avatar
Lukas Platz committed
160
    plot.output(ny=1, nx=3, xsize=24, ysize=6, name=filename_res)
Philipp Arras's avatar
Philipp Arras committed
161
    print("Saved results as '{}'.".format(filename_res))