getting_started_2.py 3.55 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 22
import nifty5 as ift
import numpy as np


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 39
if __name__ == '__main__':
    # ABOUT THIS CODE
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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
52 53
    # # Sphere with with uniform exposure
    # 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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
60 61
    domain = ift.MultiDomain.make({'xi': harmonic_space})
    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 73 74 75 76
    # Set up a sky model
    xi = ift.Variable(position)['xi']
    logsky_h = xi * A
    logsky = HT(logsky_h)
    sky = ift.PointwiseExponential(logsky)
Natalia Porqueres's avatar
Natalia Porqueres committed
77

Jakob Knollmueller's avatar
Jakob Knollmueller committed
78 79 80 81
    M = ift.DiagonalOperator(exposure)
    GR = ift.GeometryRemover(position_space)
    # Set up instrumental response
    R = GR * M
Natalia Porqueres's avatar
Natalia Porqueres committed
82

Jakob Knollmueller's avatar
Jakob Knollmueller committed
83 84 85 86
    # Generate mock data
    d_space = R.target[0]
    lamb = R(sky)
    mock_position = ift.from_random('normal', lamb.position.domain)
87 88 89
    data = lamb.at(mock_position).value
    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
90

Jakob Knollmueller's avatar
Jakob Knollmueller committed
91 92 93 94
    # Compute likelihood and Hamiltonian
    position = ift.from_random('normal', lamb.position.domain)
    likelihood = ift.PoissonianEnergy(lamb, data)
    ic_cg = ift.GradientNormController(iteration_limit=50)
95 96
    ic_newton = ift.GradientNormController(name='Newton', iteration_limit=50,
                                           tol_abs_gradnorm=1e-3)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
97
    minimizer = ift.RelaxedNewton(ic_newton)
Natalia Porqueres's avatar
Natalia Porqueres committed
98

Jakob Knollmueller's avatar
Jakob Knollmueller committed
99
    # Minimize the Hamiltonian
Martin Reinecke's avatar
Martin Reinecke committed
100 101
    H = ift.Hamiltonian(likelihood)
    H = H.makeInvertible(ic_cg)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
102
    H, convergence = minimizer(H)
Natalia Porqueres's avatar
Natalia Porqueres committed
103

Jakob Knollmueller's avatar
Jakob Knollmueller committed
104 105
    # Plot results
    result_sky = sky.at(H.position).value
106
    # PLOTTING