demo_wf2.py 4.29 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 steepest descent.

"""
from __future__ import division
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49

#from pycallgraph import PyCallGraph
#from pycallgraph import Config
#from pycallgraph import GlobbingFilter
#from pycallgraph.output import GraphvizOutput
#
#config = Config()
#config.trace_filter = GlobbingFilter(exclude=[
#    'pycallgraph.*',
#    #'*.secret_function',
#])
#
#graphviz = GraphvizOutput(output_file='steepest_profiling.png')
#
#
ultimanet's avatar
ultimanet committed
50
# comment
51

Marco Selig's avatar
Marco Selig committed
52
from nifty import *                                              # version 0.8.0
53
from nifty.operators.nifty_minimization import steepest_descent_new
Marco Selig's avatar
Marco Selig committed
54

55
if __name__ == "__main__":
Marco Selig's avatar
Marco Selig committed
56

57 58
    # some signal space; e.g., a two-dimensional regular grid
    x_space = rg_space([256, 256])                                   # define signal space
Marco Selig's avatar
Marco Selig committed
59

60
    k_space = x_space.get_codomain()                                 # get conjugate space
Marco Selig's avatar
Marco Selig committed
61

62 63
    # some power spectrum
    power = (lambda k: 42 / (k + 1) ** 3)
Marco Selig's avatar
Marco Selig committed
64

65 66
    S = power_operator(k_space, codomain=x_space, spec=power)                          # define signal covariance
    s = S.get_random_field(domain=x_space, codomain=k_space)                           # generate signal
Marco Selig's avatar
Marco Selig committed
67

68 69
    R = response_operator(x_space, codomain=k_space, sigma=0.0, mask=1.0, assign=None) # define response
    d_space = R.target                                               # get data space
Marco Selig's avatar
Marco Selig committed
70

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
    d = R(s) + n                                                     # compute data
Marco Selig's avatar
Marco Selig committed
76

77 78
    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
79 80


81

82 83 84 85
    def energy(x):
        DIx = D.inverse_times(x)
        H = 0.5 * DIx.dot(x) - j.dot(x)
        return H
86 87


88 89 90 91
    def gradient(x):
        DIx = D.inverse_times(x)
        g = DIx - j
        return g
92 93


94 95 96
    def eggs(x):
        """
            Calculation of the information Hamiltonian and its gradient.
Marco Selig's avatar
Marco Selig committed
97

98 99 100 101 102 103
        """
    #    DIx = D.inverse_times(x)
    #    H = 0.5 * DIx.dot(x) - j.dot(x)                              # compute information Hamiltonian
    #    g = DIx - j                                                  # compute its gradient
    #    return H, g
        return energy(x), gradient(x)
Marco Selig's avatar
Marco Selig committed
104 105


106
    m = Field(x_space, codomain=k_space)                               # reconstruct map
107

108 109
    #with PyCallGraph(output=graphviz, config=config):
    m, convergence = steepest_descent(eggs=eggs, note=True)(m, tol=1E-3, clevel=3)
Marco Selig's avatar
Marco Selig committed
110

111
    m = Field(x_space, codomain=k_space)
112 113
    m, convergence = steepest_descent_new(energy, gradient, note=True)(m, tol=1E-3, clevel=3)
    #s.plot(title="signal")                                           # plot signal
114
    #d_ = Field(x_space, val=d.val, target=k_space)
115 116
    #d_.plot(title="data", vmin=s.min(), vmax=s.max())                # plot data
    #m.plot(title="reconstructed map", vmin=s.min(), vmax=s.max())    # plot map
Marco Selig's avatar
Marco Selig committed
117