diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index ab7a311b9948262d91bb41e453e0acb82922114c..fb1c4c73bba77a02a634cc6ccda292dcad360fc1 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -63,133 +63,6 @@ pages: before_script: - export MPLBACKEND="agg" -run_critical_filtering: - stage: demo_runs - script: - - ls - - python setup.py install --user -f - - python3 setup.py install --user -f - - python demos/critical_filtering.py - - python3 demos/critical_filtering.py - artifacts: - paths: - - '*.png' - -run_nonlinear_critical_filter: - stage: demo_runs - script: - - python setup.py install --user -f - - python3 setup.py install --user -f - - python demos/nonlinear_critical_filter.py - - python3 demos/nonlinear_critical_filter.py - artifacts: - paths: - - '*.png' - -run_nonlinear_wiener_filter: - stage: demo_runs - script: - - python setup.py install --user -f - - python3 setup.py install --user -f - - python demos/nonlinear_wiener_filter.py - - python3 demos/nonlinear_wiener_filter.py - only: - - run_demos - artifacts: - paths: - - '*.png' - -# FIXME: disable for now. Fixing it is part of issue #244. -#run_poisson_demo: -# stage: demo_runs -# script: -# - python setup.py install --user -f -# - python3 setup.py install --user -f -# - python demos/poisson_demo.py -# - python3 demos/poisson_demo.py -# artifacts: -# paths: -# - '*.png' - -run_probing: - stage: demo_runs - script: - - python setup.py install --user -f - - python3 setup.py install --user -f - - python demos/probing.py - - python3 demos/probing.py - artifacts: - paths: - - '*.png' - -run_sampling: - stage: demo_runs - script: - - python setup.py install --user -f - - python3 setup.py install --user -f - - python demos/sampling.py - - python3 demos/sampling.py - artifacts: - paths: - - '*.png' - -run_tomography: - stage: demo_runs - script: - - python setup.py install --user -f - - python3 setup.py install --user -f - - python demos/tomography.py - - python3 demos/tomography.py - artifacts: - paths: - - '*.png' - -run_wiener_filter_data_space_noiseless: - stage: demo_runs - script: - - python setup.py install --user -f - - python3 setup.py install --user -f - - python demos/wiener_filter_data_space_noiseless.py - - python3 demos/wiener_filter_data_space_noiseless.py - artifacts: - paths: - - '*.png' - -run_wiener_filter_easy.py: - stage: demo_runs - script: - - python setup.py install --user -f - - python3 setup.py install --user -f - - python demos/wiener_filter_easy.py - - python3 demos/wiener_filter_easy.py - artifacts: - paths: - - '*.png' - -run_wiener_filter_via_curvature.py: - stage: demo_runs - script: - - pip install --user numericalunits - - pip3 install --user numericalunits - - python setup.py install --user -f - - python3 setup.py install --user -f - - python demos/wiener_filter_via_curvature.py - - python3 demos/wiener_filter_via_curvature.py - artifacts: - paths: - - '*.png' - -run_wiener_filter_via_hamiltonian.py: - stage: demo_runs - script: - - python setup.py install --user -f - - python3 setup.py install --user -f - - python demos/wiener_filter_via_hamiltonian.py - - python3 demos/wiener_filter_via_hamiltonian.py - artifacts: - paths: - - '*.png' - run_ipynb: stage: demo_runs script: diff --git a/demos/Wiener_Filter.ipynb b/demos/Wiener_Filter.ipynb index bf80aedf74362b16fb181e82542c6b941318e2b5..038efaffba09a6cf4a8a5e82e44255e832730ca0 100644 --- a/demos/Wiener_Filter.ipynb +++ b/demos/Wiener_Filter.ipynb @@ -717,21 +717,21 @@ "metadata": { "celltoolbar": "Slideshow", "kernelspec": { - "display_name": "Python 2", + "display_name": "Python 3", "language": "python", - "name": "python2" + "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.15" + "pygments_lexer": "ipython3", + "version": "3.6.5" } }, "nbformat": 4, diff --git a/demos/critical_filtering.py b/demos/critical_filtering.py deleted file mode 100644 index 0fc20a4fff09ed28dc0ae5ec79d682c44b8293a1..0000000000000000000000000000000000000000 --- a/demos/critical_filtering.py +++ /dev/null @@ -1,124 +0,0 @@ -import nifty5 as ift -import numpy as np -from nifty5 import Exponential, Linear, Tanh - -np.random.seed(42) - - -def adjust_zero_mode(m0, t0): - mtmp = m0.to_global_data().copy() - zero_position = len(m0.shape)*(0,) - zero_mode = mtmp[zero_position] - mtmp[zero_position] = zero_mode / abs(zero_mode) - ttmp = t0.to_global_data().copy() - ttmp[0] += 2 * np.log(abs(zero_mode)) - return (ift.Field.from_global_data(m0.domain, mtmp), - ift.Field.from_global_data(t0.domain, ttmp)) - - -if __name__ == "__main__": - noise_level = 1. - p_spec = (lambda k: (.5 / (k + 1) ** 3)) - - nonlinearity = Linear() - # Set up position space - s_space = ift.RGSpace((128, 128)) - h_space = s_space.get_default_codomain() - - # Define harmonic transformation and associated harmonic space - HT = ift.HarmonicTransformOperator(h_space, target=s_space) - - # Setting up power space - p_space = ift.PowerSpace(h_space, - binbounds=ift.PowerSpace.useful_binbounds( - h_space, logarithmic=True)) - - # Choosing the prior correlation structure and defining - # correlation operator - p = ift.PS_field(p_space, p_spec) - log_p = ift.log(p) - S = ift.create_power_operator(h_space, power_spectrum=lambda k: 1e-5) - - # Drawing a sample sh from the prior distribution in harmonic space - sh = S.draw_sample() - - # Choosing the measurement instrument - # Instrument = SmoothingOperator(s_space, sigma=0.01) - mask = np.ones(s_space.shape) - # mask[6000:8000] = 0. - mask[30:70, 30:70] = 0. - mask = ift.Field.from_global_data(s_space, mask) - - MaskOperator = ift.DiagonalOperator(mask) - R = ift.GeometryRemover(s_space) - R = R*MaskOperator - # R = R*HT - # R = R * ift.create_harmonic_smoothing_operator((harmonic_space,), 0, - # response_sigma) - MeasurementOperator = R - - d_space = MeasurementOperator.target - - Distributor = ift.PowerDistributor(target=h_space, power_space=p_space) - power = Distributor(ift.exp(0.5*log_p)) - # Creating the mock data - true_sky = nonlinearity(HT(power*sh)) - noiseless_data = MeasurementOperator(true_sky) - noise_amplitude = noiseless_data.val.std()*noise_level - N = ift.ScalingOperator(noise_amplitude**2, d_space) - n = N.draw_sample() - # Creating the mock data - d = noiseless_data + n - - m0 = ift.full(h_space, 1e-7) - t0 = ift.full(p_space, -4.) - power0 = Distributor.times(ift.exp(0.5 * t0)) - - plotdict = {"colormap": "Planck-like"} - zmin = true_sky.min() - zmax = true_sky.max() - ift.plot(true_sky, title="True sky", name="true_sky.png", **plotdict) - ift.plot(MeasurementOperator.adjoint_times(d), title="Data", - name="data.png", **plotdict) - - IC1 = ift.GradientNormController(name="IC1", iteration_limit=100, - tol_abs_gradnorm=1e-3) - LS = ift.LineSearchStrongWolfe(c2=0.02) - minimizer = ift.RelaxedNewton(IC1, line_searcher=LS) - - IC = ift.GradientNormController(iteration_limit=500, - tol_abs_gradnorm=1e-3) - - for i in range(20): - power0 = Distributor(ift.exp(0.5*t0)) - map0_energy = ift.library.NonlinearWienerFilterEnergy( - m0, d, MeasurementOperator, nonlinearity, HT, power0, N, S, IC) - - # Minimization with chosen minimizer - map0_energy, convergence = minimizer(map0_energy) - m0 = map0_energy.position - - # Updating parameters for correlation structure reconstruction - D0 = map0_energy.curvature - - # Initializing the power energy with updated parameters - power0_energy = ift.library.NonlinearPowerEnergy( - position=t0, d=d, N=N, xi=m0, D=D0, ht=HT, - Instrument=MeasurementOperator, nonlinearity=nonlinearity, - Distributor=Distributor, sigma=1., samples=2, - iteration_controller=IC) - - power0_energy = minimizer(power0_energy)[0] - - # Setting new power spectrum - t0 = power0_energy.position - - # break degeneracy between amplitude and excitation by setting - # excitation monopole to 1 - m0, t0 = adjust_zero_mode(m0, t0) - - ift.plot(nonlinearity(HT(power0*m0)), title="Reconstructed sky", - name="reconstructed_sky.png", zmin=zmin, zmax=zmax, **plotdict) - ymin = np.min(p.to_global_data()) - ift.plot([ift.exp(t0), p], title="Power spectra", - name="ps.png", ymin=ymin, **plotdict) diff --git a/demos/getting_started_1.py b/demos/getting_started_1.py new file mode 100644 index 0000000000000000000000000000000000000000..af5be4358e6b78c7b3c03315f85857dbf206d028 --- /dev/null +++ b/demos/getting_started_1.py @@ -0,0 +1,60 @@ +import nifty5 as ift +import numpy as np +from global_newton.models_other.apply_data import ApplyData +from global_newton.models_energy.hamiltonian import Hamiltonian +from nifty5.library.unit_log_gauss import UnitLogGauss +if __name__ == '__main__': + # s_space = ift.RGSpace([1024]) + s_space = ift.RGSpace([128,128]) + # s_space = ift.HPSpace(64) + + h_space = s_space.get_default_codomain() + total_domain = ift.MultiDomain.make({'xi': h_space}) + HT = ift.HarmonicTransformOperator(h_space, s_space) + + def sqrtpspec(k): + return 16. / (20.+k**2) + + GR = ift.GeometryRemover(s_space) + + d_space = GR.target + B = ift.FFTSmoothingOperator(s_space,0.1) + mask = np.ones(s_space.shape) + mask[64:89,76:100] = 0. + mask = ift.Field(s_space,val=mask) + Mask = ift.DiagonalOperator(mask) + R = GR * Mask * B + noise = 1. + N = ift.ScalingOperator(noise, d_space) + + p_space = ift.PowerSpace(h_space) + pd = ift.PowerDistributor(h_space, p_space) + position = ift.from_random('normal', total_domain) + xi = ift.Variable(position)['xi'] + a = ift.Constant(position, ift.PS_field(p_space, sqrtpspec)) + A = pd(a) + s_h = A * xi + s = HT(s_h) + Rs = R(s) + + + + MOCK_POSITION = ift.from_random('normal',total_domain) + data = Rs.at(MOCK_POSITION).value + N.draw_sample() + + NWR = ApplyData(data, ift.Field(d_space,val=noise), Rs) + + INITIAL_POSITION = ift.from_random('normal',total_domain) + likelihood = UnitLogGauss(INITIAL_POSITION, NWR) + + IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3) + inverter = ift.ConjugateGradient(controller=IC) + IC2 = ift.GradientNormController(name='Newton', iteration_limit=15) + minimizer = ift.RelaxedNewton(IC2) + + + H = Hamiltonian(likelihood, inverter) + H, convergence = minimizer(H) + result = s.at(H.position).value + + diff --git a/demos/nonlinear_critical_filter.py b/demos/nonlinear_critical_filter.py deleted file mode 100644 index 9c64f9a4f1542e1739c3bd241f3e2ead0347b5b5..0000000000000000000000000000000000000000 --- a/demos/nonlinear_critical_filter.py +++ /dev/null @@ -1,116 +0,0 @@ -import nifty5 as ift -from nifty5 import Exponential -import numpy as np -np.random.seed(42) - - -def adjust_zero_mode(m0, t0): - mtmp = m0.to_global_data().copy() - zero_position = len(m0.shape)*(0,) - zero_mode = mtmp[zero_position] - mtmp[zero_position] = zero_mode / abs(zero_mode) - ttmp = t0.to_global_data().copy() - ttmp[0] += 2 * np.log(abs(zero_mode)) - return (ift.Field.from_global_data(m0.domain, mtmp), - ift.Field.from_global_data(t0.domain, ttmp)) - - -if __name__ == "__main__": - noise_level = 1. - p_spec = (lambda k: (1. / (k + 1) ** 2)) - - # nonlinearity = Linear() - nonlinearity = Exponential() - # Set up position space - # s_space = ift.RGSpace([1024]) - s_space = ift.HPSpace(32) - h_space = s_space.get_default_codomain() - - # Define harmonic transformation and associated harmonic space - HT = ift.HarmonicTransformOperator(h_space, target=s_space) - - # Setting up power space - p_space = ift.PowerSpace(h_space, - binbounds=ift.PowerSpace.useful_binbounds( - h_space, logarithmic=True)) - # Choosing the prior correlation structure and defining - # correlation operator - p = ift.PS_field(p_space, p_spec) - log_p = ift.log(p) - S = ift.create_power_operator(h_space, lambda k: 1.) - - # Drawing a sample sh from the prior distribution in harmonic space - sh = S.draw_sample() - - # Choosing the measurement instrument - # Instrument = SmoothingOperator(s_space, sigma=0.01) - mask = np.ones(s_space.shape) - mask[6000:8000] = 0. - mask = ift.Field.from_global_data(s_space, mask) - - MaskOperator = ift.DiagonalOperator(mask) - R = ift.GeometryRemover(s_space) - R = R*MaskOperator - # R = R*HT - # R = R * ift.create_harmonic_smoothing_operator((harmonic_space,), 0, - # response_sigma) - MeasurementOperator = R - - d_space = MeasurementOperator.target - - Distributor = ift.PowerDistributor(target=h_space, power_space=p_space) - power = Distributor(ift.exp(0.5*log_p)) - # Creating the mock data - true_sky = nonlinearity(HT(power*sh)) - noiseless_data = MeasurementOperator(true_sky) - noise_amplitude = noiseless_data.val.std()*noise_level - N = ift.ScalingOperator(noise_amplitude**2, d_space) - n = N.draw_sample() - # Creating the mock data - d = noiseless_data + n - - m0 = ift.full(h_space, 1e-7) - t0 = ift.full(p_space, -4.) - power0 = Distributor.times(ift.exp(0.5 * t0)) - - IC1 = ift.GradientNormController(name="IC1", iteration_limit=100, - tol_abs_gradnorm=1e-3) - LS = ift.LineSearchStrongWolfe(c2=0.02) - minimizer = ift.RelaxedNewton(IC1, line_searcher=LS) - - IC = ift.GradientNormController(iteration_limit=500, - tol_abs_gradnorm=1e-3) - - for i in range(20): - power0 = Distributor(ift.exp(0.5*t0)) - map0_energy = ift.library.NonlinearWienerFilterEnergy( - m0, d, MeasurementOperator, nonlinearity, HT, power0, N, S, IC) - - # Minimization with chosen minimizer - map0_energy, convergence = minimizer(map0_energy) - m0 = map0_energy.position - - # Updating parameters for correlation structure reconstruction - D0 = map0_energy.curvature - - # Initializing the power energy with updated parameters - power0_energy = ift.library.NonlinearPowerEnergy( - position=t0, d=d, N=N, xi=m0, D=D0, ht=HT, - Instrument=MeasurementOperator, nonlinearity=nonlinearity, - Distributor=Distributor, sigma=1., samples=2, iteration_controller=IC) - - power0_energy = minimizer(power0_energy)[0] - - # Setting new power spectrum - t0 = power0_energy.position - - # break degeneracy between amplitude and excitation by setting - # excitation monopole to 1 - m0, t0 = adjust_zero_mode(m0, t0) - - plotdict = {"colormap": "Planck-like"} - ift.plot(true_sky, name="true_sky.png", **plotdict) - ift.plot(nonlinearity(HT(power0*m0)), - name="reconstructed_sky.png", **plotdict) - ift.plot(MeasurementOperator.adjoint_times(d), name="data.png", **plotdict) - ift.plot([ift.exp(t0), p], name="ps.png") diff --git a/demos/nonlinear_wiener_filter.py b/demos/nonlinear_wiener_filter.py deleted file mode 100644 index 01c2ea0f98ceba456b0421c9e7a66911801d9ae3..0000000000000000000000000000000000000000 --- a/demos/nonlinear_wiener_filter.py +++ /dev/null @@ -1,74 +0,0 @@ -import nifty5 as ift -from nifty5.library.nonlinearities import Linear, Exponential, Tanh -import numpy as np -np.random.seed(20) - -if __name__ == "__main__": - - noise_level = 0.3 - correlation_length = 0.1 - p_spec = lambda k: (1. / (k*correlation_length + 1) ** 4) - - nonlinearity = Tanh() - # nonlinearity = Linear() - # nonlinearity = Exponential() - - # Set up position space - s_space = ift.RGSpace(1024) - h_space = s_space.get_default_codomain() - - # Define harmonic transformation and associated harmonic space - HT = ift.HarmonicTransformOperator(h_space, target=s_space) - - S = ift.ScalingOperator(1., h_space) - - # Drawing a sample sh from the prior distribution in harmonic space - sh = S.draw_sample() - - # Choosing the measurement instrument - # Instrument = SmoothingOperator(s_space, sigma=0.01) - mask = np.ones(s_space.shape) - mask[600:800] = 0. - mask = ift.Field.from_global_data(s_space, mask) - - R = ift.GeometryRemover(s_space) * ift.DiagonalOperator(mask) - - d_space = R.target - - p_op = ift.create_power_operator(h_space, p_spec) - power = ift.sqrt(p_op(ift.full(h_space, 1.))) - - # Creating the mock data - true_sky = nonlinearity(HT(power*sh)) - noiseless_data = R(true_sky) - noise_amplitude = noiseless_data.val.std()*noise_level - N = ift.ScalingOperator(noise_amplitude**2, d_space) - n = N.draw_sample() - # Creating the mock data - d = noiseless_data + n - - IC1 = ift.GradientNormController(name="IC1", iteration_limit=100, - tol_abs_gradnorm=1e-4) - LS = ift.LineSearchStrongWolfe(c2=0.02) - minimizer = ift.RelaxedNewton(IC1, line_searcher=LS) - - IC = ift.GradientNormController(iteration_limit=2000, - tol_abs_gradnorm=1e-3) - - # initial guess - m = ift.full(h_space, 1e-7) - map_energy = ift.library.NonlinearWienerFilterEnergy( - m, d, R, nonlinearity, HT, power, N, S, IC) - - # Minimization with chosen minimizer - map_energy, convergence = minimizer(map_energy) - m = map_energy.position - - recsky = nonlinearity(HT(power*m)) - data = R.adjoint_times(d) - lo = np.min([true_sky.min(), recsky.min(), data.min()]) - hi = np.max([true_sky.max(), recsky.max(), data.max()]) - plotdict = {"colormap": "Planck-like", "ymin": lo, "ymax": hi} - ift.plot(true_sky, name="true_sky.png", **plotdict) - ift.plot(recsky, name="reconstructed_sky.png", **plotdict) - ift.plot(data, name="data.png", **plotdict) diff --git a/demos/paper_demos/cartesian_wiener_filter.py b/demos/paper_demos/cartesian_wiener_filter.py deleted file mode 100644 index cf851e1327910601769496ee6ffc6ee713195afd..0000000000000000000000000000000000000000 --- a/demos/paper_demos/cartesian_wiener_filter.py +++ /dev/null @@ -1,97 +0,0 @@ -import numpy as np -import nifty5 as ift - -if __name__ == "__main__": - signal_to_noise = 0.5 # The signal to noise ratio - - # Setting up parameters - L_1 = 2. # Total side-length of the domain - N_pixels_1 = 512 # Grid resolution (pixels per axis) - L_2 = 2. # Total side-length of the domain - N_pixels_2 = 512 # Grid resolution (pixels per axis) - correlation_length_1 = 1. - field_variance_1 = 2. # Variance of field in position space - response_sigma_1 = 0.05 # Smoothing length of response - correlation_length_2 = 1. - field_variance_2 = 2. # Variance of field in position space - response_sigma_2 = 0.01 # Smoothing length of response - - def power_spectrum_1(k): # note: field_variance**2 = a*k_0/4. - a = 4 * correlation_length_1 * field_variance_1**2 - return a / (1 + k * correlation_length_1) ** 4. - - def power_spectrum_2(k): # note: field_variance**2 = a*k_0/4. - a = 4 * correlation_length_2 * field_variance_2**2 - return a / (1 + k * correlation_length_2) ** 2.5 - - signal_space_1 = ift.RGSpace([N_pixels_1], distances=L_1/N_pixels_1) - harmonic_space_1 = signal_space_1.get_default_codomain() - signal_space_2 = ift.RGSpace([N_pixels_2], distances=L_2/N_pixels_2) - harmonic_space_2 = signal_space_2.get_default_codomain() - - signal_domain = ift.DomainTuple.make((signal_space_1, signal_space_2)) - harmonic_domain = ift.DomainTuple.make((harmonic_space_1, - harmonic_space_2)) - - ht_1 = ift.HarmonicTransformOperator(harmonic_domain, space=0) - ht_2 = ift.HarmonicTransformOperator(ht_1.target, space=1) - ht = ht_2*ht_1 - - S = (ift.create_power_operator(harmonic_domain, power_spectrum_1, 0) * - ift.create_power_operator(harmonic_domain, power_spectrum_2, 1)) - - np.random.seed(10) - mock_signal = S.draw_sample() - - # Setting up a exemplary response - N1_10 = int(N_pixels_1/10) - mask_1 = np.ones(signal_space_1.shape) - mask_1[N1_10*7:N1_10*9] = 0. - mask_1 = ift.Field.from_global_data(signal_space_1, mask_1) - - N2_10 = int(N_pixels_2/10) - mask_2 = np.ones(signal_space_2.shape) - mask_2[N2_10*7:N2_10*9] = 0. - mask_2 = ift.Field.from_global_data(signal_space_2, mask_2) - - R = ift.GeometryRemover(signal_domain) - R = R*ift.DiagonalOperator(mask_1, signal_domain, spaces=0) - R = R*ift.DiagonalOperator(mask_2, signal_domain, spaces=1) - R = R*ht - R = R * ift.create_harmonic_smoothing_operator(harmonic_domain, 0, - response_sigma_1) - R = R * ift.create_harmonic_smoothing_operator(harmonic_domain, 1, - response_sigma_2) - data_domain = R.target - - noiseless_data = R(mock_signal) - noise_amplitude = noiseless_data.val.std()/signal_to_noise - # Setting up the noise covariance and drawing a random noise realization - N = ift.ScalingOperator(noise_amplitude**2, data_domain) - noise = N.draw_sample() - data = noiseless_data + noise - - # Wiener filter - j = R.adjoint_times(N.inverse_times(data)) - ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1) - sampling_ctrl = ift.GradientNormController(name="sampling", - tol_abs_gradnorm=1e2) - wiener_curvature = ift.library.WienerFilterCurvature( - S=S, N=N, R=R, iteration_controller=ctrl, - iteration_controller_sampling=sampling_ctrl) - - m_k = wiener_curvature.inverse_times(j) - m = ht(m_k) - - plotdict = {"colormap": "Planck-like"} - plot_space = ift.RGSpace((N_pixels_1, N_pixels_2)) - ift.plot(ht(mock_signal).cast_domain(plot_space), - name='mock_signal.png', **plotdict) - ift.plot(data.cast_domain(plot_space), name='data.png', **plotdict) - ift.plot(m.cast_domain(plot_space), name='map.png', **plotdict) - # sampling the uncertainty map - mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 50) - ift.plot(ift.sqrt(variance).cast_domain(plot_space), - name="uncertainty.png", **plotdict) - ift.plot((mean+m).cast_domain(plot_space), - name="posterior_mean.png", **plotdict) diff --git a/demos/paper_demos/wiener_filter.py b/demos/paper_demos/wiener_filter.py deleted file mode 100644 index e9f199d2ab35acbab626904434879fd4ec9d9a33..0000000000000000000000000000000000000000 --- a/demos/paper_demos/wiener_filter.py +++ /dev/null @@ -1,67 +0,0 @@ -import nifty5 as ift -import numpy as np - - -if __name__ == "__main__": - # Setting up parameters - L = 2. # Total side-length of the domain - N_pixels = 512 # Grid resolution (pixels per axis) - correlation_length_scale = .2 # Typical distance over which the field is - # correlated - fluctuation_scale = 2. # Variance of field in position space - response_sigma = 0.05 # Smoothing length of response - signal_to_noise = 1.5 # The signal to noise ratio - np.random.seed(43) # Fixing the random seed - - def power_spectrum(k): # Defining the power spectrum - a = 4 * correlation_length_scale * fluctuation_scale**2 - return a / (1 + (k * correlation_length_scale)**2) ** 2 - - signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels) - harmonic_space = signal_space.get_default_codomain() - ht = ift.HarmonicTransformOperator(harmonic_space, target=signal_space) - - # Creating the mock signal - S = ift.create_power_operator(harmonic_space, - power_spectrum=power_spectrum) - mock_signal = S.draw_sample() - - # Setting up an exemplary response - mask = np.ones(signal_space.shape) - N10 = int(N_pixels/10) - mask[N10*5:N10*9, N10*5:N10*9] = 0. - mask = ift.Field.from_global_data(signal_space, mask).lock() - R = ift.GeometryRemover(signal_space) - R = R*ift.DiagonalOperator(mask) - R = R*ht - R = R * ift.create_harmonic_smoothing_operator((harmonic_space,), 0, - response_sigma) - data_domain = R.target[0] - - noiseless_data = R(mock_signal) - noise_amplitude = noiseless_data.val.std()/signal_to_noise - # Setting up the noise covariance and drawing a random noise realization - N = ift.ScalingOperator(noise_amplitude**2, data_domain) - noise = N.draw_sample() - data = noiseless_data + noise - - # Wiener filter - j = R.adjoint_times(N.inverse_times(data)) - ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=1e-2) - sampling_ctrl = ift.GradientNormController(name="sampling", - tol_abs_gradnorm=2e1) - wiener_curvature = ift.library.WienerFilterCurvature( - S=S, N=N, R=R, iteration_controller=ctrl, - iteration_controller_sampling=sampling_ctrl) - m_k = wiener_curvature.inverse_times(j) - m = ht(m_k) - - plotdict = {"colormap": "Planck-like"} - ift.plot(ht(mock_signal), name="mock_signal.png", **plotdict) - ift.plot(data.cast_domain(signal_space), name="data.png", **plotdict) - ift.plot(m, name="map.png", **plotdict) - - # sampling the uncertainty map - mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 50) - ift.plot(ift.sqrt(variance), name="uncertainty.png", **plotdict) - ift.plot(mean+m, name="posterior_mean.png", **plotdict) diff --git a/demos/poisson_demo.py b/demos/poisson_demo.py deleted file mode 100644 index 64a6e0ef52b41561e72a1702cd0288b372204ab3..0000000000000000000000000000000000000000 --- a/demos/poisson_demo.py +++ /dev/null @@ -1,143 +0,0 @@ -# Program to generate figures of article "Information theory for fields" -# by Torsten Ensslin, Annalen der Physik, submitted to special edition -# "Physics of Information" in April 2018, arXiv:1804.03350 -# to get exact figure of paper comment out lines marked by "COMMENT OUT" - -import numpy as np -import nifty5 as ift -import matplotlib.pyplot as plt - - -class Exp3(object): - def __call__(self, x): - return ift.exp(3*x) - - def derivative(self, x): - return 3*ift.exp(3*x) - - -if __name__ == "__main__": - np.random.seed(20) - - # Set up physical constants - nu = 1. # excitation field level - kappa = 10. # diffusion constant - eps = 1e-8 # small number to tame zero mode - sigma_n = 0.2 # noise level - sigma_n2 = sigma_n**2 - L = 1. # Total length of interval or volume the field lives on - nprobes = 1000 # Number of probes for uncertainty quantification used in paper - nprobes = 10 # COMMENT OUT TO REPRODUCE PAPER FIGURE EXACTLY - - # Define resolution (pixels per dimension) - N_pixels = 1024 - - # Define data gaps - N1a = int(0.6*N_pixels) - N1b = int(0.64*N_pixels) - N2a = int(0.67*N_pixels) - N2b = int(0.8*N_pixels) - - # Set up derived constants - amp = nu/(2*kappa) # spectral normalization - pow_spec = lambda k: amp / (eps + k**2) - lambda2 = 2*kappa*sigma_n2/nu # resulting correlation length squared - lambda1 = np.sqrt(lambda2) - pixel_width = L/N_pixels - x = np.arange(0, 1, pixel_width) - - # Set up the geometry - s_domain = ift.RGSpace([N_pixels], distances=pixel_width) - h_domain = s_domain.get_default_codomain() - HT = ift.HarmonicTransformOperator(h_domain, s_domain) - aHT = HT.adjoint - - # Create mock signal - Phi_h = ift.create_power_operator(h_domain, power_spectrum=pow_spec) - phi_h = Phi_h.draw_sample() - # remove zero mode - glob = phi_h.to_global_data() - glob[0] = 0. - phi_h = ift.Field.from_global_data(phi_h.domain, glob) - phi = HT(phi_h) - - # Setting up an exemplary response - GeoRem = ift.GeometryRemover(s_domain) - d_domain = GeoRem.target[0] - mask = np.ones(d_domain.shape) - mask[N1a:N1b] = 0. - mask[N2a:N2b] = 0. - fmask = ift.Field.from_global_data(d_domain, mask) - Mask = ift.DiagonalOperator(fmask) - R0 = Mask*GeoRem - R = R0*HT - - # Linear measurement scenario - N = ift.ScalingOperator(sigma_n2, d_domain) # Noise covariance - n = Mask(N.draw_sample()) # seting the noise to zero in masked region - d = R(phi_h) + n - - # Wiener filter - j = R.adjoint_times(N.inverse_times(d)) - IC = ift.GradientNormController(name="inverter", iteration_limit=500, - tol_abs_gradnorm=1e-3) - D = (ift.SandwichOperator.make(R, N.inverse) + Phi_h.inverse).inverse - D = ift.InversionEnabler(D, IC, approximation=Phi_h) - m = HT(D(j)) - - # Uncertainty - D = ift.SandwichOperator.make(aHT, D) # real space propagator - Dhat = ift.probe_with_posterior_samples(D.inverse, None, - nprobes=nprobes)[1] - sig = ift.sqrt(Dhat) - - # Plotting - x_mod = np.where(mask > 0, x, None) - plt.rcParams["text.usetex"] = True - c1 = (m-sig).to_global_data() - c2 = (m+sig).to_global_data() - plt.fill_between(x, c1, c2, color='pink', alpha=None) - plt.plot(x, phi.to_global_data(), label=r"$\varphi$", color='black') - plt.scatter(x_mod, d.to_global_data(), label=r'$d$', s=1, color='blue', - alpha=0.5) - plt.plot(x, m.to_global_data(), label=r'$m$', color='red') - plt.xlim([0, L]) - plt.ylim([-1, 1]) - plt.title('Wiener filter reconstruction') - plt.legend() - plt.savefig('Wiener_filter.pdf') - plt.close() - - nonlin = Exp3() - lam = R0(nonlin(HT(phi_h))) - data = ift.Field.from_local_data( - d_domain, np.random.poisson(lam.local_data).astype(np.float64)) - - # initial guess - psi0 = ift.full(h_domain, 1e-7) - energy = ift.library.PoissonEnergy(psi0, data, R0, nonlin, HT, Phi_h, IC) - IC1 = ift.GradientNormController(name="IC1", iteration_limit=200, - tol_abs_gradnorm=1e-4) - minimizer = ift.RelaxedNewton(IC1) - energy = minimizer(energy)[0] - - var = ift.probe_with_posterior_samples(energy.curvature, HT, nprobes)[1] - sig = ift.sqrt(var) - - m = HT(energy.position) - phi = HT(phi_h) - plt.rcParams["text.usetex"] = True - c1 = nonlin(m-sig).to_global_data() - c2 = nonlin(m+sig).to_global_data() - plt.fill_between(x, c1, c2, color='pink', alpha=None) - plt.plot(x, nonlin(phi).to_global_data(), label=r"$e^{3\varphi}$", - color='black') - plt.scatter(x_mod, data.to_global_data(), label=r'$d$', s=1, color='blue', - alpha=0.5) - plt.plot(x, nonlin(m).to_global_data(), - label=r'$e^{3\varphi_\mathrm{cl}}$', color='red') - plt.xlim([0, L]) - plt.ylim([-0.1, 7.5]) - plt.title('Poisson log-normal reconstruction') - plt.legend() - plt.savefig('Poisson.pdf') diff --git a/demos/probing.py b/demos/probing.py deleted file mode 100644 index 3361446b0907b35d91fac90aaee94fce580c6118..0000000000000000000000000000000000000000 --- a/demos/probing.py +++ /dev/null @@ -1,12 +0,0 @@ -import nifty5 as ift -import numpy as np - - -np.random.seed(42) -x = ift.RGSpace((8, 8)) - -f = ift.Field.from_random(domain=x, random_type='normal') -diagOp = ift.DiagonalOperator(f) - -diag = ift.probe_diagonal(diagOp, 1000) -ift.logger.info((f - diag).norm()) diff --git a/demos/sampling.py b/demos/sampling.py deleted file mode 100644 index 21f0d2e86fbf3901296070866e05bbd4ea5b4f84..0000000000000000000000000000000000000000 --- a/demos/sampling.py +++ /dev/null @@ -1,98 +0,0 @@ -import nifty5 as ift -import numpy as np -import matplotlib.pyplot as plt -from nifty5.sugar import create_power_operator - -np.random.seed(42) - -x_space = ift.RGSpace(1024) -h_space = x_space.get_default_codomain() - -d_space = x_space -N_hat = np.full(d_space.shape, 10.) -N_hat[400:450] = 0.0001 -N_hat = ift.Field.from_global_data(d_space, N_hat) -N = ift.DiagonalOperator(N_hat) - -FFT = ift.HarmonicTransformOperator(h_space, x_space) -R = ift.ScalingOperator(1., x_space) - - -def ampspec(k): return 1. / (1. + k**2.) - - -S = ift.ScalingOperator(1., h_space) -A = create_power_operator(h_space, ampspec) -s_h = S.draw_sample() -sky = FFT * A -s_x = sky(s_h) -n = N.draw_sample() -d = R(s_x) + n - -R_p = R * FFT * A -j = R_p.adjoint(N.inverse(d)) -D_inv = ift.SandwichOperator.make(R_p, N.inverse) + S.inverse - - -N_samps = 200 -N_iter = 100 - -tol = 1e-3 -IC = ift.GradientNormController(tol_abs_gradnorm=tol, iteration_limit=N_iter) -curv = ift.library.WienerFilterCurvature(S=S, N=N, R=R_p, - iteration_controller=IC, - iteration_controller_sampling=IC) -m_xi = curv.inverse_times(j) -samps_long = [curv.draw_sample(from_inverse=True) for i in range(N_samps)] - -tol = 1e2 -IC = ift.GradientNormController(tol_abs_gradnorm=tol, iteration_limit=N_iter) -curv = ift.library.WienerFilterCurvature(S=S, N=N, R=R_p, - iteration_controller=IC, - iteration_controller_sampling=IC) -samps_short = [curv.draw_sample(from_inverse=True) for i in range(N_samps)] - -# Compute mean -sc = ift.StatCalculator() -for samp in samps_long: - sc.add(samp) -m_x = sky(sc.mean + m_xi) - -plt.plot(d.to_global_data(), '+', label="data", alpha=.5) -plt.plot(s_x.to_global_data(), label="original") -plt.plot(m_x.to_global_data(), label="reconstruction") -plt.legend() -plt.savefig('reconstruction.png') -plt.close() - -pltdict = {'alpha': .3, 'linewidth': .2} -for i in range(N_samps): - if i == 0: - plt.plot(sky(samps_short[i]).to_global_data(), color='b', - label='Short samples (residuals)', - **pltdict) - plt.plot(sky(samps_long[i]).to_global_data(), color='r', - label='Long samples (residuals)', - **pltdict) - else: - plt.plot(sky(samps_short[i]).to_global_data(), color='b', **pltdict) - plt.plot(sky(samps_long[i]).to_global_data(), color='r', **pltdict) -plt.plot((s_x - m_x).to_global_data(), color='k', label='signal - mean') -plt.legend() -plt.savefig('samples_residuals.png') -plt.close() - -D_hat_old = ift.full(x_space, 0.).to_global_data() -D_hat_new = ift.full(x_space, 0.).to_global_data() -for i in range(N_samps): - D_hat_old += sky(samps_short[i]).to_global_data()**2 - D_hat_new += sky(samps_long[i]).to_global_data()**2 -plt.plot(np.sqrt(D_hat_old / N_samps), 'r--', label='Short uncertainty') -plt.plot(-np.sqrt(D_hat_old / N_samps), 'r--') -plt.fill_between(range(len(D_hat_new)), -np.sqrt(D_hat_new / N_samps), np.sqrt( - D_hat_new / N_samps), facecolor='0.5', alpha=0.5, - label='Long uncertainty') -plt.plot((s_x - m_x).to_global_data(), color='k', label='signal - mean') -plt.legend() -plt.savefig('uncertainty.png') -plt.close() diff --git a/demos/tomography.py b/demos/tomography.py deleted file mode 100644 index 3957e99c20405fa89740c7108a4022c04937a5d7..0000000000000000000000000000000000000000 --- a/demos/tomography.py +++ /dev/null @@ -1,48 +0,0 @@ -import nifty5 as ift -import numpy as np - - -if __name__ == "__main__": - np.random.seed(13) - - s_space = ift.RGSpace([512, 512]) - h_space = s_space.get_default_codomain() - ht = ift.HarmonicTransformOperator(h_space) - - pow_spec = (lambda k: 42 / (k+1)**4) - - S = ift.create_power_operator(h_space, power_spectrum=pow_spec) - - sh = S.draw_sample() - ss = ht(sh) - - def make_random_los(n): - starts = np.random.random((2, n)) - ends = np.random.random((2, n)) - return starts, ends - - nlos = 1000 - starts, ends = make_random_los(nlos) - - R = ift.library.LOSResponse( - s_space, starts=starts, ends=ends, - sigmas_up=np.full(nlos, 0.1), sigmas_low=np.full(nlos, 0.1)) - - Rh = R*ht - - signal_to_noise = 10 - noise_amplitude = Rh(sh).val.std()/signal_to_noise - N = ift.ScalingOperator(noise_amplitude**2, Rh.target) - n = N.draw_sample() - d = Rh(sh) + n - j = Rh.adjoint_times(N.inverse_times(d)) - ctrl = ift.GradientNormController(name="Iter", tol_abs_gradnorm=1e-10, - iteration_limit=300) - Di = ift.library.WienerFilterCurvature(S=S, R=Rh, N=N, - iteration_controller=ctrl) - mh = Di.inverse_times(j) - m = ht(mh) - - ift.plot(m, name="reconstruction.png") - ift.plot(ss, name="signal.png") - ift.plot(ht(j), name="j.png") diff --git a/demos/wiener_filter_data_space_noiseless.py b/demos/wiener_filter_data_space_noiseless.py deleted file mode 100644 index 9310f62660c99b631ef69e934dd5ca55f616c1a7..0000000000000000000000000000000000000000 --- a/demos/wiener_filter_data_space_noiseless.py +++ /dev/null @@ -1,115 +0,0 @@ -import numpy as np -import nifty5 as ift - - -# TODO: MAKE RESPONSE MPI COMPATIBLE OR USE LOS RESPONSE INSTEAD - -class CustomResponse(ift.LinearOperator): - """ - A custom operator that measures a specific points` - - An operator that is a delta measurement at certain points - """ - def __init__(self, domain, data_points): - self._domain = ift.DomainTuple.make(domain) - self._points = data_points - data_shape = ift.Field.full(domain, 0.).to_global_data()[data_points]\ - .shape - self._target = ift.DomainTuple.make(ift.UnstructuredDomain(data_shape)) - - def _times(self, x): - d = np.zeros(self._target.shape, dtype=np.float64) - d += x.to_global_data()[self._points] - return ift.from_global_data(self._target, d) - - def _adjoint_times(self, d): - x = np.zeros(self._domain.shape, dtype=np.float64) - x[self._points] += d.to_global_data() - return ift.from_global_data(self._domain, x) - - @property - def domain(self): - return self._domain - - @property - def target(self): - return self._target - - def apply(self, x, mode): - self._check_input(x, mode) - return self._times(x) if mode == self.TIMES else self._adjoint_times(x) - - @property - def capability(self): - return self.TIMES | self.ADJOINT_TIMES - - -if __name__ == "__main__": - np.random.seed(43) - # Set up physical constants - # Total length of interval or volume the field lives on, e.g. in meters - L = 2. - # Typical distance over which the field is correlated (in same unit as L) - correlation_length = 0.3 - # Variance of field in position space sqrt(<|s_x|^2>) (in same unit as s) - field_variance = 2. - # Smoothing length of response (in same unit as L) - response_sigma = 0.01 - # typical noise amplitude of the measurement - noise_level = 0. - - # Define resolution (pixels per dimension) - N_pixels = 256 - - # Set up derived constants - k_0 = 1./correlation_length - # defining a power spectrum with the right correlation length - # we later set the field variance to the desired value - unscaled_pow_spec = (lambda k: 1. / (1 + k/k_0) ** 4) - pixel_width = L/N_pixels - - # Set up the geometry - s_space = ift.RGSpace([N_pixels, N_pixels], distances=pixel_width) - h_space = s_space.get_default_codomain() - s_var = ift.get_signal_variance(unscaled_pow_spec, h_space) - pow_spec = (lambda k: unscaled_pow_spec(k)/s_var*field_variance**2) - - HT = ift.HarmonicTransformOperator(h_space, s_space) - - # Create mock data - - Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec) - sh = Sh.draw_sample() - - Rx = CustomResponse(s_space, [np.arange(0, N_pixels, 5)[:, np.newaxis], - np.arange(0, N_pixels, 2)[np.newaxis, :]]) - ift.extra.consistency_check(Rx) - a = ift.Field.from_random('normal', s_space) - b = ift.Field.from_random('normal', Rx.target) - R = Rx * HT - - noiseless_data = R(sh) - N = ift.ScalingOperator(noise_level**2, R.target) - n = N.draw_sample() - - d = noiseless_data + n - - # Wiener filter - - IC = ift.GradientNormController(name="inverter", iteration_limit=1000, - tol_abs_gradnorm=0.0001) - # setting up measurement precision matrix M - M = (ift.SandwichOperator.make(R.adjoint, Sh) + N) - M = ift.InversionEnabler(M, IC) - m = Sh(R.adjoint(M.inverse_times(d))) - - # Plotting - backprojection = Rx.adjoint(d) - reweighted_backprojection = (backprojection / backprojection.max() * - HT(sh).max()) - zmax = max(HT(sh).max(), reweighted_backprojection.max(), HT(m).max()) - zmin = min(HT(sh).min(), reweighted_backprojection.min(), HT(m).min()) - plotdict = {"colormap": "Planck-like", "zmax": zmax, "zmin": zmin} - ift.plot(HT(sh), name="mock_signal.png", **plotdict) - ift.plot(backprojection, name="backprojected_data.png", **plotdict) - ift.plot(HT(m), name="reconstruction.png", **plotdict) diff --git a/demos/wiener_filter_easy.py b/demos/wiener_filter_easy.py deleted file mode 100644 index 49749fe957e86e2efa9a84cfd4330999b6ad1e3a..0000000000000000000000000000000000000000 --- a/demos/wiener_filter_easy.py +++ /dev/null @@ -1,66 +0,0 @@ -import numpy as np -import nifty5 as ift - - -if __name__ == "__main__": - np.random.seed(43) - # Set up physical constants - # Total length of interval or volume the field lives on, e.g. in meters - L = 2. - # Typical distance over which the field is correlated (in same unit as L) - correlation_length = 0.1 - # Variance of field in position space sqrt(<|s_x|^2>) (in same unit as s) - field_variance = 2. - # Smoothing length of response (in same unit as L) - response_sigma = 0.01 - # typical noise amplitude of the measurement - noise_level = 1. - - # Define resolution (pixels per dimension) - N_pixels = 256 - - # Set up derived constants - k_0 = 1./correlation_length - #defining a power spectrum with the right correlation length - #we later set the field variance to the desired value - unscaled_pow_spec = (lambda k: 1. / (1 + k/k_0) ** 4) - pixel_width = L/N_pixels - - # Set up the geometry - s_space = ift.RGSpace([N_pixels, N_pixels], distances=pixel_width) - h_space = s_space.get_default_codomain() - s_var = ift.get_signal_variance(unscaled_pow_spec, h_space) - pow_spec = (lambda k: unscaled_pow_spec(k)/s_var*field_variance**2) - - HT = ift.HarmonicTransformOperator(h_space, s_space) - - # Create mock data - - Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec) - sh = Sh.draw_sample() - - R = HT*ift.create_harmonic_smoothing_operator((h_space,), 0, - response_sigma) - noiseless_data = R(sh) - N = ift.ScalingOperator(noise_level**2, s_space) - n = N.draw_sample() - - d = noiseless_data + n - - # Wiener filter - - j = R.adjoint_times(N.inverse_times(d)) - IC = ift.GradientNormController(name="inverter", iteration_limit=500, - tol_abs_gradnorm=0.1) - D = (ift.SandwichOperator.make(R, N.inverse) + Sh.inverse).inverse - D = ift.InversionEnabler(D, IC, approximation=Sh) - m = D(j) - - # Plotting - d_field = d.cast_domain(s_space) - zmax = max(HT(sh).max(), d_field.max(), HT(m).max()) - zmin = min(HT(sh).min(), d_field.min(), HT(m).min()) - plotdict = {"colormap": "Planck-like", "zmax": zmax, "zmin": zmin} - ift.plot(HT(sh), name="mock_signal.png", **plotdict) - ift.plot(d_field, name="data.png", **plotdict) - ift.plot(HT(m), name="reconstruction.png", **plotdict) diff --git a/demos/wiener_filter_via_curvature.py b/demos/wiener_filter_via_curvature.py deleted file mode 100644 index d83553e86d02e040551a4db76c97c49b3eb1439d..0000000000000000000000000000000000000000 --- a/demos/wiener_filter_via_curvature.py +++ /dev/null @@ -1,97 +0,0 @@ -import numpy as np -import nifty5 as ift -import numericalunits as nu - - -def init_nu(): - if ift.dobj.ntask > 1: - from mpi4py import MPI - comm = MPI.COMM_WORLD - data = np.random.randint(10000000, size=1) - data = comm.bcast(data, root=0) - nu.reset_units(data[0]) - - -if __name__ == "__main__": - # In MPI mode, the random seed for numericalunits must be set by hand - init_nu() - - # the number of dimensions for the problem - # For dimensionality>2, you should probably reduce N_pixels, otherwise - # the code may run out of memory - dimensionality = 2 - # Grid resolution (pixels per axis) - N_pixels = 512 - np.random.seed(43) - - # Setting up variable parameters - - # Typical distance over which the field is correlated - correlation_length = 0.05*nu.m - # sigma of field in position space sqrt(<|s_x|^2>) - field_sigma = 2. * nu.K - # smoothing length of response - response_sigma = 0.03*nu.m - # The signal to noise ratio - signal_to_noise = 1 - - # note that field_variance**2 = a*k_0/4. for this analytic form of power - # spectrum - def power_spectrum(k): - # RL FIXME: signal_amplitude is not how much signal varies - cldim = correlation_length**(2*dimensionality) - a = 4/(2*np.pi) * cldim * field_sigma**2 - # to be integrated over spherical shells later on - return a / (1 + (k*correlation_length)**(2*dimensionality)) ** 2 - - # Setting up the geometry - - # Total side-length of the domain - L = 2.*nu.m - shape = [N_pixels]*dimensionality - - signal_space = ift.RGSpace(shape, distances=L/N_pixels) - harmonic_space = signal_space.get_default_codomain() - HT = ift.HarmonicTransformOperator(harmonic_space, target=signal_space) - - # Creating the mock data - S = ift.create_power_operator(harmonic_space, - power_spectrum=power_spectrum) - np.random.seed(43) - - mock_signal = S.draw_sample() - - sensitivity = (1./nu.m)**dimensionality/nu.K - R = ift.GeometryRemover(signal_space) - R = R*ift.ScalingOperator(sensitivity, signal_space) - R = R*HT - R = R * ift.create_harmonic_smoothing_operator( - (harmonic_space,), 0, response_sigma) - data_domain = R.target[0] - - noiseless_data = R(mock_signal) - noise_amplitude = noiseless_data.val.std()/signal_to_noise - N = ift.ScalingOperator(noise_amplitude**2, data_domain) - noise = N.draw_sample() - data = noiseless_data + noise - - j = R.adjoint_times(N.inverse_times(data)) - ctrl = ift.GradientNormController( - name="inverter", tol_abs_gradnorm=1e-5/(nu.K*(nu.m**dimensionality))) - wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R, - iteration_controller=ctrl) - - m = wiener_curvature.inverse_times(j) - m_s = HT(m) - - sspace2 = ift.RGSpace(shape, distances=L/N_pixels/nu.m) - - zmax = max(m_s.max(), HT(mock_signal).max()) - zmin = min(m_s.min(), HT(mock_signal).min()) - plotdict = {"zmax": zmax/nu.K, "zmin": zmin/nu.K} - - ift.plot((HT(mock_signal)/nu.K).cast_domain(sspace2), - name="mock_signal.png", **plotdict) - ift.plot(data.cast_domain(sspace2), name="data.png") - ift.plot((m_s/nu.K).cast_domain(sspace2), name="reconstruction.png", - **plotdict) diff --git a/demos/wiener_filter_via_hamiltonian.py b/demos/wiener_filter_via_hamiltonian.py deleted file mode 100644 index bc2d2a37badddeec6d7e1cc6b75800ad35f2d679..0000000000000000000000000000000000000000 --- a/demos/wiener_filter_via_hamiltonian.py +++ /dev/null @@ -1,79 +0,0 @@ -import nifty5 as ift -import numpy as np - -np.random.seed(42) - - -if __name__ == "__main__": - # Set up position space - # s_space = ift.RGSpace([128, 128]) - s_space = ift.HPSpace(32) - - # Define associated harmonic space and harmonic transformation - h_space = s_space.get_default_codomain() - ht = ift.HarmonicTransformOperator(h_space, s_space) - - # Choose prior correlation structure and define correlation operator - p_spec = (lambda k: (42/(k+1)**3)) - - S = ift.create_power_operator(h_space, power_spectrum=p_spec) - - # Draw sample sh from the prior distribution in harmonic space - sh = S.draw_sample() - - # Choose measurement instrument - diag = np.ones(s_space.shape) - if len(s_space.shape) == 1: - diag[3000:7000] = 0 - elif len(s_space.shape) == 2: - diag[20:80, 20:80] = 0 - else: - raise NotImplementedError - - diag = ift.Field.from_global_data(s_space, diag) - Instrument = ift.DiagonalOperator(diag) - - # Add harmonic transformation to the instrument - R = Instrument*ht - noiseless_data = R(sh) - signal_to_noise = 1. - noise_amplitude = noiseless_data.val.std()/signal_to_noise - N = ift.ScalingOperator(noise_amplitude**2, s_space) - n = N.draw_sample() - - # Create mock data - d = noiseless_data + n - j = R.adjoint_times(N.inverse_times(d)) - - # Choose minimization strategy - ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1) - controller = ift.GradientNormController(name="min", tol_abs_gradnorm=0.1) - minimizer = ift.RelaxedNewton(controller=controller) - m0 = ift.full(h_space, 0.) - - # Initialize Wiener filter energy - energy = ift.library.WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S, - iteration_controller=ctrl, - iteration_controller_sampling=ctrl) - - energy, convergence = minimizer(energy) - m = energy.position - curv = energy.curvature - - # Generate plots - zmax = max(ht(sh).max(), ht(m).max()) - zmin = min(ht(sh).min(), ht(m).min()) - plotdict = {"zmax": zmax, "zmin": zmin, "colormap": "Planck-like"} - plotdict2 = {"colormap": "Planck-like"} - ift.plot(ht(sh), name="mock_signal.png", **plotdict) - ift.plot(ht(m), name="reconstruction.png", **plotdict) - - # Sample uncertainty map - mean, variance = ift.probe_with_posterior_samples(curv, ht, 50) - ift.plot(variance, name="posterior_variance.png", **plotdict2) - ift.plot(mean+ht(m), name="posterior_mean.png", **plotdict) - - # try to do the same with diagonal probing - variance = ift.probe_diagonal(ht*curv.inverse*ht.adjoint, 100) - # sm = ift.FFTSmoothingOperator(s_space, sigma=0.015) - ift.plot(variance, name="posterior_variance2.png", **plotdict)