sampling_enabler.py 3.39 KB
Newer Older
Jakob Knollmueller's avatar
Jakob Knollmueller committed
1 2 3 4 5 6 7 8 9 10 11 12 13
# 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/>.
#
14
# Copyright(C) 2013-2019 Max-Planck-Society
Jakob Knollmueller's avatar
Jakob Knollmueller committed
15
#
16
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
17

Philipp Arras's avatar
Philipp Arras committed
18
import numpy as np
19

20
from ..minimization.conjugate_gradient import ConjugateGradient
Jakob Knollmueller's avatar
Jakob Knollmueller committed
21
from ..minimization.quadratic_energy import QuadraticEnergy
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
22
from .endomorphic_operator import EndomorphicOperator
Jakob Knollmueller's avatar
Jakob Knollmueller committed
23 24


Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
25
class SamplingEnabler(EndomorphicOperator):
26
    """Class which acts as a operator object built of (`likelihood` + `prior`)
Martin Reinecke's avatar
grammar  
Martin Reinecke committed
27
    and enables sampling from its inverse even if the operator object
28 29
    itself does not support it.

Jakob Knollmueller's avatar
Jakob Knollmueller committed
30 31 32

    Parameters
    ----------
33
    likelihood : :class:`EndomorphicOperator`
Martin Reinecke's avatar
Martin Reinecke committed
34
        Metric of the likelihood
35
    prior : :class:`EndomorphicOperator`
Martin Reinecke's avatar
Martin Reinecke committed
36
        Inverse metric of the prior
37 38 39
    iteration_controller : :class:`IterationController`
        The iteration controller to use for the iterative numerical inversion
        done by a :class:`ConjugateGradient` object.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
40
    approximation : :class:`LinearOperator`, optional
Martin Reinecke's avatar
Martin Reinecke committed
41 42 43 44
        if not None, this linear operator should be an approximation to the
        operator, which supports the operation modes that the operator doesn't
        have. It is used as a preconditioner during the iterative inversion,
        to accelerate convergence.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
45 46
    """

47
    def __init__(self, likelihood, prior, iteration_controller,
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
48
                 approximation=None):
Martin Reinecke's avatar
Martin Reinecke committed
49
        self._op = likelihood + prior
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
50 51
        self._likelihood = likelihood
        self._prior = prior
52
        self._ic = iteration_controller
53
        self._approximation = approximation
Martin Reinecke's avatar
Martin Reinecke committed
54
        self._domain = self._op.domain
Martin Reinecke's avatar
Martin Reinecke committed
55
        self._capability = self._op.capability
Jakob Knollmueller's avatar
Jakob Knollmueller committed
56 57 58 59 60

    def draw_sample(self, from_inverse=False, dtype=np.float64):
        try:
            return self._op.draw_sample(from_inverse, dtype)
        except NotImplementedError:
Martin Reinecke's avatar
Martin Reinecke committed
61 62
            if not from_inverse:
                raise ValueError("from_inverse must be True here")
Jakob Knollmueller's avatar
bugfix  
Jakob Knollmueller committed
63 64
            s = self._prior.draw_sample(from_inverse=True)
            sp = self._prior(s)
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
65 66 67
            nj = self._likelihood.draw_sample()
            energy = QuadraticEnergy(s, self._op, sp + nj,
                                     _grad=self._likelihood(s) - nj)
68
            inverter = ConjugateGradient(self._ic)
69 70 71 72 73
            if self._approximation is not None:
                energy, convergence = inverter(
                    energy, preconditioner=self._approximation.inverse)
            else:
                energy, convergence = inverter(energy)
Jakob Knollmueller's avatar
Jakob Knollmueller committed
74
            return energy.position
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
75 76 77

    def apply(self, x, mode):
        return self._op.apply(x, mode)
Martin Reinecke's avatar
Martin Reinecke committed
78 79 80 81 82 83 84 85

    def __repr__(self):
        from ..utilities import indent
        return "\n".join((
            "SamplingEnabler:",
            indent("\n".join((
                "Likelihood:", self._likelihood.__repr__(),
                "Prior:", self._prior.__repr__())))))