diff --git a/demos/Wiener_Filter.ipynb b/demos/Wiener_Filter.ipynb index 19074536831d2c9912584097283a5d3beb05409f..59f046d5b8d9c91bd4430f6a8b5049dd68506d6b 100644 --- a/demos/Wiener_Filter.ipynb +++ b/demos/Wiener_Filter.ipynb @@ -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)" ] }, { diff --git a/demos/critical_filtering.py b/demos/critical_filtering.py index d859ab895223e110139e4d92d6c42c350e124092..51bcc81b4182b03feb2f3e499410ce97bda64a75 100644 --- a/demos/critical_filtering.py +++ b/demos/critical_filtering.py @@ -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] diff --git a/demos/nonlinear_critical_filter.py b/demos/nonlinear_critical_filter.py index 8df6da0aa76c204180148cdfb0bbd042aa311ac6..8aee9c6407b44aa89a0b43726b4395243f4ec49a 100644 --- a/demos/nonlinear_critical_filter.py +++ b/demos/nonlinear_critical_filter.py @@ -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] diff --git a/demos/nonlinear_wiener_filter.py b/demos/nonlinear_wiener_filter.py index 3004b5a010507fda04200db3679ea6888b3f5850..01c2ea0f98ceba456b0421c9e7a66911801d9ae3 100644 --- a/demos/nonlinear_wiener_filter.py +++ b/demos/nonlinear_wiener_filter.py @@ -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) diff --git a/demos/paper_demos/cartesian_wiener_filter.py b/demos/paper_demos/cartesian_wiener_filter.py index 2cec36eb05993dd208b188d88d0617292fd46240..cf851e1327910601769496ee6ffc6ee713195afd 100644 --- a/demos/paper_demos/cartesian_wiener_filter.py +++ b/demos/paper_demos/cartesian_wiener_filter.py @@ -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) diff --git a/demos/paper_demos/wiener_filter.py b/demos/paper_demos/wiener_filter.py index 59f29ddb320cba87e061b799a4585221b7a9d6cd..e9f199d2ab35acbab626904434879fd4ec9d9a33 100644 --- a/demos/paper_demos/wiener_filter.py +++ b/demos/paper_demos/wiener_filter.py @@ -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) diff --git a/demos/poisson_demo.py b/demos/poisson_demo.py index 26bff8186e0bda319a31554937323f450d51f66d..64a6e0ef52b41561e72a1702cd0288b372204ab3 100644 --- a/demos/poisson_demo.py +++ b/demos/poisson_demo.py @@ -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) diff --git a/demos/sampling.py b/demos/sampling.py index 43ef9d1eabe15d102935bd271cda59b53f81fa0f..21f0d2e86fbf3901296070866e05bbd4ea5b4f84 100644 --- a/demos/sampling.py +++ b/demos/sampling.py @@ -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 diff --git a/demos/tomography.py b/demos/tomography.py index b2e3cffdef65d629c9bd66194cae8fa01b7d34bc..3957e99c20405fa89740c7108a4022c04937a5d7 100644 --- a/demos/tomography.py +++ b/demos/tomography.py @@ -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) diff --git a/demos/wiener_filter_data_space_noiseless.py b/demos/wiener_filter_data_space_noiseless.py index a8911d546dd9de62b13151b47a6418e24461032d..9310f62660c99b631ef69e934dd5ca55f616c1a7 100644 --- a/demos/wiener_filter_data_space_noiseless.py +++ b/demos/wiener_filter_data_space_noiseless.py @@ -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 diff --git a/demos/wiener_filter_easy.py b/demos/wiener_filter_easy.py index 3631c29922fb9e96f3f7a4cf993e012008a3ace1..49749fe957e86e2efa9a84cfd4330999b6ad1e3a 100644 --- a/demos/wiener_filter_easy.py +++ b/demos/wiener_filter_easy.py @@ -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 diff --git a/demos/wiener_filter_via_curvature.py b/demos/wiener_filter_via_curvature.py index 5717120578a46ba9310639c67e0bd2a2fb4ff99f..d83553e86d02e040551a4db76c97c49b3eb1439d 100644 --- a/demos/wiener_filter_via_curvature.py +++ b/demos/wiener_filter_via_curvature.py @@ -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) diff --git a/demos/wiener_filter_via_hamiltonian.py b/demos/wiener_filter_via_hamiltonian.py index a63eb3635a1e2ea6eb5c151a2f8a1f66cbb385fd..bc2d2a37badddeec6d7e1cc6b75800ad35f2d679 100644 --- a/demos/wiener_filter_via_hamiltonian.py +++ b/demos/wiener_filter_via_hamiltonian.py @@ -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 diff --git a/nifty5/library/nonlinear_power_energy.py b/nifty5/library/nonlinear_power_energy.py index 5d28126d963253b1fa7e4b60d24b6a8cce68ec6c..65ca8ebec04fc7bace95692fd1981b62926416a4 100644 --- a/nifty5/library/nonlinear_power_energy.py +++ b/nifty5/library/nonlinear_power_energy.py @@ -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) diff --git a/nifty5/library/nonlinear_wiener_filter_energy.py b/nifty5/library/nonlinear_wiener_filter_energy.py index c0ecf3ef5288e43ba62a1f58d4a530a117718b75..def3a90fa23019ba141789651efa3c05479e7ad2 100644 --- a/nifty5/library/nonlinear_wiener_filter_energy.py +++ b/nifty5/library/nonlinear_wiener_filter_energy.py @@ -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) diff --git a/nifty5/library/poisson_energy.py b/nifty5/library/poisson_energy.py index f3ea413161ade66cd77ba28c2ad887c7d23f13ed..5705717461abae0fdc003644b989b486af88828a 100644 --- a/nifty5/library/poisson_energy.py +++ b/nifty5/library/poisson_energy.py @@ -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) diff --git a/nifty5/library/wiener_filter_curvature.py b/nifty5/library/wiener_filter_curvature.py index eeed8d7a7878abddbfdb4e90b548b2cac8fc27f2..4ff71df6e20f75613c067c0590051fe0d09d3457 100644 --- a/nifty5/library/wiener_filter_curvature.py +++ b/nifty5/library/wiener_filter_curvature.py @@ -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) diff --git a/nifty5/library/wiener_filter_energy.py b/nifty5/library/wiener_filter_energy.py index 276aa71e22953e908045491cf15fce2f2ec71137..c5b7dbdb19233d36bdc51a920ce35b6cde2457bc 100644 --- a/nifty5/library/wiener_filter_energy.py +++ b/nifty5/library/wiener_filter_energy.py @@ -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) diff --git a/nifty5/minimization/energy_sum.py b/nifty5/minimization/energy_sum.py index 5a4453109787e8add6b9c34785e91284ca2f1d36..5155947a79800f04e168828b1c8b0ced52558210 100644 --- a/nifty5/minimization/energy_sum.py +++ b/nifty5/minimization/energy_sum.py @@ -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) diff --git a/nifty5/operators/inversion_enabler.py b/nifty5/operators/inversion_enabler.py index e670bed811338101bcf8a88f97ef186d5adfae1c..9980d8492f0781d436232712f63d5f38ae8eeb66 100644 --- a/nifty5/operators/inversion_enabler.py +++ b/nifty5/operators/inversion_enabler.py @@ -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 diff --git a/nifty5/operators/sampling_enabler.py b/nifty5/operators/sampling_enabler.py index 4fb373e8e800bcd7a77c03f997fb5f9061e2382b..e3a66a11061d8b3bc5970ff93fac828bb5f07723 100644 --- a/nifty5/operators/sampling_enabler.py +++ b/nifty5/operators/sampling_enabler.py @@ -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 diff --git a/test/test_energies/test_map.py b/test/test_energies/test_map.py index 597f499616d50daf014e55539c2128a428bfa97c..7c7e5f5208f150dfa8968460394abb7aa6f56cf4 100644 --- a/test/test_energies/test_map.py +++ b/test/test_energies/test_map.py @@ -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) diff --git a/test/test_energies/test_noise.py b/test/test_energies/test_noise.py index d62284ce7935263e4187a0ffe74d718c4831c6a2..a762919db25b13ff86f52c5afa7f3b2596c7bd4a 100644 --- a/test/test_energies/test_noise.py +++ b/test/test_energies/test_noise.py @@ -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)] diff --git a/test/test_energies/test_power.py b/test/test_energies/test_power.py index fee623d886d306874f960b5cd9d95863fed02fc5..99669117f8a7ab42c89a49d7d43cca6cd0a0bc4f 100644 --- a/test/test_energies/test_power.py +++ b/test/test_energies/test_power.py @@ -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, diff --git a/test/test_minimization/test_minimizers.py b/test/test_minimization/test_minimizers.py index dbd029faac1497695705312f25054c79355d79fe..060fe30b35191596c7d7a9f9ffb933059e977ff4 100644 --- a/test/test_minimization/test_minimizers.py +++ b/test/test_minimization/test_minimizers.py @@ -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)