response_operator.py 2.03 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
from .. import Field, FieldArray, DomainTuple
Martin Reinecke's avatar
Martin Reinecke committed
2
3
4
from .linear_operator import LinearOperator
from .fft_smoothing_operator import FFTSmoothingOperator
from .diagonal_operator import DiagonalOperator
Martin Reinecke's avatar
Martin Reinecke committed
5
6
from .scaling_operator import ScalingOperator
import numpy as np
7

8

Martin Reinecke's avatar
Martin Reinecke committed
9
10
11
class GeometryRemover(LinearOperator):
    def __init__(self, domain):
        super(GeometryRemover, self).__init__()
Martin Reinecke's avatar
Martin Reinecke committed
12
        self._domain = DomainTuple.make(domain)
Martin Reinecke's avatar
Martin Reinecke committed
13
        target_list = [FieldArray(dom.shape) for dom in self._domain]
Martin Reinecke's avatar
Martin Reinecke committed
14
        self._target = DomainTuple.make(target_list)
15

16
17
18
19
20
21
22
23
    @property
    def domain(self):
        return self._domain

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

Martin Reinecke's avatar
Martin Reinecke committed
24
25
26
27
28
29
    @property
    def capability(self):
        return self.TIMES | self.ADJOINT_TIMES

    def apply(self, x, mode):
        self._check_input(x, mode)
Martin Reinecke's avatar
Martin Reinecke committed
30
        if mode == self.TIMES:
Martin Reinecke's avatar
Martin Reinecke committed
31
32
            return Field(self._target, val=x.weight(1).val)
        return Field(self._domain, val=x.val).weight(1)
Martin Reinecke's avatar
Martin Reinecke committed
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64


def ResponseOperator(domain, sigma, exposure):
    domain = DomainTuple.make(domain)
    ncomp = len(exposure)
    if len(sigma) != ncomp or len(domain) != ncomp:
        raise ValueError("length mismatch between sigma, exposure "
                         "and domain")
    ncomp = len(sigma)
    if ncomp == 0:
        raise ValueError("Empty response operator not allowed")

    kernel = None
    expo = None
    for i in range(ncomp):
        if sigma[i] > 0:
            op = FFTSmoothingOperator(domain, sigma[i], space=i)
            kernel = op if kernel is None else op*kernel
        if np.isscalar(exposure[i]):
            if exposure[i] != 1.:
                op = ScalingOperator(exposure[i], domain)
                expo = op if expo is None else op*expo
        elif isinstance(exposure[i], Field):
            op = DiagonalOperator(exposure[i], domain=domain, spaces=i)
            expo = op if expo is None else op*expo

    res = GeometryRemover(domain)
    if expo is not None:
        res = res * expo
    if kernel is not None:
        res = res * kernel
    return res