From 8a74ad034e4f7e3a0d18f868a70d350728dd3fa1 Mon Sep 17 00:00:00 2001 From: Torsten Ensslin Date: Sun, 6 Jan 2019 12:11:08 +0100 Subject: [PATCH] improved demos, mainly by explaining what is going on and by streamlining --- demos/bernoulli_demo.py | 33 +++++++++++++++----------- demos/getting_started_1.py | 47 +++++++++++++++++++++++++++++-------- demos/getting_started_2.py | 48 ++++++++++++++++++++++++-------------- 3 files changed, 87 insertions(+), 41 deletions(-) diff --git a/demos/bernoulli_demo.py b/demos/bernoulli_demo.py index 536e14a0..39c79156 100644 --- a/demos/bernoulli_demo.py +++ b/demos/bernoulli_demo.py @@ -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") diff --git a/demos/getting_started_1.py b/demos/getting_started_1.py index 410db8b8..37b2d11b 100644 --- a/demos/getting_started_1.py +++ b/demos/getting_started_1.py @@ -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 diff --git a/demos/getting_started_2.py b/demos/getting_started_2.py index 8c74dab8..4014720d 100644 --- a/demos/getting_started_2.py +++ b/demos/getting_started_2.py @@ -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) -- GitLab