Commit 27941a71 authored by Philipp Arras's avatar Philipp Arras
Browse files

Typos in getting_started_1

parent e9d11521
...@@ -15,36 +15,35 @@ ...@@ -15,36 +15,35 @@
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
###############################################################################
############################################################ # Compute a Wiener filter solution with NIFTy
# Wiener filter demo code # Shows how measurement gaps are filled in
# showing how measurement gaps are filled in
# 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 make_chess_mask(position_space): def make_checkerboard_mask(position_space):
# set up a checkerboard mask for 2D mode # Checkerboard mask for 2D mode
mask = np.ones(position_space.shape) mask = np.ones(position_space.shape)
for i in range(4): for i in range(4):
for j in range(4): for j in range(4):
if (i+j) % 2 == 0: if (i + j) % 2 == 0:
mask[i*128//4:(i+1)*128//4, j*128//4:(j+1)*128//4] = 0 mask[i*128//4:(i + 1)*128//4, j*128//4:(j + 1)*128//4] = 0
return mask return mask
def make_random_mask(): def make_random_mask():
# set up random mask for spherical mode # Random mask for spherical mode
mask = ift.from_random('pm1', position_space) mask = ift.from_random('pm1', position_space)
mask = (mask+1)/2 mask = (mask + 1)/2
return mask.to_global_data() return mask.to_global_data()
def mask_to_nan(mask, field): def mask_to_nan(mask, field):
# mask masked pixels for plotting data # Set masked pixels to nan for plotting
masked_data = field.local_data.copy() masked_data = field.local_data.copy()
masked_data[mask.local_data == 0] = np.nan masked_data[mask.local_data == 0] = np.nan
return ift.from_local_data(field.domain, masked_data) return ift.from_local_data(field.domain, masked_data)
...@@ -53,59 +52,61 @@ def mask_to_nan(mask, field): ...@@ -53,59 +52,61 @@ def mask_to_nan(mask, field):
if __name__ == '__main__': if __name__ == '__main__':
np.random.seed(42) np.random.seed(42)
# Choose position space of signal and masking in measurement # Choose space on which the signal field is defined
mode = 1 mode = 1
if mode == 0: if mode == 0:
# One dimensional regular grid # One-dimensional regular grid
position_space = ift.RGSpace([1024]) position_space = ift.RGSpace([1024])
mask = np.ones(position_space.shape) mask = np.ones(position_space.shape)
elif mode == 1: elif mode == 1:
# Two dimensional regular grid with chess mask # Two-dimensional regular grid with checkerboard mask
position_space = ift.RGSpace([128, 128]) position_space = ift.RGSpace([128, 128])
mask = make_chess_mask(position_space) mask = make_checkerboard_mask(position_space)
else: else:
# Sphere with half of its locations randomly masked # Sphere with half of its pixels randomly masked
position_space = ift.HPSpace(128) position_space = ift.HPSpace(128)
mask = make_random_mask() mask = make_random_mask()
# specify harmonic space corresponding to signal # Specify harmonic space corresponding to signal
# the signal degree of freedom will live here
harmonic_space = position_space.get_default_codomain() 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 # Harmonic transform from harmonic space to position space
# prior correlation covariance HT = ift.HarmonicTransformOperator(harmonic_space, target=position_space)
# define 1D isotropic power spectrum # Set prior correlation covariance with a power spectrum leading to
# homogeneous and isotropic statistics
def power_spectrum(k): def power_spectrum(k):
return 100. / (20.+k**3) return 100./(20. + k**3)
# 1D spectral space in which power spectrum lives
power_space = ift.PowerSpace(harmonic_space) # 1D spectral space on which the power spectrum is defined
# its mapping to (higher dimentional) harmonic space power_space = ift.PowerSpace(harmonic_space)
PD = ift.PowerDistributor(harmonic_space, power_space)
# applying the mapping to the 1D spectrum # Mapping to (higher dimensional) harmonic space
PD = ift.PowerDistributor(harmonic_space, power_space)
# Apply the mapping
prior_correlation_structure = PD(ift.PS_field(power_space, power_spectrum)) prior_correlation_structure = PD(ift.PS_field(power_space, power_spectrum))
# inserting the result into the diagonal of an harmonic space operator
# Insert the result into the diagonal of an harmonic space operator
S = ift.DiagonalOperator(prior_correlation_structure) S = ift.DiagonalOperator(prior_correlation_structure)
# S is now the field covariance # S is the prior field covariance
# Build instrument response consisting of a discretization, mask # Build instrument response consisting of a discretization, mask
# and harmonic transformaion # and harmonic transformaion
# data lives in geometry free space, thus we need to remove geometry # Data lives in a geometry-free space, thus the geometry is removed
GR = ift.GeometryRemover(position_space) GR = ift.GeometryRemover(position_space)
# masking operator, to model that parts of the field where not observed
# Masking operator to model that parts of the field have not been observed
mask = ift.Field.from_global_data(position_space, mask) mask = ift.Field.from_global_data(position_space, mask)
Mask = ift.DiagonalOperator(mask) Mask = ift.DiagonalOperator(mask)
# compile response operator out of
# - an harmic transform (to get to image space) # The response operator consists out of
# - an harmonic transform (to get to image space)
# - the application of the mask # - the application of the mask
# - the removal of geometric information # - the removal of geometric information
R = GR(Mask(HT)) R = GR(Mask(HT))
# find out where the data then lives
data_space = GR.target data_space = GR.target
# Set the noise covariance N # Set the noise covariance N
...@@ -117,30 +118,30 @@ if __name__ == '__main__': ...@@ -117,30 +118,30 @@ if __name__ == '__main__':
MOCK_NOISE = N.draw_sample() MOCK_NOISE = N.draw_sample()
data = R(MOCK_SIGNAL) + MOCK_NOISE data = R(MOCK_SIGNAL) + MOCK_NOISE
# Build propagator D and information source j # Build inverse propagator D and information source j
D_inv = R.adjoint(N.inverse(R)) + S.inverse D_inv = R.adjoint(N.inverse(R)) + S.inverse
j = R.adjoint_times(N.inverse_times(data)) j = R.adjoint_times(N.inverse_times(data))
# Make propagator invertible # Make D_inv invertible (via Conjugate Gradient)
IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3) IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3)
D = ift.InversionEnabler(D_inv, IC, approximation=S.inverse).inverse D = ift.InversionEnabler(D_inv, IC, approximation=S.inverse).inverse
# WIENER FILTER # Calculate WIENER FILTER solution
# m = D j
m = D(j) m = D(j)
# PLOTTING # Plotting
rg = isinstance(position_space, ift.RGSpace) rg = isinstance(position_space, ift.RGSpace)
plot = ift.Plot() plot = ift.Plot()
if rg and len(position_space.shape) == 1: if rg and len(position_space.shape) == 1:
plot.add([HT(MOCK_SIGNAL), GR.adjoint(data), HT(m)], plot.add(
label=['Mock signal', 'Data', 'Reconstruction'], [HT(MOCK_SIGNAL), GR.adjoint(data),
alpha=[1, .3, 1]) HT(m)],
plot.add(mask_to_nan(mask, HT(m-MOCK_SIGNAL)), title='Residuals') label=['Mock signal', 'Data', 'Reconstruction'],
alpha=[1, .3, 1])
plot.add(mask_to_nan(mask, HT(m - MOCK_SIGNAL)), title='Residuals')
plot.output(nx=2, ny=1, xsize=10, ysize=4, title="getting_started_1") plot.output(nx=2, ny=1, xsize=10, ysize=4, title="getting_started_1")
else: else:
plot.add(HT(MOCK_SIGNAL), title='Mock Signal') plot.add(HT(MOCK_SIGNAL), title='Mock Signal')
plot.add(mask_to_nan(mask, (GR(Mask)).adjoint(data)), plot.add(mask_to_nan(mask, (GR(Mask)).adjoint(data)), title='Data')
title='Data')
plot.add(HT(m), title='Reconstruction') plot.add(HT(m), title='Reconstruction')
plot.add(mask_to_nan(mask, HT(m-MOCK_SIGNAL)), title='Residuals') plot.add(mask_to_nan(mask, HT(m - MOCK_SIGNAL)), title='Residuals')
plot.output(nx=2, ny=2, xsize=10, ysize=10, title="getting_started_1") plot.output(nx=2, ny=2, xsize=10, ysize=10, title="getting_started_1")
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