getting_started_3.py 5.23 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
# 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/>.
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
16 17 18
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.

############################################################
Philipp Arras's avatar
Philipp Arras committed
19 20 21
# Non-linear tomography
# The data is integrated lines of sight
# Random lines (set mode=0), radial lines (mode=1)
22
#############################################################
23

Jakob Knollmueller's avatar
Jakob Knollmueller committed
24 25
import numpy as np

Philipp Arras's avatar
Philipp Arras committed
26 27
import nifty5 as ift

Jakob Knollmueller's avatar
Jakob Knollmueller committed
28

Philipp Arras's avatar
Philipp Arras committed
29
def random_los(n_los):
30
    starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
Philipp Arras's avatar
Philipp Arras committed
31 32 33 34 35 36 37
    ends = list(0.5 + 0*np.random.uniform(0, 1, (n_los, 2)).T)
    return starts, ends


def radial_los(n_los):
    starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
    ends = list(np.random.uniform(0, 1, (n_los, 2)).T)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
38 39
    return starts, ends

Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
40

Jakob Knollmueller's avatar
Jakob Knollmueller committed
41
if __name__ == '__main__':
42
    np.random.seed(420)
Philipp Arras's avatar
Philipp Arras committed
43 44 45 46 47

    # Choose between random line-of-sight response (mode=1) and radial lines
    # of sight (mode=2)
    mode = 1

48
    position_space = ift.RGSpace([128, 128])
Jakob Knollmueller's avatar
Jakob Knollmueller committed
49

Philipp Arras's avatar
Philipp Arras committed
50 51 52
    # Set up an amplitude model for the field
    # The parameters mean:
    # 64 spectral bins
53 54
    #
    # Spectral smoothness (affects Gaussian process part)
Philipp Arras's avatar
Philipp Arras committed
55
    # 3 = relatively high variance of spectral curbvature
56 57
    # 0.4 = quefrency mode below which cepstrum flattens
    #
Philipp Arras's avatar
Philipp Arras committed
58 59
    # Power-law part of spectrum:
    # -5 = preferred power-law slope
60
    # 0.5 = low variance of power-law slope
Philipp Arras's avatar
Philipp Arras committed
61 62 63 64 65
    # 0.4 = y-intercept mean
    # 0.3 = relatively high y-intercept variance
    A = ift.AmplitudeModel(position_space, 64, 3, 0.4, -5., 0.5, 0.4, 0.3)

    # Build the model for a correlated signal
Jakob Knollmueller's avatar
Jakob Knollmueller committed
66 67
    harmonic_space = position_space.get_default_codomain()
    ht = ift.HarmonicTransformOperator(harmonic_space, position_space)
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
68
    power_space = A.target[0]
Jakob Knollmueller's avatar
Jakob Knollmueller committed
69 70
    power_distributor = ift.PowerDistributor(harmonic_space, power_space)

Philipp Arras's avatar
Philipp Arras committed
71 72
    vol = ift.ScalingOperator(harmonic_space.scalar_dvol**(-0.5),
                              harmonic_space)
Martin Reinecke's avatar
Martin Reinecke committed
73 74
    correlated_field = ht(
        vol(power_distributor(A))*ift.ducktape(harmonic_space, None, 'xi'))
Philipp Arras's avatar
Philipp Arras committed
75 76
    # Alternatively, one can use:
    # correlated_field = ift.CorrelatedField(position_space, A)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
77

Philipp Arras's avatar
Philipp Arras committed
78
    # Apply a nonlinearity
Martin Reinecke's avatar
Martin Reinecke committed
79
    signal = ift.positive_tanh(correlated_field)
Martin Reinecke's avatar
Martin Reinecke committed
80

Philipp Arras's avatar
Philipp Arras committed
81 82
    # Build the line-of-sight response and define signal response
    LOS_starts, LOS_ends = random_los(100) if mode == 1 else radial_los(100)
Philipp Arras's avatar
Philipp Arras committed
83
    R = ift.LOSResponse(position_space, starts=LOS_starts, ends=LOS_ends)
Martin Reinecke's avatar
Martin Reinecke committed
84
    signal_response = R(signal)
Philipp Arras's avatar
Philipp Arras committed
85 86

    # Specify noise
Jakob Knollmueller's avatar
Jakob Knollmueller committed
87
    data_space = R.target
Jakob Knollmueller's avatar
Jakob Knollmueller committed
88
    noise = .001
89
    N = ift.ScalingOperator(noise, data_space)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
90

Philipp Arras's avatar
Philipp Arras committed
91 92 93
    # 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
