Commit a7f3e8c4 authored by Martin Reinecke's avatar Martin Reinecke

some updates to energies and curvatures

parent f16ee8d4
Pipeline #17942 passed with stage
in 3 minutes and 38 seconds
......@@ -5,7 +5,7 @@ import plotly.graph_objs as go
import plotly.offline as pl
np.random.seed(42)
#np.seterr(all="raise",under="ignore")
def plot_parameters(m, t, p, p_d):
x = ift.log(t.domain[0].kindex)
......@@ -92,14 +92,9 @@ if __name__ == "__main__":
d_data = d.val.real
pl.plot([go.Heatmap(z=d_data)], filename='data.html', auto_open=False)
# Minimization strategy
def convergence_measure(a_energy, iteration): # returns current energy
x = a_energy.value
print(x, iteration)
IC1 = ift.DefaultIterationController(verbose=True,iteration_limit=5)
IC1 = ift.DefaultIterationController(verbose=True,iteration_limit=100)
minimizer1 = ift.RelaxedNewton(IC1)
IC2 = ift.DefaultIterationController(verbose=True,iteration_limit=30)
IC2 = ift.DefaultIterationController(verbose=True,iteration_limit=100)
minimizer2 = ift.VL_BFGS(IC2, max_history_length=20)
IC3 = ift.DefaultIterationController(verbose=True,iteration_limit=100)
minimizer3 = ift.SteepestDescent(IC3)
......@@ -116,17 +111,19 @@ if __name__ == "__main__":
S0 = ift.create_power_operator(h_space, power_spectrum=ps0)
# Initialize non-linear Wiener Filter energy
ICI = ift.DefaultIterationController(verbose=True,iteration_limit=60)
inverter = ift.ConjugateGradient(controller=ICI,preconditioner=S0.times)
map_energy = ift.library.WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0, inverter=inverter)
ICI = ift.DefaultIterationController(verbose=False,iteration_limit=60)
map_inverter = ift.ConjugateGradient(controller=ICI)
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.DefaultIterationController(verbose=False,iteration_limit=60)
power_inverter = ift.ConjugateGradient(controller=ICI2)
power_energy = ift.library.CriticalPowerEnergy(position=t0, m=m0, D=D0,
smoothness_prior=10., samples=3)
smoothness_prior=10., samples=3, inverter=power_inverter)
(power_energy, convergence) = minimizer2(power_energy)
(power_energy, convergence) = minimizer3(power_energy)
# Set new power spectrum
t0.val = power_energy.position.val.real
......
......@@ -21,8 +21,7 @@ class CriticalPowerCurvature(InvertibleOperatorMixin, EndomorphicOperator):
# ---Overwritten properties and methods---
def __init__(self, theta, T, inverter=None, preconditioner=None, **kwargs):
def __init__(self, theta, T, inverter, preconditioner=None, **kwargs):
self.theta = DiagonalOperator(theta.domain, diagonal=theta)
self.T = T
if preconditioner is None:
......
......@@ -54,7 +54,8 @@ class CriticalPowerEnergy(Energy):
# ---Overwritten properties and methods---
def __init__(self, position, m, D=None, alpha=1.0, q=0.,
smoothness_prior=0., logarithmic=True, samples=3, w=None):
smoothness_prior=0., logarithmic=True, samples=3, w=None,
inverter=None):
super(CriticalPowerEnergy, self).__init__(position=position)
self.m = m
self.D = D
......@@ -65,7 +66,8 @@ class CriticalPowerEnergy(Energy):
strength=smoothness_prior,
logarithmic=logarithmic)
self.rho = self.position.domain[0].rho
self._w = w if w is not None else None
self._w = w
self._inverter = inverter
# ---Mandatory properties and methods---
......@@ -73,7 +75,8 @@ class CriticalPowerEnergy(Energy):
return self.__class__(position, self.m, D=self.D, alpha=self.alpha,
q=self.q, smoothness_prior=self.smoothness_prior,
logarithmic=self.logarithmic,
w=self.w, samples=self.samples)
w=self.w, samples=self.samples,
inverter=self._inverter)
@property
def value(self):
......@@ -93,7 +96,7 @@ class CriticalPowerEnergy(Energy):
@property
def curvature(self):
curvature = CriticalPowerCurvature(theta=self._theta.weight(-1),
T=self.T)
T=self.T, inverter=self._inverter)
return curvature
# ---Added properties and methods---
......
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