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
ultimanet's avatar
ultimanet committed
39 40 41
#import imp
#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
    # some signal space; e.g., a two-dimensional regular grid
theos's avatar
theos committed
49 50
    #shape = [1024, 1024]
    #x_space = rg_space(shape)
51
    #y_space = point_space(1280*1280)
theos's avatar
theos committed
52
    x_space = hp_space(32)
53
    #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
    # some noise variance; e.g., signal-to-noise ratio of 1
    N = diagonal_operator(d_space, diag=s.var(), bare=True)          # define noise covariance
75
    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
    n_power = x_space.enforce_power(s.var()/x_space.dim)
97
    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
    d_ = Field(x_space, val=d.val, target=k_space)
105
    d_.plot(title="data", vmin=s.min(), vmax=s.max(), save = 'plot_d.png')
Ultima's avatar
Ultima committed
106 107


108
    n_ = Field(x_space, val=n.val, target=k_space)
109
    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)

    #