scaling_operator.py 4.29 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 .endomorphic_operator import EndomorphicOperator
Martin Reinecke's avatar
Martin Reinecke committed
21
22
23


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

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

    Notes
    -----
Philipp Arras's avatar
Docs    
Philipp Arras committed
35
36
    :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
37
38
39
40
41
42
43

    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
44
45
46
47
48
49
    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
50
51
    """

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

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

    def apply(self, x, mode):
62
63
        from ..sugar import full

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

        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
74
            fct = np.conj(fct)
Lukas Platz's avatar
Lukas Platz committed
75
        if (mode & MODES_WITH_INVERSE) != 0:
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
76
77
            fct = 1./fct
        return x*fct
Martin Reinecke's avatar
Martin Reinecke committed
78

79
    def _flip_modes(self, trafo):
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
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
85
        return ScalingOperator(self._domain, fct)
Martin Reinecke's avatar
Martin Reinecke committed
86

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

    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
98

99
    def __call__(self, other):
Vincent Eberle's avatar
Vincent Eberle committed
100
        res = EndomorphicOperator.__call__(self, other)
Martin Reinecke's avatar
more    
Martin Reinecke committed
101
102
103
        from .operator import Operator
        if isinstance(res, Operator):
            return res
104
        if np.isreal(self._factor) and self._factor >= 0:
Martin Reinecke's avatar
more    
Martin Reinecke committed
105
            if other.jac is not None:
106
107
108
109
110
111
                if other.metric is not None:
                    from .sandwich_operator import SandwichOperator
                    sqrt_fac = np.sqrt(self._factor)
                    newop = ScalingOperator(other.metric.domain, sqrt_fac)
                    met = SandwichOperator.make(newop, other.metric)
                    res = res.add_metric(met)
Vincent Eberle's avatar
Vincent Eberle committed
112
        return res
Reimar Leike's avatar
Reimar Leike committed
113

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