Commit f5fae368 authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'tweak_demos' into 'master'

cosmetics and adjustments

See merge request !182
parents ec5edd53 6d16a1f2
Pipeline #16563 passed with stages
in 32 minutes and 40 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