response_operator.py 4.54 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
from builtins import range
2
import numpy as np
3 4 5 6 7 8
from ... import Field,\
                FieldArray
from ..linear_operator import LinearOperator
from ..smoothing_operator import SmoothingOperator
from ..composed_operator import ComposedOperator
from ..diagonal_operator import DiagonalOperator
9

10

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

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

    Attributes
    ----------
34 35 36
    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
37
        The domain in which the outcome of the operator lives.
38 39
    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, sigma=[1.], exposure=[1.],
71 72
                 default_spaces=None):
        super(ResponseOperator, self).__init__(default_spaces)
73 74

        self._domain = self._parse_domain(domain)
75

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

79
        if len(sigma) != len(exposure):
80 81 82
            raise ValueError("Length of smoothing kernel and length of"
                             "exposure do not match")

Martin Reinecke's avatar
Martin Reinecke committed
83 84
        for ii in range(len(kernel_smoothing)):
            kernel_smoothing[ii] = SmoothingOperator.make(self._domain[ii],
85
                                                          sigma=sigma[ii])
86 87 88 89 90
            kernel_exposure[ii] = DiagonalOperator(self._domain[ii],
                                                   diagonal=exposure[ii])

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

92 93 94 95 96 97
        target_list = []
        for space in self.domain:
            target_list += [FieldArray(space.shape)]

        self._target = self._parse_domain(target_list)

98 99 100 101 102
    def _add_attributes_to_copy(self, copy, **kwargs):
        copy._domain = self._domain
        copy._target = self._target
        copy._composed_kernel = self._composed_kernel.copy()
        copy._composed_exposure = self._composed_exposure.copy()
Martin Reinecke's avatar
bug fix  
Martin Reinecke committed
103
        copy = super(ResponseOperator, self)._add_attributes_to_copy(copy,
104 105 106
                                                                     **kwargs)
        return copy

107 108 109 110 111 112 113 114 115 116
    @property
    def domain(self):
        return self._domain

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

    @property
    def unitary(self):
117
        return False
118 119

    def _times(self, x, spaces):
120 121
        res = self._composed_kernel.times(x, spaces)
        res = self._composed_exposure.times(res, spaces)
122 123 124
        # res = res.weight(power=1)
        # removing geometric information
        return Field(self._target, val=res.val)
125 126 127

    def _adjoint_times(self, x, spaces):
        # setting correct spaces
128 129
        res = Field(self.domain, val=x.val)
        res = self._composed_exposure.adjoint_times(res, spaces)
130
        res = res.weight(power=-1)
131
        res = self._composed_kernel.adjoint_times(res, spaces)
132
        return res