getting_started_3.py 5.28 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 19 20 21 22 23
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.

############################################################
# Non-linear tomography 
# data is line of sight (LOS) field
# random lines (set mode=0), radial lines (mode=1)
#############################################################
mode = 0
24

Jakob Knollmueller's avatar
Jakob Knollmueller committed
25 26 27 28 29
import nifty5 as ift
import numpy as np


def get_random_LOS(n_los):
30
    # Setting up LOS
31
    starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
32 33 34 35
    if mode == 0:
        ends = list(np.random.uniform(0, 1, (n_los, 2)).T)
    else:
        ends = list(0.5+0*np.random.uniform(0, 1, (n_los, 2)).T)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
36 37
    return starts, ends

Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
38

Jakob Knollmueller's avatar
Jakob Knollmueller committed
39
if __name__ == '__main__':
Philipp Arras's avatar
Philipp Arras committed
40
    # FIXME description of the tutorial
41
    np.random.seed(420)
Philipp Arras's avatar
Philipp Arras committed
42
    np.seterr(all='raise')
43
    position_space = ift.RGSpace([128, 128])
Jakob Knollmueller's avatar
Jakob Knollmueller committed
44

45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
    # Setting up an amplitude model for the field
    A = ift.AmplitudeModel(position_space, 64, 3, 0.4, -5., 0.5, 0.4, 0.3)
    # made choices:
    # 64 spectral bins 
    #
    # Spectral smoothness (affects Gaussian process part)
    # 3 = relatively high variance of spectral curbvature 
    # 0.4 = quefrency mode below which cepstrum flattens
    #
    # power law part of spectrum:
    # -5= preferred power-law slope
    # 0.5 = low variance of power-law slope
    #
    # Gaussian process part of log-spectrum
    # 0.4 = y-intercept mean of additional power 
    # 0.3 = y-intercept variance of additional power
    
Jakob Knollmueller's avatar
Jakob Knollmueller committed
62 63 64
    # Building the model for a correlated signal
    harmonic_space = position_space.get_default_codomain()
    ht = ift.HarmonicTransformOperator(harmonic_space, position_space)
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
65
    power_space = A.target[0]
Jakob Knollmueller's avatar
Jakob Knollmueller committed
66 67
    power_distributor = ift.PowerDistributor(harmonic_space, power_space)

68
    vol = harmonic_space.scalar_dvol
Martin Reinecke's avatar
Martin Reinecke committed
69
    vol = ift.ScalingOperator(vol**(-0.5), harmonic_space)
Martin Reinecke's avatar
Martin Reinecke committed
70 71
    correlated_field = ht(
        vol(power_distributor(A))*ift.ducktape(harmonic_space, None, 'xi'))
72
    # alternatively to the block above one can do:
73
    #correlated_field = ift.CorrelatedField(position_space, A)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
74 75

    # apply some nonlinearity
Martin Reinecke's avatar
Martin Reinecke committed
76
    signal = ift.positive_tanh(correlated_field)
Martin Reinecke's avatar
Martin Reinecke committed
77

Jakob Knollmueller's avatar
Jakob Knollmueller committed
78
    # Building the Line of Sight response
79
    LOS_starts, LOS_ends = get_random_LOS(100)
Philipp Arras's avatar
Philipp Arras committed
80
    R = ift.LOSResponse(position_space, starts=LOS_starts, ends=LOS_ends)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
81
    # build signal response model and model likelihood
Martin Reinecke's avatar
Martin Reinecke committed
82
    signal_response = R(signal)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
83 84
    # specify noise
    data_space = R.target
Jakob Knollmueller's avatar
Jakob Knollmueller committed
85
    noise = .001
86
    N = ift.ScalingOperator(noise, data_space)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
87

88
    # generate mock signal and data
Philipp Arras's avatar
Philipp Arras committed
89
    MOCK_POSITION = ift.from_random('normal', signal_response.domain)
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
90
    data = signal_response(MOCK_POSITION) + N.draw_sample()
Jakob Knollmueller's avatar
Jakob Knollmueller committed
91 92

    # set up model likelihood
Philipp Arras's avatar
Philipp Arras committed
93
    likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
94 95

    # set up minimization and inversion schemes
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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
101
    # build model Hamiltonian
Martin Reinecke's avatar
Martin Reinecke committed
102
    H = ift.Hamiltonian(likelihood, ic_sampling)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
103

104
    INITIAL_POSITION = ift.MultiField.full(H.domain, 0.)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
105
    position = INITIAL_POSITION
Jakob Knollmueller's avatar
Jakob Knollmueller committed
106

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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
113
    # number of samples used to estimate the KL
114
    N_samples = 20
115 116
    
    # five intermediate steps in minimization to illustrate progress
117
    for i in range(5):
118
        # set up KL
Martin Reinecke's avatar
Martin Reinecke committed
119
        KL = ift.KL_Energy(position, H, N_samples)
120
        # minimize KL until iteration limit reached
Jakob Knollmueller's avatar
Jakob Knollmueller committed
121
        KL, convergence = minimizer(KL)
122
        # read out position
Jakob Knollmueller's avatar
Jakob Knollmueller committed
123
        position = KL.position
124
        # plot momentariy field
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

130
    # final plot 
Philipp Arras's avatar
Philipp Arras committed
131
    KL = ift.KL_Energy(position, H, N_samples)
132
    plot = ift.Plot()
133
    # do statistics
Martin Reinecke's avatar
Martin Reinecke committed
134
    sc = ift.StatCalculator()
135
    for sample in KL.samples:
Philipp Arras's avatar
Philipp Arras committed
136
        sc.add(signal(sample + KL.position))
Martin Reinecke's avatar
merge  
Martin Reinecke committed
137 138
    plot.add(sc.mean, title="Posterior Mean")
    plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation")
Martin Reinecke's avatar
Martin Reinecke committed
139

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