MultiCorrelatedField-getting_startet_2.py 7.79 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
# 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-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.

###############################################################################
# Log-normal field reconstruction from Poissonian data with inhomogenous
# exposure (in case for 2D mode)
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
###############################################################################

import sys
import numpy as np
import nifty5 as ift


def exposure_2d():
    # Structured exposure for 2D mode
31
    x_shape, y_shape = sky_target.shape
32

33
    exposure = np.ones(sky_target.shape)
34
35
36
37
38
39
40
    exposure[x_shape//3:x_shape//2, :] *= 2.
    exposure[x_shape*4//5:x_shape, :] *= .1
    exposure[x_shape//2:x_shape*3//2, :] *= 3.
    exposure[:, x_shape//3:x_shape//2] *= 2.
    exposure[:, x_shape*4//5:x_shape] *= .1
    exposure[:, x_shape//2:x_shape*3//2] *= 3.

41
    return ift.Field.from_global_data(sky_target, exposure)
42
43


44
45
46
47
def cd(input):
    return input.cast_domain(plot_domain)


48
49
50
if __name__ == '__main__':
    np.random.seed(42)

Lukas Platz's avatar
Lukas Platz committed
51
52
53
54
55
56
57
    if len(sys.argv) == 2:
        mode = int(sys.argv[1])
    else:
        mode = 1
    if not mode in (0, 1):
        raise ValueError("Only modes 0 and 1 defined")

58
59
    filename = f"test-MultiCorrelatedField-{mode}" + "-{}.png"

60
    # Two-dimensional signal on RGs with inhomogeneous exposure
61
    N = 128
62
63
    position_space = ift.RGSpace(N)
    energy_space = ift.RGSpace(N)
64
    sky_target = ift.DomainTuple.make((position_space, energy_space))
65
    plot_domain = ift.RGSpace((N, N))
66
67

    exposure = exposure_2d()
Lukas Platz's avatar
Lukas Platz committed
68
69
70
    if mode == 0:
        # need a much higher exposure be able to reconstruct the spectrum
        exposure = exposure * 100
71
72
73
74
75
76
77
78

    # Define harmonic space and harmonic transform
    harmonic_space_p = position_space.get_default_codomain()
    harmonic_space_e = energy_space.get_default_codomain()

    power_space_p = ift.PowerSpace(harmonic_space_p)
    power_space_e = ift.PowerSpace(harmonic_space_e)

79
80
81
82
83
84
85
86
87
88
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
116
117
118
119
120
    amp_p_dict = {
        'target': power_space_p,
        'n_pix': 64,  # 64 spectral bins
        'keys': ['tau_p', 'phi_p', 'zm_p'],

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

        # Power-law part of spectrum
        'sm': -4.,  # Expected exponent of power law
        'sv': 1.,  # Prior standard deviation of -||-
        'im': 0.,  # Expected y-intercept of power law
        'iv': 0.,  # Prior standard deviation of -||-

        # zero-mode
        'za': 5.,  # IG alpha parameter
        'zq': 0.001  # IG q parameter
    }

    amp_e_dict = {
        'target': power_space_e,
        'n_pix': 64,  # 64 spectral bins
        'keys': ['tau_e', 'phi_e', 'zm_e'],

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

        # Power-law part of spectrum
        'sm': -3.,  # Expected exponent of power law
        'sv': 1.,  # Prior standard deviation of -||-
        'im': 0.,  # Expected y-intercept of power law
        'iv': 0.,  # Prior standard deviation of -||-

        # zero-mode
        'za': 6.,  # IG alpha parameter
        'zq': 0.002  # IG q parameter
    }

    amp_p = ift.SLAmplitude(**amp_p_dict)
    amp_e = ift.SLAmplitude(**amp_e_dict)
121
122

    # Define sky operator
123
    # za=4, zq=5 leads to an amplitude median of one
Lukas Platz's avatar
Lukas Platz committed
124
125
126
127
    if mode == 1:
        amplitudes = [amp_p, amp_e]
    elif mode == 0:
        amplitudes = [None, amp_e]
128

129
130
131
    diff_zm_params = {'gs_a': 11., 'gs_q': 50.}
    log_diffuse = ift.MultiCorrelatedField(sky_target, amplitudes,
                                           **diff_zm_params)
Lukas Platz's avatar
Lukas Platz committed
132
    sky = log_diffuse.exp()
133
134
135
136
137
138
139
140

    M = ift.DiagonalOperator(exposure)
    GR = ift.GeometryRemover(M.target)
    # Define instrumental response
    R = GR(M)

    # Generate mock data and define likelihood operator
    lamb = R(sky)
141
    mock_position = ift.from_random('normal', sky.domain)
142
    data = lamb(mock_position)
143
144
    signal = sky(mock_position)

145
    # DEBUG
146
    mock_log_diffuse_vals = log_diffuse(mock_position).to_global_data()
Lukas Platz's avatar
Lukas Platz committed
147
    print("std(mock_log_diffuse):", mock_log_diffuse_vals.std())
148

149
    data = np.random.poisson(data.to_global_data().astype(np.float64))
150
    data = ift.Field.from_global_data(R.target, data)
151
152
153
    likelihood = ift.PoissonianEnergy(data)(lamb)

    # Settings for minimization
154
155
156
157
158
159
160
161
    ic_sampling = ift.AbsDeltaEnergyController(name='Sampling',
                                               iteration_limit=200,
                                               deltaE=0.5)
    H = ift.StandardHamiltonian(likelihood, ic_sampling)

    ic_newton = ift.AbsDeltaEnergyController(name='Newton',
                                             iteration_limit=30,
                                             deltaE=0.5)
162
    minimizer = ift.NewtonCG(ic_newton)
163
164
165
166
    pos = ift.from_random('normal', sky.domain) * 0.01
    iteration = 0

    # Minimize
167
    for i in range(10):
168
169
170
171
172
173
174
175
176
177
        KL = ift.MetricGaussianKL(pos, H, 10)
        KL, convergence = minimizer(KL)
        pos = KL.position
        iteration += 1

        # Plotting
        sc = ift.StatCalculator()
        for sample in KL.samples:
            sc.add(sky(pos + sample))
        reconst = sc.mean
178

Lukas Platz's avatar
Lukas Platz committed
179
        plot = ift.Plot()
180
181
182
183
        plot.add(cd(signal), title='Signal')
        plot.add(cd(GR.adjoint(data)), title='Data')
        plot.add(cd(reconst), title='Posterior Mean')
        plot.add(cd(reconst - signal), title='Error')
184
185
186
        plot.output(xsize=12,
                    ysize=10,
                    name=filename.format(f'{iteration:02d}-reconstr'))
187
188
189
190
191
192
193
194
195
196
197
198
199
200

        p_p = ift.PowerDistributor(harmonic_space_p, power_space_p) @ amp_p
        p_e = ift.PowerDistributor(harmonic_space_e, power_space_e) @ amp_e

        if mode == 1:
            amplitude_specs = [(p_p, 'position'), (p_e, 'energy')]
        elif mode == 0:
            amplitude_specs = [(p_e, 'energy')]

        for p, n in amplitude_specs:
            sc = ift.StatCalculator()
            for sample in KL.samples:
                sc.add(ift.power_analyze(p.force(pos + sample)))
            plot = ift.Plot()
201
202
203
204
205
206
207
208
            plot.add(
                [sc.mean, ift.power_analyze(p.force(mock_position))],
                linewidth=[1., 3.],
                label=['Posterior Mean', 'Ground Truth'],
                title=f"Posterior {n} spectrum")
            plot.output(xsize=12,
                        ysize=10,
                        name=filename.format(f'{iteration:02d}-power_{n}'))
209
210

        print("Saved results as '{}'.".format(filename.format('*')))
211
212
213
214
215
216
217
218
219
220
221
222
223

    # Calculate InverseGamma prior misfits
    keys = ['zm_e', 'zm_p', 'xi_pspec_scaling']
    hyperparams = [(amp_e_dict['za'], amp_e_dict['zq']),
                   (amp_p_dict['za'], amp_p_dict['zq']),
                   (diff_zm_params['gs_a'], diff_zm_params['gs_q'])]

    print("InverseGamma missfits:")
    for key, hpar in zip(keys, hyperparams):
        op = ift.InverseGammaOperator(sky.domain[key], *hpar)
        print(f"{key}:")
        print(f"\tmock: {float(op(mock_position[key]).local_data):+2.3f}")
        print(f"\t rec: {float(op(pos[key]).local_data):+2.3f}\n")