bernoulli_demo.py 3.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
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
17

18
19
#####################################################################
# Bernoulli reconstruction
Philipp Arras's avatar
Philipp Arras committed
20
# Reconstruct an event probability field with values between 0 and 1
21
22
23
# from the observed events
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
#####################################################################
Philipp Arras's avatar
Philipp Arras committed
24

Jakob Knollmueller's avatar
Jakob Knollmueller committed
25
26
import numpy as np

Philipp Arras's avatar
Philipp Arras committed
27
import nifty5 as ift
Jakob Knollmueller's avatar
Jakob Knollmueller committed
28
29
30
31
32

if __name__ == '__main__':
    np.random.seed(41)

    # Set up the position space of the signal
33
34
    mode = 2
    if mode == 0:
Philipp Arras's avatar
Philipp Arras committed
35
        # One-dimensional regular grid
36
37
        position_space = ift.RGSpace(1024)
    elif mode == 1:
Philipp Arras's avatar
Philipp Arras committed
38
        # Two-dimensional regular grid
39
40
41
42
        position_space = ift.RGSpace([512, 512])
    else:
        # Sphere
        position_space = ift.HPSpace(128)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
43

Philipp Arras's avatar
Philipp Arras committed
44
    # Define harmonic space and transform
Jakob Knollmueller's avatar
Jakob Knollmueller committed
45
46
47
    harmonic_space = position_space.get_default_codomain()
    HT = ift.HarmonicTransformOperator(harmonic_space, position_space)

Martin Reinecke's avatar
Martin Reinecke committed
48
    position = ift.from_random('normal', harmonic_space)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
49
50
51

    # Define power spectrum and amplitudes
    def sqrtpspec(k):
Philipp Arras's avatar
Philipp Arras committed
52
        return 1./(20. + k**2)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
53

Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
54
    A = ift.create_power_operator(harmonic_space, sqrtpspec)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
55

Philipp Arras's avatar
Philipp Arras committed
56
    # Set up a sky operator and instrumental response
Jakob Knollmueller's avatar
Jakob Knollmueller committed
57
    sky = ift.sigmoid(HT(A))
Jakob Knollmueller's avatar
Jakob Knollmueller committed
58
59
60
61
    GR = ift.GeometryRemover(position_space)
    R = GR

    # Generate mock data
Martin Reinecke's avatar
Martin Reinecke committed
62
    p = R(sky)
Martin Reinecke's avatar
Martin Reinecke committed
63
    mock_position = ift.from_random('normal', harmonic_space)
Martin Reinecke's avatar
Martin Reinecke committed
64
65
66
    tmp = p(mock_position).to_global_data().astype(np.float64)
    data = np.random.binomial(1, tmp)
    data = ift.Field.from_global_data(R.target, data)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
67
68

    # Compute likelihood and Hamiltonian
Martin Reinecke's avatar
Martin Reinecke committed
69
    position = ift.from_random('normal', harmonic_space)
Philipp Arras's avatar
Philipp Arras committed
70
    likelihood = ift.BernoulliEnergy(data)(p)
Philipp Arras's avatar
Philipp Arras committed
71
72
    ic_newton = ift.DeltaEnergyController(
        name='Newton', iteration_limit=100, tol_rel_deltaE=1e-8)
73
    minimizer = ift.NewtonCG(ic_newton)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
74
75
76
    ic_sampling = ift.GradientNormController(iteration_limit=100)

    # Minimize the Hamiltonian
77
    H = ift.StandardHamiltonian(likelihood, ic_sampling)
Martin Reinecke's avatar
fixes  
Martin Reinecke committed
78
    H = ift.EnergyAdapter(position, H, want_metric=True)
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
79
    # minimizer = ift.L_BFGS(ic_newton)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
80
    H, convergence = minimizer(H)
Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
81

Martin Reinecke's avatar
Martin Reinecke committed
82
    reconstruction = sky(H.position)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
83

84
85
86
87
    plot = ift.Plot()
    plot.add(reconstruction, title='reconstruction')
    plot.add(GR.adjoint_times(data), title='data')
    plot.add(sky(mock_position), title='truth')
Philipp Arras's avatar
Philipp Arras committed
88
    plot.output(nx=3, xsize=16, ysize=9, title="results", name="bernoulli.png")