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)