getting_started_2.py 4.24 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
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

class EnergyAdapter(ift.Energy):
    def __init__(self, position, op):
        super(EnergyAdapter, self).__init__(position)
        self._op = op
        pvar = ift.Linearization.make_var(position)
        self._res = op(pvar)

    def at(self, position):
        return EnergyAdapter(position, self._op)

    @property
    def value(self):
        return self._res.val.local_data[()]

    @property
    def gradient(self):
        return self._res.gradient

    @property
    def metric(self):
        return self._res.metric

Natalia Porqueres's avatar
Natalia Porqueres committed
45

Jakob Knollmueller's avatar
Jakob Knollmueller committed
46
47
def get_2D_exposure():
    x_shape, y_shape = position_space.shape
Natalia Porqueres's avatar
Natalia Porqueres committed
48

Jakob Knollmueller's avatar
Jakob Knollmueller committed
49
    exposure = np.ones(position_space.shape)
Martin Reinecke's avatar
Martin Reinecke committed
50
51
52
53
54
55
    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.
56
57

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

Natalia Porqueres's avatar
Natalia Porqueres committed
60

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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
65
    # Set up the position space of the signal
66
    #
Jakob Knollmueller's avatar
Jakob Knollmueller committed
67
68
69
    # # 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
70

71
    # Two-dimensional regular grid with inhomogeneous exposure
Jakob Knollmueller's avatar
Jakob Knollmueller committed
72
73
    position_space = ift.RGSpace([512, 512])
    exposure = get_2D_exposure()
Natalia Porqueres's avatar
Natalia Porqueres committed
74

Jakob Knollmueller's avatar
Jakob Knollmueller committed
75
76
    # # Sphere with with uniform exposure
    # position_space = ift.HPSpace(128)
77
    # exposure = ift.Field.full(position_space, 1.)
Natalia Porqueres's avatar
Natalia Porqueres committed
78

Jakob Knollmueller's avatar
Jakob Knollmueller committed
79
80
81
    # 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
82

Martin Reinecke's avatar
Martin Reinecke committed
83
    domain = ift.DomainTuple.make(harmonic_space)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
84
    position = ift.from_random('normal', domain)
Natalia Porqueres's avatar
Natalia Porqueres committed
85

Jakob Knollmueller's avatar
Jakob Knollmueller committed
86
87
88
    # Define power spectrum and amplitudes
    def sqrtpspec(k):
        return 1. / (20. + k**2)
Natalia Porqueres's avatar
Natalia Porqueres committed
89

Jakob Knollmueller's avatar
Jakob Knollmueller committed
90
91
92
93
    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
94

Jakob Knollmueller's avatar
Jakob Knollmueller committed
95
    # Set up a sky model
Martin Reinecke's avatar
Martin Reinecke committed
96
    sky = lambda inp: (HT(inp*A)).exp()
Natalia Porqueres's avatar
Natalia Porqueres committed
97

Jakob Knollmueller's avatar
Jakob Knollmueller committed
98
99
100
101
    M = ift.DiagonalOperator(exposure)
    GR = ift.GeometryRemover(position_space)
    # Set up instrumental response
    R = GR * M
Natalia Porqueres's avatar
Natalia Porqueres committed
102

Jakob Knollmueller's avatar
Jakob Knollmueller committed
103
104
    # Generate mock data
    d_space = R.target[0]
Martin Reinecke's avatar
Martin Reinecke committed
105
106
    lamb = lambda inp: R(sky(inp))
    mock_position = ift.from_random('normal', domain)
Martin Reinecke's avatar
Martin Reinecke committed
107
108
109
    #ift.extra.check_value_gradient_consistency2(lamb, mock_position)
    #testl = GaussianEnergy2(None, M)
    #ift.extra.check_value_gradient_metric_consistency2(testl, sky(mock_position))
Martin Reinecke's avatar
Martin Reinecke committed
110
    data = lamb(mock_position)
111
112
    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
113

Jakob Knollmueller's avatar
Jakob Knollmueller committed
114
    # Compute likelihood and Hamiltonian
Martin Reinecke's avatar
Martin Reinecke committed
115
    likelihood = ift.PoissonianEnergy(lamb, data)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
116
    ic_cg = ift.GradientNormController(iteration_limit=50)
117
118
    ic_newton = ift.GradientNormController(name='Newton', iteration_limit=50,
                                           tol_abs_gradnorm=1e-3)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
119
    minimizer = ift.RelaxedNewton(ic_newton)
Natalia Porqueres's avatar
Natalia Porqueres committed
120

Jakob Knollmueller's avatar
Jakob Knollmueller committed
121
    # Minimize the Hamiltonian
Martin Reinecke's avatar
Martin Reinecke committed
122
    H = ift.Hamiltonian(likelihood)
Martin Reinecke's avatar
Martin Reinecke committed
123
124
    H = EnergyAdapter(position, H)
    #ift.extra.check_value_gradient_consistency(H)
125
    H = H.make_invertible(ic_cg)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
126
    H, convergence = minimizer(H)
Natalia Porqueres's avatar
Natalia Porqueres committed
127

Jakob Knollmueller's avatar
Jakob Knollmueller committed
128
    # Plot results
Martin Reinecke's avatar
Martin Reinecke committed
129
    result_sky = sky(H.position)
Martin Reinecke's avatar
Martin Reinecke committed
130
131
    ift.plot(result_sky)
    ift.plot_finish()
Martin Reinecke's avatar
Martin Reinecke committed
132
    # FIXME PLOTTING