log_normal_wiener_filter.py 3.93 KB
Newer Older
1
2
3
# -*- coding: utf-8 -*-

from nifty import *
Martin Reinecke's avatar
Martin Reinecke committed
4
import nifty as ift
5

Theo Steininger's avatar
Theo Steininger committed
6

7
8
9
10
if __name__ == "__main__":
    # Setting up parameters    |\label{code:wf_parameters}|
    correlation_length = 1.     # Typical distance over which the field is correlated
    field_variance = 2.         # Variance of field in position space
Theo Steininger's avatar
Theo Steininger committed
11
12
    response_sigma = 0.02       # Smoothing length of response (in same unit as L)
    signal_to_noise = 100       # The signal to noise ratio
13
14
15
16
17
18
19
    np.random.seed(43)          # Fixing the random seed
    def power_spectrum(k):      # Defining the power spectrum
        a = 4 * correlation_length * field_variance**2
        return a / (1 + k * correlation_length) ** 4

    # Setting up the geometry |\label{code:wf_geometry}|
    L = 2. # Total side-length of the domain
Theo Steininger's avatar
Theo Steininger committed
20
    N_pixels = 128 # Grid resolution (pixels per axis)
21
22
    #signal_space = RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
    signal_space = HPSpace(16)
23
    harmonic_space = FFTOperator.get_default_codomain(signal_space)
Martin Reinecke's avatar
Martin Reinecke committed
24
    fft = FFTOperator(harmonic_space, target=signal_space)
25
26
27
28
29
30
31
32
33
34
    power_space = PowerSpace(harmonic_space)

    # Creating the mock signal |\label{code:wf_mock_signal}|
    S = create_power_operator(harmonic_space, power_spectrum=power_spectrum)
    mock_power = Field(power_space, val=power_spectrum)
    mock_signal = fft(mock_power.power_synthesize(real_signal=True))

    # Setting up an exemplary response
    mask = Field(signal_space, val=1.)
    N10 = int(N_pixels/10)
Theo Steininger's avatar
Theo Steininger committed
35
    #mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
36
37
38
39
40
41
42
43
44
45
46
47
    R = ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}|
    data_domain = R.target[0]
    R_harmonic = ComposedOperator([fft, R], default_spaces=[0, 0])

    # Setting up the noise covariance and drawing a random noise realization
    N = DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise, bare=True)
    noise = Field.from_random(domain=data_domain, random_type='normal',
                              std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
    data = R(exp(mock_signal)) + noise #|\label{code:wf_mock_data}|

    # Wiener filter
    m0 = Field(harmonic_space, val=0.j)
Martin Reinecke's avatar
Martin Reinecke committed
48
49
50
51
52
53
54
    ctrl = ift.DefaultIterationController(verbose=False,tol_abs_gradnorm=1)
    ctrl2 = ift.DefaultIterationController(verbose=True,tol_abs_gradnorm=0.1, name="outer")
    inverter = ift.ConjugateGradient(controller=ctrl, preconditioner=S.times)
    energy = library.LogNormalWienerFilterEnergy(m0, data, R_harmonic, N, S, inverter=inverter)
    minimizer1 = VL_BFGS(controller=ctrl2,max_history_length=20)
    minimizer2 = RelaxedNewton(controller=ctrl2)
    minimizer3 = SteepestDescent(controller=ctrl2)
Theo Steininger's avatar
Theo Steininger committed
55
56


57
58
#    me1 = minimizer1(energy)
#    me2 = minimizer2(energy)
Martin Reinecke's avatar
Martin Reinecke committed
59
    me3 = minimizer3(energy)
Theo Steininger's avatar
Theo Steininger committed
60

61
62
63
64
#    m1 = fft(me1[0].position)
#    m2 = fft(me2[0].position)
#    m3 = fft(me3[0].position)
#
65
66
67
68
69
70
71
72
73
74
75


#    # Probing the variance
#    class Proby(DiagonalProberMixin, Prober): pass
#    proby = Proby(signal_space, probe_count=100)
#    proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z))))
#
#    sm = SmoothingOperator(signal_space, sigma=0.02)
#    variance = sm(proby.diagonal.weight(-1))

    #Plotting #|\label{code:wf_plotting}|
76
77
    #plotter = plotting.RG2DPlotter(color_map=plotting.colormaps.PlankCmap())
    plotter = plotting.HealpixPlotter(color_map=plotting.colormaps.PlankCmap())
Theo Steininger's avatar
Theo Steininger committed
78

79
80
81
    plotter.figure.xaxis = plotting.Axis(label='Pixel Index')
    plotter.figure.yaxis = plotting.Axis(label='Pixel Index')

Theo Steininger's avatar
Theo Steininger committed
82
    plotter.plot.zmax = 5; plotter.plot.zmin = -5
83
84
85
86
87
88
89
90
91
92
##    plotter(variance, path = 'variance.html')
#    #plotter.plot.zmin = exp(mock_signal.min());
#    plotter(mock_signal.real, path='mock_signal.html')
#    plotter(Field(signal_space, val=np.log(data.val.get_full_data().real).reshape(signal_space.shape)),
#            path = 'log_of_data.html')
#
#    plotter(m1.real, path='m_LBFGS.html')
#    plotter(m2.real, path='m_Newton.html')
#    plotter(m3.real, path='m_SteepestDescent.html')
#