harmonic_propagator_operator.py 4.98 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# NIFTy
# Copyright (C) 2017  Theo Steininger
#
# Author: Theo Steininger
#
# 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/>.

from nifty.operators import EndomorphicOperator,\
                            FFTOperator,\
                            InvertibleOperatorMixin


class HarmonicPropagatorOperator(InvertibleOperatorMixin, EndomorphicOperator):
Theo Steininger's avatar
Theo Steininger committed
25
    """ NIFTY Harmonic Propagator Operator D.
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

    The propagator operator D, is known from the Wiener Filter.
    Its inverse functional form might look like:
    D = (S^(-1) + M)^(-1)
    D = (S^(-1) + N^(-1))^(-1)
    D = (S^(-1) + R^(\dagger) N^(-1) R)^(-1)
    In contrast to the PropagatorOperator the inference is done in the
    harmonic space.

    Parameters
    ----------
        S : LinearOperator
            Covariance of the signal prior.
        M : LinearOperator
            Likelihood contribution.
        R : LinearOperator
            Response operator translating signal to (noiseless) data.
        N : LinearOperator
            Covariance of the noise prior or the likelihood, respectively.
        inverter : class to invert explicitly defined operators
            (default:ConjugateGradient)
        preconditioner : Field
            numerical preconditioner to speed up convergence
49
50
51
        default_spaces : tuple of ints *optional*
            Defines on which space(s) of a given field the Operator acts by
            default (default: None)
52
53
54

    Attributes
    ----------
55
56
57
58
59
60
61
62
63
    domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
        The domain on which the Operator's input Field lives.
    target : tuple of DomainObjects, i.e. Spaces and FieldTypes
        The domain in which the outcome of the operator lives. As the Operator
        is endomorphic this is the same as its domain.
    unitary : boolean
        Indicates whether the Operator is unitary or not.
    self_adjoint : boolean
        Indicates whether the operator is self_adjoint or not.
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

    Raises
    ------
    ValueError
        is raised if
            * neither N nor M is given

    Notes
    -----

    Examples
    --------

    See Also
    --------
    Scientific reference
    https://arxiv.org/abs/0806.3474
Theo Steininger's avatar
Theo Steininger committed
81

82
83
84
85
    """

    # ---Overwritten properties and methods---

86
    def __init__(self, S, M=None, R=None, N=None, inverter=None,
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
                 preconditioner=None):
        """
            Sets the standard operator properties and `codomain`, `_A1`, `_A2`,
            and `RN` if required.

            Parameters
            ----------
            S : operator
                Covariance of the signal prior.
            M : operator
                Likelihood contribution.
            R : operator
                Response operator translating signal to (noiseless) data.
            N : operator
                Covariance of the noise prior or the likelihood, respectively.

        """
        # infer domain, and target
        if M is not None:
            self._codomain = M.domain
            self._likelihood = M.times

        elif N is None:
            raise ValueError("Either M or N must be given!")

        elif R is not None:
            self._codomain = R.domain
            self._likelihood = \
                lambda z: R.adjoint_times(N.inverse_times(R.times(z)))
        else:
            self._codomain = N.domain
            self._likelihood = lambda z: N.inverse_times(z)

        self._domain = S.domain
        self._S = S
        self._fft_S = FFTOperator(self._domain, target=self._codomain)

        super(HarmonicPropagatorOperator, self).__init__(inverter=inverter,
                                                 preconditioner=preconditioner)

    # ---Mandatory properties and methods---

    @property
    def domain(self):
        return self._domain

    @property
    def self_adjoint(self):
        return True

    @property
    def unitary(self):
        return False

    # ---Added properties and methods---
    def _likelihood_times(self, x, spaces=None):
        transformed_x = self._fft_S.times(x, spaces=spaces)
        y = self._likelihood(transformed_x)
145
        transformed_y = self._fft_S.adjoint_times(y, spaces=spaces)
146
147
148
149
150
        result = x.copy_empty()
        result.set_val(transformed_y, copy=False)
        return result

    def _inverse_times(self, x, spaces):
151
        pre_result = self._S.inverse_times(x, spaces)
152
153
154
155
        pre_result += self._likelihood_times(x)
        result = x.copy_empty()
        result.set_val(pre_result, copy=False)
        return result