Skip to content

GitLab

Commit 011fe3fa authored by Martin Reinecke's avatar Martin Reinecke
Browse files

adjust

parent 6b18bbdf
Hide whitespace changes
Inline Side-by-side
demos/critical_filtering.py
View file @ 011fe3fa
......@@ -2,17 +2,15 @@ import numpy as np
import nifty2go as ift
np.random.seed(42)
#np.seterr(all="raise",under="ignore")
# np.seterr(all="raise",under="ignore")
def plot_parameters(m, t, p, p_d):
x = np.log(t.domain[0].k_lengths)
m = fft.adjoint_times(m)
t = t.val.real
p = p.val.real
p_d = p_d.val.real
ift.plotting.plot(m.real, name='map.pdf')
#pl.plot([go.Scatter(x=x, y=t), go.Scatter(x=x, y=p),
# go.Scatter(x=x, y=p_d)], filename="t.html", auto_open=False)
class AdjointFFTResponse(ift.LinearOperator):
......@@ -52,7 +50,9 @@ if __name__ == "__main__":
h_space = fft.target[0]
# Set up power space
p_space = ift.PowerSpace(h_space,binbounds=ift.PowerSpace.useful_binbounds(h_space,logarithmic=True))
p_space = ift.PowerSpace(h_space,
binbounds=ift.PowerSpace.useful_binbounds(
h_space, logarithmic=True))
# Choose the prior correlation structure and defining correlation operator
p_spec = (lambda k: (.5 / (k + 1) ** 3))
......@@ -60,7 +60,7 @@ if __name__ == "__main__":
# Draw a sample sh from the prior distribution in harmonic space
sp = ift.Field(p_space, val=p_spec(p_space.k_lengths))
sh = sp.power_synthesize(real_signal=True)
sh = ift.power_synthesize(sp, real_signal=True)
# Choose the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.01)
......@@ -72,32 +72,37 @@ if __name__ == "__main__":
R = AdjointFFTResponse(fft, Instrument)
noise = 1.
N = ift.DiagonalOperator(ift.Field.full(s_space,noise))
N = ift.DiagonalOperator(ift.Field.full(s_space, noise))
n = ift.Field.from_random(domain=s_space,
random_type='normal',
std=np.sqrt(noise),
mean=0)
random_type='normal',
std=np.sqrt(noise),
mean=0)
# Create mock data
d = R(sh) + n
# The information source
j = R.adjoint_times(N.inverse_times(d))
realized_power = ift.log(sh.power_analyze(binbounds=p_space.binbounds))
data_power = ift.log(fft(d).power_analyze(binbounds=p_space.binbounds))
realized_power = ift.log(ift.power_analyze(sh,
binbounds=p_space.binbounds))
data_power = ift.log(ift.power_analyze(fft(d),
binbounds=p_space.binbounds))
d_data = d.val.real
ift.plotting.plot(d.real, name="data.pdf")
IC1 = ift.GradientNormController(verbose=True,iteration_limit=100,tol_abs_gradnorm=0.1)
IC1 = ift.GradientNormController(verbose=True, iteration_limit=100,
tol_abs_gradnorm=0.1)
minimizer1 = ift.RelaxedNewton(IC1)
IC2 = ift.GradientNormController(verbose=True,iteration_limit=100,tol_abs_gradnorm=0.1)
IC2 = ift.GradientNormController(verbose=True, iteration_limit=100,
tol_abs_gradnorm=0.1)
minimizer2 = ift.VL_BFGS(IC2, max_history_length=20)
IC3 = ift.GradientNormController(verbose=True,iteration_limit=100,tol_abs_gradnorm=0.1)
IC3 = ift.GradientNormController(verbose=True, iteration_limit=100,
tol_abs_gradnorm=0.1)
minimizer3 = ift.SteepestDescent(IC3)
# Set starting position
flat_power = ift.Field.full(p_space, 1e-8)
m0 = flat_power.power_synthesize(real_signal=True)
m0 = ift.power_synthesize(flat_power, real_signal=True)
def ps0(k):
return (1./(1.+k)**2)
......@@ -107,17 +112,25 @@ if __name__ == "__main__":
S0 = ift.create_power_operator(h_space, power_spectrum=ps0)
# Initialize non-linear Wiener Filter energy
ICI = ift.GradientNormController(verbose=False,iteration_limit=500,tol_abs_gradnorm=0.1)
ICI = ift.GradientNormController(verbose=False, name="ICI",
iteration_limit=500,
tol_abs_gradnorm=0.1)
map_inverter = ift.ConjugateGradient(controller=ICI)
map_energy = ift.library.WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0, inverter=map_inverter)
map_energy = ift.library.WienerFilterEnergy(position=m0,
d=d, R=R, N=N, S=S0,
inverter=map_inverter)
# Solve the Wiener Filter analytically
D0 = map_energy.curvature
m0 = D0.inverse_times(j)
# Initialize power energy with updated parameters
ICI2 = ift.GradientNormController(name="powI",verbose=True,iteration_limit=200,tol_abs_gradnorm=1e-5)
ICI2 = ift.GradientNormController(name="powI",
verbose=False,
iteration_limit=200,
tol_abs_gradnorm=1e-5)
power_inverter = ift.ConjugateGradient(controller=ICI2)
power_energy = ift.library.CriticalPowerEnergy(position=t0, m=m0, D=D0,
smoothness_prior=10., samples=3, inverter=power_inverter)
power_energy = ift.library.CriticalPowerEnergy(
position=t0, m=m0, D=D0, smoothness_prior=10., samples=3,
inverter=power_inverter)
(power_energy, convergence) = minimizer1(power_energy)
......
demos/log_normal_wiener_filter.py
View file @ 011fe3fa
......@@ -25,7 +25,7 @@ if __name__ == "__main__":
# Creating the mock signal |\label{code:wf_mock_signal}|
S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum)
mock_power = ift.Field(power_space, val=power_spectrum(power_space.k_lengths))
mock_signal = fft(mock_power.power_synthesize(real_signal=True))
mock_signal = fft(ift.power_synthesize(mock_power, real_signal=True))
# Setting up an exemplary response
mask = ift.Field.ones(signal_space)
......
demos/paper_demos/cartesian_wiener_filter.py
View file @ 011fe3fa
......@@ -57,14 +57,14 @@ if __name__ == "__main__":
mock_power = ift.Field(domain=(power_space_1, power_space_2),
val=np.outer(mock_power_1.val, mock_power_2.val))
diagonal = mock_power.power_synthesize_special(spaces=(0, 1))**2
diagonal = ift.power_synthesize_special(mock_power, spaces=(0, 1))**2
diagonal = diagonal.real
S = ift.DiagonalOperator(diagonal.weight(-1))
np.random.seed(10)
mock_signal = fft(mock_power.power_synthesize(real_signal=True))
mock_signal = fft(ift.power_synthesize(mock_power, real_signal=True))
# Setting up a exemplary response
N1_10 = int(N_pixels_1/10)
......@@ -93,7 +93,8 @@ if __name__ == "__main__":
j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic, inverter=inverter)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic)
wiener_curvature = ift.InversionEnabler(wiener_curvature, inverter)
m_k = wiener_curvature.inverse_times(j) #|\label{code:wf_wiener_filter}|
m = fft(m_k)
......
demos/paper_demos/wiener_filter.py
View file @ 011fe3fa
......@@ -24,7 +24,7 @@ if __name__ == "__main__":
# Creating the mock signal |\label{code:wf_mock_signal}|
S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum)
mock_power = ift.Field(power_space, val=power_spectrum(power_space.k_lengths))
mock_signal = fft(mock_power.power_synthesize(real_signal=True))
mock_signal = fft(ift.power_synthesize(mock_power, real_signal=True))
# Setting up an exemplary response
mask = ift.Field.ones(signal_space)
......@@ -45,7 +45,8 @@ if __name__ == "__main__":
j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(verbose=True,tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic,inverter=inverter)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic)
wiener_curvature =ift.InversionEnabler(wiener_curvature, inverter)
m_k = wiener_curvature.inverse_times(j) #|\label{code:wf_wiener_filter}|
m = fft(m_k)
......
demos/wiener_filter_via_curvature.py
View file @ 011fe3fa
......@@ -11,7 +11,7 @@ if __name__ == "__main__":
# Typical distance over which the field is correlated
correlation_length = 0.05*nu.m
# sigma of field in position space sqrt(<|s_x|^2>)
field_sigma = 2.* nu.K
field_sigma = 2. * nu.K
# smoothing length of response
response_sigma = 0.01*nu.m
# The signal to noise ratio
......@@ -22,7 +22,8 @@ if __name__ == "__main__":
def power_spectrum(k):
cldim = correlation_length**(2*dimensionality)
a = 4/(2*np.pi) * cldim * field_sigma**2
return a / (1 + (k*correlation_length)**(2*dimensionality)) ** 2 # to be integrated over spherical shells later on
# to be integrated over spherical shells later on
return a / (1 + (k*correlation_length)**(2*dimensionality)) ** 2
# Setting up the geometry
......@@ -38,39 +39,47 @@ if __name__ == "__main__":
power_space = ift.PowerSpace(harmonic_space)
# Creating the mock data
S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum)
S = ift.create_power_operator(harmonic_space,
power_spectrum=power_spectrum)
np.random.seed(43)
mock_power = ift.Field(power_space, val=power_spectrum(power_space.k_lengths))
mock_harmonic = mock_power.power_synthesize(real_signal=True)
mock_power = ift.Field(power_space,
val=power_spectrum(power_space.k_lengths))
mock_harmonic = ift.power_synthesize(mock_power, real_signal=True)
mock_harmonic = mock_harmonic.real
mock_signal = fft(mock_harmonic)
exposure = 1.
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),exposure=(exposure,))
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
exposure=(exposure,))
data_domain = R.target[0]
R_harmonic = ift.ComposedOperator([fft, R])
N = ift.DiagonalOperator(ift.Field.full(data_domain,mock_signal.var()/signal_to_noise))
noise = ift.Field.from_random(domain=data_domain,
random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise),
mean=0)
N = ift.DiagonalOperator(ift.Field.full(data_domain,
mock_signal.var()/signal_to_noise))
noise = ift.Field.from_random(
domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(mock_signal) + noise
# Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(verbose=True,tol_abs_gradnorm=1e-4/nu.K)
ctrl = ift.GradientNormController(
verbose=True, tol_abs_gradnorm=1e-4/nu.K/(nu.m**(0.5*dimensionality)))
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N,
R=R_harmonic)
inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R_harmonic, inverter=inverter)
wiener_curvature = ift.InversionEnabler(wiener_curvature, inverter)
m = wiener_curvature.inverse_times(j)
m_s = fft(m)
sspace2=ift.RGSpace(shape, distances=L/N_pixels/nu.m)
sspace2 = ift.RGSpace(shape, distances=L/N_pixels/nu.m)
ift.plotting.plot(ift.Field(sspace2,mock_signal.real.val)/nu.K,name="mock_signal.pdf")
ift.plotting.plot(ift.Field(sspace2,
val=data.val.real.reshape(signal_space.shape))/nu.K, name="data.pdf")
ift.plotting.plot(ift.Field(sspace2,m_s.real.val)/nu.K, name="map.pdf")
ift.plotting.plot(ift.Field(sspace2, mock_signal.real.val)/nu.K,
name="mock_signal.pdf")
ift.plotting.plot(ift.Field(
sspace2, val=data.val.real.reshape(signal_space.shape))/nu.K,
name="data.pdf")
ift.plotting.plot(ift.Field(sspace2, m_s.real.val)/nu.K, name="map.pdf")
demos/wiener_filter_via_hamiltonian.py
View file @ 011fe3fa
......@@ -51,7 +51,7 @@ if __name__ == "__main__":
# Drawing a sample sh from the prior distribution in harmonic space
sp = ift.Field(p_space, val=p_spec(p_space.k_lengths))
sh = sp.power_synthesize(real_signal=True)
sh = ift.power_synthesize(sp, real_signal=True)
ss = fft.adjoint_times(sh)
# Choosing the measurement instrument
......
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