demo_wf1.py 3.52 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
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
x_space = rg_space([128, 128])                                   # define signal space
Ultima's avatar
Ultima committed
48
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
52

#x_space = gl_space(800)
Marco Selig's avatar
Marco Selig committed
53

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

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

Ultima's avatar
Ultima committed
59
S = power_operator(k_space, codomain=x_space, spec=power)                          # define signal covariance
Marco Selig's avatar
Marco Selig committed
60
s = S.get_random_field(domain=x_space)                           # generate signal
Ultima's avatar
Ultima committed
61
62
63
#my_mask = x_space.cast(1)
#my_mask[400:900,400:900] = 0
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
Marco Selig's avatar
Marco Selig committed
66
d_space = R.target                                               # get data space
Marco Selig's avatar
Marco Selig committed
67

Marco Selig's avatar
Marco Selig committed
68
69
70
# 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
71

72

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

Marco Selig's avatar
Marco Selig committed
75
76
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
77

Ultima's avatar
Ultima committed
78
79
80
#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
81

Ultima's avatar
Ultima committed
82
83
84
85
86
87
88
89
90
91
xi = field(x_space, random='gau', target=k_space)
m = D(xi, W=ident, tol=1E-8, limii=10, note=True, force=True)
temp_result = (D.inverse_times(m)-xi)
print (temp_result.dot(temp_result))
print (temp_result.val)

s.plot(title="signal", save = 'plot_s.png')                                           # plot signal
d_ = field(x_space, val=d.val, target=k_space)
d_.plot(title="data", vmin=s.min(), vmax=s.max(), save = 'plot_d.png')                # plot data
m.plot(title="reconstructed map", vmin=s.min(), vmax=s.max(), save = 'plot_m.png')    # plot map
Ultima's avatar
Ultima committed
92
#