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

remove ResponseOperator

parent 6ad2e1ec
......@@ -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
......
......@@ -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))
......@@ -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"]
......@@ -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
# 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)
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