Commit b33dac3d authored by Martin Reinecke's avatar Martin Reinecke

rewrite ResponseOperator

parent d9d0e0b4
Pipeline #23389 passed with stage
in 4 minutes and 36 seconds
from builtins import range
from .. import Field, FieldArray, DomainTuple
from .linear_operator import LinearOperator
from .fft_smoothing_operator import FFTSmoothingOperator
......@@ -7,58 +6,11 @@ from .scaling_operator import ScalingOperator
import numpy as np
class ResponseOperator(LinearOperator):
""" NIFTy ResponseOperator (example)
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
----------
domain : tuple of DomainObjects
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
spaces : tuple of ints *optional*
Defines on which subdomain the individual components act
(default: None)
"""
def __init__(self, domain, sigma, exposure, spaces=None):
super(ResponseOperator, self).__init__()
class GeometryRemover(LinearOperator):
def __init__(self, domain):
super(GeometryRemover, self).__init__()
self._domain = DomainTuple.make(domain)
if spaces is None:
spaces = tuple(range(len(self._domain)))
if len(sigma) != len(exposure) or len(sigma) != len(spaces):
raise ValueError("length mismatch between sigma, exposure, "
"and spaces ")
ncomp = len(sigma)
if ncomp == 0:
raise ValueError("Empty response operator not allowed")
self._kernel = None
for i in range(ncomp):
op = FFTSmoothingOperator(self._domain, sigma[i], space=spaces[i])
self._kernel = op if self._kernel is None else op*self._kernel
self._exposure = None
for i in range(ncomp):
if np.isscalar(exposure[i]):
op = ScalingOperator(exposure[i], self._domain)
elif isinstance(exposure[i], Field):
op = DiagonalOperator(exposure[i], domain=self._domain,
spaces=(spaces[i],))
self._exposure = op if self._exposure is None\
else op*self._exposure
target_list = [FieldArray(self._domain[i].shape) for i in spaces]
target_list = [FieldArray(dom.shape) for dom in self._domain]
self._target = DomainTuple.make(target_list)
@property
......@@ -73,18 +25,40 @@ class ResponseOperator(LinearOperator):
def capability(self):
return self.TIMES | self.ADJOINT_TIMES
def _times(self, x):
res = self._exposure(self._kernel(x))
# removing geometric information
return Field(self._target, val=res.val)
def _adjoint_times(self, x):
# setting correct spaces
res = Field(self.domain, val=x.val)
res = self._exposure.adjoint_times(res)
res = res.weight(power=-1)
return self._kernel.adjoint_times(res)
def apply(self, x, mode):
self._check_input(x, mode)
return self._times(x) if mode == self.TIMES else self._adjoint_times(x)
if mode == self.TIMES:
return Field(self._target, val=x.val)
return Field(self._domain, val=x.val).weight(power=-1)
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
......@@ -6,7 +6,7 @@ from test.common import expand
class ResponseOperator_Tests(unittest.TestCase):
spaces = [ift.RGSpace(128)]
spaces = [ift.RGSpace(128), ift.GLSpace(nlat=37)]
@expand(product(spaces, [0., 5., 1.], [0., 1., .33]))
def test_property(self, space, sigma, exposure):
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment