diff --git a/demos/log_normal_wiener_filter.py b/demos/log_normal_wiener_filter.py
index 5a4eabc8906ed862686ca628d64854e4543ceea4..f5bcccf64e4bd7b2c8c00f493d4c7452b58398b0 100644
--- a/demos/log_normal_wiener_filter.py
+++ b/demos/log_normal_wiener_filter.py
@@ -37,8 +37,6 @@ if __name__ == "__main__":
R = R*ift.DiagonalOperator(mask)
R = R*ht
R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma)
- #R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
- # sensitivity=(mask,))
data_domain = R.target[0]
# Setting up the noise covariance and drawing a random noise realization
diff --git a/demos/nonlinear_critical_filter.py b/demos/nonlinear_critical_filter.py
index 9a67ca09ab4ff326f14e3407e63933bb8da889da..501f6f0ef1f059114eca0b1c415ebf38eee30ea9 100644
--- a/demos/nonlinear_critical_filter.py
+++ b/demos/nonlinear_critical_filter.py
@@ -24,10 +24,10 @@ if __name__ == "__main__":
# Set up position space
# s_space = ift.RGSpace([1024])
s_space = ift.HPSpace(32)
+ h_space = s_space.get_default_codomain()
# Define harmonic transformation and associated harmonic space
- FFT = ift.FFTOperator(s_space)
- h_space = FFT.target[0]
+ HT = ift.HarmonicTransformOperator(h_space, target=s_space)
# Setting up power space
p_space = ift.PowerSpace(h_space,
@@ -51,22 +51,28 @@ if __name__ == "__main__":
mask = ift.Field(s_space, val=ift.dobj.from_global_data(mask))
MaskOperator = ift.DiagonalOperator(mask)
- InstrumentResponse = ift.ResponseOperator(s_space, sigma=[0.0],
- exposure=[1.0])
- MeasurementOperator = InstrumentResponse*MaskOperator
+ R = ift.GeometryRemover(s_space)
+ R = R*MaskOperator
+ #R = R*HT
+ #R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma)
+ MeasurementOperator = R
d_space = MeasurementOperator.target
- noise_covariance = ift.Field(d_space, val=noise_level**2).weight()
- N = ift.DiagonalOperator(noise_covariance)
- n = ift.Field.from_random(domain=d_space, random_type='normal',
- std=noise_level)
Projection = ift.PowerProjectionOperator(domain=h_space,
power_space=p_space)
power = Projection.adjoint_times(ift.exp(0.5*log_p))
# Creating the mock data
- true_sky = nonlinearity(FFT.adjoint_times(power*sh))
- d = MeasurementOperator(true_sky) + n
+ true_sky = nonlinearity(HT(power*sh))
+ noiseless_data = MeasurementOperator(true_sky)
+ noise_amplitude = noiseless_data.val.std()*noise_level
+ N = ift.DiagonalOperator(
+ ift.Field.full(d_space, noise_amplitude**2))
+ n = ift.Field.from_random(
+ domain=d_space, random_type='normal',
+ std=noise_amplitude, mean=0)
+ # Creating the mock data
+ d = noiseless_data + n
m0 = ift.power_synthesize(ift.Field(p_space, val=1e-7))
t0 = ift.Field(p_space, val=-4.)
@@ -84,7 +90,7 @@ if __name__ == "__main__":
for i in range(20):
power0 = Projection.adjoint_times(ift.exp(0.5*t0))
map0_energy = ift.library.NonlinearWienerFilterEnergy(
- m0, d, MeasurementOperator, nonlinearity, FFT, power0, N, S,
+ m0, d, MeasurementOperator, nonlinearity, HT, power0, N, S,
inverter=inverter)
# Minimization with chosen minimizer
@@ -96,7 +102,7 @@ if __name__ == "__main__":
# Initializing the power energy with updated parameters
power0_energy = ift.library.NonlinearPowerEnergy(
- position=t0, d=d, N=N, m=m0, D=D0, FFT=FFT,
+ position=t0, d=d, N=N, xi=m0, D=D0, ht=HT,
Instrument=MeasurementOperator, nonlinearity=nonlinearity,
Projection=Projection, sigma=1., samples=2, inverter=inverter)
@@ -110,6 +116,6 @@ if __name__ == "__main__":
m0, t0 = adjust_zero_mode(m0, t0)
ift.plot(true_sky)
- ift.plot(nonlinearity(FFT.adjoint_times(power0*m0)),
+ ift.plot(nonlinearity(HT(power0*m0)),
title='reconstructed_sky')
ift.plot(MeasurementOperator.adjoint_times(d))
diff --git a/nifty4/__init__.py b/nifty4/__init__.py
index a2fa13eae989212bd8f6c1ad464d00d16f162331..6d7457bddcf9602a34695924d549742abdcdcd11 100644
--- a/nifty4/__init__.py
+++ b/nifty4/__init__.py
@@ -22,8 +22,7 @@ from .operators.harmonic_transform_operator import HarmonicTransformOperator
from .operators.fft_operator import FFTOperator
from .operators.fft_smoothing_operator import FFTSmoothingOperator
from .operators.direct_smoothing_operator import DirectSmoothingOperator
-from .operators.response_operator import ResponseOperator
-from .operators.response_operator import GeometryRemover
+from .operators.geometry_remover import GeometryRemover
from .operators.laplace_operator import LaplaceOperator
from .operators.power_projection_operator import PowerProjectionOperator
from .operators.inversion_enabler import InversionEnabler
@@ -58,7 +57,7 @@ __all__ = ["Domain", "UnstructuredDomain", "StructuredDomain", "RGSpace", "LMSpa
"HPSpace", "GLSpace", "DOFSpace", "PowerSpace", "DomainTuple",
"LinearOperator", "EndomorphicOperator", "ScalingOperator",
"DiagonalOperator", "FFTOperator", "FFTSmoothingOperator",
- "DirectSmoothingOperator", "ResponseOperator", "LaplaceOperator",
+ "DirectSmoothingOperator", "LaplaceOperator",
"PowerProjectionOperator", "InversionEnabler",
"Field", "sqrt", "exp", "log",
"Prober", "DiagonalProberMixin", "TraceProberMixin"]
diff --git a/nifty4/operators/response_operator.py b/nifty4/operators/geometry_remover.py
similarity index 54%
rename from nifty4/operators/response_operator.py
rename to nifty4/operators/geometry_remover.py
index 1b980641df8b5509352b18bd73b8c17c6e88a6d7..f001517b7ccfe471c4f1ab320c10bae65e164443 100644
--- a/nifty4/operators/response_operator.py
+++ b/nifty4/operators/geometry_remover.py
@@ -20,10 +20,6 @@ from ..field import Field
from ..spaces.unstructured_domain import UnstructuredDomain
from ..domain_tuple import DomainTuple
from .linear_operator import LinearOperator
-from .fft_smoothing_operator import FFTSmoothingOperator
-from .diagonal_operator import DiagonalOperator
-from .scaling_operator import ScalingOperator
-import numpy as np
class GeometryRemover(LinearOperator):
@@ -50,38 +46,3 @@ class GeometryRemover(LinearOperator):
if mode == self.TIMES:
return Field(self._target, val=x.weight(1).val)
return Field(self._domain, val=x.val).weight(1)
-
-
-def ResponseOperator(domain, sigma, sensitivity):
- # sensitivity has units 1/field/volume and gives a measure of how much
- # the instrument will excited when it is exposed to a certain field
- # volume amplitude
- domain = DomainTuple.make(domain)
- ncomp = len(sensitivity)
- if len(sigma) != ncomp or len(domain) != ncomp:
- raise ValueError("length mismatch between sigma, sensitivity "
- "and domain")
- ncomp = len(sigma)
- if ncomp == 0:
- raise ValueError("Empty response operator not allowed")
-
- kernel = None
- sensi = 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(sensitivity[i]):
- if sensitivity[i] != 1.:
- op = ScalingOperator(sensitivity[i], domain)
- sensi = op if sensi is None else op*sensi
- elif isinstance(sensitivity[i], Field):
- op = DiagonalOperator(sensitivity[i], domain=domain, spaces=i)
- sensi = op if sensi is None else op*sensi
-
- res = GeometryRemover(domain)
- if sensi is not None:
- res = res * sensi
- if kernel is not None:
- res = res * kernel
- return res
diff --git a/test/test_operators/test_response_operator.py b/test/test_operators/test_response_operator.py
deleted file mode 100644
index 1d9a23170d2a0a2b1ec2ece57ef27d0ac5bb3302..0000000000000000000000000000000000000000
--- a/test/test_operators/test_response_operator.py
+++ /dev/null
@@ -1,48 +0,0 @@
-# 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.
-
-import unittest
-from numpy.testing import assert_allclose
-import nifty4 as ift
-from itertools import product
-from test.common import expand
-
-
-class ResponseOperator_Tests(unittest.TestCase):
- spaces = [ift.RGSpace(128), ift.GLSpace(nlat=37)]
-
- @expand(product(spaces, [0., 5., 1.], [0., 1., .33]))
- def test_property(self, space, sigma, sensitivity):
- if not isinstance(space, ift.RGSpace): # no smoothing supported
- sigma = 0.
- op = ift.ResponseOperator(space, sigma=[sigma],
- sensitivity=[sensitivity])
- if op.domain[0] != space:
- raise TypeError
-
- @expand(product(spaces, [0., 5., 1.], [0., 1., .33]))
- def test_times_adjoint_times(self, space, sigma, sensitivity):
- if not isinstance(space, ift.RGSpace): # no smoothing supported
- sigma = 0.
- op = ift.ResponseOperator(space, sigma=[sigma],
- sensitivity=[sensitivity])
- rand1 = ift.Field.from_random('normal', domain=space)
- rand2 = ift.Field.from_random('normal', domain=op.target[0])
- tt1 = rand2.vdot(op.times(rand1))
- tt2 = rand1.vdot(op.adjoint_times(rand2))
- assert_allclose(tt1, tt2)