getting_started_3.py 5.28 KB
Newer Older
 Philipp Arras committed Jul 10, 2018 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ``````# This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . # # Copyright(C) 2013-2018 Max-Planck-Society # `````` Torsten Ensslin committed Jan 06, 2019 16 17 18 19 20 21 22 23 ``````# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. ############################################################ # Non-linear tomography # data is line of sight (LOS) field # random lines (set mode=0), radial lines (mode=1) ############################################################# mode = 0 `````` Philipp Arras committed Jul 10, 2018 24 `````` `````` Jakob Knollmueller committed Jul 02, 2018 25 26 27 28 29 ``````import nifty5 as ift import numpy as np def get_random_LOS(n_los): `````` Torsten Ensslin committed Jan 06, 2019 30 `````` # Setting up LOS `````` 31 `````` starts = list(np.random.uniform(0, 1, (n_los, 2)).T) `````` Torsten Ensslin committed Jan 06, 2019 32 33 34 35 `````` if mode == 0: ends = list(np.random.uniform(0, 1, (n_los, 2)).T) else: ends = list(0.5+0*np.random.uniform(0, 1, (n_los, 2)).T) `````` Jakob Knollmueller committed Jul 02, 2018 36 37 `````` return starts, ends `````` Martin Reinecke committed Aug 02, 2018 38 `````` `````` Jakob Knollmueller committed Jul 02, 2018 39 ``````if __name__ == '__main__': `````` Philipp Arras committed Jul 10, 2018 40 `````` # FIXME description of the tutorial `````` Torsten Ensslin committed Jan 06, 2019 41 `````` np.random.seed(420) `````` Philipp Arras committed Oct 02, 2018 42 `````` np.seterr(all='raise') `````` 43 `````` position_space = ift.RGSpace([128, 128]) `````` Jakob Knollmueller committed Jul 02, 2018 44 `````` `````` Torsten Ensslin committed Jan 06, 2019 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 `````` # Setting up an amplitude model for the field A = ift.AmplitudeModel(position_space, 64, 3, 0.4, -5., 0.5, 0.4, 0.3) # made choices: # 64 spectral bins # # Spectral smoothness (affects Gaussian process part) # 3 = relatively high variance of spectral curbvature # 0.4 = quefrency mode below which cepstrum flattens # # power law part of spectrum: # -5= preferred power-law slope # 0.5 = low variance of power-law slope # # Gaussian process part of log-spectrum # 0.4 = y-intercept mean of additional power # 0.3 = y-intercept variance of additional power `````` Jakob Knollmueller committed Jul 02, 2018 62 63 64 `````` # Building the model for a correlated signal harmonic_space = position_space.get_default_codomain() ht = ift.HarmonicTransformOperator(harmonic_space, position_space) `````` Martin Reinecke committed Jul 27, 2018 65 `````` power_space = A.target[0] `````` Jakob Knollmueller committed Jul 02, 2018 66 67 `````` power_distributor = ift.PowerDistributor(harmonic_space, power_space) `````` Philipp Frank committed Sep 24, 2018 68 `````` vol = harmonic_space.scalar_dvol `````` Martin Reinecke committed Nov 05, 2018 69 `````` vol = ift.ScalingOperator(vol**(-0.5), harmonic_space) `````` Martin Reinecke committed Nov 26, 2018 70 71 `````` correlated_field = ht( vol(power_distributor(A))*ift.ducktape(harmonic_space, None, 'xi')) `````` Martin Reinecke committed Jul 03, 2018 72 `````` # alternatively to the block above one can do: `````` Philipp Frank committed Sep 24, 2018 73 `````` #correlated_field = ift.CorrelatedField(position_space, A) `````` Jakob Knollmueller committed Jul 02, 2018 74 75 `````` # apply some nonlinearity `````` Martin Reinecke committed Aug 11, 2018 76 `````` signal = ift.positive_tanh(correlated_field) `````` Martin Reinecke committed Jul 08, 2018 77 `````` `````` Jakob Knollmueller committed Jul 02, 2018 78 `````` # Building the Line of Sight response `````` Jakob Knollmueller committed Jul 02, 2018 79 `````` LOS_starts, LOS_ends = get_random_LOS(100) `````` Philipp Arras committed Aug 24, 2018 80 `````` R = ift.LOSResponse(position_space, starts=LOS_starts, ends=LOS_ends) `````` Jakob Knollmueller committed Jul 02, 2018 81 `````` # build signal response model and model likelihood `````` Martin Reinecke committed Aug 05, 2018 82 `````` signal_response = R(signal) `````` Jakob Knollmueller committed Jul 02, 2018 83 84 `````` # specify noise data_space = R.target `````` Jakob Knollmueller committed Jul 02, 2018 85 `````` noise = .001 `````` 86 `````` N = ift.ScalingOperator(noise, data_space) `````` Jakob Knollmueller committed Jul 02, 2018 87 `````` `````` Torsten Ensslin committed Jan 06, 2019 88 `````` # generate mock signal and data `````` Philipp Arras committed Sep 27, 2018 89 `````` MOCK_POSITION = ift.from_random('normal', signal_response.domain) `````` Martin Reinecke committed Jul 27, 2018 90 `````` data = signal_response(MOCK_POSITION) + N.draw_sample() `````` Jakob Knollmueller committed Jul 02, 2018 91 92 `````` # set up model likelihood `````` Philipp Arras committed Aug 24, 2018 93 `````` likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response) `````` Jakob Knollmueller committed Jul 02, 2018 94 95 `````` # set up minimization and inversion schemes `````` Jakob Knollmueller committed Jul 02, 2018 96 `````` ic_sampling = ift.GradientNormController(iteration_limit=100) `````` Martin Reinecke committed Aug 29, 2018 97 `````` ic_newton = ift.GradInfNormController( `````` Lukas Platz committed Oct 29, 2018 98 `````` name='Newton', tol=1e-7, iteration_limit=35) `````` Martin Reinecke committed Aug 18, 2018 99 `````` minimizer = ift.NewtonCG(ic_newton) `````` Jakob Knollmueller committed Jul 02, 2018 100 `````` `````` Jakob Knollmueller committed Jul 02, 2018 101 `````` # build model Hamiltonian `````` Martin Reinecke committed Jul 02, 2018 102 `````` H = ift.Hamiltonian(likelihood, ic_sampling) `````` Jakob Knollmueller committed Jul 02, 2018 103 `````` `````` Lukas Platz committed Oct 29, 2018 104 `````` INITIAL_POSITION = ift.MultiField.full(H.domain, 0.) `````` Jakob Knollmueller committed Jul 02, 2018 105 `````` position = INITIAL_POSITION `````` Jakob Knollmueller committed Jul 02, 2018 106 `````` `````` Martin Reinecke committed Aug 21, 2018 107 `````` plot = ift.Plot() `````` Martin Reinecke committed Aug 24, 2018 108 109 `````` plot.add(signal(MOCK_POSITION), title='Ground Truth') plot.add(R.adjoint_times(data), title='Data') `````` Philipp Arras committed Oct 02, 2018 110 `````` plot.add([A.force(MOCK_POSITION)], title='Power Spectrum') `````` Martin Reinecke committed Aug 24, 2018 111 `````` plot.output(ny=1, nx=3, xsize=24, ysize=6, name="setup.png") `````` Jakob Knollmueller committed Jul 02, 2018 112 `````` `````` Jakob Knollmueller committed Jul 02, 2018 113 `````` # number of samples used to estimate the KL `````` Jakob Knollmueller committed Jul 02, 2018 114 `````` N_samples = 20 `````` Torsten Ensslin committed Jan 06, 2019 115 116 `````` # five intermediate steps in minimization to illustrate progress `````` Lukas Platz committed Oct 29, 2018 117 `````` for i in range(5): `````` Torsten Ensslin committed Jan 06, 2019 118 `````` # set up KL `````` Martin Reinecke committed Aug 29, 2018 119 `````` KL = ift.KL_Energy(position, H, N_samples) `````` Torsten Ensslin committed Jan 06, 2019 120 `````` # minimize KL until iteration limit reached `````` Jakob Knollmueller committed Jul 02, 2018 121 `````` KL, convergence = minimizer(KL) `````` Torsten Ensslin committed Jan 06, 2019 122 `````` # read out position `````` Jakob Knollmueller committed Jul 02, 2018 123 `````` position = KL.position `````` Torsten Ensslin committed Jan 06, 2019 124 `````` # plot momentariy field `````` Martin Reinecke committed Aug 21, 2018 125 `````` plot = ift.Plot() `````` Martin Reinecke committed Aug 24, 2018 126 `````` plot.add(signal(KL.position), title="reconstruction") `````` Philipp Arras committed Oct 02, 2018 127 `````` plot.add([A.force(KL.position), A.force(MOCK_POSITION)], title="power") `````` Lukas Platz committed Oct 29, 2018 128 `````` plot.output(ny=1, ysize=6, xsize=16, name="loop-{:02}.png".format(i)) `````` Jakob Knollmueller committed Jul 02, 2018 129 `````` `````` Torsten Ensslin committed Jan 06, 2019 130 `````` # final plot `````` Philipp Arras committed Sep 21, 2018 131 `````` KL = ift.KL_Energy(position, H, N_samples) `````` Martin Reinecke committed Aug 21, 2018 132 `````` plot = ift.Plot() `````` Torsten Ensslin committed Jan 06, 2019 133 `````` # do statistics `````` Martin Reinecke committed Jul 04, 2018 134 `````` sc = ift.StatCalculator() `````` Martin Reinecke committed Aug 20, 2018 135 `````` for sample in KL.samples: `````` Philipp Arras committed Sep 11, 2018 136 `````` sc.add(signal(sample + KL.position)) `````` Martin Reinecke committed Aug 24, 2018 137 138 `````` plot.add(sc.mean, title="Posterior Mean") plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation") `````` Martin Reinecke committed Jul 04, 2018 139 `````` `````` Philipp Arras committed Oct 02, 2018 140 `````` powers = [A.force(s + KL.position) for s in KL.samples] `````` Martin Reinecke committed Aug 24, 2018 141 `````` plot.add( `````` Lukas Platz committed Oct 29, 2018 142 143 144 `````` powers + [A.force(KL.position), A.force(MOCK_POSITION)], title="Sampled Posterior Power Spectrum", linewidth=[1.]*len(powers) + [3., 3.]) `````` Martin Reinecke committed Aug 24, 2018 145 `` plot.output(ny=1, nx=3, xsize=24, ysize=6, name="results.png")``