Commit a5859e54 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

introduce better response operator

parent 7fe10569
Pipeline #24148 passed with stage
in 4 minutes and 29 seconds
......@@ -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")
......@@ -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
......
......@@ -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)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment