getting_started_3.py 3.59 KB
Newer Older
Jakob Knollmueller's avatar
Jakob Knollmueller committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import nifty5 as ift
from nifty5.library.los_response import LOSResponse
from nifty5.library.amplitude_model import make_amplitude_model
from nifty5.library.smooth_sky import make_correlated_field
import numpy as np


def get_random_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)

    return starts, ends


if __name__ == '__main__':
    np.random.seed(41)
    position_space = ift.RGSpace([128,128])

Jakob Knollmueller's avatar
Jakob Knollmueller committed
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
    # Setting up an amplitude model
    A, amplitude_internals = make_amplitude_model(position_space,16, 1, 10, -4., 1, 0., 1.)

    # Building the model for a correlated signal
    harmonic_space = position_space.get_default_codomain()
    ht = ift.HarmonicTransformOperator(harmonic_space, position_space)
    power_space = A.value.domain[0]
    power_distributor = ift.PowerDistributor(harmonic_space, power_space)
    position = {}
    position['xi'] = ift.Field.from_random('normal', harmonic_space)
    position = ift.MultiField(position)

    xi = ift.Variable(position)['xi']
    Amp = power_distributor(A)
    correlated_field_h = Amp * xi
    correlated_field = ht(correlated_field_h)
    # # alternatively to the block above one can do:
    # correlated_field, _ = make_correlated_field(position_space,A)

    # apply some nonlinearity
    signal = ift.PointwisePositiveTanh(correlated_field)
    # Building the Line of Sight response
    LOS_starts, LOS_ends = get_random_LOS(1000)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
42
    R = LOSResponse(position_space, starts=LOS_starts, ends=LOS_ends)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
43
44

    # build signal response model and model likelihood
Jakob Knollmueller's avatar
Jakob Knollmueller committed
45
    signal_response = R(signal)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
46
47
    # specify noise
    data_space = R.target
Jakob Knollmueller's avatar
Jakob Knollmueller committed
48
49
50
    noise = .001
    N = ift.ScalingOperator(noise,data_space)

Jakob Knollmueller's avatar
Jakob Knollmueller committed
51
52
53

    # generate mock data
    MOCK_POSITION = ift.from_random('normal', signal.position.domain)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
54
    data = signal_response.at(MOCK_POSITION).value + N.draw_sample()
Jakob Knollmueller's avatar
Jakob Knollmueller committed
55
56
57
58
59

    # set up model likelihood
    likelihood = ift.GaussianEnergy(signal_response, mean=data, covariance=N)

    # set up minimization and inversion schemes
Jakob Knollmueller's avatar
Jakob Knollmueller committed
60
61
    ic_cg = ift.GradientNormController(iteration_limit=10)
    ic_sampling = ift.GradientNormController(iteration_limit=100)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
62
63
    ic_newton = ift.GradientNormController(name='Newton', iteration_limit=100)
    minimizer = ift.RelaxedNewton(ic_newton)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
64

Jakob Knollmueller's avatar
Jakob Knollmueller committed
65
    # build model Hamiltonian
Jakob Knollmueller's avatar
Jakob Knollmueller committed
66
    H = ift.Hamiltonian(likelihood,ic_cg,iteration_controller_sampling=ic_sampling)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
67

Jakob Knollmueller's avatar
Jakob Knollmueller committed
68
69
    INITIAL_POSITION = ift.from_random('normal',H.position.domain)
    position = INITIAL_POSITION
Jakob Knollmueller's avatar
Jakob Knollmueller committed
70

Jakob Knollmueller's avatar
Jakob Knollmueller committed
71
72
73
74
    ift.plot(signal.at(MOCK_POSITION).value,name='truth.pdf')
    ift.plot(R.adjoint_times(data),name='data.pdf')
    ift.plot([ A.at(MOCK_POSITION).value], name='power.pdf')

Jakob Knollmueller's avatar
Jakob Knollmueller committed
75
76
    # number of samples used to estimate the KL
    N_samples = 10
Jakob Knollmueller's avatar
Jakob Knollmueller committed
77
78
79
80
81
82
83
84
85
86
87
    for i in range(5):
        H = H.at(position)
        samples = [H.curvature.draw_sample(from_inverse=True) for _ in range(N_samples)]

        KL = ift.SampledKullbachLeiblerDivergence(H, samples, ic_cg)
        KL, convergence = minimizer(KL)
        position = KL.position

        ift.plot(signal.at(position).value,name='reconstruction.pdf')

        ift.plot([A.at(position).value, A.at(MOCK_POSITION).value],name='power.pdf')
Jakob Knollmueller's avatar
Jakob Knollmueller committed
88

Jakob Knollmueller's avatar
Jakob Knollmueller committed
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
    avrg = 0.
    va = 0.
    powers = []
    for sample in samples:
        sam = signal.at(sample + position).value
        powers.append(A.at(sample+position).value)
        avrg += sam
        va += sam**2

    avrg /= len(samples)
    va /= len(samples)
    va -= avrg**2
    std = ift.sqrt(va)
    ift.plot(avrg, name='avrg.pdf')
    ift.plot(std, name='std.pdf')
    ift.plot([A.at(position).value, A.at(MOCK_POSITION).value]+powers, name='power.pdf')