response_operator.py 2.42 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
11

class ResponseOperator(LinearOperator):

    def __init__(self, domain,
12
                 sigma=[1.], exposure=[1.],
13
14
15
                 unitary=False):

        self._domain = self._parse_domain(domain)
16
17
18
19
20
21
22

        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)

23
        self._target = self._parse_domain(FieldArray(shape_target))
24
        self._sigma = sigma
25
        self._exposure = exposure
26
27
        self._unitary = unitary

28
29
30
31
32
33
34
35

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

        for ii in xrange(len(self._kernel)):
            self._kernel[ii] = SmoothingOperator(self._domain[ii],
                                        sigma=self._sigma[ii])

        self._composed_kernel = ComposedOperator(self._kernel)
36

37
38
39
40
41
42
43
44
45
46

        self._exposure_op = len(self._domain)*[None]
        if len(self._exposure_op)!= len(self._kernel):
            raise ValueError("Definition of kernel and exposure do not suit each other")
        else:
            for ii in xrange(len(self._exposure_op)):
                self._exposure_op[ii] = DiagonalOperator(self._domain[ii],
                                                      diagonal=self._exposure[ii])
            self._composed_exposure = ComposedOperator(self._exposure_op)

47
48
49
50
51
52
53
54
55
56
57
58
59
60

    @property
    def domain(self):
        return self._domain

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

    @property
    def unitary(self):
        return self._unitary

    def _times(self, x, spaces):
61
62
        res = self._composed_kernel.times(x, spaces)
        res = self._composed_exposure.times(res, spaces)
63
64
65
        # res = res.weight(power=1)
        # removing geometric information
        return Field(self._target, val=res.val)
66
67
68

    def _adjoint_times(self, x, spaces):
        # setting correct spaces
69
70
        res = Field(self.domain, val=x.val)
        res = self._composed_exposure.adjoint_times(res, spaces)
71
        res = res.weight(power=-1)
72
        res = self._composed_kernel.adjoint_times(res, spaces)
73
        return res