Commit a3f357c3 authored by Martin Reinecke's avatar Martin Reinecke

tweaks

parent fe1ce4a4
Pipeline #26452 passed with stage
in 5 minutes and 23 seconds
......@@ -35,7 +35,8 @@ if __name__ == "__main__":
d_space = R.target
power = ift.sqrt(ift.create_power_operator(h_space, p_spec).diagonal)
p_op = ift.create_power_operator(h_space, p_spec)
power = ift.sqrt(p_op(ift.Field.full(h_space, 1.)))
# Creating the mock data
true_sky = nonlinearity(HT(power*sh))
......
......@@ -36,12 +36,17 @@ if __name__ == "__main__":
ht_1 = ift.HarmonicTransformOperator(harmonic_domain, space=0)
ht_2 = ift.HarmonicTransformOperator(ht_1.target, space=1)
ht = ht_2*ht_1
del ht_1, ht_2
S = (ift.create_power_operator(harmonic_domain, power_spectrum_1, 0) *
ift.create_power_operator(harmonic_domain, power_spectrum_2, 1))
np.random.seed(10)
mock_signal = S.draw_sample()
plotdict = {"colormap": "Planck-like"}
plot_space = ift.RGSpace((N_pixels_1, N_pixels_2))
#ift.plot(ht(mock_signal).cast_domain(plot_space),
# name='mock_signal.png', **plotdict)
# Setting up a exemplary response
N1_10 = int(N_pixels_1/10)
......@@ -54,9 +59,10 @@ if __name__ == "__main__":
mask_2[N2_10*7:N2_10*9] = 0.
mask_2 = ift.Field.from_global_data(signal_space_2, mask_2)
R = ift.GeometryRemover(signal_domain)
R = R*ift.DiagonalOperator(mask_1, signal_domain, spaces=0)
#R = ift.GeometryRemover(signal_domain)
R = ift.DiagonalOperator(mask_1, signal_domain, spaces=0)
R = R*ift.DiagonalOperator(mask_2, signal_domain, spaces=1)
del mask_1, mask_2
R = R*ht
R = R * ift.create_harmonic_smoothing_operator(harmonic_domain, 0,
response_sigma_1)
......@@ -65,31 +71,31 @@ if __name__ == "__main__":
data_domain = R.target
noiseless_data = R(mock_signal)
del mock_signal
noise_amplitude = noiseless_data.val.std()/signal_to_noise
# Setting up the noise covariance and drawing a random noise realization
N = ift.ScalingOperator(noise_amplitude**2, data_domain)
noise = N.draw_sample()
data = noiseless_data + noise
#ift.plot(data.cast_domain(plot_space), name='data.png', **plotdict)
del noiseless_data, noise
# Wiener filter
j = R.adjoint_times(N.inverse_times(data))
del data
ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R, inverter=inverter)
del S, N, R, inverter
m_k = wiener_curvature.inverse_times(j)
m = ht(m_k)
del j
plotdict = {"colormap": "Planck-like"}
plot_space = ift.RGSpace((N_pixels_1, N_pixels_2))
ift.plot(ht(mock_signal).cast_domain(plot_space),
name='mock_signal.png', **plotdict)
ift.plot(data.cast_domain(plot_space), name='data.png', **plotdict)
ift.plot(m.cast_domain(plot_space), name='map.png', **plotdict)
#ift.plot(ht(m_k).cast_domain(plot_space), name='map.png', **plotdict)
# sampling the uncertainty map
mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 10)
ift.plot(ift.sqrt(variance).cast_domain(plot_space),
name="uncertainty.png", **plotdict)
ift.plot((mean+m).cast_domain(plot_space),
ift.plot((mean+ht(m_k)).cast_domain(plot_space),
name="posterior_mean.png", **plotdict)
......@@ -59,7 +59,7 @@ class NonlinearPowerEnergy(Energy):
self.D = D
self.d = d
self.N = N
self.T = SmoothnessOperator(domain=self.position.domain[0],
self.T = SmoothnessOperator(domain=position.domain[0],
strength=sigma, logarithmic=True)
self.ht = ht
self.Instrument = Instrument
......@@ -76,19 +76,15 @@ class NonlinearPowerEnergy(Energy):
self.inverter = inverter
A = Distributor(exp(.5 * position))
map_s = self.ht(A * xi)
Tpos = self.T(position)
self._gradient = None
for xi_sample in self.xi_sample_list:
map_s = self.ht(A * xi_sample)
LinR = LinearizedPowerResponse(
self.Instrument, self.nonlinearity, self.ht, self.Distributor,
self.position, xi_sample)
map_s = ht(A*xi_sample)
LinR = LinearizedPowerResponse(Instrument, nonlinearity, ht,
Distributor, position, xi_sample)
residual = self.d - \
self.Instrument(self.nonlinearity(map_s))
tmp = self.N.inverse_times(residual)
residual = self.d - Instrument(nonlinearity(map_s))
tmp = N.inverse_times(residual)
lh = 0.5 * residual.vdot(tmp)
grad = LinR.adjoint_times(tmp)
......@@ -100,7 +96,8 @@ class NonlinearPowerEnergy(Energy):
self._gradient += grad
self._value *= 1. / len(self.xi_sample_list)
self._value += 0.5 * self.position.vdot(Tpos)
Tpos = self.T(position)
self._value += 0.5 * position.vdot(Tpos)
self._gradient *= -1. / len(self.xi_sample_list)
self._gradient += Tpos
self._gradient.lock()
......
......@@ -31,20 +31,18 @@ class NonlinearWienerFilterEnergy(Energy):
self.nonlinearity = nonlinearity
self.ht = ht
self.power = power
m = self.ht(self.power*self.position)
self.LinearizedResponse = LinearizedSignalResponse(
Instrument, nonlinearity, ht, power, m)
m = ht(power*position)
residual = d - Instrument(nonlinearity(m))
self.N = N
self.S = S
self.inverter = inverter
t1 = self.S.inverse_times(self.position)
t2 = self.N.inverse_times(residual)
tmp = self.position.vdot(t1) + residual.vdot(t2)
self._value = 0.5 * tmp.real
self._gradient = t1 - self.LinearizedResponse.adjoint_times(t2)
self._gradient.lock()
t1 = S.inverse_times(position)
t2 = N.inverse_times(residual)
self._value = 0.5 * (position.vdot(t1) + residual.vdot(t2)).real
self.R = LinearizedSignalResponse(Instrument, nonlinearity, ht, power,
m)
self._gradient = (t1 - self.R.adjoint_times(t2)).lock()
def at(self, position):
return self.__class__(position, self.d, self.Instrument,
......@@ -62,5 +60,4 @@ class NonlinearWienerFilterEnergy(Energy):
@property
@memo
def curvature(self):
return WienerFilterCurvature(R=self.LinearizedResponse, N=self.N,
S=self.S, inverter=self.inverter)
return WienerFilterCurvature(self.R, self.N, self.S, self.inverter)
......@@ -51,12 +51,11 @@ class WienerFilterEnergy(Energy):
self._curvature = WienerFilterCurvature(R, N, S, inverter)
self._inverter = inverter
if _j is None:
_j = self.R.adjoint_times(self.N.inverse_times(d))
_j = R.adjoint_times(N.inverse_times(d))
self._j = _j
Dx = self._curvature(self.position)
self._value = 0.5*self.position.vdot(Dx) - self._j.vdot(self.position)
self._gradient = Dx - self._j
self._gradient.lock()
self._value = 0.5*position.vdot(Dx) - self._j.vdot(position)
self._gradient = (Dx - self._j).lock()
def at(self, position):
return self.__class__(position=position, d=None, R=self.R, N=self.N,
......
......@@ -71,8 +71,6 @@ class DiagonalOperator(EndomorphicOperator):
self._spaces = utilities.parse_spaces(spaces, len(self._domain))
if len(self._spaces) != len(diagonal.domain):
raise ValueError("spaces and domain must have the same length")
# if nspc==len(self.diagonal.domain),
# we could do some optimization
for i, j in enumerate(self._spaces):
if diagonal.domain[i] != self._domain[j]:
raise ValueError("domain mismatch")
......@@ -168,6 +166,8 @@ class DiagonalOperator(EndomorphicOperator):
@property
def adjoint(self):
if np.issubdtype(self._ldiag.dtype, np.floating):
return self
res = self._skeleton(())
res._ldiag = self._ldiag.conjugate()
return res
......
......@@ -61,7 +61,7 @@ class ScalingOperator(EndomorphicOperator):
if self._factor == 1.:
return x.copy()
if self._factor == 0.:
return Field.zeros_like(x, dtype=x.dtype)
return Field.zeros_like(x)
if mode == self.TIMES:
return x*self._factor
......@@ -81,6 +81,8 @@ class ScalingOperator(EndomorphicOperator):
@property
def adjoint(self):
if np.issubdtype(type(self._factor), np.floating):
return self
return ScalingOperator(np.conj(self._factor), self._domain)
@property
......
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