getting_started_2.py 4.17 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 18 19 20 21 22
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.

############################################################
# Lognormal field reconstruction from Poisson data  
# from inhomogenous exposure (in case for 2D mode)
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
#############################################################
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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
28
def get_2D_exposure():
29
    # set up a 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 40

    exposure = ift.Field.from_global_data(position_space, exposure)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
41
    return 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 description of the tutorial
46
    np.random.seed(41)
Natalia Porqueres's avatar
Natalia Porqueres committed
47 48


49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
    # Choose problem geometry and masking
    mode = 2
    # 1D (mode=0), 2D (mode=1), or on the sphere (mode=2)

    # Set up the position space of the signal depending on mode
    if mode == 0:
        # One dimensional regular grid with uniform exposure
        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])
        exposure = get_2D_exposure()
    else:
        #  Sphere with uniform exposure
        position_space = ift.HPSpace(128)
        exposure = ift.Field.full(position_space, 1.)
        
    # Defining harmonic space and harmonic transform
Jakob Knollmueller's avatar
Jakob Knollmueller committed
68 69
    harmonic_space = position_space.get_default_codomain()
    HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
Natalia Porqueres's avatar
Natalia Porqueres committed
70

71
    # domain over which the field's degrees of freedom live
Martin Reinecke's avatar
Martin Reinecke committed
72
    domain = ift.DomainTuple.make(harmonic_space)
73 74
    
    # true value of the of those
Jakob Knollmueller's avatar
Jakob Knollmueller committed
75
    position = ift.from_random('normal', domain)
Natalia Porqueres's avatar
Natalia Porqueres committed
76

Jakob Knollmueller's avatar
Jakob Knollmueller committed
77 78 79
    # Define power spectrum and amplitudes
    def sqrtpspec(k):
        return 1. / (20. + k**2)
Natalia Porqueres's avatar
Natalia Porqueres committed
80

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

Jakob Knollmueller's avatar
Jakob Knollmueller committed
86
    # Set up a sky model
Martin Reinecke's avatar
Martin Reinecke committed
87
    sky = ift.exp(HT(ift.makeOp(A)))
Natalia Porqueres's avatar
Natalia Porqueres committed
88

Jakob Knollmueller's avatar
Jakob Knollmueller committed
89 90 91
    M = ift.DiagonalOperator(exposure)
    GR = ift.GeometryRemover(position_space)
    # Set up instrumental response
Martin Reinecke's avatar
Martin Reinecke committed
92
    R = GR(M)
Natalia Porqueres's avatar
Natalia Porqueres committed
93

Jakob Knollmueller's avatar
Jakob Knollmueller committed
94 95
    # Generate mock data
    d_space = R.target[0]
Martin Reinecke's avatar
Martin Reinecke committed
96
    lamb = R(sky)
Martin Reinecke's avatar
Martin Reinecke committed
97 98
    mock_position = ift.from_random('normal', domain)
    data = lamb(mock_position)
99 100
    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
101

Jakob Knollmueller's avatar
Jakob Knollmueller committed
102
    # Compute likelihood and Hamiltonian
103 104
    ic_newton = ift.DeltaEnergyController(name='Newton', iteration_limit=100,
                                          tol_rel_deltaE=1e-8)
105
    likelihood = ift.PoissonianEnergy(data)(lamb)
106
    minimizer = ift.NewtonCG(ic_newton)
Natalia Porqueres's avatar
Natalia Porqueres committed
107

Jakob Knollmueller's avatar
Jakob Knollmueller committed
108
    # Minimize the Hamiltonian
Martin Reinecke's avatar
Martin Reinecke committed
109
    H = ift.Hamiltonian(likelihood)
Martin Reinecke's avatar
fixes  
Martin Reinecke committed
110
    H = ift.EnergyAdapter(position, H, want_metric=True)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
111
    H, convergence = minimizer(H)
Natalia Porqueres's avatar
Natalia Porqueres committed
112

Jakob Knollmueller's avatar
Jakob Knollmueller committed
113
    # Plot results
Lukas Platz's avatar
Lukas Platz committed
114 115
    signal = sky(mock_position)
    reconst = sky(H.position)
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')
121
    plot.output(name='getting_started_2.pdf', xsize=16, ysize=16)