wiener_filter_curvature.py 3.66 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
# 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/>.
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.

19 20 21
from ..operators import EndomorphicOperator, InversionEnabler
from ..field import Field, sqrt
from ..sugar import power_analyze, power_synthesize
22

Martin Reinecke's avatar
PEP8  
Martin Reinecke committed
23

Martin Reinecke's avatar
Martin Reinecke committed
24
class WienerFilterCurvature(EndomorphicOperator):
Jakob Knollmueller's avatar
Jakob Knollmueller committed
25 26 27 28
    """The curvature of the WienerFilterEnergy.

    This operator implements the second derivative of the
    WienerFilterEnergy used in some minimization algorithms or
29 30
    for error estimates of the posterior maps. It is the
    inverse of the propagator operator.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
31 32 33 34

    Parameters
    ----------
    R: LinearOperator,
Martin Reinecke's avatar
Martin Reinecke committed
35 36 37
       The response operator of the Wiener filter measurement.
    N: EndomorphicOperator
       The noise covariance.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
38
    S: DiagonalOperator,
Martin Reinecke's avatar
Martin Reinecke committed
39
       The prior signal covariance
Jakob Knollmueller's avatar
Jakob Knollmueller committed
40 41
    """

42
    def __init__(self, R, N, S, inverter):
Martin Reinecke's avatar
Martin Reinecke committed
43
        super(WienerFilterCurvature, self).__init__()
44 45 46
        self.R = R
        self.N = N
        self.S = S
Martin Reinecke's avatar
Martin Reinecke committed
47 48
        op = R.adjoint*N.inverse*R + S.inverse
        self._op = InversionEnabler(op, inverter, S.times)
49

50 51
    @property
    def domain(self):
Martin Reinecke's avatar
Martin Reinecke committed
52
        return self._op.domain
53 54

    @property
Martin Reinecke's avatar
Martin Reinecke committed
55 56
    def capability(self):
        return self._op.capability
57

Martin Reinecke's avatar
Martin Reinecke committed
58 59
    def apply(self, x, mode):
        return self._op.apply(x, mode)
60

Martin Reinecke's avatar
Martin Reinecke committed
61
    def generate_posterior_sample(self):
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
        """ Generates a posterior sample from a Gaussian distribution with
        given mean and covariance.

        This method generates samples by setting up the observation and
        reconstruction of a mock signal in order to obtain residuals of the
        right correlation which are added to the given mean.

        Returns
        -------
        sample : Field
            Returns the a sample from the Gaussian of given mean and
            covariance.
        """

        power = sqrt(power_analyze(self.S.diagonal()))
        mock_signal = power_synthesize(power, real_signal=True)

        noise = self.N.diagonal().weight(-1)

        mock_noise = Field.from_random(random_type="normal",
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
82
                                       domain=self.N.domain, dtype=noise.dtype)
83 84 85 86 87 88
        mock_noise *= sqrt(noise)

        mock_data = self.R(mock_signal) + mock_noise

        mock_j = self.R.adjoint_times(self.N.inverse_times(mock_data))
        mock_m = self.inverse_times(mock_j)
Martin Reinecke's avatar
Martin Reinecke committed
89
        return mock_signal - mock_m
90 91 92 93

    def generate_posterior_sample2(self):
        power = self.S.diagonal()
        mock_signal = Field.from_random(random_type="normal",
94
                                        domain=self.S.domain, dtype=power.dtype)
95 96 97 98 99 100 101 102 103 104 105 106
        mock_signal *= sqrt(power)

        noise = self.N.diagonal()
        mock_noise = Field.from_random(random_type="normal",
                                       domain=self.N.domain, dtype=noise.dtype)
        mock_noise *= sqrt(noise)

        mock_data = self.R(mock_signal) + mock_noise

        mock_j = self.R.adjoint_times(self.N.inverse_times(mock_data))
        mock_m = self.inverse_times(mock_j)
        return mock_signal - mock_m