getting_started_3.py 6.12 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

Philipp Frank's avatar
Philipp Frank committed
34
ift.random.push_sseq_from_seed(27)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
35

Philipp Arras's avatar
Philipp Arras committed
36
def random_los(n_los):
Martin Reinecke's avatar
Martin Reinecke committed
37
38
    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
39
40
41
42
    return starts, ends


def radial_los(n_los):
Martin Reinecke's avatar
Martin Reinecke committed
43
44
    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
45
46
    return starts, ends

Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
47

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

62
63
    args = {
        'offset_mean': 0,
64
        'offset_std': (1e-3, 1e-6),
Lukas Platz's avatar
Lukas Platz committed
65

Lukas Platz's avatar
Lukas Platz committed
66
        # Amplitude of field fluctuations
Philipp Frank's avatar
Philipp Frank committed
67
        'fluctuations': (1., 0.8),  # 1.0, 1e-2
Lukas Platz's avatar
Lukas Platz committed
68

Lukas Platz's avatar
Lukas Platz committed
69
        # Exponent of power law power spectrum component
Philipp Frank's avatar
Philipp Frank committed
70
        'loglogavgslope': (-3., 1),  # -6.0, 1
Lukas Platz's avatar
Lukas Platz committed
71
72

        # Amplitude of integrated Wiener process power spectrum component
Philipp Frank's avatar
Philipp Frank committed
73
        'flexibility': (2, 1.),  # 1.0, 0.5
Lukas Platz's avatar
Lukas Platz committed
74

Lukas Platz's avatar
Lukas Platz committed
75
        # How ragged the integrated Wiener process component is
Philipp Frank's avatar
Philipp Frank committed
76
        'asperity': (0.5, 0.4)  # 0.1, 0.5
Lukas Platz's avatar
Lukas Platz committed
77
78
    }

79
    correlated_field = ift.SimpleCorrelatedField(position_space, **args)
80
    pspec = correlated_field.power_spectrum
Jakob Knollmueller's avatar
Jakob Knollmueller committed
81

Philipp Arras's avatar
Philipp Arras committed
82
    # Apply a nonlinearity
Jakob Knollmueller's avatar
Jakob Knollmueller committed
83
    signal = ift.sigmoid(correlated_field)
Martin Reinecke's avatar
Martin Reinecke committed
84

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

    # Specify noise
Jakob Knollmueller's avatar
Jakob Knollmueller committed
91
    data_space = R.target
Jakob Knollmueller's avatar
Jakob Knollmueller committed
92
    noise = .001
93
    N = ift.ScalingOperator(data_space, noise)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
94

Philipp Arras's avatar
Philipp Arras committed
95
    # Generate mock signal and data
96
    mock_position = ift.from_random(signal_response.domain, 'normal')
Philipp Arras's avatar
Philipp Arras committed
97
    data = signal_response(mock_position) + N.draw_sample_with_dtype(dtype=np.float64)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
98

Philipp Arras's avatar
Philipp Arras committed
99
    # Minimization parameters
100
    ic_sampling = ift.AbsDeltaEnergyController(
Philipp Arras's avatar
Philipp Arras committed
101
        deltaE=0.05, iteration_limit=100)
102
103
    ic_newton = ift.AbsDeltaEnergyController(
        name='Newton', deltaE=0.5, iteration_limit=35)
104
    minimizer = ift.NewtonCG(ic_newton)
Philipp Frank's avatar
Philipp Frank committed
105
106
107
    ic_sampling_nl = ift.AbsDeltaEnergyController(
        name='Sampling', deltaE=0.5, iteration_limit=15, convergence_level=2)
    minimizer_sampling = ift.NewtonCG(ic_sampling_nl)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
108

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

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

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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
123
    # number of samples used to estimate the KL
Philipp Frank's avatar
Philipp Frank committed
124
    N_samples = 10
Philipp Arras's avatar
Philipp Arras committed
125

Philipp Frank's avatar
Philipp Frank committed
126
127
128
129
130
    # Draw new samples to approximate the KL six times
    for i in range(6):
        if i==5:
            # Double the number of samples in the last step for better statistics
            N_samples = 2*N_samples
Philipp Arras's avatar
Philipp Arras committed
131
        # Draw new samples and minimize KL
Philipp Frank's avatar
Philipp Frank committed
132
        KL = ift.GeoMetricKL(mean, H, N_samples, minimizer_sampling, True)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
133
        KL, convergence = minimizer(KL)
134
        mean = KL.position
Philipp Arras's avatar
Philipp Arras committed
135
        ift.extra.minisanity(data, lambda x: N.inverse, signal_response,
136
                             KL.position, KL.samples)
Philipp Arras's avatar
Philipp Arras committed
137
138

        # Plot current reconstruction
139
        plot = ift.Plot()
Philipp Frank's avatar
Philipp Frank committed
140
        plot.add(signal(KL.position), title="Latent mean", zmin = 0, zmax = 1)
141
        plot.add([pspec.force(KL.position + ss) for ss in KL.samples],
Philipp Arras's avatar
Philipp Arras committed
142
                 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

Martin Reinecke's avatar
Martin Reinecke committed
146
    sc = ift.StatCalculator()
147
    for sample in KL.samples:
Philipp Arras's avatar
Philipp Arras committed
148
        sc.add(signal(sample + KL.position))
Philipp Arras's avatar
Philipp Arras committed
149
150

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

156
    powers = [pspec.force(s + KL.position) for s in KL.samples]
157
158
    sc = ift.StatCalculator()
    for pp in powers:
Philipp Frank's avatar
Philipp Frank committed
159
        sc.add(pp.log())
Martin Reinecke's avatar
merge    
Martin Reinecke committed
160
    plot.add(
161
        powers + [pspec.force(mock_position),
Philipp Frank's avatar
Philipp Frank committed
162
                  pspec.force(KL.position), sc.mean.exp()],
Lukas Platz's avatar
Lukas Platz committed
163
        title="Sampled Posterior Power Spectrum",
164
165
        linewidth=[1.]*len(powers) + [3., 3., 3.],
        label=[None]*len(powers) + ['Ground truth', 'Posterior latent mean', 'Posterior mean'])
Lukas Platz's avatar
Lukas Platz committed
166
    plot.output(ny=1, nx=3, xsize=24, ysize=6, name=filename_res)
Philipp Arras's avatar
Philipp Arras committed
167
    print("Saved results as '{}'.".format(filename_res))
Philipp Arras's avatar
Philipp Arras committed
168
169
170
171


if __name__ == '__main__':
    main()