response_operator.py 4.36 KB
Newer Older
1
import numpy as np
2
3
4
5
6
from nifty import Field,\
                  FieldArray
from nifty.operators.linear_operator import LinearOperator
from nifty.operators.smoothing_operator import SmoothingOperator
from nifty.operators.composed_operator import ComposedOperator
7
from nifty.operators.diagonal_operator import DiagonalOperator
8

9

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

13
14
15
16
17
18
19
    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
20
    domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
21
22
23
24
25
26
        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
27
28
29
    default_spaces : tuple of ints *optional*
        Defines on which space(s) of a given field the Operator acts by
        default (default: None)
30
31
32

    Attributes
    ----------
33
34
35
36
37
38
39
    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
        The domain in which the outcome of the operator lives. As the Operator
        is endomorphic this is the same as its domain.
    unitary : boolean
        Indicates whether the Operator is unitary or not.
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

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

    Notes
    -----

    Examples
    --------
    >>> x1 = RGSpace(5)
    >>> x2 = RGSpace(10)
    >>> R = ResponseOperator(domain=(x1,x2), sigma=[.5, .25],
                             exposure=[2.,3.])
    >>> f = Field((x1,x2), val=4.)
    >>> R.times(f)
58
59
60
61
62
63
    <distributed_data_object>
    array([[ 24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.],
           [ 24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.],
           [ 24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.],
           [ 24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.],
           [ 24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.,  24.]])
64
65
66
67
68
69

    See Also
    --------

    """

70
    def __init__(self, domain,
71
72
73
                 sigma=[1.], exposure=[1.],
                 default_spaces=None):
        super(ResponseOperator, self).__init__(default_spaces)
74
75

        self._domain = self._parse_domain(domain)
76
77
78
79
80
81
82

        shapes = len(self._domain)*[None]
        shape_target = []
        for ii in xrange(len(shapes)):
            shapes[ii] = self._domain[ii].shape
            shape_target = np.append(shape_target, self._domain[ii].shape)

83
        self._target = self._parse_domain(FieldArray(shape_target))
84
85
86
87

        kernel_smoothing = len(self._domain)*[None]
        kernel_exposure = len(self._domain)*[None]

88
        if len(sigma) != len(exposure):
89
90
91
92
93
94
95
96
97
98
99
            raise ValueError("Length of smoothing kernel and length of"
                             "exposure do not match")

        for ii in xrange(len(kernel_smoothing)):
            kernel_smoothing[ii] = SmoothingOperator(self._domain[ii],
                                                     sigma=sigma[ii])
            kernel_exposure[ii] = DiagonalOperator(self._domain[ii],
                                                   diagonal=exposure[ii])

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

101
102
103
104
105
106
107
108
109
110
    @property
    def domain(self):
        return self._domain

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

    @property
    def unitary(self):
111
        return False
112
113

    def _times(self, x, spaces):
114
115
        res = self._composed_kernel.times(x, spaces)
        res = self._composed_exposure.times(res, spaces)
116
117
118
        # res = res.weight(power=1)
        # removing geometric information
        return Field(self._target, val=res.val)
119
120
121

    def _adjoint_times(self, x, spaces):
        # setting correct spaces
122
123
        res = Field(self.domain, val=x.val)
        res = self._composed_exposure.adjoint_times(res, spaces)
124
        res = res.weight(power=-1)
125
        res = self._composed_kernel.adjoint_times(res, spaces)
126
        return res