getting_started_2.py 4.06 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
# 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
#
16 17
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.

Philipp Arras's avatar
Philipp Arras committed
18 19 20
###############################################################################
# Log-normal field reconstruction from Poissonian data with inhomogenous
# exposure (in case for 2D mode)
21
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
Philipp Arras's avatar
Philipp Arras committed
22
###############################################################################
23

Natalia Porqueres's avatar
Natalia Porqueres committed
24 25 26
import nifty5 as ift
import numpy as np

Martin Reinecke's avatar
Martin Reinecke committed
27

Philipp Arras's avatar
Philipp Arras committed
28 29
def exposure_2d():
    # Structured exposure for 2D mode
Jakob Knollmueller's avatar
Jakob Knollmueller committed
30
    x_shape, y_shape = position_space.shape
Natalia Porqueres's avatar
Natalia Porqueres committed
31

Jakob Knollmueller's avatar
Jakob Knollmueller committed
32
    exposure = np.ones(position_space.shape)
Martin Reinecke's avatar
Martin Reinecke committed
33 34 35 36 37 38
    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.
39

Philipp Arras's avatar
Philipp Arras committed
40
    return ift.Field.from_global_data(position_space, exposure)
Natalia Porqueres's avatar
Natalia Porqueres committed
41

Natalia Porqueres's avatar
Natalia Porqueres committed
42

Jakob Knollmueller's avatar
Jakob Knollmueller committed
43
if __name__ == '__main__':
Philipp Arras's avatar
Philipp Arras committed
44
    # FIXME All random seeds to 42
45
    np.random.seed(41)
Natalia Porqueres's avatar
Natalia Porqueres committed
46

Philipp Arras's avatar
Philipp Arras committed
47
    # Choose space on which the signal field is defined
48 49
    mode = 2
    if mode == 0:
Philipp Arras's avatar
Philipp Arras committed
50
        # One-dimensional regular grid with uniform exposure
51 52 53 54 55
        position_space = ift.RGSpace(1024)
        exposure = ift.Field.full(position_space, 1.)
    elif mode == 1:
        # Two-dimensional regular grid with inhomogeneous exposure
        position_space = ift.RGSpace([512, 512])
Philipp Arras's avatar
Philipp Arras committed
56
        exposure = exposure_2d()
57
    else:
Philipp Arras's avatar
Philipp Arras committed
58
        # Sphere with uniform exposure
59 60
        position_space = ift.HPSpace(128)
        exposure = ift.Field.full(position_space, 1.)
Philipp Arras's avatar
Philipp Arras committed
61 62

    # Define harmonic space and harmonic transform
Jakob Knollmueller's avatar
Jakob Knollmueller committed
63 64
    harmonic_space = position_space.get_default_codomain()
    HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
Natalia Porqueres's avatar
Natalia Porqueres committed
65

Philipp Arras's avatar
Philipp Arras committed
66
    # Domain on which the field's degrees of freedom are defined
Martin Reinecke's avatar
Martin Reinecke committed
67
    domain = ift.DomainTuple.make(harmonic_space)
Natalia Porqueres's avatar
Natalia Porqueres committed
68

Philipp Arras's avatar
Philipp Arras committed
69
    # Define amplitude (square root of power spectrum)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
70
    def sqrtpspec(k):
Philipp Arras's avatar
Philipp Arras committed
71
        return 1./(20. + k**2)
Natalia Porqueres's avatar
Natalia Porqueres committed
72

Jakob Knollmueller's avatar
Jakob Knollmueller committed
73 74 75 76
    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
77

Philipp Arras's avatar
Philipp Arras committed
78
    # Define sky model
Martin Reinecke's avatar
Martin Reinecke committed
79
    sky = ift.exp(HT(ift.makeOp(A)))
Natalia Porqueres's avatar
Natalia Porqueres committed
80

Jakob Knollmueller's avatar
Jakob Knollmueller committed
81 82
    M = ift.DiagonalOperator(exposure)
    GR = ift.GeometryRemover(position_space)
Philipp Arras's avatar
Philipp Arras committed
83
    # Define instrumental response
Martin Reinecke's avatar
Martin Reinecke committed
84
    R = GR(M)
Natalia Porqueres's avatar
Natalia Porqueres committed
85

Philipp Arras's avatar
Philipp Arras committed
86
    # Generate mock data and define likelihood operator
Jakob Knollmueller's avatar
Jakob Knollmueller committed
87
    d_space = R.target[0]
Martin Reinecke's avatar
Martin Reinecke committed
88
    lamb = R(sky)
Martin Reinecke's avatar
Martin Reinecke committed
89 90
    mock_position = ift.from_random('normal', domain)
    data = lamb(mock_position)
91 92
    data = np.random.poisson(data.to_global_data().astype(np.float64))
    data = ift.Field.from_global_data(d_space, data)
93
    likelihood = ift.PoissonianEnergy(data)(lamb)
Philipp Arras's avatar
Philipp Arras committed
94 95 96 97

    # Settings for minimization
    ic_newton = ift.DeltaEnergyController(
        name='Newton', iteration_limit=100, tol_rel_deltaE=1e-8)
98
    minimizer = ift.NewtonCG(ic_newton)
Natalia Porqueres's avatar
Natalia Porqueres committed
99

Philipp Arras's avatar
Philipp Arras committed
100
    # Compute MAP solution by minimizing the information Hamiltonian
Martin Reinecke's avatar
Martin Reinecke committed
101
    H = ift.Hamiltonian(likelihood)
Philipp Arras's avatar
Philipp Arras committed
102 103
    initial_position = ift.from_random('normal', domain)
    H = ift.EnergyAdapter(initial_position, H, want_metric=True)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
104
    H, convergence = minimizer(H)
Natalia Porqueres's avatar
Natalia Porqueres committed
105

Philipp Arras's avatar
Philipp Arras committed
106
    # Plotting
Lukas Platz's avatar
Lukas Platz committed
107 108
    signal = sky(mock_position)
    reconst = sky(H.position)
109 110 111 112 113
    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')
114
    plot.output(name='getting_started_2.pdf', xsize=16, ysize=16)