response_operator.py 3.42 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
from builtins import range
Martin Reinecke's avatar
Martin Reinecke committed
2
from .. import Field, FieldArray, DomainTuple
Martin Reinecke's avatar
Martin Reinecke committed
3
4
5
6
from .linear_operator import LinearOperator
from .fft_smoothing_operator import FFTSmoothingOperator
from .composed_operator import ComposedOperator
from .diagonal_operator import DiagonalOperator
7

8

9
class ResponseOperator(LinearOperator):
Theo Steininger's avatar
Theo Steininger committed
10
    """ NIFTy ResponseOperator (example)
11

12
13
14
15
16
17
18
    This NIFTy ResponseOperator provides the user with an example how a
    ResponseOperator can look like. It smoothes and exposes a field. The
    outcome of the Operator is geometrically not ordered as typical data
    set are.

    Parameters
    ----------
Theo Steininger's avatar
Theo Steininger committed
19
    domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
20
21
22
23
24
25
        The domains on which the operator lives. Either one space or a list
        of spaces
    sigma : list(np.float)
        Defines the smoothing length of the operator for each space it lives on
    exposure : list(np.float)
        Defines the exposure of the operator for each space it lives on
26
27
28
    default_spaces : tuple of ints *optional*
        Defines on which space(s) of a given field the Operator acts by
        default (default: None)
29
30
31

    Attributes
    ----------
32
33
34
    domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
        The domain on which the Operator's input Field lives.
    target : tuple of DomainObjects, i.e. Spaces and FieldTypes
Martin Reinecke's avatar
Martin Reinecke committed
35
        The domain in which the outcome of the operator lives.
36
37
    unitary : boolean
        Indicates whether the Operator is unitary or not.
38
39
40
41
42
43
44
45

    Raises
    ------
    ValueError:
        raised if:
            * len of sigma-list and exposure-list are not equal
    """

46
47
    def __init__(self, domain, sigma=[1.], exposure=[1.], spaces=None):
        super(ResponseOperator, self).__init__()
48

49
        if len(sigma) != len(exposure):
50
51
            raise ValueError("Length of smoothing kernel and length of"
                             "exposure do not match")
Martin Reinecke's avatar
updates    
Martin Reinecke committed
52
        nsigma = len(sigma)
53

Martin Reinecke's avatar
Martin Reinecke committed
54
        self._domain = DomainTuple.make(domain)
Martin Reinecke's avatar
updates    
Martin Reinecke committed
55

56
57
58
59
60
        if spaces is None:
            spaces = range(len(self._domain))

        kernel_smoothing = [FFTSmoothingOperator(self._domain, sigma[x],
                                                 space=spaces[x])
Martin Reinecke's avatar
updates    
Martin Reinecke committed
61
                            for x in range(nsigma)]
62
        kernel_exposure = [DiagonalOperator(Field(self._domain[spaces[x]],
63
                                                  exposure[x]),
64
65
                                            domain=self._domain,
                                            spaces=(spaces[x],))
Martin Reinecke's avatar
Martin Reinecke committed
66
                           for x in range(nsigma)]
67
68
69

        self._composed_kernel = ComposedOperator(kernel_smoothing)
        self._composed_exposure = ComposedOperator(kernel_exposure)
70

71
        target_list = [FieldArray(self._domain[i].shape) for i in spaces]
Martin Reinecke's avatar
Martin Reinecke committed
72
        self._target = DomainTuple.make(target_list)
73

74
75
76
77
78
79
80
81
82
83
    @property
    def domain(self):
        return self._domain

    @property
    def target(self):
        return self._target

    @property
    def unitary(self):
84
        return False
85

86
87
88
    def _times(self, x):
        res = self._composed_kernel.times(x)
        res = self._composed_exposure.times(res)
89
90
        # removing geometric information
        return Field(self._target, val=res.val)
91

92
    def _adjoint_times(self, x):
93
        # setting correct spaces
94
        res = Field(self.domain, val=x.val)
95
        res = self._composed_exposure.adjoint_times(res)
96
        res = res.weight(power=-1)
97
        return self._composed_kernel.adjoint_times(res)