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

Natalia Porqueres's avatar
Natalia Porqueres committed
18
19
20
import nifty5 as ift
import numpy as np

Martin Reinecke's avatar
Martin Reinecke committed
21

Jakob Knollmueller's avatar
Jakob Knollmueller committed
22
23
def get_2D_exposure():
    x_shape, y_shape = position_space.shape
Natalia Porqueres's avatar
Natalia Porqueres committed
24

Jakob Knollmueller's avatar
Jakob Knollmueller committed
25
    exposure = np.ones(position_space.shape)
Martin Reinecke's avatar
Martin Reinecke committed
26
27
28
29
30
31
    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.
32
33

    exposure = ift.Field.from_global_data(position_space, exposure)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
34
    return exposure
Natalia Porqueres's avatar
Natalia Porqueres committed
35

Natalia Porqueres's avatar
Natalia Porqueres committed
36

Jakob Knollmueller's avatar
Jakob Knollmueller committed
37
if __name__ == '__main__':
Philipp Arras's avatar
Philipp Arras committed
38
    # FIXME description of the tutorial
39
    np.random.seed(41)
Natalia Porqueres's avatar
Natalia Porqueres committed
40

Jakob Knollmueller's avatar
Jakob Knollmueller committed
41
    # Set up the position space of the signal
42
    #
Jakob Knollmueller's avatar
Jakob Knollmueller committed
43
44
    # # One dimensional regular grid with uniform exposure
    # position_space = ift.RGSpace(1024)
Lukas Platz's avatar
Lukas Platz committed
45
    # exposure = ift.Field.full(position_space, 1.)
Natalia Porqueres's avatar
Natalia Porqueres committed
46

47
    # Two-dimensional regular grid with inhomogeneous exposure
Jakob Knollmueller's avatar
Jakob Knollmueller committed
48
49
    position_space = ift.RGSpace([512, 512])
    exposure = get_2D_exposure()
Natalia Porqueres's avatar
Natalia Porqueres committed
50

Martin Reinecke's avatar
Martin Reinecke committed
51
    #  Sphere with uniform exposure
Jakob Knollmueller's avatar
Jakob Knollmueller committed
52
    # position_space = ift.HPSpace(128)
53
    # exposure = ift.Field.full(position_space, 1.)
Natalia Porqueres's avatar
Natalia Porqueres committed
54

Jakob Knollmueller's avatar
Jakob Knollmueller committed
55
56
57
    # Defining harmonic space and transform
    harmonic_space = position_space.get_default_codomain()
    HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
Natalia Porqueres's avatar
Natalia Porqueres committed
58

Martin Reinecke's avatar
Martin Reinecke committed
59
    domain = ift.DomainTuple.make(harmonic_space)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
60
    position = ift.from_random('normal', domain)
Natalia Porqueres's avatar
Natalia Porqueres committed
61

Jakob Knollmueller's avatar
Jakob Knollmueller committed
62
63
64
    # Define power spectrum and amplitudes
    def sqrtpspec(k):
        return 1. / (20. + k**2)
Natalia Porqueres's avatar
Natalia Porqueres committed
65

Jakob Knollmueller's avatar
Jakob Knollmueller committed
66
67
68
69
    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
70

Jakob Knollmueller's avatar
Jakob Knollmueller committed
71
    # Set up a sky model
Martin Reinecke's avatar
Martin Reinecke committed
72
    sky = ift.exp(HT(ift.makeOp(A)))
Natalia Porqueres's avatar
Natalia Porqueres committed
73

Jakob Knollmueller's avatar
Jakob Knollmueller committed
74
75
76
    M = ift.DiagonalOperator(exposure)
    GR = ift.GeometryRemover(position_space)
    # Set up instrumental response
Martin Reinecke's avatar
Martin Reinecke committed
77
    R = GR(M)
Natalia Porqueres's avatar
Natalia Porqueres committed
78

Jakob Knollmueller's avatar
Jakob Knollmueller committed
79
80
    # Generate mock data
    d_space = R.target[0]
Martin Reinecke's avatar
Martin Reinecke committed
81
    lamb = R(sky)
Martin Reinecke's avatar
Martin Reinecke committed
82
83
    mock_position = ift.from_random('normal', domain)
    data = lamb(mock_position)
84
85
    data = np.random.poisson(data.to_global_data().astype(np.float64))
    data = ift.Field.from_global_data(d_space, data)
Natalia Porqueres's avatar
Natalia Porqueres committed
86

Jakob Knollmueller's avatar
Jakob Knollmueller committed
87
    # Compute likelihood and Hamiltonian
88
89
    ic_newton = ift.DeltaEnergyController(name='Newton', iteration_limit=100,
                                          tol_rel_deltaE=1e-8)
90
    likelihood = ift.PoissonianEnergy(data)(lamb)
91
    minimizer = ift.NewtonCG(ic_newton)
Natalia Porqueres's avatar
Natalia Porqueres committed
92

Jakob Knollmueller's avatar
Jakob Knollmueller committed
93
    # Minimize the Hamiltonian
Martin Reinecke's avatar
Martin Reinecke committed
94
    H = ift.Hamiltonian(likelihood)
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
95
    H = ift.EnergyAdapter(position, H, want_metric=True)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
96
    H, convergence = minimizer(H)
Natalia Porqueres's avatar
Natalia Porqueres committed
97

Jakob Knollmueller's avatar
Jakob Knollmueller committed
98
    # Plot results
Lukas Platz's avatar
Lukas Platz committed
99
100
    signal = sky(mock_position)
    reconst = sky(H.position)
101
102
103
104
105
106
    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')
    plot.output(name='getting_started_2.png', xsize=16, ysize=16)