getting_started_2.py 2.71 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)
9 10 11 12 13 14 15 16
    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.

    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 82 83
    # Minimize the Hamiltonian
    H = ift.Hamiltonian(likelihood, ic_cg)
    H, convergence = minimizer(H)
Natalia Porqueres's avatar
Natalia Porqueres committed
84

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