demo_wf2.py 3.35 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
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
## 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 steepest descent.

"""
from __future__ import division
from nifty import *                                                   # version 0.6.0
from nifty.nifty_tools import *



# some signal space; e.g., a two-dimensional regular grid
x_space = rg_space([256, 256])                                        # define signal space

k_space = x_space.get_codomain()                                      # get conjugate space

# some power spectrum
power = (lambda k: 42 / (k + 1) ** 3)

S = power_operator(k_space, spec=power)                               # define signal covariance
s = S.get_random_field(domain=x_space)                                # generate signal

R = response_operator(x_space, sigma=0.0, mask=1.0, assign=None)      # define response
d_space = R.target                                                    # get data space

# some noise variance; e.g., 100
N = diagonal_operator(d_space, diag=100, bare=True)                   # define noise covariance
n = N.get_random_field(domain=d_space)                                # generate noise

d = R(s) + n                                                          # compute data



j = R.adjoint_times(N.inverse_times(d))                               # define information source
D = propagator_operator(S=S,N=N,R=R)                                  # define information propagator


def eggs(x):
    """
        Calculation of the information Hamiltonian and its gradient.

    """
    Dx = D.inverse_times(x)
    H = 0.5 * Dx.dot(x) - j.dot(x)                                    # compute information Hamiltonian
    g = Dx - j                                                        # compute its gradient
    return H,g


m = field(x_space, target=k_space)                                    # reconstruct map
m,convergence = steepest_descent(eggs=eggs, note=True)(m, tol=1E-4, clevel=3)



s.plot(title="signal")                                                # plot signal
d_ = field(x_space, val=d.val, target=k_space)
d_.plot(title="data", vmin=s.val.min(), vmax=s.val.max())             # plot data
m.plot(title="reconstructed map", vmin=s.val.min(), vmax=s.val.max()) # plot map