94

Philipp Arras's avatar
Philipp Arras committed
95
    # Minimization parameters
Jakob Knollmueller's avatar
Jakob Knollmueller committed
96
    ic_sampling = ift.GradientNormController(iteration_limit=100)
Martin Reinecke's avatar
Martin Reinecke committed
97
    ic_newton = ift.GradInfNormController(
98
        name='Newton', tol=1e-7, iteration_limit=35)
99
    minimizer = ift.NewtonCG(ic_newton)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
100

Philipp Arras's avatar
Philipp Arras committed
101 102
    # Set up model likelihood and information Hamiltonian
    likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response)
Martin Reinecke's avatar
Martin Reinecke committed
103
    H = ift.Hamiltonian(likelihood, ic_sampling)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
104

Philipp Arras's avatar
Philipp Arras committed
105 106
    initial_position = ift.MultiField.full(H.domain, 0.)
    position = initial_position
Jakob Knollmueller's avatar
Jakob Knollmueller committed
107

108
    plot = ift.Plot()
Philipp Arras's avatar
Philipp Arras committed
109
    plot.add(signal(mock_position), title='Ground Truth')
Martin Reinecke's avatar
merge  
Martin Reinecke committed
110
    plot.add(R.adjoint_times(data), title='Data')
Philipp Arras's avatar
Philipp Arras committed
111
    plot.add([A.force(mock_position)], title='Power Spectrum')
Martin Reinecke's avatar
merge  
Martin Reinecke committed
112
    plot.output(ny=1, nx=3, xsize=24, ysize=6, name="setup.png")
Jakob Knollmueller's avatar
Jakob Knollmueller committed
113

Jakob Knollmueller's avatar
Jakob Knollmueller committed
114
    # number of samples used to estimate the KL
115
    N_samples = 20
Philipp Arras's avatar
Philipp Arras committed
116 117

    # Draw new samples to approximate the KL five times
118
    for i in range(5):
Martin Reinecke's avatar
Martin Reinecke committed
119
        KL = ift.KL_Energy(position, H, N_samples)
Philipp Arras's avatar
Philipp Arras committed
120
        # Minimize KL
Jakob Knollmueller's avatar
Jakob Knollmueller committed
121
        KL, convergence = minimizer(KL)
Philipp Arras's avatar
Philipp Arras committed
122
        # Update position
Jakob Knollmueller's avatar
Jakob Knollmueller committed
123
        position = KL.position
Philipp Arras's avatar
Philipp Arras committed
124 125

        # Plot current reconstruction
126
        plot = ift.Plot()
Martin Reinecke's avatar
merge  
Martin Reinecke committed
127
        plot.add(signal(KL.position), title="reconstruction")
Philipp Arras's avatar
Philipp Arras committed
128
        plot.add([A.force(KL.position), A.force(mock_position)], title="power")
Lukas Platz's avatar
Lukas Platz committed
129
        plot.output(ny=1, ysize=6, xsize=16, name="loop-{:02}.png".format(i))
Jakob Knollmueller's avatar
Jakob Knollmueller committed
130

Philipp Arras's avatar
Philipp Arras committed
131
    # Draw posterior samples
Philipp Arras's avatar
Philipp Arras committed
132
    KL = ift.KL_Energy(position, H, N_samples)
Martin Reinecke's avatar
Martin Reinecke committed
133
    sc = ift.StatCalculator()
134
    for sample in KL.samples:
Philipp Arras's avatar
Philipp Arras committed
135
        sc.add(signal(sample + KL.position))
Philipp Arras's avatar
Philipp Arras committed
136 137 138

    # Plotting
    plot = ift.Plot()
Martin Reinecke's avatar
merge  
Martin Reinecke committed
139 140
    plot.add(sc.mean, title="Posterior Mean")
    plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation")
Martin Reinecke's avatar
Martin Reinecke committed
141

Philipp Arras's avatar
Philipp Arras committed
142
    powers = [A.force(s + KL.position) for s in KL.samples]
Martin Reinecke's avatar
merge  
Martin Reinecke committed
143
    plot.add(
Philipp Arras's avatar
Philipp Arras committed
144 145
        powers + [A.force(KL.position),
                  A.force(mock_position)],
Lukas Platz's avatar
Lukas Platz committed
146 147
        title="Sampled Posterior Power Spectrum",
        linewidth=[1.]*len(powers) + [3., 3.])
Martin Reinecke's avatar
merge  
Martin Reinecke committed
148
    plot.output(ny=1, nx=3, xsize=24, ysize=6, name="results.png")