scaling_operator.py 3.87 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-2018 Max-Planck-Society
Martin Reinecke's avatar
Martin Reinecke committed
15
16
17
18
19
20
21
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.

from __future__ import division
import numpy as np
from ..field import Field
22
from ..multi.multi_field import MultiField
Martin Reinecke's avatar
Martin Reinecke committed
23
from .endomorphic_operator import EndomorphicOperator
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
24
from ..domain_tuple import DomainTuple
Martin Reinecke's avatar
Martin Reinecke committed
25
26
27


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

Martin Reinecke's avatar
Martin Reinecke committed
30
    The NIFTy ScalingOperator class is a subclass derived from the
Martin Reinecke's avatar
Martin Reinecke committed
31
32
33
34
35
36
    EndomorphicOperator. It multiplies an input field with a given factor.

    Parameters
    ----------
    factor : scalar
        The multiplication factor
Martin Reinecke's avatar
docs    
Martin Reinecke committed
37
    domain : Domain or tuple of Domain or DomainTuple
Martin Reinecke's avatar
Martin Reinecke committed
38
        The domain on which the Operator's input Field lives.
Martin Reinecke's avatar
Martin Reinecke committed
39
40
41
42
43
44
45
46
47
48

    Notes
    -----
    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).

    This shortcoming will hopefully be fixed in the future.
Martin Reinecke's avatar
Martin Reinecke committed
49
50
51
    """

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

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

    def apply(self, x, mode):
        self._check_input(x, mode)

Martin Reinecke's avatar
Martin Reinecke committed
63
64
65
        if self._factor == 1.:
            return x.copy()
        if self._factor == 0.:
66
            return x.empty_copy().fill(0.)
Martin Reinecke's avatar
Martin Reinecke committed
67

Martin Reinecke's avatar
Martin Reinecke committed
68
69
70
71
72
73
74
75
76
        if mode == self.TIMES:
            return x*self._factor
        elif mode == self.ADJOINT_TIMES:
            return x*np.conj(self._factor)
        elif mode == self.INVERSE_TIMES:
            return x*(1./self._factor)
        else:
            return x*(1./np.conj(self._factor))

77
78
79
80
81
    def _flip_modes(self, trafo):
        ADJ = self.ADJOINT_BIT
        INV = self.INVERSE_BIT

        if trafo == 0:
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
82
            return self
83
        if trafo == ADJ and np.issubdtype(type(self._factor), np.floating):
Martin Reinecke's avatar
Martin Reinecke committed
84
            return self
85
        if trafo == ADJ:
Martin Reinecke's avatar
Martin Reinecke committed
86
            return ScalingOperator(np.conj(self._factor), self._domain)
87
        elif trafo == INV:
Martin Reinecke's avatar
Martin Reinecke committed
88
            return ScalingOperator(1./self._factor, self._domain)
89
        elif trafo == ADJ | INV:
Martin Reinecke's avatar
Martin Reinecke committed
90
            return ScalingOperator(1./np.conj(self._factor), self._domain)
91
        raise ValueError("invalid operator transformation")
Martin Reinecke's avatar
Martin Reinecke committed
92

Martin Reinecke's avatar
Martin Reinecke committed
93
94
95
96
97
98
    @property
    def domain(self):
        return self._domain

    @property
    def capability(self):
Martin Reinecke's avatar
Martin Reinecke committed
99
        return self._all_ops
100

101
102
    def draw_sample(self, from_inverse=False, dtype=np.float64):
        fct = self._factor
103
104
105
        if fct.imag != 0. or fct.real < 0.:
            raise ValueError("operator not positive definite")
        if fct.real == 0. and from_inverse:
clienhar's avatar
clienhar committed
106
            raise ValueError("operator not positive definite")
107
        fct = 1./np.sqrt(fct) if from_inverse else np.sqrt(fct)
108
109
        cls = Field if isinstance(self._domain, DomainTuple) else MultiField
        return cls.from_random(
Martin Reinecke's avatar
Martin Reinecke committed
110
           random_type="normal", domain=self._domain, std=fct, dtype=dtype)