getting_started_1.py 2.35 KB
 Jakob Knollmueller committed Jun 27, 2018 1 2 ``````import nifty5 as ift import numpy as np `````` Jakob Knollmueller committed Jun 28, 2018 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 `````` def make_chess_mask(): mask = np.ones(position_space.shape) for i in range(4): for j in range(4): if (i+j)%2 == 0: mask[i*128/4:(i+1)*128/4, j*128/4:(j+1)*128/4] = 0 return mask def make_random_mask(): mask = ift.from_random('pm1',position_space) mask = (mask+1)/2 return mask.val `````` Jakob Knollmueller committed Jun 27, 2018 18 ``````if __name__ == '__main__': `````` Jakob Knollmueller committed Jun 28, 2018 19 `````` ## describtion of the tutorial ### `````` Jakob Knollmueller committed Jun 27, 2018 20 `````` `````` Jakob Knollmueller committed Jun 28, 2018 21 `````` # Choose problem geometry and masking `````` Jakob Knollmueller committed Jun 27, 2018 22 `````` `````` Jakob Knollmueller committed Jun 28, 2018 23 24 25 `````` # # One dimensional regular grid # position_space = ift.RGSpace([1024]) # mask = np.ones(position_space.shape) `````` Jakob Knollmueller committed Jun 27, 2018 26 `````` `````` Jakob Knollmueller committed Jun 28, 2018 27 28 29 `````` # # Two dimensional regular grid with chess mask # position_space = ift.RGSpace([128,128]) # mask = make_chess_mask() `````` Jakob Knollmueller committed Jun 27, 2018 30 `````` `````` Jakob Knollmueller committed Jun 28, 2018 31 32 33 `````` # # Sphere with half of its locations randomly masked position_space = ift.HPSpace(128) mask = make_random_mask() `````` Jakob Knollmueller committed Jun 27, 2018 34 `````` `````` Jakob Knollmueller committed Jun 28, 2018 35 36 37 `````` # set up corresponding harmonic transform and space harmonic_space = position_space.get_default_codomain() HT = ift.HarmonicTransformOperator(harmonic_space, target=position_space) `````` Jakob Knollmueller committed Jun 27, 2018 38 `````` `````` Jakob Knollmueller committed Jun 28, 2018 39 40 41 42 43 44 `````` # set correlation structure with a power spectrum and build prior correlation covariance def power_spectrum(k): return 100. / (20.+k**3) power_space = ift.PowerSpace(harmonic_space) PD = ift.PowerDistributor(harmonic_space, power_space) prior_correlation_structure = PD(ift.PS_field(power_space, power_spectrum)) `````` Jakob Knollmueller committed Jun 27, 2018 45 `````` `````` Jakob Knollmueller committed Jun 28, 2018 46 `````` S = ift.DiagonalOperator(prior_correlation_structure) `````` Jakob Knollmueller committed Jun 27, 2018 47 `````` `````` Jakob Knollmueller committed Jun 28, 2018 48 49 50 51 52 53 54 55 56 57 58 `````` # build instrument response consisting of a discretization, mask and harmonic transformaion GR = ift.GeometryRemover(position_space) mask = ift.Field(position_space,val=mask) Mask = ift.DiagonalOperator(mask) R = GR * Mask * HT data_space = GR.target # setting the noise covariance noise = 5. N = ift.ScalingOperator(noise, data_space) `````` Jakob Knollmueller committed Jun 27, 2018 59 60 `````` `````` Jakob Knollmueller committed Jun 28, 2018 61 62 63 64 `````` # creating mock data MOCK_SIGNAL = S.draw_sample() MOCK_NOISE = N.draw_sample() data = R(MOCK_SIGNAL) + MOCK_NOISE `````` Jakob Knollmueller committed Jun 27, 2018 65 `````` `````` Jakob Knollmueller committed Jun 28, 2018 66 67 68 69 `````` # building 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 `````` Jakob Knollmueller committed Jun 27, 2018 70 `````` IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3) `````` Jakob Knollmueller committed Jun 28, 2018 71 72 73 74 75 76 77 `````` D = ift.InversionEnabler(D_inv,IC,approximation=S.inverse).inverse # WIENER FILTER m = D(j) ##PLOTTING #Truth, data, reconstruction, residuals `````` Jakob Knollmueller committed Jun 27, 2018 78 79 80 81 `````` ``````