getting_started_2.py 3.53 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 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
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.

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

Martin Reinecke's avatar
Martin Reinecke committed
22

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

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

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

Natalia Porqueres's avatar
Natalia Porqueres committed
37

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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
42
    # Set up the position space of the signal
43
    #
Jakob Knollmueller's avatar
Jakob Knollmueller committed
44
45
46
    # # One dimensional regular grid with uniform exposure
    # position_space = ift.RGSpace(1024)
    # exposure = np.ones(position_space.shape)
Natalia Porqueres's avatar
Natalia Porqueres committed
47

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

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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
56
57
58
    # 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
59

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

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

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

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

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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
80
81
    # Generate mock data
    d_space = R.target[0]
Martin Reinecke's avatar
Martin Reinecke committed
82
    lamb = R(sky)
Martin Reinecke's avatar
Martin Reinecke committed
83
84
    mock_position = ift.from_random('normal', domain)
    data = lamb(mock_position)
85
86
    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
87

Jakob Knollmueller's avatar
Jakob Knollmueller committed
88
    # Compute likelihood and Hamiltonian
Martin Reinecke's avatar
Martin Reinecke committed
89
    likelihood = ift.PoissonianEnergy(lamb, data)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
90
    ic_cg = ift.GradientNormController(iteration_limit=50)
91
92
    ic_newton = ift.GradientNormController(name='Newton', iteration_limit=50,
                                           tol_abs_gradnorm=1e-3)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
93
    minimizer = ift.RelaxedNewton(ic_newton)
Natalia Porqueres's avatar
Natalia Porqueres committed
94

Jakob Knollmueller's avatar
Jakob Knollmueller committed
95
    # Minimize the Hamiltonian
Martin Reinecke's avatar
Martin Reinecke committed
96
    H = ift.Hamiltonian(likelihood)
Martin Reinecke's avatar
Martin Reinecke committed
97
    H = ift.EnergyAdapter(position, H)
98
    H = H.make_invertible(ic_cg)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
99
    H, convergence = minimizer(H)
Natalia Porqueres's avatar
Natalia Porqueres committed
100

Jakob Knollmueller's avatar
Jakob Knollmueller committed
101
    # Plot results
Martin Reinecke's avatar
merge  
Martin Reinecke committed
102
    ift.plot(sky(H.position), title='Reconstruction')
103
    ift.plot(GR.adjoint(data), title='Data')
Martin Reinecke's avatar
typo  
Martin Reinecke committed
104
    ift.plot_finish(name='getting_started_2.png', xsize=16, ysize=16)