sampling_enabler.py 3.82 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`
36
        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.
45
46
47
48
    start_from_zero : boolean
        If true, the conjugate gradient algorithm starts from a field filled
        with zeros. Otherwise, it starts from a prior samples. Default is
        False.
Jakob Knollmueller's avatar
Jakob Knollmueller committed
49
50
    """

51
    def __init__(self, likelihood, prior, iteration_controller,
52
                 approximation=None, start_from_zero=False):
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
53
54
        self._likelihood = likelihood
        self._prior = prior
55
        self._ic = iteration_controller
56
        self._approximation = approximation
57
58
        self._start_from_zero = bool(start_from_zero)
        self._op = likelihood + prior
Martin Reinecke's avatar
Martin Reinecke committed
59
        self._domain = self._op.domain
Martin Reinecke's avatar
Martin Reinecke committed
60
        self._capability = self._op.capability
Jakob Knollmueller's avatar
Jakob Knollmueller committed
61
62
63
64
65

    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
66
67
            if not from_inverse:
                raise ValueError("from_inverse must be True here")
68
69
70
71
72
73
74
75
76
            if self._start_from_zero:
                b = self._op.draw_sample()
                energy = QuadraticEnergy(0*b, self._op, b)
            else:
                s = self._prior.draw_sample(from_inverse=True)
                sp = self._prior(s)
                nj = self._likelihood.draw_sample()
                energy = QuadraticEnergy(s, self._op, sp + nj,
                                         _grad=self._likelihood(s) - nj)
77
            inverter = ConjugateGradient(self._ic)
78
79
80
81
82
            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
83
            return energy.position
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
84
85
86

    def apply(self, x, mode):
        return self._op.apply(x, mode)
Martin Reinecke's avatar
Martin Reinecke committed
87
88
89
90
91
92
93
94

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