demo_wf1.py 3.95 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
44
about.warnings.off()
Marco Selig's avatar
Marco Selig committed
45
46

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

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

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

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

Ultima's avatar
Ultima committed
65
R = response_operator(x_space, sigma=0.01, mask=my_mask, assign=None) # define response
Ultima's avatar
Ultima committed
66
67
68
R = response_operator(x_space, assign=None) #
#R = identity_operator(x_space)

Marco Selig's avatar
Marco Selig committed
69
d_space = R.target                                               # get data space
Marco Selig's avatar
Marco Selig committed
70

Marco Selig's avatar
Marco Selig committed
71
72
73
# 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
Marco Selig's avatar
Marco Selig committed
74

75

Ultima's avatar
Ultima committed
76
d = R(s) + n                                                     # compute data
Marco Selig's avatar
Marco Selig committed
77

Marco Selig's avatar
Marco Selig committed
78
79
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
80

Ultima's avatar
Ultima committed
81
82
83
#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)
Marco Selig's avatar
Marco Selig committed
84

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


m = D(j, W=S, tol=1E-8, limii=100, note=True)


#temp_result = (D.inverse_times(m)-xi)

Ultima's avatar
Ultima committed
93

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

Ultima's avatar
Ultima committed
97
98
99
100
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
101

Ultima's avatar
Ultima committed
102
103
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
104
105


Ultima's avatar
Ultima committed
106
107
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
108
109


Ultima's avatar
Ultima committed
110
111
112
113
114

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)

Ultima's avatar
Ultima committed
115
#