getting_started_2.py 2.75 KB
Newer Older
Natalia Porqueres's avatar
Natalia Porqueres committed
1
2
3
4
import nifty5 as ift
import numpy as np


Jakob Knollmueller's avatar
Jakob Knollmueller committed
5
6
def get_2D_exposure():
    x_shape, y_shape = position_space.shape
Natalia Porqueres's avatar
Natalia Porqueres committed
7

Jakob Knollmueller's avatar
Jakob Knollmueller committed
8
    exposure = np.ones(position_space.shape)
Martin Reinecke's avatar
Martin Reinecke committed
9
10
11
12
13
14
    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.
15
16

    exposure = ift.Field.from_global_data(position_space, exposure)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
17
    return exposure
Natalia Porqueres's avatar
Natalia Porqueres committed
18

Natalia Porqueres's avatar
Natalia Porqueres committed
19

Jakob Knollmueller's avatar
Jakob Knollmueller committed
20
21
if __name__ == '__main__':
    # ABOUT THIS CODE
22
    np.random.seed(41)
Natalia Porqueres's avatar
Natalia Porqueres committed
23

Jakob Knollmueller's avatar
Jakob Knollmueller committed
24
    # Set up the position space of the signal
25
    #
Jakob Knollmueller's avatar
Jakob Knollmueller committed
26
27
28
    # # One dimensional regular grid with uniform exposure
    # position_space = ift.RGSpace(1024)
    # exposure = np.ones(position_space.shape)
Natalia Porqueres's avatar
Natalia Porqueres committed
29

30
    # Two-dimensional regular grid with inhomogeneous exposure
Jakob Knollmueller's avatar
Jakob Knollmueller committed
31
32
    position_space = ift.RGSpace([512, 512])
    exposure = get_2D_exposure()
Natalia Porqueres's avatar
Natalia Porqueres committed
33

Jakob Knollmueller's avatar
Jakob Knollmueller committed
34
35
    # # Sphere with with uniform exposure
    # position_space = ift.HPSpace(128)
36
    # exposure = ift.Field.full(position_space, 1.)
Natalia Porqueres's avatar
Natalia Porqueres committed
37

Jakob Knollmueller's avatar
Jakob Knollmueller committed
38
39
40
    # Defining harmonic space and transform
    harmonic_space = position_space.get_default_codomain()
    HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
Natalia Porqueres's avatar
Natalia Porqueres committed
41

Jakob Knollmueller's avatar
Jakob Knollmueller committed
42
43
    domain = ift.MultiDomain.make({'xi': harmonic_space})
    position = ift.from_random('normal', domain)
Natalia Porqueres's avatar
Natalia Porqueres committed
44

Jakob Knollmueller's avatar
Jakob Knollmueller committed
45
46
47
    # Define power spectrum and amplitudes
    def sqrtpspec(k):
        return 1. / (20. + k**2)
Natalia Porqueres's avatar
Natalia Porqueres committed
48

Jakob Knollmueller's avatar
Jakob Knollmueller committed
49
50
51
52
    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
53

Jakob Knollmueller's avatar
Jakob Knollmueller committed
54
55
56
57
58
    # Set up a sky model
    xi = ift.Variable(position)['xi']
    logsky_h = xi * A
    logsky = HT(logsky_h)
    sky = ift.PointwiseExponential(logsky)
Natalia Porqueres's avatar
Natalia Porqueres committed
59

Jakob Knollmueller's avatar
Jakob Knollmueller committed
60
61
62
63
    M = ift.DiagonalOperator(exposure)
    GR = ift.GeometryRemover(position_space)
    # Set up instrumental response
    R = GR * M
Natalia Porqueres's avatar
Natalia Porqueres committed
64

Jakob Knollmueller's avatar
Jakob Knollmueller committed
65
66
67
68
    # Generate mock data
    d_space = R.target[0]
    lamb = R(sky)
    mock_position = ift.from_random('normal', lamb.position.domain)
69
70
71
    data = lamb.at(mock_position).value
    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
72

Jakob Knollmueller's avatar
Jakob Knollmueller committed
73
74
75
76
    # Compute likelihood and Hamiltonian
    position = ift.from_random('normal', lamb.position.domain)
    likelihood = ift.PoissonianEnergy(lamb, data)
    ic_cg = ift.GradientNormController(iteration_limit=50)
77
78
    ic_newton = ift.GradientNormController(name='Newton', iteration_limit=50,
                                           tol_abs_gradnorm=1e-3)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
79
    minimizer = ift.RelaxedNewton(ic_newton)
Natalia Porqueres's avatar
Natalia Porqueres committed
80

Jakob Knollmueller's avatar
Jakob Knollmueller committed
81
    # Minimize the Hamiltonian
Martin Reinecke's avatar
Martin Reinecke committed
82
    H = ift.Hamiltonian(likelihood)
83
    H = H.make_invertible(ic_cg)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
84
    H, convergence = minimizer(H)
Natalia Porqueres's avatar
Natalia Porqueres committed
85

Jakob Knollmueller's avatar
Jakob Knollmueller committed
86
87
    # Plot results
    result_sky = sky.at(H.position).value
88
    # PLOTTING