bernoulli_demo.py 2.51 KB
Newer Older
Jakob Knollmueller's avatar
Jakob Knollmueller committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
import nifty5 as ift
import numpy as np


if __name__ == '__main__':
    # ABOUT THIS CODE
    np.random.seed(41)

    # Set up the position space of the signal
    #
    # # One dimensional regular grid with uniform exposure
    # position_space = ift.RGSpace(1024)
    # exposure = np.ones(position_space.shape)

    # Two-dimensional regular grid with inhomogeneous exposure
    position_space = ift.RGSpace([512, 512])

    # # Sphere with with uniform exposure
    # position_space = ift.HPSpace(128)
    # exposure = ift.Field.full(position_space, 1.)

    # Defining harmonic space and transform
    harmonic_space = position_space.get_default_codomain()
    HT = ift.HarmonicTransformOperator(harmonic_space, position_space)

    domain = ift.MultiDomain.make({'xi': harmonic_space})
    position = ift.from_random('normal', domain)

    # Define power spectrum and amplitudes
    def sqrtpspec(k):
        return 1. / (20. + k**2)

    p_space = ift.PowerSpace(harmonic_space)
    pd = ift.PowerDistributor(harmonic_space, p_space)
    a = ift.PS_field(p_space, sqrtpspec)
    A = pd(a)

    # Set up a sky model
    xi = ift.Variable(position)['xi']
    logsky_h = xi * A
    logsky = HT(logsky_h)
    sky = ift.PointwisePositiveTanh(logsky)

    GR = ift.GeometryRemover(position_space)
    # Set up instrumental response
    R = GR

    # Generate mock data
    d_space = R.target[0]
    p = R(sky)
    mock_position = ift.from_random('normal', p.position.domain)
    pp = p.at(mock_position).value
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
53
    data = np.random.binomial(1, pp.to_global_data().astype(np.float64))
Jakob Knollmueller's avatar
Jakob Knollmueller committed
54
55
56
57
58
59
60
61
62
63
64
65
66
    data = ift.Field.from_global_data(d_space, data)

    # Compute likelihood and Hamiltonian
    position = ift.from_random('normal', p.position.domain)
    likelihood = ift.BernoulliEnergy(p, data)
    ic_cg = ift.GradientNormController(iteration_limit=50)
    ic_newton = ift.GradientNormController(name='Newton', iteration_limit=30,
                                           tol_abs_gradnorm=1e-3)
    minimizer = ift.RelaxedNewton(ic_newton)
    ic_sampling = ift.GradientNormController(iteration_limit=100)

    # Minimize the Hamiltonian
    H = ift.Hamiltonian(likelihood, ic_sampling)
67
    H = H.make_invertible(ic_cg)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
68
69
    # minimizer = ift.SteepestDescent(ic_newton)
    H, convergence = minimizer(H)
Martin Reinecke's avatar
PEP8    
Martin Reinecke committed
70

Jakob Knollmueller's avatar
Jakob Knollmueller committed
71
72
73
74
75
    reconstruction = sky.at(H.position).value

    ift.plot(reconstruction, title='reconstruction', name='reconstruction.pdf')
    ift.plot(GR.adjoint_times(data), title='data', name='data.pdf')
    ift.plot(sky.at(mock_position).value, title='truth', name='truth.pdf')