response_operator.py 4.01 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

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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
        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


    Attributes
    ----------

    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)
48
49
50
51
52
53
    <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.]])
54
55
56
57
58
59

    See Also
    --------

    """

60
    def __init__(self, domain,
61
62
63
                 sigma=[1.], exposure=[1.],
                 default_spaces=None):
        super(ResponseOperator, self).__init__(default_spaces)
64
65

        self._domain = self._parse_domain(domain)
66
67
68
69
70
71
72

        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)

73
        self._target = self._parse_domain(FieldArray(shape_target))
74
        self._sigma = sigma
75
        self._exposure = exposure
76
77
78
79
80

        self._kernel = len(self._domain)*[None]

        for ii in xrange(len(self._kernel)):
            self._kernel[ii] = SmoothingOperator(self._domain[ii],
Theo Steininger's avatar
Theo Steininger committed
81
                                                 sigma=self._sigma[ii])
82
83

        self._composed_kernel = ComposedOperator(self._kernel)
84

85
        self._exposure_op = len(self._domain)*[None]
Theo Steininger's avatar
Theo Steininger committed
86
87
88
        if len(self._exposure_op) != len(self._kernel):
            raise ValueError("Definition of kernel and exposure do not suit "
                             "each other")
89
90
        else:
            for ii in xrange(len(self._exposure_op)):
Theo Steininger's avatar
Theo Steininger committed
91
92
93
                self._exposure_op[ii] = DiagonalOperator(
                                                self._domain[ii],
                                                diagonal=self._exposure[ii])
94
95
            self._composed_exposure = ComposedOperator(self._exposure_op)

96
97
98
99
100
101
102
103
104
105
    @property
    def domain(self):
        return self._domain

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

    @property
    def unitary(self):
106
        return False
107
108

    def _times(self, x, spaces):
109
110
        res = self._composed_kernel.times(x, spaces)
        res = self._composed_exposure.times(res, spaces)
111
112
113
        # res = res.weight(power=1)
        # removing geometric information
        return Field(self._target, val=res.val)
114
115
116

    def _adjoint_times(self, x, spaces):
        # setting correct spaces
117
118
        res = Field(self.domain, val=x.val)
        res = self._composed_exposure.adjoint_times(res, spaces)
119
        res = res.weight(power=-1)
120
        res = self._composed_kernel.adjoint_times(res, spaces)
121
        return res