getting_started_2.py 4.09 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13
# 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/>.
#
14
# Copyright(C) 2013-2019 Max-Planck-Society
15
#
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
import numpy as np

Philipp Arras's avatar
Philipp Arras committed
26 27
import nifty5 as ift

Martin Reinecke's avatar
Martin Reinecke committed
28

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

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

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

Natalia Porqueres's avatar
Natalia Porqueres committed
43

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

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

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

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

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

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

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

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

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

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

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

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