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