demo_wf1.py 4.16 KB
Newer Older
Marco Selig's avatar
Marco Selig committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
## NIFTY (Numerical Information Field Theory) has been developed at the
## Max-Planck-Institute for Astrophysics.
##
## Copyright (C) 2013 Max-Planck-Society
##
## Author: Marco Selig
## Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
##
## 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 <http://www.gnu.org/licenses/>.

"""
    ..                  __   ____   __
    ..                /__/ /   _/ /  /_
    ..      __ ___    __  /  /_  /   _/  __   __
    ..    /   _   | /  / /   _/ /  /   /  / /  /
    ..   /  / /  / /  / /  /   /  /_  /  /_/  /
    ..  /__/ /__/ /__/ /__/    \___/  \___   /  demo
    ..                               /______/

    NIFTY demo applying a Wiener filter using conjugate gradient.

"""
from __future__ import division
Ultima's avatar
Ultima committed
35
36
37

import matplotlib as mpl
mpl.use('Agg')
Ultima's avatar
Ultima committed
38
import gc
Ultima's avatar
Ultima committed
39
import imp
Ultima's avatar
Ultima committed
40
41
nifty = imp.load_module('nifty', None,
                        '/home/steininger/Downloads/nifty', ('','',5))
Ultima's avatar
Ultima committed
42

Marco Selig's avatar
Marco Selig committed
43
from nifty import *                                              # version 0.8.0
Marco Selig's avatar
Marco Selig committed
44

45
46
if __name__ == "__main__":
    about.warnings.off()
Marco Selig's avatar
Marco Selig committed
47

48
49
50
51
52
53
    # some signal space; e.g., a two-dimensional regular grid
    shape = [1024]
    x_space = rg_space(shape)
    #y_space = point_space(1280*1280)
    #x_space = hp_space(32)
    #x_space = gl_space(800)
Marco Selig's avatar
Marco Selig committed
54

55
    k_space = x_space.get_codomain()                                 # get conjugate space
Marco Selig's avatar
Marco Selig committed
56

57
58
    # some power spectrum
    power = (lambda k: 42 / (k + 1) ** 4)
Marco Selig's avatar
Marco Selig committed
59

60
61
62
63
64
65
    S = power_operator(k_space, codomain=x_space, spec=power)                          # define signal covariance
    s = S.get_random_field(domain=x_space)                           # generate signal
    #my_mask = x_space.cast(1)
    #stretch = 0.6
    #my_mask[shape[0]/2*stretch:shape[0]/2/stretch, shape[1]/2*stretch:shape[1]/2/stretch] = 0
    my_mask = 1
Ultima's avatar
Ultima committed
66

67
68
69
    R = response_operator(x_space, sigma=0.01, mask=my_mask, assign=None) # define response
    R = response_operator(x_space, assign=None) #
    #R = identity_operator(x_space)
Marco Selig's avatar
Marco Selig committed
70

71
    d_space = R.target                                               # get data space
Marco Selig's avatar
Marco Selig committed
72

73
74
75
    # some noise variance; e.g., signal-to-noise ratio of 1
    N = diagonal_operator(d_space, diag=s.var(), bare=True)          # define noise covariance
    n = N.get_random_field(domain=d_space)                           # generate noise
76

Marco Selig's avatar
Marco Selig committed
77

78
    d = R(s) + n                                                     # compute data
Marco Selig's avatar
Marco Selig committed
79

80
81
    j = R.adjoint_times(N.inverse_times(d))                          # define information source
    D = propagator_operator(S=S, N=N, R=R)                           # define information propagator
Marco Selig's avatar
Marco Selig committed
82

83
84
85
    #m = D(j, W=S, tol=1E-8, limii=100, note=True)
    #m = D(j, tol=1E-8, limii=20, note=True, force=True)
    ident = identity(x_space)
Ultima's avatar
Ultima committed
86

87
    #xi = field(x_space, random='gau', target=k_space)
Ultima's avatar
Ultima committed
88
89


90
    m = D(j, W=S, tol=1E-8, limii=100, note=True)
Ultima's avatar
Ultima committed
91
92


93
    #temp_result = (D.inverse_times(m)-xi)
Ultima's avatar
Ultima committed
94

Ultima's avatar
Ultima committed
95

96
97
    n_power = x_space.enforce_power(s.var()/np.prod(shape))
    s_power = S.get_power()
Ultima's avatar
Ultima committed
98

99
100
101
102
    s.plot(title="signal", save = 'plot_s.png')
    s.plot(title="signal power", power=True, other=power,
           mono=False, save = 'power_plot_s.png', nbin=1000, log=True,
           vmax = 100, vmin=10e-7)
Ultima's avatar
Ultima committed
103

104
105
    d_ = field(x_space, val=d.val, target=k_space)
    d_.plot(title="data", vmin=s.min(), vmax=s.max(), save = 'plot_d.png')
Ultima's avatar
Ultima committed
106
107


108
109
    n_ = field(x_space, val=n.val, target=k_space)
    n_.plot(title="data", vmin=s.min(), vmax=s.max(), save = 'plot_n.png')
Ultima's avatar
Ultima committed
110

Ultima's avatar
Ultima committed
111
112


113
114
115
116
117
    m.plot(title="reconstructed map", vmin=s.min(), vmax=s.max(), save = 'plot_m.png')
    m.plot(title="reconstructed power", power=True, other=(n_power, s_power),
           save = 'power_plot_m.png', vmin=0.001, vmax=10, mono=False)

    #