getting_started_3.py 5.23 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 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])
49 50 51
    harmonic_space = position_space.get_default_codomain()
    ht = ift.HarmonicTransformOperator(harmonic_space, position_space)
    power_space = ift.PowerSpace(harmonic_space)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
52

Philipp Arras's avatar
Philipp Arras committed
53
    # Set up an amplitude operator for the field
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
    dct = {
        'target': power_space,
        'n_pix': 64,  # 64 spectral bins

        # Spectral smoothness (affects Gaussian process part)
        'a': 3,  # relatively high variance of spectral curbvature
        'k0': .4,  # quefrency mode below which cepstrum flattens

        # Power-law part of spectrum:
        'sm': -5,  # preferred power-law slope
        'sv': .5,  # low variance of power-law slope
        'im': .4,  # y-intercept mean
        'iv': .3  # relatively high y-intercept variance
    }
    A = ift.AmplitudeOperator(**dct)
Philipp Arras's avatar
Philipp Arras committed
69

Philipp Arras's avatar
Philipp Arras committed
70
    # Build the operator for a correlated signal
Jakob Knollmueller's avatar
Jakob Knollmueller committed
71
    power_distributor = ift.PowerDistributor(harmonic_space, power_space)
72 73 74
    vol = harmonic_space.scalar_dvol**-0.5
    xi = ift.ducktape(harmonic_space, None, 'xi')
    correlated_field = ht(vol*power_distributor(A)*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
Jakob Knollmueller's avatar
Jakob Knollmueller committed
79
    signal = ift.sigmoid(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
    # Set up likelihood and information Hamiltonian
Philipp Arras's avatar
Philipp Arras committed
102
    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):
Philipp Arras's avatar
Philipp Arras committed
119
        # Draw new samples and minimize KL
Martin Reinecke's avatar
Martin Reinecke committed
120
        KL = ift.KL_Energy(position, H, N_samples)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
121 122
        KL, convergence = minimizer(KL)
        position = KL.position
Philipp Arras's avatar
Philipp Arras committed
123 124

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

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

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

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