Commit 8a74ad03 authored by Torsten Ensslin's avatar Torsten Ensslin

improved demos, mainly by explaining what is going on and by streamlining

parent 452c8ba0
......@@ -13,9 +13,15 @@
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
#####################################################################
# Bernoulli reconstruction
# reconstruct an event probability field with values between 0 and 1
# from the observed events
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
#####################################################################
import nifty5 as ift
import numpy as np
......@@ -25,17 +31,16 @@ if __name__ == '__main__':
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 uniform exposure
# position_space = ift.HPSpace(128)
# exposure = ift.Field.full(position_space, 1.)
mode = 2
if mode == 0:
# One dimensional regular grid
position_space = ift.RGSpace(1024)
elif mode == 1:
# Two dimensional regular grid
position_space = ift.RGSpace([512, 512])
else:
# Sphere
position_space = ift.HPSpace(128)
# Defining harmonic space and transform
harmonic_space = position_space.get_default_codomain()
......@@ -83,5 +88,5 @@ if __name__ == '__main__':
plot.add(reconstruction, title='reconstruction')
plot.add(GR.adjoint_times(data), title='data')
plot.add(sky(mock_position), title='truth')
plot.output(nx=3, xsize=16, ysize=5, title="results",
plot.output(nx=3, xsize=16, ysize=9, title="results",
name="bernoulli.png")
......@@ -13,14 +13,21 @@
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
############################################################
# Wiener filter demo code
# showing how measurement gaps are filled in
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
#############################################################
import nifty5 as ift
import numpy as np
def make_chess_mask(position_space):
# set up a checkerboard mask for 2D mode
mask = np.ones(position_space.shape)
for i in range(4):
for j in range(4):
......@@ -30,12 +37,14 @@ def make_chess_mask(position_space):
def make_random_mask():
# set up random mask for spherical mode
mask = ift.from_random('pm1', position_space)
mask = (mask+1)/2
return mask.to_global_data()
def mask_to_nan(mask, field):
# mask masked pixels for plotting data
masked_data = field.local_data.copy()
masked_data[mask.local_data == 0] = np.nan
return ift.from_local_data(field.domain, masked_data)
......@@ -43,9 +52,8 @@ def mask_to_nan(mask, field):
if __name__ == '__main__':
np.random.seed(42)
# FIXME description of the tutorial
# Choose problem geometry and masking
# Choose position space of signal and masking in measurement
mode = 1
if mode == 0:
# One dimensional regular grid
......@@ -60,29 +68,47 @@ if __name__ == '__main__':
position_space = ift.HPSpace(128)
mask = make_random_mask()
# specify harmonic space corresponding to signal
# the signal degree of freedom will live here
harmonic_space = position_space.get_default_codomain()
# harmonic transform from harmonic space to position space
HT = ift.HarmonicTransformOperator(harmonic_space, target=position_space)
# Set correlation structure with a power spectrum and build
# prior correlation covariance
# define 1D isotropic power spectrum
def power_spectrum(k):
return 100. / (20.+k**3)
power_space = ift.PowerSpace(harmonic_space)
PD = ift.PowerDistributor(harmonic_space, power_space)
# 1D spectral space in which power spectrum lives
power_space = ift.PowerSpace(harmonic_space)
# its mapping to (higher dimentional) harmonic space
PD = ift.PowerDistributor(harmonic_space, power_space)
# applying the mapping to the 1D spectrum
prior_correlation_structure = PD(ift.PS_field(power_space, power_spectrum))
# inserting the result into the diagonal of an harmonic space operator
S = ift.DiagonalOperator(prior_correlation_structure)
# S is now the field covariance
# Build instrument response consisting of a discretization, mask
# and harmonic transformaion
# data lives in geometry free space, thus we need to remove geometry
GR = ift.GeometryRemover(position_space)
# masking operator, to model that parts of the field where not observed
mask = ift.Field.from_global_data(position_space, mask)
Mask = ift.DiagonalOperator(mask)
# compile response operator out of
# - an harmic transform (to get to image space)
# - the application of the mask
# - the removal of geometric information
R = GR(Mask(HT))
# find out where the data then lives
data_space = GR.target
# Set the noise covariance
# Set the noise covariance N
noise = 5.
N = ift.ScalingOperator(noise, data_space)
......@@ -92,13 +118,14 @@ if __name__ == '__main__':
data = R(MOCK_SIGNAL) + MOCK_NOISE
# Build propagator D and information source j
j = R.adjoint_times(N.inverse_times(data))
D_inv = R.adjoint(N.inverse(R)) + S.inverse
# Make it invertible
j = R.adjoint_times(N.inverse_times(data))
# Make propagator invertible
IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3)
D = ift.InversionEnabler(D_inv, IC, approximation=S.inverse).inverse
# WIENER FILTER
# m = D j
m = D(j)
# PLOTTING
......
......@@ -13,14 +13,20 @@
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
############################################################
# Lognormal field reconstruction from Poisson data
# from inhomogenous exposure (in case for 2D mode)
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
#############################################################
import nifty5 as ift
import numpy as np
def get_2D_exposure():
# set up a structured exposure for 2D mode
x_shape, y_shape = position_space.shape
exposure = np.ones(position_space.shape)
......@@ -39,25 +45,33 @@ if __name__ == '__main__':
# FIXME description of the tutorial
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 = ift.Field.full(position_space, 1.)
# Two-dimensional regular grid with inhomogeneous exposure
position_space = ift.RGSpace([512, 512])
exposure = get_2D_exposure()
# Sphere with uniform exposure
# position_space = ift.HPSpace(128)
# exposure = ift.Field.full(position_space, 1.)
# Defining harmonic space and transform
# Choose problem geometry and masking
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:
# One dimensional regular grid with uniform exposure
position_space = ift.RGSpace(1024)
exposure = ift.Field.full(position_space, 1.)
elif mode == 1:
# Two-dimensional regular grid with inhomogeneous exposure
position_space = ift.RGSpace([512, 512])
exposure = get_2D_exposure()
else:
# Sphere with uniform exposure
position_space = ift.HPSpace(128)
exposure = ift.Field.full(position_space, 1.)
# Defining harmonic space and harmonic transform
harmonic_space = position_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
# domain over which the field's degrees of freedom live
domain = ift.DomainTuple.make(harmonic_space)
# true value of the of those
position = ift.from_random('normal', domain)
# Define power spectrum and amplitudes
......@@ -104,4 +118,4 @@ if __name__ == '__main__':
plot.add(GR.adjoint(data), title='Data')
plot.add(reconst, title='Reconstruction')
plot.add(reconst - signal, title='Residuals')
plot.output(name='getting_started_2.png', xsize=16, ysize=16)
plot.output(name='getting_started_2.pdf', xsize=16, ysize=16)
Markdown is supported
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