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
41
42
        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
Jakob Knollmueller's avatar
Jakob Knollmueller committed
43
44
45
46
        preconditioner during the iterative inversion, to accelerate
        convergence.
    """

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__())))))