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

fix 2D part of notebook

parent 28deed25
Pipeline #24370 passed with stage
in 6 minutes and 33 seconds
...@@ -291,16 +291,17 @@ ...@@ -291,16 +291,17 @@
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
sc = ift.probing.utils.StatCalculator() sc = ift.probing.utils.StatCalculator()
IC = ift.GradientNormController(name="inverter", iteration_limit=50000, IC = ift.GradientNormController(iteration_limit=50000,
tol_abs_gradnorm=0.1) tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=IC) inverter = ift.ConjugateGradient(controller=IC)
curv = ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,Sh,inverter) curv = ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,Sh,inverter)
for i in range(200): for i in range(200):
print i
sc.add(HT(curv.generate_posterior_sample())) sc.add(HT(curv.generate_posterior_sample()))
m_var = sc.var m_var = sc.var
``` ```
...@@ -357,37 +358,36 @@ ...@@ -357,37 +358,36 @@
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
N_pixels = 256 # Number of pixels N_pixels = 256 # Number of pixels
sigma2 = 1000 # Noise variance sigma2 = 10. # Noise variance
def pow_spec(k): def pow_spec(k):
P0, k0, gamma = [.2, 20, 4] P0, k0, gamma = [.2, 5, 4]
return P0 * (1. + (k/k0)**2)**(- gamma / 2) return P0 * (1. + (k/k0)**2)**(- gamma / 2)
s_space = ift.RGSpace([N_pixels, N_pixels]) s_space = ift.RGSpace([N_pixels, N_pixels])
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
h_space = s_space.get_default_codomain() h_space = s_space.get_default_codomain()
fft = ift.FFTOperator(s_space,h_space) HT = ift.HarmonicTransformOperator(h_space,s_space)
p_space = ift.PowerSpace(h_space) p_space = ift.PowerSpace(h_space)
# Operators # Operators
Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec) Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)
N = ift.ScalingOperator(sigma2,s_space) N = ift.ScalingOperator(sigma2,s_space)
R = ift.FFTSmoothingOperator(s_space, sigma=.01) R = ift.FFTSmoothingOperator(s_space, sigma=.01)
#D = PropagatorOperator(R=R, N=N, Sh=Sh) #D = PropagatorOperator(R=R, N=N, Sh=Sh)
# Fields and data # Fields and data
sh = ift.power_synthesize(ift.PS_field(p_space,pow_spec),real_signal=True) sh = ift.power_synthesize(ift.PS_field(p_space,pow_spec),real_signal=True)
s = fft.adjoint_times(sh)
n = ift.Field.from_random(domain=s_space, random_type='normal', n = ift.Field.from_random(domain=s_space, random_type='normal',
std=np.sqrt(sigma2), mean=0) std=np.sqrt(sigma2), mean=0)
# Lose some data # Lose some data
...@@ -395,45 +395,43 @@ ...@@ -395,45 +395,43 @@
h = int(N_pixels * 0.2 * 2) h = int(N_pixels * 0.2 * 2)
mask = ift.Field(s_space, val=1) mask = ift.Field(s_space, val=1)
mask.val[l:h,l:h] = 0 mask.val[l:h,l:h] = 0
R = ift.DiagonalOperator(mask) R = ift.DiagonalOperator(mask)*HT
n.val[l:h, l:h] = 0 n.val[l:h, l:h] = 0
D = PropagatorOperator(R=R, N=N, Sh=fft.inverse*Sh*fft) D = PropagatorOperator(R=R, N=N, Sh=Sh)
d = R(s) + n d = R(sh) + n
j = R.adjoint_times(N.inverse_times(d)) j = R.adjoint_times(N.inverse_times(d))
# Run Wiener filter # Run Wiener filter
m = D(j) m = D(j)
# Uncertainty # Uncertainty
sc = ift.probing.utils.StatCalculator() sc = ift.probing.utils.StatCalculator()
IC = ift.GradientNormController(name="inverter", iteration_limit=50000, IC = ift.GradientNormController(iteration_limit=50000,
tol_abs_gradnorm=0.1) tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=IC) inverter = ift.ConjugateGradient(controller=IC)
curv = ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,fft.inverse*Sh*fft,inverter) curv = ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,Sh,inverter)
for i in range(20): for i in range(20):
print i
sc.add(HT(curv.generate_posterior_sample())) sc.add(HT(curv.generate_posterior_sample()))
m_var = sc.var m_var = sc.var
diagProber = DiagonalProber(domain=s_space, probe_dtype=np.complex, probe_count=10)
diagProber(D)
m_var = Field(s_space, val=diagProber.diagonal.val).weight(-1)
# Get data # Get data
s_power = sh.power_analyze() s_power = ift.power_analyze(sh)
m_power = fft(m).power_analyze() m_power = ift.power_analyze(m)
s_power_data = s_power.val.get_full_data().real s_power_data = s_power.val.real
m_power_data = m_power.val.get_full_data().real m_power_data = m_power.val.real
s_data = s.val.get_full_data().real s_data = HT(sh).val.real
m_data = m.val.get_full_data().real m_data = HT(m).val.real
m_var_data = m_var.val.get_full_data().real m_var_data = m_var.val.real
d_data = d.val.get_full_data().real d_data = d.val.real
uncertainty = np.sqrt(np.abs(m_var_data)) uncertainty = np.sqrt(np.abs(m_var_data))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
...@@ -460,16 +458,10 @@ ...@@ -460,16 +458,10 @@
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig
```
%% Cell type:code id: tags:
``` python
mi = np.min(s_data) mi = np.min(s_data)
ma = np.max(s_data) ma = np.max(s_data)
fig, axes = plt.subplots(2, 2, figsize=(15, 15)) fig, axes = plt.subplots(2, 2, figsize=(15, 15))
...@@ -483,16 +475,10 @@ ...@@ -483,16 +475,10 @@
fig.subplots_adjust(right=0.8) fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([.85, 0.15, 0.05, 0.7]) cbar_ax = fig.add_axes([.85, 0.15, 0.05, 0.7])
fig.colorbar(im, cax=cbar_ax) fig.colorbar(im, cax=cbar_ax)
``` ```
%% Cell type:code id: tags:
``` python
fig
```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Is the uncertainty map reliable? ### Is the uncertainty map reliable?
%% Cell type:code id: tags: %% Cell type:code id: tags:
...@@ -502,11 +488,10 @@ ...@@ -502,11 +488,10 @@
print("Error within uncertainty map bounds: " + str(np.sum(precise) * 100 / N_pixels**2) + "%") print("Error within uncertainty map bounds: " + str(np.sum(precise) * 100 / N_pixels**2) + "%")
fig = plt.figure() fig = plt.figure()
plt.imshow(precise.astype(float), cmap="brg") plt.imshow(precise.astype(float), cmap="brg")
plt.colorbar() plt.colorbar()
fig
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Start Coding # Start Coding
......
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