getting_started_3.py 5.75 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/>.
#
Philipp Arras's avatar
Philipp Arras committed
14
# Copyright(C) 2013-2020 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
Martin Reinecke committed
32
import nifty7 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

Philipp Arras's avatar
Philipp Arras committed
47
def 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
59
    #  For a detailed showcase of the effects the parameters
    #  of the CorrelatedField model have on the generated fields,
Lukas Platz's avatar
Lukas Platz committed
60
    #  see 'getting_started_4_CorrelatedFields.ipynb'.
Lukas Platz's avatar
Lukas Platz committed
61

Lukas Platz's avatar
Lukas Platz committed
62
    cfmaker = ift.CorrelatedFieldMaker.make(
Philipp Arras's avatar
Philipp Arras committed
63
64
65
66
        offset_mean=      0.0,
        offset_std_mean= 1e-3,
        offset_std_std=  1e-6,
        prefix='')
Lukas Platz's avatar
Lukas Platz committed
67
68

    fluctuations_dict = {
Lukas Platz's avatar
Lukas Platz committed
69
        # Amplitude of field fluctuations
Lukas Platz's avatar
Lukas Platz committed
70
        'fluctuations_mean':   2.0,  # 1.0
Lukas Platz's avatar
Lukas Platz committed
71
72
        'fluctuations_stddev': 1.0,  # 1e-2

Lukas Platz's avatar
Lukas Platz committed
73
74
        # Exponent of power law power spectrum component
        'loglogavgslope_mean': -2.0,  # -3.0
Philipp Arras's avatar
Philipp Arras committed
75
        'loglogavgslope_stddev': 0.5,  # 0.5
Lukas Platz's avatar
Lukas Platz committed
76
77

        # Amplitude of integrated Wiener process power spectrum component
Lukas Platz's avatar
Lukas Platz committed
78
79
80
        'flexibility_mean':   2.5,  # 1.0
        'flexibility_stddev': 1.0,  # 0.5

Lukas Platz's avatar
Lukas Platz committed
81
        # How ragged the integrated Wiener process component is
Lukas Platz's avatar
Lukas Platz committed
82
        'asperity_mean':   0.5,  # 0.1
Lukas Platz's avatar
Lukas Platz committed
83
        'asperity_stddev': 0.5  # 0.5
Lukas Platz's avatar
Lukas Platz committed
84
85
86
    }
    cfmaker.add_fluctuations(position_space, **fluctuations_dict)

87
    correlated_field = cfmaker.finalize()
88
    A = cfmaker.amplitude
Jakob Knollmueller's avatar
Jakob Knollmueller committed
89

Philipp Arras's avatar
Philipp Arras committed
90
    # Apply a nonlinearity
Jakob Knollmueller's avatar
Jakob Knollmueller committed
91
    signal = ift.sigmoid(correlated_field)
Martin Reinecke's avatar
Martin Reinecke committed
92

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

    # Specify noise
Jakob Knollmueller's avatar
Jakob Knollmueller committed
99
    data_space = R.target
Jakob Knollmueller's avatar
Jakob Knollmueller committed
100
    noise = .001
101
    N = ift.ScalingOperator(data_space, noise)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
102

Philipp Arras's avatar
Philipp Arras committed
103
    # Generate mock signal and data
104
    mock_position = ift.from_random(signal_response.domain, 'normal')
Philipp Arras's avatar
Philipp Arras committed
105
    data = signal_response(mock_position) + N.draw_sample_with_dtype(dtype=np.float64)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
106

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

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

119
120
    initial_mean = ift.MultiField.full(H.domain, 0.)
    mean = initial_mean
Jakob Knollmueller's avatar
Jakob Knollmueller committed
121

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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
128
    # number of samples used to estimate the KL
129
    N_samples = 20
Philipp Arras's avatar
Philipp Arras committed
130
131

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

        # Plot current reconstruction
139
        plot = ift.Plot()
Philipp Arras's avatar
Philipp Arras committed
140
141
142
        plot.add(signal(KL.position), title="Latent mean")
        plot.add([A.force(KL.position + ss) for ss in KL.samples],
                 title="Samples power spectrum")
Lukas Platz's avatar
Lukas Platz committed
143
        plot.output(ny=1, ysize=6, xsize=16,
Philipp Arras's avatar
Philipp Arras committed
144
                    name=filename.format("loop_{:02d}".format(i)))
Jakob Knollmueller's avatar
Jakob Knollmueller committed
145

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

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

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


if __name__ == '__main__':
    main()