Commit b3b3a72a authored by Philipp Arras's avatar Philipp Arras

Do not pass inverters around anymore where it is not necessary

parent 0ac48f64
Pipeline #31375 failed with stages
in 4 minutes and 14 seconds
......@@ -169,10 +169,9 @@
"def Curvature(R, N, Sh):\n",
" IC = ift.GradientNormController(iteration_limit=50000,\n",
" tol_abs_gradnorm=0.1)\n",
" inverter = ift.ConjugateGradient(controller=IC)\n",
" # WienerFilterCurvature is (R.adjoint*N.inverse*R + Sh.inverse) plus some handy\n",
" # helper methods.\n",
" return ift.library.WienerFilterCurvature(R,N,Sh,inverter, sampling_inverter=inverter)"
" return ift.library.WienerFilterCurvature(R,N,Sh,iteration_controller=IC)"
]
},
{
......
......@@ -85,15 +85,13 @@ if __name__ == "__main__":
LS = ift.LineSearchStrongWolfe(c2=0.02)
minimizer = ift.RelaxedNewton(IC1, line_searcher=LS)
ICI = ift.GradientNormController(iteration_limit=500,
tol_abs_gradnorm=1e-3)
inverter = ift.ConjugateGradient(controller=ICI)
IC = ift.GradientNormController(iteration_limit=500,
tol_abs_gradnorm=1e-3)
for i in range(20):
power0 = Distributor(ift.exp(0.5*t0))
map0_energy = ift.library.NonlinearWienerFilterEnergy(
m0, d, MeasurementOperator, nonlinearity, HT, power0, N, S,
inverter=inverter)
m0, d, MeasurementOperator, nonlinearity, HT, power0, N, S, IC)
# Minimization with chosen minimizer
map0_energy, convergence = minimizer(map0_energy)
......@@ -106,7 +104,8 @@ if __name__ == "__main__":
power0_energy = ift.library.NonlinearPowerEnergy(
position=t0, d=d, N=N, xi=m0, D=D0, ht=HT,
Instrument=MeasurementOperator, nonlinearity=nonlinearity,
Distributor=Distributor, sigma=1., samples=2, inverter=inverter)
Distributor=Distributor, sigma=1., samples=2,
iteration_controller=IC)
power0_energy = minimizer(power0_energy)[0]
......
......@@ -78,15 +78,13 @@ if __name__ == "__main__":
LS = ift.LineSearchStrongWolfe(c2=0.02)
minimizer = ift.RelaxedNewton(IC1, line_searcher=LS)
ICI = ift.GradientNormController(iteration_limit=500,
tol_abs_gradnorm=1e-3)
inverter = ift.ConjugateGradient(controller=ICI)
IC = ift.GradientNormController(iteration_limit=500,
tol_abs_gradnorm=1e-3)
for i in range(20):
power0 = Distributor(ift.exp(0.5*t0))
map0_energy = ift.library.NonlinearWienerFilterEnergy(
m0, d, MeasurementOperator, nonlinearity, HT, power0, N, S,
inverter=inverter)
m0, d, MeasurementOperator, nonlinearity, HT, power0, N, S, IC)
# Minimization with chosen minimizer
map0_energy, convergence = minimizer(map0_energy)
......@@ -99,7 +97,7 @@ if __name__ == "__main__":
power0_energy = ift.library.NonlinearPowerEnergy(
position=t0, d=d, N=N, xi=m0, D=D0, ht=HT,
Instrument=MeasurementOperator, nonlinearity=nonlinearity,
Distributor=Distributor, sigma=1., samples=2, inverter=inverter)
Distributor=Distributor, sigma=1., samples=2, iteration_controller=IC)
power0_energy = minimizer(power0_energy)[0]
......
......@@ -52,14 +52,13 @@ if __name__ == "__main__":
LS = ift.LineSearchStrongWolfe(c2=0.02)
minimizer = ift.RelaxedNewton(IC1, line_searcher=LS)
ICI = ift.GradientNormController(iteration_limit=2000,
tol_abs_gradnorm=1e-3)
inverter = ift.ConjugateGradient(controller=ICI)
IC = ift.GradientNormController(iteration_limit=2000,
tol_abs_gradnorm=1e-3)
# initial guess
m = ift.full(h_space, 1e-7)
map_energy = ift.library.NonlinearWienerFilterEnergy(
m, d, R, nonlinearity, HT, power, N, S, inverter=inverter)
m, d, R, nonlinearity, HT, power, N, S, IC)
# Minimization with chosen minimizer
map_energy, convergence = minimizer(map_energy)
......
......@@ -76,10 +76,9 @@ if __name__ == "__main__":
ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1)
sampling_ctrl = ift.GradientNormController(name="sampling",
tol_abs_gradnorm=1e2)
inverter = ift.ConjugateGradient(controller=ctrl)
sampling_inverter = ift.ConjugateGradient(controller=sampling_ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R, inverter=inverter, sampling_inverter=sampling_inverter)
S=S, N=N, R=R, iteration_controller=ctrl,
iteration_controller_sampling=sampling_ctrl)
m_k = wiener_curvature.inverse_times(j)
m = ht(m_k)
......
......@@ -50,10 +50,9 @@ if __name__ == "__main__":
ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=1e-2)
sampling_ctrl = ift.GradientNormController(name="sampling",
tol_abs_gradnorm=2e1)
inverter = ift.ConjugateGradient(controller=ctrl)
sampling_inverter = ift.ConjugateGradient(controller=sampling_ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R, inverter=inverter, sampling_inverter=sampling_inverter)
S=S, N=N, R=R, iteration_controller=ctrl,
iteration_controller_sampling=sampling_ctrl)
m_k = wiener_curvature.inverse_times(j)
m = ht(m_k)
......
......@@ -81,9 +81,8 @@ if __name__ == "__main__":
j = R.adjoint_times(N.inverse_times(d))
IC = ift.GradientNormController(name="inverter", iteration_limit=500,
tol_abs_gradnorm=1e-3)
inverter = ift.ConjugateGradient(controller=IC)
D = (ift.SandwichOperator.make(R, N.inverse) + Phi_h.inverse).inverse
D = ift.InversionEnabler(D, inverter, approximation=Phi_h)
D = ift.InversionEnabler(D, IC, approximation=Phi_h)
m = HT(D(j))
# Uncertainty
......@@ -116,8 +115,7 @@ if __name__ == "__main__":
# initial guess
psi0 = ift.full(h_domain, 1e-7)
energy = ift.library.PoissonEnergy(psi0, data, R0, nonlin, HT, Phi_h,
inverter)
energy = ift.library.PoissonEnergy(psi0, data, R0, nonlin, HT, Phi_h, IC)
IC1 = ift.GradientNormController(name="IC1", iteration_limit=200,
tol_abs_gradnorm=1e-4)
minimizer = ift.RelaxedNewton(IC1)
......
......@@ -39,17 +39,17 @@ N_iter = 100
tol = 1e-3
IC = ift.GradientNormController(tol_abs_gradnorm=tol, iteration_limit=N_iter)
inverter = ift.ConjugateGradient(IC)
curv = ift.library.WienerFilterCurvature(S=S, N=N, R=R_p, inverter=inverter,
sampling_inverter=inverter)
curv = ift.library.WienerFilterCurvature(S=S, N=N, R=R_p,
iteration_controller=IC,
iteration_controller_sampling=IC)
m_xi = curv.inverse_times(j)
samps_long = [curv.draw_sample(from_inverse=True) for i in range(N_samps)]
tol = 1e2
IC = ift.GradientNormController(tol_abs_gradnorm=tol, iteration_limit=N_iter)
inverter = ift.ConjugateGradient(IC)
curv = ift.library.WienerFilterCurvature(S=S, N=N, R=R_p, inverter=inverter,
sampling_inverter=inverter)
curv = ift.library.WienerFilterCurvature(S=S, N=N, R=R_p,
iteration_controller=IC,
iteration_controller_sampling=IC)
samps_short = [curv.draw_sample(from_inverse=True) for i in range(N_samps)]
# Compute mean
......
......@@ -38,8 +38,8 @@ if __name__ == "__main__":
j = Rh.adjoint_times(N.inverse_times(d))
ctrl = ift.GradientNormController(name="Iter", tol_abs_gradnorm=1e-10,
iteration_limit=300)
inverter = ift.ConjugateGradient(controller=ctrl)
Di = ift.library.WienerFilterCurvature(S=S, R=Rh, N=N, inverter=inverter)
Di = ift.library.WienerFilterCurvature(S=S, R=Rh, N=N,
iteration_controller=ctrl)
mh = Di.inverse_times(j)
m = ht(mh)
......
......@@ -98,10 +98,9 @@ if __name__ == "__main__":
IC = ift.GradientNormController(name="inverter", iteration_limit=1000,
tol_abs_gradnorm=0.0001)
inverter = ift.ConjugateGradient(controller=IC)
# setting up measurement precision matrix M
M = (ift.SandwichOperator.make(R.adjoint, Sh) + N)
M = ift.InversionEnabler(M, inverter)
M = ift.InversionEnabler(M, IC)
m = Sh(R.adjoint(M.inverse_times(d)))
# Plotting
......
......@@ -52,9 +52,8 @@ if __name__ == "__main__":
j = R.adjoint_times(N.inverse_times(d))
IC = ift.GradientNormController(name="inverter", iteration_limit=500,
tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=IC)
D = (ift.SandwichOperator.make(R, N.inverse) + Sh.inverse).inverse
D = ift.InversionEnabler(D, inverter, approximation=Sh)
D = ift.InversionEnabler(D, IC, approximation=Sh)
m = D(j)
# Plotting
......
......@@ -78,9 +78,8 @@ if __name__ == "__main__":
j = R.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(
name="inverter", tol_abs_gradnorm=1e-5/(nu.K*(nu.m**dimensionality)))
inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R, inverter=inverter)
wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, R=R,
iteration_controller=ctrl)
m = wiener_curvature.inverse_times(j)
m_s = HT(m)
......
......@@ -47,15 +47,14 @@ if __name__ == "__main__":
# Choose minimization strategy
ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl)
controller = ift.GradientNormController(name="min", tol_abs_gradnorm=0.1)
minimizer = ift.RelaxedNewton(controller=controller)
m0 = ift.full(h_space, 0.)
# Initialize Wiener filter energy
energy = ift.library.WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S,
inverter=inverter,
sampling_inverter=inverter)
iteration_controller=ctrl,
iteration_controller_sampling=ctrl)
energy, convergence = minimizer(energy)
m = energy.position
......
......@@ -63,7 +63,7 @@ class NonlinearPowerEnergy(Energy):
# MR FIXME: docstring incomplete and outdated
def __init__(self, position, d, N, xi, D, ht, Instrument, nonlinearity,
Distributor, sigma=0., samples=3, xi_sample_list=None,
inverter=None):
iteration_controller=None):
super(NonlinearPowerEnergy, self).__init__(position)
self.xi = xi
self.D = D
......@@ -83,7 +83,7 @@ class NonlinearPowerEnergy(Energy):
xi_sample_list = [D.draw_sample(from_inverse=True) + xi
for _ in range(samples)]
self.xi_sample_list = xi_sample_list
self.inverter = inverter
self._ic = iteration_controller
A = Distributor(exp(.5 * position))
......@@ -118,7 +118,7 @@ class NonlinearPowerEnergy(Energy):
self.Distributor, sigma=self.sigma,
samples=len(self.xi_sample_list),
xi_sample_list=self.xi_sample_list,
inverter=self.inverter)
iteration_controller=self._ic)
@property
def value(self):
......@@ -139,4 +139,4 @@ class NonlinearPowerEnergy(Energy):
op = LinearizedResponse.adjoint*self.N.inverse*LinearizedResponse
result = op if result is None else result + op
result = result*(1./len(self.xi_sample_list)) + self.T
return InversionEnabler(result, self.inverter)
return InversionEnabler(result, self._ic)
......@@ -24,8 +24,7 @@ from ..sugar import makeOp
class NonlinearWienerFilterEnergy(Energy):
def __init__(self, position, d, Instrument, nonlinearity, ht, power, N, S,
inverter=None,
sampling_inverter=None):
iteration_controller=None, iteration_controller_sampling=None):
super(NonlinearWienerFilterEnergy, self).__init__(position=position)
self.d = d.lock()
self.Instrument = Instrument
......@@ -37,10 +36,10 @@ class NonlinearWienerFilterEnergy(Energy):
residual = d - Instrument(nonlinearity(m))
self.N = N
self.S = S
self.inverter = inverter
if sampling_inverter is None:
sampling_inverter = inverter
self.sampling_inverter = sampling_inverter
self._ic = iteration_controller
if iteration_controller_sampling is None:
iteration_controller_sampling = self._ic
self._ic_samp = iteration_controller_sampling
t1 = S.inverse_times(position)
t2 = N.inverse_times(residual)
self._value = 0.5 * (position.vdot(t1) + residual.vdot(t2)).real
......@@ -51,7 +50,7 @@ class NonlinearWienerFilterEnergy(Energy):
def at(self, position):
return self.__class__(position, self.d, self.Instrument,
self.nonlinearity, self.ht, self.power, self.N,
self.S, self.inverter)
self.S, self._ic, self._ic_samp)
@property
def value(self):
......@@ -64,5 +63,5 @@ class NonlinearWienerFilterEnergy(Energy):
@property
@memo
def curvature(self):
return WienerFilterCurvature(self.R, self.N, self.S, self.inverter,
self.sampling_inverter)
return WienerFilterCurvature(self.R, self.N, self.S, self._ic,
self._ic_samp)
......@@ -25,9 +25,9 @@ from ..sugar import log
class PoissonEnergy(Energy):
def __init__(self, position, d, Instrument, nonlinearity, ht, Phi_h,
inverter=None):
iteration_controller=None):
super(PoissonEnergy, self).__init__(position=position)
self._inverter = inverter
self._ic = iteration_controller
self._d = d
self._Instrument = Instrument
self._nonlinearity = nonlinearity
......@@ -51,7 +51,7 @@ class PoissonEnergy(Energy):
def at(self, position):
return self.__class__(position, self._d, self._Instrument,
self._nonlinearity, self._ht, self._Phi_h,
self._inverter)
self._ic)
@property
def value(self):
......@@ -63,5 +63,4 @@ class PoissonEnergy(Energy):
@property
def curvature(self):
return InversionEnabler(self._curv, self._inverter,
approximation=self._Phi_h.inverse)
return InversionEnabler(self._curv, self._ic, self._Phi_h.inverse)
......@@ -21,7 +21,8 @@ from ..operators.inversion_enabler import InversionEnabler
from ..operators.sampling_enabler import SamplingEnabler
def WienerFilterCurvature(R, N, S, inverter, sampling_inverter=None):
def WienerFilterCurvature(R, N, S, iteration_controller=None,
iteration_controller_sampling=None):
"""The curvature of the WienerFilterEnergy.
This operator implements the second derivative of the
......@@ -37,16 +38,16 @@ def WienerFilterCurvature(R, N, S, inverter, sampling_inverter=None):
The noise covariance.
S : DiagonalOperator
The prior signal covariance
inverter : Minimizer
The minimizer to use during numerical inversion
sampling_inverter : Minimizer
The minimizer to use during numerical sampling
if None, it is not possible to draw inverse samples
default: None
iteration_controller : IterationController
The iteration controller to use during numerical inversion via
ConjugateGradient.
iteration_controller_sampling : IterationController
The iteration controller to use for sampling.
"""
M = SandwichOperator.make(R, N.inverse)
if sampling_inverter is not None:
op = SamplingEnabler(M, S.inverse, sampling_inverter, S.inverse)
if iteration_controller is not None:
op = SamplingEnabler(M, S.inverse, iteration_controller_sampling,
S.inverse)
else:
op = M + S.inverse
return InversionEnabler(op, inverter, S.inverse)
return InversionEnabler(op, iteration_controller, S.inverse)
......@@ -20,8 +20,8 @@ from ..minimization.quadratic_energy import QuadraticEnergy
from .wiener_filter_curvature import WienerFilterCurvature
def WienerFilterEnergy(position, d, R, N, S, inverter=None,
sampling_inverter=None):
def WienerFilterEnergy(position, d, R, N, S, iteration_controller=None,
iteration_controller_sampling=None):
"""The Energy for the Wiener filter.
It covers the case of linear measurement with
......@@ -48,6 +48,7 @@ def WienerFilterEnergy(position, d, R, N, S, inverter=None,
if None, it is not possible to draw inverse samples
default: None
"""
op = WienerFilterCurvature(R, N, S, inverter, sampling_inverter)
op = WienerFilterCurvature(R, N, S, iteration_controller,
iteration_controller_sampling)
vec = R.adjoint_times(N.inverse_times(d))
return QuadraticEnergy(position, op, vec)
......@@ -61,6 +61,4 @@ class EnergySum(Energy):
if precon is None and self._precon_idx is not None:
precon = self._energies[self._precon_idx].curvature
from ..operators.inversion_enabler import InversionEnabler
from .conjugate_gradient import ConjugateGradient
return InversionEnabler(
res, ConjugateGradient(self._min_controller), precon)
return InversionEnabler(res, self._min_controller, precon)
......@@ -16,11 +16,13 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
from ..minimization.quadratic_energy import QuadraticEnergy
from ..minimization.iteration_controller import IterationController
import numpy as np
from ..logger import logger
from ..minimization.conjugate_gradient import ConjugateGradient
from ..minimization.iteration_controller import IterationController
from ..minimization.quadratic_energy import QuadraticEnergy
from .endomorphic_operator import EndomorphicOperator
import numpy as np
class InversionEnabler(EndomorphicOperator):
......@@ -34,9 +36,9 @@ class InversionEnabler(EndomorphicOperator):
The InversionEnabler object will support the same operation modes as
`op`, and additionally the inverse set. The newly-added modes will
be computed by iterative inversion.
inverter : :class:`Minimizer`
The minimizer to use for the iterative numerical inversion.
Typically, this is a :class:`ConjugateGradient` object.
iteration_controller : :class:`IterationController`
The iteration controller to use for the iterative numerical inversion
done by a :class:`ConjugateGradient` object.
approximation : :class:`LinearOperator`, optional
if not None, this operator should be an approximation to `op`, which
supports the operation modes that `op` doesn't have. It is used as a
......@@ -44,10 +46,10 @@ class InversionEnabler(EndomorphicOperator):
convergence.
"""
def __init__(self, op, inverter, approximation=None):
def __init__(self, op, iteration_controller, approximation=None):
super(InversionEnabler, self).__init__()
self._op = op
self._inverter = inverter
self._ic = iteration_controller
self._approximation = approximation
@property
......@@ -70,7 +72,8 @@ class InversionEnabler(EndomorphicOperator):
if prec is not None:
prec = prec._flip_modes(self._ilog[mode])
energy = QuadraticEnergy(x0, invop, x)
r, stat = self._inverter(energy, preconditioner=prec)
inverter = ConjugateGradient(self._ic)
r, stat = inverter(energy, preconditioner=prec)
if stat != IterationController.CONVERGED:
logger.warning("Error detected during operator inversion")
return r.position
......
......@@ -18,6 +18,7 @@
import numpy as np
from ..minimization.conjugate_gradient import ConjugateGradient
from ..minimization.quadratic_energy import QuadraticEnergy
from .endomorphic_operator import EndomorphicOperator
......@@ -33,9 +34,9 @@ class SamplingEnabler(EndomorphicOperator):
The InversionEnabler object will support the same operation modes as
`op`, and additionally the inverse set. The newly-added modes will
be computed by iterative inversion.
inverter : :class:`Minimizer`
The minimizer to use for the iterative numerical inversion.
Typically, this is a :class:`ConjugateGradient` object.
iteration_controller : :class:`IterationController`
The iteration controller to use for the iterative numerical inversion
done by a :class:`ConjugateGradient` object.
approximation : :class:`LinearOperator`, optional
if not None, this operator should be an approximation to `op`, which
supports the operation modes that `op` doesn't have. It is used as a
......@@ -43,13 +44,13 @@ class SamplingEnabler(EndomorphicOperator):
convergence.
"""
def __init__(self, likelihood, prior, sampling_inverter,
def __init__(self, likelihood, prior, iteration_controller,
approximation=None):
self._op = likelihood + prior
super(SamplingEnabler, self).__init__()
self._likelihood = likelihood
self._prior = prior
self._sampling_inverter = sampling_inverter
self._ic = iteration_controller
self._approximation = approximation
def draw_sample(self, from_inverse=False, dtype=np.float64):
......@@ -61,8 +62,9 @@ class SamplingEnabler(EndomorphicOperator):
nj = self._likelihood.draw_sample()
energy = QuadraticEnergy(s, self._op, sp + nj,
_grad=self._likelihood(s) - nj)
energy, convergence = self._sampling_inverter(
energy, preconditioner=self._approximation.inverse)
inverter = ConjugateGradient(self._ic)
energy, convergence = inverter(energy,
preconditioner=self._approximation.inverse)
return energy.position
@property
......
......@@ -56,11 +56,10 @@ class Energy_Tests(unittest.TestCase):
IC = ift.GradientNormController(
iteration_limit=100,
tol_abs_gradnorm=1e-5)
inverter = ift.ConjugateGradient(IC)
S = ift.create_power_operator(hspace, power_spectrum=_flat_PS)
energy = ift.library.WienerFilterEnergy(
position=s0, d=d, R=R, N=N, S=S, inverter=inverter)
position=s0, d=d, R=R, N=N, S=S, iteration_controller=IC)
ift.extra.check_value_gradient_curvature_consistency(
energy, tol=1e-4, ntries=10)
......
......@@ -65,7 +65,6 @@ class Noise_Energy_Tests(unittest.TestCase):
IC = ift.GradientNormController(
iteration_limit=100,
tol_abs_gradnorm=1e-5)
inverter = ift.ConjugateGradient(IC)
S = ift.create_power_operator(hspace, power_spectrum=_flat_PS)
C = ift.library.NonlinearWienerFilterEnergy(
......@@ -77,7 +76,7 @@ class Noise_Energy_Tests(unittest.TestCase):
power=A,
N=N,
S=S,
inverter=inverter).curvature
iteration_controller=IC).curvature
res_sample_list = [d - R(f(ht(C.draw_sample(from_inverse=True) + xi)))
for _ in range(10)]
......
......@@ -57,7 +57,6 @@ class Energy_Tests(unittest.TestCase):
IC = ift.GradientNormController(
iteration_limit=100,
tol_abs_gradnorm=1e-5)
inverter = ift.ConjugateGradient(IC)
S = ift.create_power_operator(hspace, power_spectrum=_flat_PS)
D = ift.library.NonlinearWienerFilterEnergy(
......@@ -69,7 +68,7 @@ class Energy_Tests(unittest.TestCase):
N=N,
S=S,
ht=ht,
inverter=inverter).curvature
iteration_controller=IC).curvature
energy = ift.library.NonlinearPowerEnergy(
position=tau0,
......
......@@ -116,9 +116,7 @@ class Test_Minimizers(unittest.TestCase):
t1 = ift.GradientNormController(tol_abs_gradnorm=1e-5,
iteration_limit=1000)
t2 = ift.ConjugateGradient(controller=t1)
return ift.InversionEnabler(RBCurv(self._position),
inverter=t2)
return ift.InversionEnabler(RBCurv(self._position), t1)
try:
minimizer = eval(minimizer)
......
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