scaling_operator.py 3.74 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke 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
Martin Reinecke's avatar
Martin Reinecke committed
15
#
16
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
17

Martin Reinecke's avatar
Martin Reinecke committed
18
import numpy as np
19

20
from ..sugar import full
21
from .endomorphic_operator import EndomorphicOperator
Martin Reinecke's avatar
Martin Reinecke committed
22 23 24


class ScalingOperator(EndomorphicOperator):
Martin Reinecke's avatar
Martin Reinecke committed
25
    """Operator which multiplies a Field with a scalar.
Martin Reinecke's avatar
Martin Reinecke committed
26 27 28 29 30

    Parameters
    ----------
    factor : scalar
        The multiplication factor
Martin Reinecke's avatar
docs  
Martin Reinecke committed
31
    domain : Domain or tuple of Domain or DomainTuple
Philipp Arras's avatar
Philipp Arras committed
32
        The domain on which the Operator's input Field is defined.
Martin Reinecke's avatar
Martin Reinecke committed
33 34 35

    Notes
    -----
Philipp Arras's avatar
Docs  
Philipp Arras committed
36 37
    :class:`Operator` supports the multiplication with a scalar. So one does
    not need instantiate :class:`ScalingOperator` explicitly in most cases.
Philipp Arras's avatar
Philipp Arras committed
38 39 40 41 42 43 44

    Formally, this operator always supports all operation modes (times,
    adjoint_times, inverse_times and inverse_adjoint_times), even if `factor`
    is 0 or infinity. It is the user's responsibility to apply the operator
    only in appropriate ways (e.g. call inverse_times only if `factor` is
    nonzero).

Philipp Arras's avatar
Philipp Arras committed
45 46 47 48 49 50
    Along with this behaviour comes the feature that it is possible to draw an
    inverse sample from a :class:`ScalingOperator` (which is a zero-field).
    This occurs if one draws an inverse sample of a positive definite sum of
    two operators each of which are only positive semi-definite. However, it
    is unclear whether this beviour does not lead to unwanted effects
    somewhere else.
Martin Reinecke's avatar
Martin Reinecke committed
51 52 53
    """

    def __init__(self, factor, domain):
Martin Reinecke's avatar
fix  
Martin Reinecke committed
54
        from ..sugar import makeDomain
Martin Reinecke's avatar
Martin Reinecke committed
55 56 57 58

        if not np.isscalar(factor):
            raise TypeError("Scalar required")
        self._factor = factor
Martin Reinecke's avatar
fix  
Martin Reinecke committed
59
        self._domain = makeDomain(domain)
Philipp Arras's avatar
Philipp Arras committed
60
        self._capability = self._all_ops
Martin Reinecke's avatar
Martin Reinecke committed
61 62 63

    def apply(self, x, mode):
        self._check_input(x, mode)
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
64 65
        fct = self._factor
        if fct == 1.:
66
            return x
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
67
        if fct == 0.:
68
            return full(self.domain, 0.)
Lukas Platz's avatar
Lukas Platz committed
69 70 71 72

        MODES_WITH_ADJOINT = self.ADJOINT_TIMES | self.ADJOINT_INVERSE_TIMES
        MODES_WITH_INVERSE = self.INVERSE_TIMES | self.ADJOINT_INVERSE_TIMES
        if (mode & MODES_WITH_ADJOINT) != 0:
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
73
            fct = np.conj(fct)
Lukas Platz's avatar
Lukas Platz committed
74
        if (mode & MODES_WITH_INVERSE) != 0:
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
75 76
            fct = 1./fct
        return x*fct
Martin Reinecke's avatar
Martin Reinecke committed
77

78
    def _flip_modes(self, trafo):
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
79 80 81 82 83 84
        fct = self._factor
        if trafo & self.ADJOINT_BIT:
            fct = np.conj(fct)
        if trafo & self.INVERSE_BIT:
            fct = 1./fct
        return ScalingOperator(fct, self._domain)
Martin Reinecke's avatar
Martin Reinecke committed
85

86
    def _get_fct(self, from_inverse):
87
        fct = self._factor
88 89
        if (fct.imag != 0. or fct.real < 0. or
                (fct.real == 0. and from_inverse)):
Martin Reinecke's avatar
Martin Reinecke committed
90
            raise ValueError("operator not positive definite")
91 92
        return 1./np.sqrt(fct) if from_inverse else np.sqrt(fct)

Martin Reinecke's avatar
Martin Reinecke committed
93 94
#     def process_sample(self, samp, from_inverse):
#         return samp*self._get_fct(from_inverse)
95 96 97 98 99

    def draw_sample(self, from_inverse=False, dtype=np.float64):
        from ..sugar import from_random
        return from_random(random_type="normal", domain=self._domain,
                           std=self._get_fct(from_inverse), dtype=dtype)
Martin Reinecke's avatar
Martin Reinecke committed
100 101

    def __repr__(self):
Martin Reinecke's avatar
Martin Reinecke committed
102
        return "ScalingOperator ({})".format(self._factor)