Commit 6d16a1f2 authored by Martin Reinecke's avatar Martin Reinecke

cosmetics and adjustments

parent 89ff8afe
Pipeline #16533 passed with stage
in 19 minutes and 15 seconds
...@@ -23,9 +23,9 @@ def plot_parameters(m, t, p, p_d): ...@@ -23,9 +23,9 @@ def plot_parameters(m, t, p, p_d):
t = t.val.get_full_data().real t = t.val.get_full_data().real
p = p.val.get_full_data().real p = p.val.get_full_data().real
p_d = p_d.val.get_full_data().real p_d = p_d.val.get_full_data().real
pl.plot([go.Heatmap(z=m)], filename='map.html') pl.plot([go.Heatmap(z=m)], filename='map.html', auto_open=False)
pl.plot([go.Scatter(x=x, y=t), go.Scatter(x=x, y=p), pl.plot([go.Scatter(x=x, y=t), go.Scatter(x=x, y=p),
go.Scatter(x=x, y=p_d)], filename="t.html") go.Scatter(x=x, y=p_d)], filename="t.html", auto_open=False)
class AdjointFFTResponse(LinearOperator): class AdjointFFTResponse(LinearOperator):
...@@ -106,7 +106,7 @@ if __name__ == "__main__": ...@@ -106,7 +106,7 @@ if __name__ == "__main__":
data_power = log(fft(d).power_analyze(binbounds=p_space.binbounds)) data_power = log(fft(d).power_analyze(binbounds=p_space.binbounds))
d_data = d.val.get_full_data().real d_data = d.val.get_full_data().real
if rank == 0: if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html') pl.plot([go.Heatmap(z=d_data)], filename='data.html', auto_open=False)
# Minimization strategy # Minimization strategy
def convergence_measure(a_energy, iteration): # returns current energy def convergence_measure(a_energy, iteration): # returns current energy
......
...@@ -10,6 +10,7 @@ rank = comm.rank ...@@ -10,6 +10,7 @@ rank = comm.rank
np.random.seed(42) np.random.seed(42)
class AdjointFFTResponse(LinearOperator): class AdjointFFTResponse(LinearOperator):
def __init__(self, FFT, R, default_spaces=None): def __init__(self, FFT, R, default_spaces=None):
super(AdjointFFTResponse, self).__init__(default_spaces) super(AdjointFFTResponse, self).__init__(default_spaces)
...@@ -23,6 +24,7 @@ class AdjointFFTResponse(LinearOperator): ...@@ -23,6 +24,7 @@ class AdjointFFTResponse(LinearOperator):
def _adjoint_times(self, x, spaces=None): def _adjoint_times(self, x, spaces=None):
return self.FFT(self.R.adjoint_times(x)) return self.FFT(self.R.adjoint_times(x))
@property @property
def domain(self): def domain(self):
return self._domain return self._domain
...@@ -36,13 +38,12 @@ class AdjointFFTResponse(LinearOperator): ...@@ -36,13 +38,12 @@ class AdjointFFTResponse(LinearOperator):
return False return False
if __name__ == "__main__": if __name__ == "__main__":
distribution_strategy = 'not' distribution_strategy = 'not'
# Set up position space # Set up position space
s_space = RGSpace([128,128]) s_space = RGSpace([128, 128])
# s_space = HPSpace(32) # s_space = HPSpace(32)
# Define harmonic transformation and associated harmonic space # Define harmonic transformation and associated harmonic space
...@@ -52,7 +53,8 @@ if __name__ == "__main__": ...@@ -52,7 +53,8 @@ if __name__ == "__main__":
# Setting up power space # Setting up power space
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy) p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
# Choosing the prior correlation structure and defining correlation operator # 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 = create_power_operator(h_space, power_spectrum=p_spec, S = create_power_operator(h_space, power_spectrum=p_spec,
...@@ -69,7 +71,7 @@ if __name__ == "__main__": ...@@ -69,7 +71,7 @@ if __name__ == "__main__":
Instrument = DiagonalOperator(s_space, diagonal=1.) Instrument = DiagonalOperator(s_space, diagonal=1.)
# Instrument._diagonal.val[200:400, 200:400] = 0 # Instrument._diagonal.val[200:400, 200:400] = 0
#Adding a harmonic transformation to the instrument # Adding a harmonic transformation to the instrument
R = AdjointFFTResponse(fft, Instrument) R = AdjointFFTResponse(fft, Instrument)
signal_to_noise = 1. signal_to_noise = 1.
N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True) N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
...@@ -84,9 +86,9 @@ if __name__ == "__main__": ...@@ -84,9 +86,9 @@ if __name__ == "__main__":
# Choosing the minimization strategy # Choosing the minimization strategy
def convergence_measure(energy, iteration): # returns current energy def convergence_measure(energy, iteration): # returns current energy
x = energy.value x = energy.value
print (x, iteration) print(x, iteration)
# minimizer = SteepestDescent(convergence_tolerance=0, # minimizer = SteepestDescent(convergence_tolerance=0,
# iteration_limit=50, # iteration_limit=50,
...@@ -109,20 +111,19 @@ if __name__ == "__main__": ...@@ -109,20 +111,19 @@ if __name__ == "__main__":
m0 = Field(h_space, val=.0) m0 = Field(h_space, val=.0)
# Initializing the Wiener Filter energy # Initializing the Wiener Filter energy
energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S, inverter=inverter) energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S)
D0 = energy.curvature D0 = energy.curvature
# Solving the problem analytically # Solving the problem analytically
m0 = D0.inverse_times(j) m0 = D0.inverse_times(j)
sample_variance = Field(sh.domain,val=0. + 0j) sample_variance = Field(sh.domain, val=0. + 0j)
sample_mean = Field(sh.domain,val=0. + 0j) sample_mean = Field(sh.domain, val=0. + 0j)
# sampling the uncertainty map # sampling the uncertainty map
n_samples = 1 n_samples = 1
for i in range(n_samples): for i in range(n_samples):
sample = sugar.generate_posterior_sample(m0,D0) sample = sugar.generate_posterior_sample(m0, D0)
sample_variance += sample**2 sample_variance += sample**2
sample_mean += sample sample_mean += sample
variance = sample_variance/n_samples - (sample_mean/n_samples) variance = sample_variance/n_samples - (sample_mean/n_samples)
...@@ -7,7 +7,7 @@ class WienerFilterEnergy(Energy): ...@@ -7,7 +7,7 @@ class WienerFilterEnergy(Energy):
"""The Energy for the Wiener filter. """The Energy for the Wiener filter.
It covers the case of linear measurement with It covers the case of linear measurement with
Gaussian noise and Gaussain signal prior with known covariance. Gaussian noise and Gaussian signal prior with known covariance.
Parameters Parameters
---------- ----------
......
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