Skip to content
Snippets Groups Projects
Commit b7f06ed2 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

evolutionary changes

parent 51f49206
No related branches found
No related tags found
No related merge requests found
Pipeline #
...@@ -52,7 +52,7 @@ if __name__ == "__main__": ...@@ -52,7 +52,7 @@ if __name__ == "__main__":
h_space = fft.target[0] h_space = fft.target[0]
# Set up power space # Set up power space
p_space = ift.PowerSpace(h_space) p_space = ift.PowerSpace(h_space,binbounds=ift.PowerSpace.useful_binbounds(h_space,logarithmic=True))
# Choose the prior correlation structure and defining correlation operator # Choose the prior correlation structure and defining correlation operator
p_spec = (lambda k: (.5 / (k + 1) ** 3)) p_spec = (lambda k: (.5 / (k + 1) ** 3))
...@@ -88,11 +88,11 @@ if __name__ == "__main__": ...@@ -88,11 +88,11 @@ if __name__ == "__main__":
d_data = d.val.real d_data = d.val.real
ift.plotting.plot(d.real, name="data.pdf") ift.plotting.plot(d.real, name="data.pdf")
IC1 = ift.DefaultIterationController(verbose=True,iteration_limit=100) IC1 = ift.DefaultIterationController(verbose=True,iteration_limit=100,tol_abs_gradnorm=0.1)
minimizer1 = ift.RelaxedNewton(IC1) minimizer1 = ift.RelaxedNewton(IC1)
IC2 = ift.DefaultIterationController(verbose=True,iteration_limit=100) IC2 = ift.DefaultIterationController(verbose=True,iteration_limit=100,tol_abs_gradnorm=0.1)
minimizer2 = ift.VL_BFGS(IC2, max_history_length=20) minimizer2 = ift.VL_BFGS(IC2, max_history_length=20)
IC3 = ift.DefaultIterationController(verbose=True,iteration_limit=100) IC3 = ift.DefaultIterationController(verbose=True,iteration_limit=100,tol_abs_gradnorm=0.1)
minimizer3 = ift.SteepestDescent(IC3) minimizer3 = ift.SteepestDescent(IC3)
# Set starting position # Set starting position
...@@ -107,22 +107,22 @@ if __name__ == "__main__": ...@@ -107,22 +107,22 @@ if __name__ == "__main__":
S0 = ift.create_power_operator(h_space, power_spectrum=ps0) S0 = ift.create_power_operator(h_space, power_spectrum=ps0)
# Initialize non-linear Wiener Filter energy # Initialize non-linear Wiener Filter energy
ICI = ift.DefaultIterationController(verbose=False,iteration_limit=60) ICI = ift.DefaultIterationController(verbose=False,iteration_limit=500,tol_abs_gradnorm=0.1)
map_inverter = ift.ConjugateGradient(controller=ICI) 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 # Solve the Wiener Filter analytically
D0 = map_energy.curvature D0 = map_energy.curvature
m0 = D0.inverse_times(j) m0 = D0.inverse_times(j)
# Initialize power energy with updated parameters # Initialize power energy with updated parameters
ICI2 = ift.DefaultIterationController(verbose=False,iteration_limit=60) ICI2 = ift.DefaultIterationController(name="powI",verbose=True,iteration_limit=200,tol_abs_gradnorm=1e-5)
power_inverter = ift.ConjugateGradient(controller=ICI2) power_inverter = ift.ConjugateGradient(controller=ICI2)
power_energy = ift.library.CriticalPowerEnergy(position=t0, m=m0, D=D0, power_energy = ift.library.CriticalPowerEnergy(position=t0, m=m0, D=D0,
smoothness_prior=10., samples=3, inverter=power_inverter) smoothness_prior=10., samples=3, inverter=power_inverter)
(power_energy, convergence) = minimizer3(power_energy) (power_energy, convergence) = minimizer1(power_energy)
# Set new power spectrum # Set new power spectrum
t0.val = power_energy.position.val.real t0 = power_energy.position.real
# Plot current estimate # Plot current estimate
print(i) print(i)
......
...@@ -90,7 +90,7 @@ class CriticalPowerEnergy(Energy): ...@@ -90,7 +90,7 @@ class CriticalPowerEnergy(Energy):
gradient = -self._theta.weight(-1) gradient = -self._theta.weight(-1)
gradient += (self._rho_prime).weight(-1) gradient += (self._rho_prime).weight(-1)
gradient += self._Tt gradient += self._Tt
gradient.val = gradient.val.real gradient = gradient.real
return gradient return gradient
@property @property
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment