getting_started_2.py 4.07 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
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.

Philipp Arras's avatar
Philipp Arras committed
18
19
20
###############################################################################
# Log-normal field reconstruction from Poissonian data with inhomogenous
# exposure (in case for 2D mode)
21
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
Philipp Arras's avatar
Philipp Arras committed
22
###############################################################################
23

Natalia Porqueres's avatar
Natalia Porqueres committed
24
25
import numpy as np

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

Martin Reinecke's avatar
Martin Reinecke committed
28

Philipp Arras's avatar
Philipp Arras committed
29
30
def exposure_2d():
    # Structured exposure for 2D mode
Jakob Knollmueller's avatar
Jakob Knollmueller committed
31
    x_shape, y_shape = position_space.shape
Natalia Porqueres's avatar
Natalia Porqueres committed
32

Jakob Knollmueller's avatar
Jakob Knollmueller committed
33
    exposure = np.ones(position_space.shape)
Martin Reinecke's avatar
Martin Reinecke committed
34
35
36
37
38
39
    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.
40

Philipp Arras's avatar
Philipp Arras committed
41
    return ift.Field.from_global_data(position_space, exposure)
Natalia Porqueres's avatar
Natalia Porqueres committed
42

Natalia Porqueres's avatar
Natalia Porqueres committed
43

Jakob Knollmueller's avatar
Jakob Knollmueller committed
44
if __name__ == '__main__':
Philipp Arras's avatar
Philipp Arras committed
45
    # FIXME All random seeds to 42
46
    np.random.seed(41)
Natalia Porqueres's avatar
Natalia Porqueres committed
47

Philipp Arras's avatar
Philipp Arras committed
48
    # Choose space on which the signal field is defined
49
50
    mode = 2
    if mode == 0:
Philipp Arras's avatar
Philipp Arras committed
51
        # One-dimensional regular grid with uniform exposure
52
53
54
55
56
        position_space = ift.RGSpace(1024)
        exposure = ift.Field.full(position_space, 1.)
    elif mode == 1:
        # Two-dimensional regular grid with inhomogeneous exposure
        position_space = ift.RGSpace([512, 512])
Philipp Arras's avatar
Philipp Arras committed
57
        exposure = exposure_2d()
58
    else:
Philipp Arras's avatar
Philipp Arras committed
59
        # Sphere with uniform exposure
60
61
        position_space = ift.HPSpace(128)
        exposure = ift.Field.full(position_space, 1.)
Philipp Arras's avatar
Philipp Arras committed
62
63

    # Define harmonic space and harmonic transform
Jakob Knollmueller's avatar
Jakob Knollmueller committed
64
65
    harmonic_space = position_space.get_default_codomain()
    HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
Natalia Porqueres's avatar
Natalia Porqueres committed
66

Philipp Arras's avatar
Philipp Arras committed
67
    # Domain on which the field's degrees of freedom are defined
Martin Reinecke's avatar
Martin Reinecke committed
68
    domain = ift.DomainTuple.make(harmonic_space)
Natalia Porqueres's avatar
Natalia Porqueres committed
69

Philipp Arras's avatar
Philipp Arras committed
70
    # Define amplitude (square root of power spectrum)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
71
    def sqrtpspec(k):
Philipp Arras's avatar
Philipp Arras committed
72
        return 1./(20. + k**2)
Natalia Porqueres's avatar
Natalia Porqueres committed
73

Jakob Knollmueller's avatar
Jakob Knollmueller committed
74
75
76
77
    p_space = ift.PowerSpace(harmonic_space)
    pd = ift.PowerDistributor(harmonic_space, p_space)
    a = ift.PS_field(p_space, sqrtpspec)
    A = pd(a)
Natalia Porqueres's avatar
Natalia Porqueres committed
78

Philipp Arras's avatar
Philipp Arras committed
79
    # Define sky operator
Martin Reinecke's avatar
Martin Reinecke committed
80
    sky = ift.exp(HT(ift.makeOp(A)))
Natalia Porqueres's avatar
Natalia Porqueres committed
81

Jakob Knollmueller's avatar
Jakob Knollmueller committed
82
83
    M = ift.DiagonalOperator(exposure)
    GR = ift.GeometryRemover(position_space)
Philipp Arras's avatar
Philipp Arras committed
84
    # Define instrumental response
Martin Reinecke's avatar
Martin Reinecke committed
85
    R = GR(M)
Natalia Porqueres's avatar
Natalia Porqueres committed
86

Philipp Arras's avatar
Philipp Arras committed
87
    # Generate mock data and define likelihood operator
Jakob Knollmueller's avatar
Jakob Knollmueller committed
88
    d_space = R.target[0]
Martin Reinecke's avatar
Martin Reinecke committed
89
    lamb = R(sky)
Martin Reinecke's avatar
Martin Reinecke committed
90
91
    mock_position = ift.from_random('normal', domain)
    data = lamb(mock_position)
92
93
    data = np.random.poisson(data.to_global_data().astype(np.float64))
    data = ift.Field.from_global_data(d_space, data)
94
    likelihood = ift.PoissonianEnergy(data)(lamb)
Philipp Arras's avatar
Philipp Arras committed
95
96
97
98

    # Settings for minimization
    ic_newton = ift.DeltaEnergyController(
        name='Newton', iteration_limit=100, tol_rel_deltaE=1e-8)
99
    minimizer = ift.NewtonCG(ic_newton)
Natalia Porqueres's avatar
Natalia Porqueres committed
100

Philipp Arras's avatar
Philipp Arras committed
101
    # Compute MAP solution by minimizing the information Hamiltonian
102
    H = ift.StandardHamiltonian(likelihood)
Philipp Arras's avatar
Philipp Arras committed
103
104
    initial_position = ift.from_random('normal', domain)
    H = ift.EnergyAdapter(initial_position, H, want_metric=True)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
105
    H, convergence = minimizer(H)
Natalia Porqueres's avatar
Natalia Porqueres committed
106

Philipp Arras's avatar
Philipp Arras committed
107
    # Plotting
Lukas Platz's avatar
Lukas Platz committed
108
109
    signal = sky(mock_position)
    reconst = sky(H.position)
110
111
112
113
114
    plot = ift.Plot()
    plot.add(signal, title='Signal')
    plot.add(GR.adjoint(data), title='Data')
    plot.add(reconst, title='Reconstruction')
    plot.add(reconst - signal, title='Residuals')
115
    plot.output(name='getting_started_2.pdf', xsize=16, ysize=16)