Commit 09456c35 authored by Philipp Arras's avatar Philipp Arras
Browse files

Improvements demo 2

parent 27941a71
...@@ -15,18 +15,18 @@ ...@@ -15,18 +15,18 @@
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
############################################################ ###############################################################################
# Lognormal field reconstruction from Poisson data # Log-normal field reconstruction from Poissonian data with inhomogenous
# from inhomogenous exposure (in case for 2D mode) # exposure (in case for 2D mode)
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2) # 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
############################################################# ###############################################################################
import nifty5 as ift import nifty5 as ift
import numpy as np import numpy as np
def get_2D_exposure(): def exposure_2d():
# set up a structured exposure for 2D mode # Structured exposure for 2D mode
x_shape, y_shape = position_space.shape x_shape, y_shape = position_space.shape
exposure = np.ones(position_space.shape) exposure = np.ones(position_space.shape)
...@@ -37,80 +37,73 @@ def get_2D_exposure(): ...@@ -37,80 +37,73 @@ def get_2D_exposure():
exposure[:, x_shape*4//5:x_shape] *= .1 exposure[:, x_shape*4//5:x_shape] *= .1
exposure[:, x_shape//2:x_shape*3//2] *= 3. exposure[:, x_shape//2:x_shape*3//2] *= 3.
exposure = ift.Field.from_global_data(position_space, exposure) return ift.Field.from_global_data(position_space, exposure)
return exposure
if __name__ == '__main__': if __name__ == '__main__':
# FIXME description of the tutorial # FIXME All random seeds to 42
np.random.seed(41) np.random.seed(41)
# Choose space on which the signal field is defined
# Choose problem geometry and masking
mode = 2 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: if mode == 0:
# One dimensional regular grid with uniform exposure # One-dimensional regular grid with uniform exposure
position_space = ift.RGSpace(1024) position_space = ift.RGSpace(1024)
exposure = ift.Field.full(position_space, 1.) exposure = ift.Field.full(position_space, 1.)
elif mode == 1: elif mode == 1:
# Two-dimensional regular grid with inhomogeneous exposure # Two-dimensional regular grid with inhomogeneous exposure
position_space = ift.RGSpace([512, 512]) position_space = ift.RGSpace([512, 512])
exposure = get_2D_exposure() exposure = exposure_2d()
else: else:
# Sphere with uniform exposure # Sphere with uniform exposure
position_space = ift.HPSpace(128) position_space = ift.HPSpace(128)
exposure = ift.Field.full(position_space, 1.) exposure = ift.Field.full(position_space, 1.)
# Defining harmonic space and harmonic transform # Define harmonic space and harmonic transform
harmonic_space = position_space.get_default_codomain() harmonic_space = position_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(harmonic_space, position_space) HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
# domain over which the field's degrees of freedom live # Domain on which the field's degrees of freedom are defined
domain = ift.DomainTuple.make(harmonic_space) domain = ift.DomainTuple.make(harmonic_space)
# true value of the of those
position = ift.from_random('normal', domain)
# Define power spectrum and amplitudes # Define amplitude (square root of power spectrum)
def sqrtpspec(k): def sqrtpspec(k):
return 1. / (20. + k**2) return 1./(20. + k**2)
p_space = ift.PowerSpace(harmonic_space) p_space = ift.PowerSpace(harmonic_space)
pd = ift.PowerDistributor(harmonic_space, p_space) pd = ift.PowerDistributor(harmonic_space, p_space)
a = ift.PS_field(p_space, sqrtpspec) a = ift.PS_field(p_space, sqrtpspec)
A = pd(a) A = pd(a)
# Set up a sky model # Define sky model
sky = ift.exp(HT(ift.makeOp(A))) sky = ift.exp(HT(ift.makeOp(A)))
M = ift.DiagonalOperator(exposure) M = ift.DiagonalOperator(exposure)
GR = ift.GeometryRemover(position_space) GR = ift.GeometryRemover(position_space)
# Set up instrumental response # Define instrumental response
R = GR(M) R = GR(M)
# Generate mock data # Generate mock data and define likelihood operator
d_space = R.target[0] d_space = R.target[0]
lamb = R(sky) lamb = R(sky)
mock_position = ift.from_random('normal', domain) mock_position = ift.from_random('normal', domain)
data = lamb(mock_position) data = lamb(mock_position)
data = np.random.poisson(data.to_global_data().astype(np.float64)) data = np.random.poisson(data.to_global_data().astype(np.float64))
data = ift.Field.from_global_data(d_space, data) data = ift.Field.from_global_data(d_space, data)
# Compute likelihood and Hamiltonian
ic_newton = ift.DeltaEnergyController(name='Newton', iteration_limit=100,
tol_rel_deltaE=1e-8)
likelihood = ift.PoissonianEnergy(data)(lamb) likelihood = ift.PoissonianEnergy(data)(lamb)
# Settings for minimization
ic_newton = ift.DeltaEnergyController(
name='Newton', iteration_limit=100, tol_rel_deltaE=1e-8)
minimizer = ift.NewtonCG(ic_newton) minimizer = ift.NewtonCG(ic_newton)
# Minimize the Hamiltonian # Compute MAP solution by minimizing the information Hamiltonian
H = ift.Hamiltonian(likelihood) H = ift.Hamiltonian(likelihood)
H = ift.EnergyAdapter(position, H, want_metric=True) initial_position = ift.from_random('normal', domain)
H = ift.EnergyAdapter(initial_position, H, want_metric=True)
H, convergence = minimizer(H) H, convergence = minimizer(H)
# Plot results # Plotting
signal = sky(mock_position) signal = sky(mock_position)
reconst = sky(H.position) reconst = sky(H.position)
plot = ift.Plot() plot = ift.Plot()
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment