demo_wf1.py 3.21 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
38

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

Marco Selig's avatar
Marco Selig committed
42
from nifty import *                                              # version 0.8.0
43
about.warnings.off()
Marco Selig's avatar
Marco Selig committed
44
45

# some signal space; e.g., a two-dimensional regular grid
Ultima's avatar
Ultima committed
46
47
48
x_space = rg_space([4280, 1280])                                   # define signal space
#y_space = point_space(1280*1280)

Marco Selig's avatar
Marco Selig committed
49
#x_space = hp_space(32)
Ultima's avatar
Ultima committed
50
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
Marco Selig's avatar
Marco Selig committed
60

Ultima's avatar
Ultima committed
61
R = response_operator(x_space, sigma=0.0, mask=1.0, assign=None) # define response
Marco Selig's avatar
Marco Selig committed
62
d_space = R.target                                               # get data space
Marco Selig's avatar
Marco Selig committed
63

Marco Selig's avatar
Marco Selig committed
64
65
66
# 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
67

68

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

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

Ultima's avatar
Ultima committed
74
m = D(j, W=S, tol=1E-8, limii=10, note=True)                               # reconstruct map
Marco Selig's avatar
Marco Selig committed
75

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