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

adjust for new minimizers

parent 6e37062c
import numpy as np
from nifty import (VL_BFGS, DiagonalOperator, FFTOperator, Field,
LinearOperator, PowerSpace, RelaxedNewton, RGSpace,
SteepestDescent, create_power_operator, exp, log, sqrt)
import nifty as ift
from nifty.library.critical_filter import CriticalPowerEnergy
from nifty.library.wiener_filter import WienerFilterEnergy
......@@ -17,7 +15,7 @@ np.random.seed(42)
def plot_parameters(m, t, p, p_d):
x = log(t.domain[0].kindex)
x = np.log(t.domain[0].kindex)
m = fft.adjoint_times(m)
m = m.val.get_full_data().real
t = t.val.get_full_data().real
......@@ -28,7 +26,7 @@ def plot_parameters(m, t, p, p_d):
go.Scatter(x=x, y=p_d)], filename="t.html", auto_open=False)
class AdjointFFTResponse(LinearOperator):
class AdjointFFTResponse(ift.LinearOperator):
def __init__(self, FFT, R, default_spaces=None):
super(AdjointFFTResponse, self).__init__(default_spaces)
self._domain = FFT.target
......@@ -60,30 +58,32 @@ if __name__ == "__main__":
distribution_strategy = 'not'
# Set up position space
s_space = RGSpace([128, 128])
# s_space = HPSpace(32)
s_space = ift.RGSpace([128, 128])
# s_space = ift.HPSpace(32)
# Define harmonic transformation and associated harmonic space
fft = FFTOperator(s_space)
fft = ift.FFTOperator(s_space)
h_space = fft.target[0]
# Set up power space
p_space = PowerSpace(h_space, binbounds=PowerSpace.useful_binbounds(h_space,logarithmic=True),
distribution_strategy=distribution_strategy)
p_space = ift.PowerSpace(h_space,
binbounds=ift.PowerSpace.useful_binbounds(
h_space, logarithmic=True),
distribution_strategy=distribution_strategy)
# Choose the prior correlation structure and defining correlation operator
p_spec = (lambda k: (.5 / (k + 1) ** 3))
S = create_power_operator(h_space, power_spectrum=p_spec,
distribution_strategy=distribution_strategy)
S = ift.create_power_operator(h_space, power_spectrum=p_spec,
distribution_strategy=distribution_strategy)
# Draw a sample sh from the prior distribution in harmonic space
sp = Field(p_space, val=p_spec,
distribution_strategy=distribution_strategy)
sp = ift.Field(p_space, val=p_spec,
distribution_strategy=distribution_strategy)
sh = sp.power_synthesize(real_signal=True)
# Choose the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.01)
Instrument = DiagonalOperator(s_space, diagonal=1.)
# Instrument = ift.SmoothingOperator(s_space, sigma=0.01)
Instrument = ift.DiagonalOperator(s_space, diagonal=1.)
# Instrument._diagonal.val[200:400, 200:400] = 0
# Instrument._diagonal.val[64:512-64, 64:512-64] = 0
......@@ -91,20 +91,20 @@ if __name__ == "__main__":
R = AdjointFFTResponse(fft, Instrument)
noise = 1.
ndiag = Field(s_space, noise).weight(1)
N = DiagonalOperator(s_space, ndiag)
n = Field.from_random(domain=s_space,
random_type='normal',
std=sqrt(noise),
mean=0)
ndiag = ift.Field(s_space, noise).weight(1)
N = ift.DiagonalOperator(s_space, ndiag)
n = ift.Field.from_random(domain=s_space,
random_type='normal',
std=ift.sqrt(noise),
mean=0)
# Create mock data
d = R(sh) + n
# The information source
j = R.adjoint_times(N.inverse_times(d))
realized_power = log(sh.power_analyze(binbounds=p_space.binbounds))
data_power = log(fft(d).power_analyze(binbounds=p_space.binbounds))
realized_power = ift.log(sh.power_analyze(binbounds=p_space.binbounds))
data_power = ift.log(fft(d).power_analyze(binbounds=p_space.binbounds))
d_data = d.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html', auto_open=False)
......@@ -114,28 +114,25 @@ if __name__ == "__main__":
x = a_energy.value
print(x, iteration)
minimizer1 = RelaxedNewton(convergence_tolerance=1e-4,
convergence_level=1,
iteration_limit=5,
callback=convergence_measure)
minimizer2 = VL_BFGS(convergence_tolerance=1e-10,
convergence_level=1,
iteration_limit=30,
callback=convergence_measure,
max_history_length=20)
minimizer3 = SteepestDescent(convergence_tolerance=1e-4,
iteration_limit=100,
callback=convergence_measure)
IC1 = ift.GradientNormController(iteration_limit=5,
tol_abs_gradnorm=0.1)
minimizer1 = ift.RelaxedNewton(IC1)
IC2 = ift.GradientNormController(iteration_limit=30,
tol_abs_gradnorm=0.1)
minimizer2 = ift.VL_BFGS(IC2, max_history_length=20)
IC3 = ift.GradientNormController(iteration_limit=100,
tol_abs_gradnorm=0.1)
minimizer3 = ift.SteepestDescent(IC3)
# Set starting position
flat_power = Field(p_space, val=1e-8)
flat_power = ift.Field(p_space, val=1e-8)
m0 = flat_power.power_synthesize(real_signal=True)
t0 = Field(p_space, val=log(1./(1+p_space.kindex)**2))
t0 = ift.Field(p_space, val=ift.log(1./(1+p_space.kindex)**2))
for i in range(500):
S0 = create_power_operator(h_space, power_spectrum=exp(t0),
distribution_strategy=distribution_strategy)
S0 = ift.create_power_operator(h_space, power_spectrum=ift.exp(t0),
distribution_strategy=distribution_strategy)
# Initialize non-linear Wiener Filter energy
map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0)
......@@ -154,4 +151,4 @@ if __name__ == "__main__":
# Plot current estimate
print(i)
if i % 5 == 0:
plot_parameters(m0, t0, log(sp), data_power)
plot_parameters(m0, t0, ift.log(sp), data_power)
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