getting_started_3.py 3.39 KB
 Jakob Knollmueller committed Jul 02, 2018 1 2 3 4 5 ``````import nifty5 as ift import numpy as np def get_random_LOS(n_los): `````` 6 7 `````` starts = list(np.random.uniform(0, 1, (n_los, 2)).T) ends = list(np.random.uniform(0, 1, (n_los, 2)).T) `````` Jakob Knollmueller committed Jul 02, 2018 8 9 `````` return starts, ends `````` 10 `````` `````` Jakob Knollmueller committed Jul 02, 2018 11 ``````if __name__ == '__main__': `````` 12 `````` # ## ABOUT THIS TUTORIAL `````` Jakob Knollmueller committed Jul 02, 2018 13 `````` np.random.seed(42) `````` 14 `````` position_space = ift.RGSpace([128, 128]) `````` Jakob Knollmueller committed Jul 02, 2018 15 `````` `````` Jakob Knollmueller committed Jul 02, 2018 16 `````` # Setting up an amplitude model `````` Martin Reinecke committed Jul 03, 2018 17 `````` A, amplitude_internals = ift.make_amplitude_model( `````` 18 `````` position_space, 16, 1, 10, -4., 1, 0., 1.) `````` Jakob Knollmueller committed Jul 02, 2018 19 20 21 22 23 24 `````` # Building the model for a correlated signal harmonic_space = position_space.get_default_codomain() ht = ift.HarmonicTransformOperator(harmonic_space, position_space) power_space = A.value.domain[0] power_distributor = ift.PowerDistributor(harmonic_space, power_space) `````` Martin Reinecke committed Jul 08, 2018 25 26 `````` position = ift.MultiField.from_dict( {'xi': ift.Field.from_random('normal', harmonic_space)}) `````` Jakob Knollmueller committed Jul 02, 2018 27 28 29 30 31 `````` xi = ift.Variable(position)['xi'] Amp = power_distributor(A) correlated_field_h = Amp * xi correlated_field = ht(correlated_field_h) `````` Martin Reinecke committed Jul 03, 2018 32 `````` # alternatively to the block above one can do: `````` Martin Reinecke committed Jul 03, 2018 33 `````` # correlated_field,_ = ift.make_correlated_field(position_space, A) `````` Jakob Knollmueller committed Jul 02, 2018 34 35 36 `````` # apply some nonlinearity signal = ift.PointwisePositiveTanh(correlated_field) `````` Martin Reinecke committed Jul 08, 2018 37 `````` `````` Jakob Knollmueller committed Jul 02, 2018 38 `````` # Building the Line of Sight response `````` Jakob Knollmueller committed Jul 02, 2018 39 `````` LOS_starts, LOS_ends = get_random_LOS(100) `````` Martin Reinecke committed Jul 03, 2018 40 41 `````` R = ift.LOSResponse(position_space, starts=LOS_starts, ends=LOS_ends) `````` Jakob Knollmueller committed Jul 02, 2018 42 `````` # build signal response model and model likelihood `````` Jakob Knollmueller committed Jul 02, 2018 43 `````` signal_response = R(signal) `````` Jakob Knollmueller committed Jul 02, 2018 44 45 `````` # specify noise data_space = R.target `````` Jakob Knollmueller committed Jul 02, 2018 46 `````` noise = .001 `````` 47 `````` N = ift.ScalingOperator(noise, data_space) `````` Jakob Knollmueller committed Jul 02, 2018 48 49 50 `````` # generate mock data MOCK_POSITION = ift.from_random('normal', signal.position.domain) `````` Jakob Knollmueller committed Jul 02, 2018 51 `````` data = signal_response.at(MOCK_POSITION).value + N.draw_sample() `````` Jakob Knollmueller committed Jul 02, 2018 52 53 54 55 56 `````` # set up model likelihood likelihood = ift.GaussianEnergy(signal_response, mean=data, covariance=N) # set up minimization and inversion schemes `````` Jakob Knollmueller committed Jul 02, 2018 57 58 `````` ic_cg = ift.GradientNormController(iteration_limit=10) ic_sampling = ift.GradientNormController(iteration_limit=100) `````` Jakob Knollmueller committed Jul 02, 2018 59 60 `````` ic_newton = ift.GradientNormController(name='Newton', iteration_limit=100) minimizer = ift.RelaxedNewton(ic_newton) `````` Jakob Knollmueller committed Jul 02, 2018 61 `````` `````` Jakob Knollmueller committed Jul 02, 2018 62 `````` # build model Hamiltonian `````` Martin Reinecke committed Jul 02, 2018 63 `````` H = ift.Hamiltonian(likelihood, ic_sampling) `````` Jakob Knollmueller committed Jul 02, 2018 64 `````` `````` 65 `````` INITIAL_POSITION = ift.from_random('normal', H.position.domain) `````` Jakob Knollmueller committed Jul 02, 2018 66 `````` position = INITIAL_POSITION `````` Jakob Knollmueller committed Jul 02, 2018 67 `````` `````` 68 69 70 `````` ift.plot(signal.at(MOCK_POSITION).value, name='truth.pdf') ift.plot(R.adjoint_times(data), name='data.pdf') ift.plot([A.at(MOCK_POSITION).value], name='power.pdf') `````` Jakob Knollmueller committed Jul 02, 2018 71 `````` `````` Jakob Knollmueller committed Jul 02, 2018 72 `````` # number of samples used to estimate the KL `````` Jakob Knollmueller committed Jul 02, 2018 73 `````` N_samples = 20 `````` Jakob Knollmueller committed Jul 02, 2018 74 75 `````` for i in range(5): H = H.at(position) `````` Martin Reinecke committed Jul 03, 2018 76 `````` samples = [H.metric.draw_sample(from_inverse=True) `````` 77 `````` for _ in range(N_samples)] `````` Jakob Knollmueller committed Jul 02, 2018 78 `````` `````` Martin Reinecke committed Jul 02, 2018 79 80 `````` KL = ift.SampledKullbachLeiblerDivergence(H, samples) KL = KL.makeInvertible(ic_cg) `````` Jakob Knollmueller committed Jul 02, 2018 81 82 83 `````` KL, convergence = minimizer(KL) position = KL.position `````` 84 `````` ift.plot(signal.at(position).value, name='reconstruction.pdf') `````` Jakob Knollmueller committed Jul 02, 2018 85 `````` `````` 86 87 `````` ift.plot([A.at(position).value, A.at(MOCK_POSITION).value], name='power.pdf') `````` Jakob Knollmueller committed Jul 02, 2018 88 `````` `````` Martin Reinecke committed Jul 04, 2018 89 `````` sc = ift.StatCalculator() `````` Jakob Knollmueller committed Jul 02, 2018 90 `````` for sample in samples: `````` Martin Reinecke committed Jul 04, 2018 91 92 93 94 95 `````` sc.add(signal.at(sample+position).value) ift.plot(sc.mean, name='avrg.pdf') ift.plot(ift.sqrt(sc.var), name='std.pdf') powers = [A.at(s+position).value for s in samples] `````` 96 97 `````` ift.plot([A.at(position).value, A.at(MOCK_POSITION).value]+powers, name='power.pdf')``````