From a5859e542d205c19d36cc2ed4c79ccdc500850bb Mon Sep 17 00:00:00 2001 From: Martin Reinecke <martin@mpa-garching.mpg.de> Date: Mon, 29 Jan 2018 14:21:48 +0100 Subject: [PATCH] introduce better response operator --- demos/wiener_filter_via_curvature.py | 24 ++++++++++++------------ nifty4/__init__.py | 1 + nifty4/sugar.py | 9 ++++++++- 3 files changed, 21 insertions(+), 13 deletions(-) diff --git a/demos/wiener_filter_via_curvature.py b/demos/wiener_filter_via_curvature.py index d885f94d5..16ee1acb4 100644 --- a/demos/wiener_filter_via_curvature.py +++ b/demos/wiener_filter_via_curvature.py @@ -6,7 +6,7 @@ import numericalunits as nu if __name__ == "__main__": # In MPI mode, the random seed for numericalunits must be set by hand #nu.reset_units(43) - dimensionality = 1 + dimensionality = 2 np.random.seed(43) # Setting up variable parameters @@ -48,14 +48,14 @@ if __name__ == "__main__": np.random.seed(43) mock_power = ift.PS_field(power_space, power_spectrum) - mock_harmonic = ift.power_synthesize(mock_power, real_signal=True) - mock_signal = fft(mock_harmonic) + mock_signal = ift.power_synthesize(mock_power, real_signal=True) sensitivity = (1./nu.m)**dimensionality/nu.K - R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), - sensitivity=(sensitivity,)) + R = ift.GeometryRemover(signal_space) + R = R*ift.ScalingOperator(sensitivity, signal_space) + R = R*fft + R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma) data_domain = R.target[0] - R_harmonic = R*fft noise_amplitude = 1./signal_to_noise*field_sigma*sensitivity*((L/N_pixels)**dimensionality) print "noise amplitude:", noise_amplitude @@ -67,22 +67,22 @@ if __name__ == "__main__": data = R(mock_signal) + noise # Wiener filter - j = R_harmonic.adjoint_times(N.inverse_times(data)) + j = R.adjoint_times(N.inverse_times(data)) ctrl = ift.GradientNormController( verbose=True, tol_abs_gradnorm=1e-5/(nu.K*(nu.m**dimensionality))) inverter = ift.ConjugateGradient(controller=ctrl) wiener_curvature = ift.library.WienerFilterCurvature( - S=S, N=N, R=R_harmonic, inverter=inverter) + S=S, N=N, R=R, inverter=inverter) m = wiener_curvature.inverse_times(j) m_s = fft(m) sspace2 = ift.RGSpace(shape, distances=L/N_pixels/nu.m) - ift.plot(ift.Field(sspace2, mock_signal.val)/nu.K, title="mock_signal.png") + ift.plot(ift.Field(sspace2, fft(mock_signal).val)/nu.K, name="mock_signal.png") #data = ift.dobj.to_global_data(data.val).reshape(sspace2.shape) #data = ift.Field(sspace2, val=ift.dobj.from_global_data(data)) - ift.plot(ift.Field(sspace2, val=R.adjoint_times(data).val), title="data.png") - print "msig",np.min(mock_signal.val)/nu.K, np.max(mock_signal.val)/nu.K + ift.plot(ift.Field(sspace2, val=fft(R.adjoint_times(data)).val), name="data.png") + print "msig",np.min(fft(mock_signal).val)/nu.K, np.max(fft(mock_signal).val)/nu.K print "map",np.min(m_s.val)/nu.K, np.max(m_s.val)/nu.K - ift.plot(ift.Field(sspace2, m_s.val)/nu.K, title="map.png") + ift.plot(ift.Field(sspace2, m_s.val)/nu.K, name="map.png") diff --git a/nifty4/__init__.py b/nifty4/__init__.py index 8171219ec..b2d198f9f 100644 --- a/nifty4/__init__.py +++ b/nifty4/__init__.py @@ -24,6 +24,7 @@ from .operators.fft_operator import FFTOperator from .operators.fft_smoothing_operator import FFTSmoothingOperator from .operators.direct_smoothing_operator import DirectSmoothingOperator from .operators.response_operator import ResponseOperator +from .operators.response_operator import GeometryRemover from .operators.laplace_operator import LaplaceOperator from .operators.power_projection_operator import PowerProjectionOperator from .operators.inversion_enabler import InversionEnabler diff --git a/nifty4/sugar.py b/nifty4/sugar.py index 4a10a45c5..c4550d0d2 100644 --- a/nifty4/sugar.py +++ b/nifty4/sugar.py @@ -32,7 +32,8 @@ __all__ = ['PS_field', 'power_synthesize_nonrandom', 'create_power_field', 'create_power_operator', - 'create_composed_fft_operator'] + 'create_composed_fft_operator', + 'create_harmonic_smoothing_operator'] def PS_field(pspace, func, dtype=None): @@ -256,3 +257,9 @@ def create_composed_fft_operator(domain, codomain=None, all_to='other'): if res is None: raise ValueError("empty operator") return res + + +def create_harmonic_smoothing_operator(domain, space, sigma): + kfunc = domain[space].get_fft_smoothing_kernel_function(sigma) + return DiagonalOperator(kfunc(domain[space].get_k_length_array()), domain, + space) -- GitLab