Skip to content

GitLab

Commit af9bdeb9 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

remove ComposedOperator

parent 32a08573
Pipeline #23263 failed with stage
in 4 minutes and 1 second
Hide whitespace changes
Inline Side-by-side
demos/critical_filtering.py
View file @ af9bdeb9
......@@ -34,7 +34,7 @@ if __name__ == "__main__":
# Instrument._diagonal.val[64:512-64, 64:512-64] = 0
# Add a harmonic transformation to the instrument
R = ift.ComposedOperator([fft, Instrument])
R = Instrument*fft
noise = 1.
N = ift.DiagonalOperator(ift.Field.full(s_space, noise).weight(1))
......
demos/log_normal_wiener_filter.py
View file @ af9bdeb9
......@@ -36,7 +36,7 @@ if __name__ == "__main__":
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
exposure=(mask,))
data_domain = R.target[0]
R_harmonic = ift.ComposedOperator([fft, R])
R_harmonic = R*fft
# Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise)
......@@ -66,20 +66,16 @@ if __name__ == "__main__":
# m3 = fft(me3[0].position)
# Plotting
ift.plotting.plot(mock_signal, name='mock_signal.png', colormap="plasma",
xlabel="Pixel Index", ylabel="Pixel Index")
plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
"colormap": "Planck-like"}
ift.plotting.plot(mock_signal, name="mock_signal.png", **plotdict)
logdata = np.log(ift.dobj.to_global_data(data.val)).reshape(signal_space.shape)
ift.plotting.plot(ift.Field(signal_space,
val=ift.dobj.from_global_data(logdata)),
name="log_of_data.png", colormap="plasma",
xlabel="Pixel Index", ylabel="Pixel Index")
# ift.plotting.plot(m1,name='m_LBFGS.png', colormap="plasma",
# xlabel="Pixel Index", ylabel="Pixel Index")
ift.plotting.plot(m2, name='m_Newton.png', colormap="plasma",
xlabel="Pixel Index", ylabel="Pixel Index")
# ift.plotting.plot(m3, name='m_SteepestDescent.png',
# colormap="plasma", xlabel="Pixel Index",
# ylabel="Pixel Index")
name="log_of_data.png", **plotdict)
# ift.plotting.plot(m1,name='m_LBFGS.png', **plotdict)
ift.plotting.plot(m2, name='m_Newton.png', **plotdict)
# ift.plotting.plot(m3, name='m_SteepestDescent.png', **plotdict)
# Probing the variance
class Proby(ift.DiagonalProberMixin, ift.Prober):
......@@ -89,4 +85,4 @@ if __name__ == "__main__":
sm = ift.FFTSmoothingOperator(signal_space, sigma=0.02)
variance = sm(proby.diagonal.weight(-1))
ift.plotting.plot(variance, name='variance.png')
ift.plotting.plot(variance, name='variance.png', **plotdict)
demos/nonlinear_critical_filter.py
View file @ af9bdeb9
......@@ -50,8 +50,8 @@ if __name__ == "__main__":
MaskOperator = ift.DiagonalOperator(mask)
InstrumentResponse = ift.ResponseOperator(s_space, sigma=[0.0])
MeasurementOperator = ift.ComposedOperator([MaskOperator,
InstrumentResponse])
MeasurementOperator = InstrumentResponse*MaskOperator
d_space = MeasurementOperator.target
noise_covariance = ift.Field(d_space, val=noise_level**2).weight()
......
demos/paper_demos/cartesian_wiener_filter.py
View file @ af9bdeb9
......@@ -44,7 +44,7 @@ if __name__ == "__main__":
mock_power_2 = ift.PS_field(power_space_2, power_spectrum_2)
fft = ift.ComposedOperator((fft_1, fft_2))
fft = fft_2*fft_1
mock_power = ift.Field(
(power_space_1, power_space_2),
......@@ -74,7 +74,7 @@ if __name__ == "__main__":
sigma=(response_sigma_1, response_sigma_2),
exposure=(mask_1, mask_2))
data_domain = R.target
R_harmonic = ift.ComposedOperator([fft, R])
R_harmonic = R*fft
# Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise)
......@@ -105,13 +105,15 @@ if __name__ == "__main__":
plot_space = ift.RGSpace((N_pixels_1, N_pixels_2))
sm = ift.FFTSmoothingOperator(plot_space, sigma=0.03)
plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
"colormap": "Planck-like"}
ift.plotting.plot(
ift.log(ift.sqrt(sm(ift.Field(plot_space, val=variance.val.real)))),
name='uncertainty.png', zmin=0., zmax=3., title="Uncertainty map",
colormap="Planck-like")
**plotdict)
ift.plotting.plot(ift.Field(plot_space, val=mock_signal.val.real),
name='mock_signal.png', colormap="Planck-like")
name='mock_signal.png', **plotdict)
ift.plotting.plot(ift.Field(plot_space, val=data.val.real),
name='data.png', colormap="Planck-like")
name='data.png', **plotdict)
ift.plotting.plot(ift.Field(plot_space, val=m.val.real),
name='map.png', colormap="Planck-like")
name='map.png', **plotdict)
demos/paper_demos/wiener_filter.py
View file @ af9bdeb9
......@@ -38,7 +38,7 @@ if __name__ == "__main__":
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
exposure=(mask,))
data_domain = R.target[0]
R_harmonic = ift.ComposedOperator([fft, R])
R_harmonic = R * fft
# Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise)
......@@ -67,12 +67,10 @@ if __name__ == "__main__":
variance = ift.sqrt(sm(proby.diagonal.weight(-1)))
# Plotting
ift.plotting.plot(variance, name="uncertainty.png", xlabel='Pixel index',
ylabel='Pixel index')
ift.plotting.plot(mock_signal, name="mock_signal.png",
xlabel='Pixel index', ylabel='Pixel index')
plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
"colormap": "Planck-like"}
ift.plotting.plot(variance, name="uncertainty.png", **plotdict)
ift.plotting.plot(mock_signal, name="mock_signal.png", **plotdict)
ift.plotting.plot(ift.Field(signal_space, val=data.val),
name="data.png", xlabel='Pixel index',
ylabel='Pixel index')
ift.plotting.plot(m, name="map.png", xlabel='Pixel index',
ylabel='Pixel index')
name="data.png", **plotdict)
ift.plotting.plot(m, name="map.png", **plotdict)
demos/wiener_filter_via_curvature.py
View file @ af9bdeb9
......@@ -53,7 +53,7 @@ if __name__ == "__main__":
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
exposure=(exposure,))
data_domain = R.target[0]
R_harmonic = ift.ComposedOperator([fft, R])
R_harmonic = R*fft
N = ift.DiagonalOperator(
ift.Field.full(data_domain,
......
demos/wiener_filter_via_hamiltonian.py
View file @ af9bdeb9
......@@ -18,7 +18,7 @@ if __name__ == "__main__":
# Choosing the prior correlation structure and defining
# correlation operator
p_spec = (lambda k: (42 / (k + 1) ** 3))
p_spec = (lambda k: (42/(k+1)**3))
S = ift.create_power_operator(h_space, power_spectrum=p_spec)
......@@ -35,7 +35,7 @@ if __name__ == "__main__":
Instrument = ift.DiagonalOperator(diag)
# Adding a harmonic transformation to the instrument
R = ift.ComposedOperator([fft, Instrument])
R = Instrument*fft
signal_to_noise = 1.
ndiag = ift.Field.full(s_space, ss.var()/signal_to_noise)
N = ift.DiagonalOperator(ndiag.weight(1))
......
nifty/library/wiener_filter_curvature.py
View file @ af9bdeb9
......@@ -31,7 +31,7 @@ class WienerFilterCurvature(EndomorphicOperator):
@property
def domain(self):
return self.S.domain
return self._op.domain
@property
def capability(self):
......
nifty/operators/__init__.py
View file @ af9bdeb9
......@@ -6,7 +6,6 @@ from .fft_smoothing_operator import FFTSmoothingOperator
from .direct_smoothing_operator import DirectSmoothingOperator
from .fft_operator import FFTOperator
from .inversion_enabler import InversionEnabler
from .composed_operator import ComposedOperator
from .response_operator import ResponseOperator
from .laplace_operator import LaplaceOperator
from .smoothness_operator import SmoothnessOperator
......
nifty/operators/composed_operator.py deleted 100644 → 0
View file @ 32a08573
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from .linear_operator import LinearOperator
class ComposedOperator(LinearOperator):
""" NIFTY class for composed operators.
The NIFTY composed operator class combines multiple linear operators.
Parameters
----------
operators : tuple of NIFTy Operators
The tuple of LinearOperators.
Attributes
----------
domain : DomainTuple
The NIFTy.space in which the operator is defined.
target : DomainTuple
The NIFTy.space in which the outcome of the operator lives
"""
def __init__(self, operators):
super(ComposedOperator, self).__init__()
self._operator_store = ()
old_op = None
self._capability = operators[0].capability
for op in operators:
if not isinstance(op, LinearOperator):
raise TypeError("The elements of the operator list must be"
"instances of the LinearOperator base class")
if old_op is not None and op.domain != old_op.target:
raise ValueError("incompatible domains")
self._operator_store += (op,)
self._capability &= op.capability
old_op = op
if self._capability == 0:
raise ValueError("composed operator does not support any mode")
@property
def domain(self):
return self._operator_store[0].domain
@property
def target(self):
return self._operator_store[-1].target
@property
def capability(self):
return self._capability
def apply(self, x, mode):
self._check_mode(mode)
if mode == self.TIMES or mode == self.ADJOINT_INVERSE_TIMES:
for op in self._operator_store:
x = op.apply(x, mode)
else:
for op in reversed(self._operator_store):
x = op.apply(x, mode)
return x
nifty/operators/fft_operator.py
View file @ af9bdeb9
......@@ -54,20 +54,6 @@ class FFTOperator(LinearOperator):
For GLSpace, HPSpace, and LMSpace, a sensible (but not unique)
co-domain is chosen that should work satisfactorily in most situations,
but for full control, the user should explicitly specify a codomain.
Attributes
----------
domain: Tuple of Spaces
The domain of the data that is input by "times" and output by
"adjoint_times".
target: Tuple of Spaces
The domain of the data that is output by "times" and input by
"adjoint_times".
Raises
------
ValueError:
if "domain" or "target" are not of the proper type.
"""
def __init__(self, domain, target=None, space=None):
......@@ -96,7 +82,8 @@ class FFTOperator(LinearOperator):
else:
self._trafo = SphericalTransformation(pdom, hdom, self._space)
def _times_helper(self, x):
def apply(self, x, mode):
self._check_input(x, mode)
if np.issubdtype(x.dtype, np.complexfloating):
res = (self._trafo.transform(x.real) +
1j * self._trafo.transform(x.imag))
......@@ -104,10 +91,6 @@ class FFTOperator(LinearOperator):
res = self._trafo.transform(x)
return res
def apply(self, x, mode):
self._check_input(x, mode)
return self._times_helper(x)
@property
def domain(self):
return self._domain
......
nifty/operators/response_operator.py
View file @ af9bdeb9
......@@ -2,8 +2,9 @@ from builtins import range
from .. import Field, FieldArray, DomainTuple
from .linear_operator import LinearOperator
from .fft_smoothing_operator import FFTSmoothingOperator
from .composed_operator import ComposedOperator
from .diagonal_operator import DiagonalOperator
from .scaling_operator import ScalingOperator
import numpy as np
class ResponseOperator(LinearOperator):
......@@ -23,47 +24,39 @@ class ResponseOperator(LinearOperator):
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
default_spaces : tuple of ints *optional*
Defines on which space(s) of a given field the Operator acts by
default (default: None)
Attributes
----------
domain : DomainTuple
The domain on which the Operator's input Field lives.
target : DomainTuple
The domain in which the outcome of the operator lives.
Raises
------
ValueError:
raised if:
* len of sigma-list and exposure-list are not equal
spaces : tuple of ints *optional*
Defines on which subdomain the individual components act
(default: None)
"""
def __init__(self, domain, sigma=[1.], exposure=[1.], spaces=None):
def __init__(self, domain, sigma, exposure, spaces=None):
super(ResponseOperator, self).__init__()
if len(sigma) != len(exposure):
raise ValueError("Length of smoothing kernel and length of"
"exposure do not match")
nsigma = len(sigma)
self._domain = DomainTuple.make(domain)
if spaces is None:
spaces = range(len(self._domain))
kernel_smoothing = [
FFTSmoothingOperator(self._domain, sigma[x], space=spaces[x])
for x in range(nsigma)]
self._composed_kernel = ComposedOperator(kernel_smoothing)
kernel_exposure = [
DiagonalOperator(Field(self._domain[spaces[x]], exposure[x]),
domain=self._domain, spaces=(spaces[x],))
for x in range(nsigma)]
self._composed_exposure = ComposedOperator(kernel_exposure)
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]
self._target = DomainTuple.make(target_list)
......@@ -81,19 +74,17 @@ class ResponseOperator(LinearOperator):
return self.TIMES | self.ADJOINT_TIMES
def _times(self, x):
res = self._composed_kernel.times(x)
res = self._composed_exposure.times(res)
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._composed_exposure.adjoint_times(res)
res = self._exposure.adjoint_times(res)
res = res.weight(power=-1)
return self._composed_kernel.adjoint_times(res)
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)
nifty/operators/smoothness_operator.py
View file @ af9bdeb9
......@@ -38,6 +38,7 @@ class SmoothnessOperator(EndomorphicOperator):
def domain(self):
return self._laplace._domain
# MR FIXME: shouldn't this operator actually be self-adjoint?
@property
def capability(self):
return self.TIMES
......
nifty/sugar.py
View file @ af9bdeb9
......@@ -17,7 +17,7 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
import numpy as np
from . import Space, PowerSpace, Field, ComposedOperator, DiagonalOperator,\
from . import Space, PowerSpace, Field, DiagonalOperator,\
PowerProjectionOperator, FFTOperator, sqrt, DomainTuple, dobj,\
utilities
......@@ -241,11 +241,10 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None):
def create_composed_fft_operator(domain, codomain=None, all_to='other'):
fft_op_list = []
if codomain is None:
codomain = [None]*len(domain)
tdom = list(domain)
res = None
for i, space in enumerate(domain):
if not isinstance(space, Space):
continue
......@@ -256,6 +255,9 @@ def create_composed_fft_operator(domain, codomain=None, all_to='other'):
tgt = domain[i].get_default_codomain()
else:
tgt = codomain[i]
fft_op_list += [FFTOperator(domain=tdom, target=tgt, space=i)]
op = FFTOperator(domain=tdom, target=tgt, space=i)
res = op if res is None else op*res
tdom[i] = tgt
return ComposedOperator(fft_op_list)
if res is None:
raise ValueError("empty operator")
return res
test/test_operators/test_composed_operator.py
View file @ af9bdeb9
......@@ -17,7 +17,7 @@ class ComposedOperator_Tests(unittest.TestCase):
op1 = ift.DiagonalOperator(diag1, cspace, spaces=(0,))
op2 = ift.DiagonalOperator(diag2, cspace, spaces=(1,))
op = ift.ComposedOperator((op1, op2))
op = op2*op1
rand1 = ift.Field.from_random('normal', domain=(space1, space2))
rand2 = ift.Field.from_random('normal', domain=(space1, space2))
......@@ -34,7 +34,7 @@ class ComposedOperator_Tests(unittest.TestCase):
op1 = ift.DiagonalOperator(diag1, cspace, spaces=(0,))
op2 = ift.DiagonalOperator(diag2, cspace, spaces=(1,))
op = ift.ComposedOperator((op1, op2))
op = op2*op1
rand1 = ift.Field.from_random('normal', domain=(space1, space2))
tt1 = op.inverse_times(op.times(rand1))
......
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