getting_started_2.py 4.23 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

Philipp Arras's avatar
Philipp Arras committed
24 25
import sys

Natalia Porqueres's avatar
Natalia Porqueres committed
26 27
import numpy as np

Philipp Arras's avatar
Philipp Arras committed
28 29
import nifty5 as ift

Martin Reinecke's avatar
Martin Reinecke committed
30

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

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

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

Natalia Porqueres's avatar
Natalia Porqueres committed
45

Jakob Knollmueller's avatar
Jakob Knollmueller committed
46
if __name__ == '__main__':
Torsten Ensslin's avatar
Torsten Ensslin committed
47
    np.random.seed(42)
Natalia Porqueres's avatar
Natalia Porqueres committed
48

Philipp Arras's avatar
Philipp Arras committed
49
    # Choose space on which the signal field is defined
Lukas Platz's avatar
Lukas Platz committed
50 51 52 53 54
    if len(sys.argv) == 2:
        mode = int(sys.argv[1])
    else:
        mode = 1

55
    if mode == 0:
Torsten Ensslin's avatar
Torsten Ensslin committed
56
        # One-dimensional regular grid with uniform exposure of 10
57
        position_space = ift.RGSpace(1024)
Torsten Ensslin's avatar
Torsten Ensslin committed
58
        exposure = ift.Field.full(position_space, 10.)
59 60 61
    elif mode == 1:
        # Two-dimensional regular grid with inhomogeneous exposure
        position_space = ift.RGSpace([512, 512])
Philipp Arras's avatar
Philipp Arras committed
62
        exposure = exposure_2d()
63
    else:
Torsten Ensslin's avatar
Torsten Ensslin committed
64
        # Sphere with uniform exposure of 100
65
        position_space = ift.HPSpace(128)
Torsten Ensslin's avatar
Torsten Ensslin committed
66
        exposure = ift.Field.full(position_space, 100.)
Philipp Arras's avatar
Philipp Arras committed
67 68

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

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

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

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

Philipp Arras's avatar
Philipp Arras committed
84
    # Define sky operator
Martin Reinecke's avatar
Martin Reinecke committed
85
    sky = ift.exp(HT(ift.makeOp(A)))
Natalia Porqueres's avatar
Natalia Porqueres committed
86

Jakob Knollmueller's avatar
Jakob Knollmueller committed
87 88
    M = ift.DiagonalOperator(exposure)
    GR = ift.GeometryRemover(position_space)
Philipp Arras's avatar
Philipp Arras committed
89
    # Define instrumental response
Martin Reinecke's avatar
Martin Reinecke committed
90
    R = GR(M)
Natalia Porqueres's avatar
Natalia Porqueres committed
91

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

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

Philipp Arras's avatar
Philipp Arras committed
106
    # Compute MAP solution by minimizing the information Hamiltonian
107
    H = ift.StandardHamiltonian(likelihood)
Philipp Arras's avatar
Philipp Arras committed
108 109
    initial_position = ift.from_random('normal', domain)
    H = ift.EnergyAdapter(initial_position, H, want_metric=True)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
110
    H, convergence = minimizer(H)
Natalia Porqueres's avatar
Natalia Porqueres committed
111

Philipp Arras's avatar
Philipp Arras committed
112
    # Plotting
113 114
    signal = sky(mock_position)
    reconst = sky(H.position)
Philipp Arras's avatar
Philipp Arras committed
115
    filename = "getting_started_2_mode_{}.png".format(mode)
116 117 118 119 120
    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')
Lukas Platz's avatar
Lukas Platz committed
121
    plot.output(xsize=12, ysize=10, name=filename)
Philipp Arras's avatar
Philipp Arras committed
122
    print("Saved results as '{}'.".format(filename))