diff --git a/demos/nonlinear_critical_filter.py b/demos/nonlinear_critical_filter.py index 4e8f0e825a678b1c71cd7126bf11b90e8f2b2fe5..afa14cd55477cf76435c66abd854b899a06886c9 100644 --- a/demos/nonlinear_critical_filter.py +++ b/demos/nonlinear_critical_filter.py @@ -1,6 +1,5 @@ import nifty2go as ift -from nifty2go.library import NonlinearWienerFilterEnergy, NonlinearPowerEnergy -from nifty2go.library.nonlinearities import * +from nifty2go.library.nonlinearities import Exponential import numpy as np np.random.seed(42) @@ -32,9 +31,10 @@ if __name__ == "__main__": p_space = ift.PowerSpace(h_space, binbounds=ift.PowerSpace.useful_binbounds( h_space, logarithmic=True)) - s_spec = ift.Field(p_space,val=1.) - # Choosing the prior correlation structure and defining correlation operator - p = ift.Field(p_space,val = p_spec(p_space.k_lengths)) + s_spec = ift.Field(p_space, val=1.) + # Choosing the prior correlation structure and defining + # correlation operator + p = ift.Field(p_space, val=p_spec(p_space.k_lengths)) log_p = ift.log(p) S = ift.create_power_operator(h_space, power_spectrum=s_spec) @@ -42,40 +42,40 @@ if __name__ == "__main__": sp = ift.Field(p_space, val=s_spec) sh = ift.power_synthesize(sp) - # Choosing the measurement instrument # Instrument = SmoothingOperator(s_space, sigma=0.01) mask = np.ones(s_space.shape) mask[6000:8000] = 0. - mask = ift.Field(s_space,val=mask) + mask = ift.Field(s_space, val=mask) MaskOperator = ift.DiagonalOperator(mask) InstrumentResponse = ift.ResponseOperator(s_space, sigma=[0.0]) - MeasurementOperator = ift.ComposedOperator([MaskOperator, InstrumentResponse]) + MeasurementOperator = ift.ComposedOperator([MaskOperator, + InstrumentResponse]) d_space = MeasurementOperator.target noise_covariance = ift.Field(d_space, val=noise_level**2).weight() N = ift.DiagonalOperator(noise_covariance) n = ift.Field.from_random(domain=d_space, - random_type='normal', - std=noise_level, - mean=0) - Projection = ift.PowerProjectionOperator(domain= h_space, power_space=p_space) + random_type='normal', + std=noise_level, + mean=0) + Projection = ift.PowerProjectionOperator(domain=h_space, + power_space=p_space) power = Projection.adjoint_times(ift.exp(0.5 * log_p)) # Creating the mock data true_sky = nonlinearity(FFT.adjoint_times(power * sh)) d = MeasurementOperator(true_sky) + n - - - flat_power = ift.Field(p_space,val=10e-8) + flat_power = ift.Field(p_space, val=1e-7) m0 = ift.power_synthesize(flat_power) - t0 = ift.Field(p_space,val=-4.) + t0 = ift.Field(p_space, val=-4.) power0 = Projection.adjoint_times(ift.exp(0.5 * t0)) IC1 = ift.GradientNormController(verbose=True, iteration_limit=100, tol_abs_gradnorm=1e-3) - minimizer = ift.RelaxedNewton(IC1) + LS = ift.LineSearchStrongWolfe(c2=0.02) + minimizer = ift.RelaxedNewton(IC1, line_searcher=LS) ICI = ift.GradientNormController(verbose=False, name="ICI", iteration_limit=500, @@ -85,22 +85,22 @@ if __name__ == "__main__": for i in range(20): power0 = Projection.adjoint_times(ift.exp(0.5*t0)) - map0_energy = NonlinearWienerFilterEnergy(m0, d, MeasurementOperator, - nonlinearity, FFT, power0, N, S, inverter=inverter) + map0_energy = ift.library.NonlinearWienerFilterEnergy( + m0, d, MeasurementOperator, nonlinearity, FFT, power0, N, S, + inverter=inverter) # Minimization with chosen minimizer map0_energy, convergence = minimizer(map0_energy) m0 = map0_energy.position - # Updating parameters for correlation structure reconstruction D0 = map0_energy.curvature # Initializing the power energy with updated parameters - power0_energy = NonlinearPowerEnergy(position=t0, d=d, N=N, m=m0, - D=D0, FFT=FFT, Instrument=MeasurementOperator, - nonlinearity=nonlinearity,Projection=Projection, - sigma=1., samples=2, inverter=inverter) + power0_energy = ift.library.NonlinearPowerEnergy( + position=t0, d=d, N=N, m=m0, D=D0, FFT=FFT, + Instrument=MeasurementOperator, nonlinearity=nonlinearity, + Projection=Projection, sigma=1., samples=2, inverter=inverter) (power0_energy, convergence) = minimizer(power0_energy) @@ -108,9 +108,11 @@ if __name__ == "__main__": t0_ = t0.copy() t0 = power0_energy.position - # break degeneracy between amplitude and excitation by setting excitation monopole to 1 - m0, t0 = adjust_zero_mode(m0,t0) + # break degeneracy between amplitude and excitation by setting + # excitation monopole to 1 + m0, t0 = adjust_zero_mode(m0, t0) ift.plotting.plot(true_sky) - ift.plotting.plot(nonlinearity(FFT.adjoint_times(power0*m0)),title='reconstructed_sky') + ift.plotting.plot(nonlinearity(FFT.adjoint_times(power0*m0)), + title='reconstructed_sky') ift.plotting.plot(MeasurementOperator.adjoint_times(d)) diff --git a/demos/paper_demos/cartesian_wiener_filter.py b/demos/paper_demos/cartesian_wiener_filter.py index a4dcfda4436e611fddeb34430421642077ab0d89..d64c6ad76e5dad2a32a93eb8815e919354683ea5 100644 --- a/demos/paper_demos/cartesian_wiener_filter.py +++ b/demos/paper_demos/cartesian_wiener_filter.py @@ -93,9 +93,8 @@ if __name__ == "__main__": j = R_harmonic.adjoint_times(N.inverse_times(data)) ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1) inverter = ift.ConjugateGradient(controller=ctrl) - wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, - R=R_harmonic) - wiener_curvature = ift.InversionEnabler(wiener_curvature, inverter) + wiener_curvature = ift.library.WienerFilterCurvature( + S=S, N=N, R=R_harmonic, inverter=inverter) m_k = wiener_curvature.inverse_times(j) m = fft(m_k) diff --git a/demos/paper_demos/wiener_filter.py b/demos/paper_demos/wiener_filter.py index 9f414b9eb545b31ecb99c99d038bde81f8c21e26..306605c2bf978bd237a0c97546069c9d122f0d0d 100644 --- a/demos/paper_demos/wiener_filter.py +++ b/demos/paper_demos/wiener_filter.py @@ -49,9 +49,8 @@ if __name__ == "__main__": j = R_harmonic.adjoint_times(N.inverse_times(data)) ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1) inverter = ift.ConjugateGradient(controller=ctrl) - wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, - R=R_harmonic) - wiener_curvature = ift.InversionEnabler(wiener_curvature, inverter) + wiener_curvature = ift.library.WienerFilterCurvature( + S=S, N=N, R=R_harmonic, inverter=inverter) m_k = wiener_curvature.inverse_times(j) m = fft(m_k) diff --git a/demos/wiener_filter_easy.py b/demos/wiener_filter_easy.py index f7f6251b32187c4645db5830324d19572e0db26c..2ff3709e729c87a7b2867f4c020a016999743faa 100644 --- a/demos/wiener_filter_easy.py +++ b/demos/wiener_filter_easy.py @@ -5,9 +5,10 @@ import nifty2go as ift # Note that the constructor of PropagatorOperator takes as arguments the # response R and noise covariance N operating on signal space and signal # covariance operating on harmonic space. -class PropagatorOperator(ift.EndomorphicOperator): - def __init__(self, R, N, Sh): - super(PropagatorOperator, self).__init__() +class PropagatorOperator(ift.InversionEnabler, ift.EndomorphicOperator): + def __init__(self, R, N, Sh, inverter): + ift.InversionEnabler.__init__(self, inverter) + ift.EndomorphicOperator.__init__(self) self.R = R self.N = N @@ -84,7 +85,6 @@ if __name__ == "__main__": j = R.adjoint_times(N.inverse_times(d)) IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=0.1) inverter = ift.ConjugateGradient(controller=IC) - D = ift.InversionEnabler(PropagatorOperator(Sh=Sh, N=N, R=R), - inverter=inverter) + D = PropagatorOperator(Sh=Sh, N=N, R=R, inverter=inverter) m = D(j) diff --git a/demos/wiener_filter_via_curvature.py b/demos/wiener_filter_via_curvature.py index fb2340e7f001a079d6188f4b12953dc8d6996e8f..3758085d47a13b05720d79ad5a263c75cb367d8e 100644 --- a/demos/wiener_filter_via_curvature.py +++ b/demos/wiener_filter_via_curvature.py @@ -69,10 +69,9 @@ if __name__ == "__main__": j = R_harmonic.adjoint_times(N.inverse_times(data)) ctrl = ift.GradientNormController( verbose=True, tol_abs_gradnorm=1e-4/nu.K/(nu.m**(0.5*dimensionality))) - wiener_curvature = ift.library.WienerFilterCurvature(S=S, N=N, - R=R_harmonic) inverter = ift.ConjugateGradient(controller=ctrl) - wiener_curvature = ift.InversionEnabler(wiener_curvature, inverter) + wiener_curvature = ift.library.WienerFilterCurvature( + S=S, N=N, R=R_harmonic, inverter=inverter) m = wiener_curvature.inverse_times(j) m_s = fft(m) diff --git a/demos/wiener_filter_via_hamiltonian.py b/demos/wiener_filter_via_hamiltonian.py index 3bf23165b515478c6377cd7a60257787eba3a893..354cbb76aa8dfaaff156bafbf634348c857da39f 100644 --- a/demos/wiener_filter_via_hamiltonian.py +++ b/demos/wiener_filter_via_hamiltonian.py @@ -71,7 +71,7 @@ if __name__ == "__main__": n_samples = 50 for i in range(n_samples): - sample = fft(ift.sugar.generate_posterior_sample(m, D)) + sample = fft(D.generate_posterior_sample(m)) sample_variance += sample**2 sample_mean += sample sample_mean /= n_samples diff --git a/nifty/library/critical_power_curvature.py b/nifty/library/critical_power_curvature.py index 9efe32706b7e855bca91d5bd2730e187f4f7c3d1..d84b3f91c7c45667d8fb3162c78fc28c037d9d0c 100644 --- a/nifty/library/critical_power_curvature.py +++ b/nifty/library/critical_power_curvature.py @@ -1,8 +1,7 @@ -from ..operators.endomorphic_operator import EndomorphicOperator -from ..operators.diagonal_operator import DiagonalOperator +from ..operators import EndomorphicOperator, InversionEnabler, DiagonalOperator -class CriticalPowerCurvature(EndomorphicOperator): +class CriticalPowerCurvature(InversionEnabler, EndomorphicOperator): """The curvature of the CriticalPowerEnergy. This operator implements the second derivative of the @@ -17,15 +16,12 @@ class CriticalPowerCurvature(EndomorphicOperator): The smoothness prior contribution to the curvature. """ - def __init__(self, theta, T): - super(CriticalPowerCurvature, self).__init__() + def __init__(self, theta, T, inverter): + EndomorphicOperator.__init__(self) self._theta = DiagonalOperator(theta) + InversionEnabler.__init__(self, inverter, self._theta.inverse_times) self._T = T - @property - def preconditioner(self): - return self._theta.inverse_times - def _times(self, x): return self._T(x) + self._theta(x) diff --git a/nifty/library/critical_power_energy.py b/nifty/library/critical_power_energy.py index 528ae5fa8438ff5294168f355178dadfa48b67ff..d38739780854ddde7426d6752eb8e0a0c4d321d3 100644 --- a/nifty/library/critical_power_energy.py +++ b/nifty/library/critical_power_energy.py @@ -1,11 +1,9 @@ from ..minimization.energy import Energy from ..operators.smoothness_operator import SmoothnessOperator from ..operators.power_projection_operator import PowerProjectionOperator -from ..operators.inversion_enabler import InversionEnabler from .critical_power_curvature import CriticalPowerCurvature from ..utilities import memo from .. import Field, exp -from ..sugar import generate_posterior_sample class CriticalPowerEnergy(Energy): @@ -67,14 +65,12 @@ class CriticalPowerEnergy(Energy): self._inverter = inverter if w is None: - # self.logger.info("Initializing w") P = PowerProjectionOperator(domain=self.m.domain, power_space=self.position.domain[0]) if self.D is not None: w = Field.zeros(self.position.domain, dtype=self.m.dtype) for i in range(self.samples): - # self.logger.info("Drawing sample %i" % i) - sample = generate_posterior_sample(self.m, self.D) + sample = self.D.generate_posterior_sample(self.m) w += P(abs(sample)**2) w *= 1./self.samples @@ -82,15 +78,15 @@ class CriticalPowerEnergy(Energy): w = P(abs(self.m)**2) self._w = w - theta = exp(-self.position) * (self.q + self._w*0.5) + self._theta = exp(-self.position) * (self.q + self._w*0.5) Tt = self.T(self.position) - energy = theta.integrate() + energy = self._theta.integrate() energy += self.position.integrate()*(self.alpha-0.5) energy += 0.5*self.position.vdot(Tt) self._value = energy.real - gradient = -theta + gradient = -self._theta gradient += self.alpha-0.5 gradient += Tt self._gradient = gradient.real @@ -99,7 +95,7 @@ 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, - samples=self.samples, w=self.w, + samples=self.samples, w=self._w, inverter=self._inverter) @property @@ -113,8 +109,8 @@ class CriticalPowerEnergy(Energy): @property @memo def curvature(self): - curv = CriticalPowerCurvature(theta=self._theta, T=self.T) - return InversionEnabler(curv, inverter=self._inverter) + return CriticalPowerCurvature(theta=self._theta, T=self.T, + inverter=self._inverter) @property def logarithmic(self): diff --git a/nifty/library/log_normal_wiener_filter_curvature.py b/nifty/library/log_normal_wiener_filter_curvature.py index c7b0e4bd592971eeb51671e98747a6aead518444..068dba17a5b3ec1128e56f195155710a3445255b 100644 --- a/nifty/library/log_normal_wiener_filter_curvature.py +++ b/nifty/library/log_normal_wiener_filter_curvature.py @@ -1,9 +1,9 @@ -from ..operators import EndomorphicOperator +from ..operators import EndomorphicOperator, InversionEnabler from ..utilities import memo from ..field import exp -class LogNormalWienerFilterCurvature(EndomorphicOperator): +class LogNormalWienerFilterCurvature(InversionEnabler, EndomorphicOperator): """The curvature of the LogNormalWienerFilterEnergy. This operator implements the second derivative of the @@ -21,8 +21,9 @@ class LogNormalWienerFilterCurvature(EndomorphicOperator): The prior signal covariance """ - def __init__(self, R, N, S, position, fft4exp): - super(LogNormalWienerFilterCurvature, self).__init__() + def __init__(self, R, N, S, position, fft4exp, inverter): + InversionEnabler.__init__(self, inverter) + EndomorphicOperator.__init__(self) self.R = R self.N = N self.S = S diff --git a/nifty/library/log_normal_wiener_filter_energy.py b/nifty/library/log_normal_wiener_filter_energy.py index f1ecc3725644ea9ec59d82159295d6d49394b179..5d1c760be0a8c4d40f2a4dc323f294c69c5b05ba 100644 --- a/nifty/library/log_normal_wiener_filter_energy.py +++ b/nifty/library/log_normal_wiener_filter_energy.py @@ -2,7 +2,6 @@ from ..minimization.energy import Energy from ..utilities import memo from .log_normal_wiener_filter_curvature import LogNormalWienerFilterCurvature from ..sugar import create_composed_fft_operator -from ..operators.inversion_enabler import InversionEnabler class LogNormalWienerFilterEnergy(Energy): @@ -58,11 +57,9 @@ class LogNormalWienerFilterEnergy(Energy): @property @memo def curvature(self): - return InversionEnabler( - LogNormalWienerFilterCurvature(R=self.R, N=self.N, S=self.S, - position=self.position, - fft4exp=self._fft), - inverter=self._inverter) + return LogNormalWienerFilterCurvature( + R=self.R, N=self.N, S=self.S, position=self.position, + fft4exp=self._fft, inverter=self._inverter) @property @memo @@ -72,7 +69,7 @@ class LogNormalWienerFilterEnergy(Energy): @property @memo def _Rexppd(self): - expp = self._fft.adjoint_times(self.curvature.op._expp_sspace) + expp = self._fft.adjoint_times(self.curvature._expp_sspace) return self.R(expp) - self.d @property @@ -84,5 +81,5 @@ class LogNormalWienerFilterEnergy(Energy): @memo def _exppRNRexppd(self): return self._fft.adjoint_times( - self.curvature.op._expp_sspace * + self.curvature._expp_sspace * self._fft(self.R.adjoint_times(self._NRexppd))) diff --git a/nifty/library/noise_energy.py b/nifty/library/noise_energy.py index 66acb8f5de26360850d6053dc8b417446d969814..137bc8b6c1d2cfc678095b2f3e5ff0b6f1600de3 100644 --- a/nifty/library/noise_energy.py +++ b/nifty/library/noise_energy.py @@ -1,6 +1,5 @@ from .. import Field, exp from ..operators.diagonal_operator import DiagonalOperator -from ..sugar import generate_posterior_sample from ..minimization.energy import Energy from ..utilities import memo @@ -31,7 +30,7 @@ class NoiseEnergy(Energy): sample_list.append(self.m) else: for i in range(samples): - sample = generate_posterior_sample(m, D) + sample = D.generate_posterior_sample(m) sample = FFT(Field(FFT.domain, val=( FFT.adjoint_times(sample).val))) sample_list.append(sample) diff --git a/nifty/library/nonlinear_power_curvature.py b/nifty/library/nonlinear_power_curvature.py index 6c51404853611a42fd7f1260040d482624cebc58..d2845013690edc09efff284c1e18c977e08bea61 100644 --- a/nifty/library/nonlinear_power_curvature.py +++ b/nifty/library/nonlinear_power_curvature.py @@ -1,12 +1,13 @@ -from ..operators.endomorphic_operator import EndomorphicOperator +from ..operators import EndomorphicOperator, InversionEnabler from .response_operators import LinearizedPowerResponse -class NonlinearPowerCurvature(EndomorphicOperator): +class NonlinearPowerCurvature(InversionEnabler, EndomorphicOperator): def __init__(self, position, FFT, Instrument, nonlinearity, - Projection, N, T, sample_list): - super(NonlinearPowerCurvature, self).__init__() + Projection, N, T, sample_list, inverter): + InversionEnabler.__init__(self, inverter) + EndomorphicOperator.__init__(self) self.N = N self.FFT = FFT self.Instrument = Instrument diff --git a/nifty/library/nonlinear_power_energy.py b/nifty/library/nonlinear_power_energy.py index f774d5745a0c9663bbe3f603ed54267397f1db16..a79f4c825b0c17f38037088b56901f0864484db2 100644 --- a/nifty/library/nonlinear_power_energy.py +++ b/nifty/library/nonlinear_power_energy.py @@ -1,11 +1,9 @@ from .. import exp from ..utilities import memo from ..operators.smoothness_operator import SmoothnessOperator -from ..sugar import generate_posterior_sample from .nonlinear_power_curvature import NonlinearPowerCurvature from .response_operators import LinearizedPowerResponse from ..minimization.energy import Energy -from ..operators.inversion_enabler import InversionEnabler class NonlinearPowerEnergy(Energy): @@ -56,7 +54,7 @@ class NonlinearPowerEnergy(Energy): if samples is None or samples == 0: sample_list = [m] else: - sample_list = [generate_posterior_sample(m, D) + sample_list = [D.generate_posterior_sample(m) for _ in range(samples)] self.sample_list = sample_list self.inverter = inverter @@ -105,7 +103,7 @@ class NonlinearPowerEnergy(Energy): @property @memo def curvature(self): - curvature = NonlinearPowerCurvature( + return NonlinearPowerCurvature( self.position, self.FFT, self.Instrument, self.nonlinearity, - self.Projection, self.N, self.T, self.sample_list) - return InversionEnabler(curvature, inverter=self.inverter) + self.Projection, self.N, self.T, self.sample_list, + inverter=self.inverter) diff --git a/nifty/library/nonlinear_wiener_filter_energy.py b/nifty/library/nonlinear_wiener_filter_energy.py index 1781332e7f7931649716b6fb977cf9e45ab40906..065ecb6b8a4edc17428c4b75fc4f89837d1b0a62 100644 --- a/nifty/library/nonlinear_wiener_filter_energy.py +++ b/nifty/library/nonlinear_wiener_filter_energy.py @@ -1,7 +1,6 @@ from .wiener_filter_curvature import WienerFilterCurvature from ..utilities import memo from ..minimization.energy import Energy -from ..operators.inversion_enabler import InversionEnabler from .response_operators import LinearizedSignalResponse @@ -45,6 +44,5 @@ class NonlinearWienerFilterEnergy(Energy): @property @memo def curvature(self): - curvature = WienerFilterCurvature(R=self.LinearizedResponse, - N=self.N, S=self.S) - return InversionEnabler(curvature, inverter=self.inverter) + return WienerFilterCurvature(R=self.LinearizedResponse, N=self.N, + S=self.S, inverter=self.inverter) diff --git a/nifty/library/response_operators.py b/nifty/library/response_operators.py index f2f05c129672310d93da2b31e9206ed707342ffe..7b698c9b25640765ccf3a1fd7ca733dfdd2a58a7 100644 --- a/nifty/library/response_operators.py +++ b/nifty/library/response_operators.py @@ -42,7 +42,7 @@ class LinearizedPowerResponse(LinearOperator): self.Instrument = Instrument self.FFT = FFT self.Projection = Projection - self.power = exp(0.5 * t) + self.power = exp(0.5*t) self.m = m position = FFT.adjoint_times( self.Projection.adjoint_times(self.power) * self.m) diff --git a/nifty/library/wiener_filter_curvature.py b/nifty/library/wiener_filter_curvature.py index 20f55637375ddd20cf3d78c56d719c580e2072c4..50933501acbbce8a08580062b6e36396355474c7 100644 --- a/nifty/library/wiener_filter_curvature.py +++ b/nifty/library/wiener_filter_curvature.py @@ -1,7 +1,9 @@ -from ..operators import EndomorphicOperator +from ..operators import EndomorphicOperator, InversionEnabler +from ..field import Field, sqrt +from ..sugar import power_analyze, power_synthesize -class WienerFilterCurvature(EndomorphicOperator): +class WienerFilterCurvature(InversionEnabler, EndomorphicOperator): """The curvature of the WienerFilterEnergy. This operator implements the second derivative of the @@ -19,16 +21,13 @@ class WienerFilterCurvature(EndomorphicOperator): The prior signal covariance """ - def __init__(self, R, N, S): - super(WienerFilterCurvature, self).__init__() + def __init__(self, R, N, S, inverter): + EndomorphicOperator.__init__(self) + InversionEnabler.__init__(self, inverter, S.times) self.R = R self.N = N self.S = S - @property - def preconditioner(self): - return self.S.times - @property def domain(self): return self.S.domain @@ -45,3 +44,40 @@ class WienerFilterCurvature(EndomorphicOperator): res = self.R.adjoint_times(self.N.inverse_times(self.R(x))) res += self.S.inverse_times(x) return res + + def generate_posterior_sample(self, mean): + """ Generates a posterior sample from a Gaussian distribution with + given mean and covariance. + + This method generates samples by setting up the observation and + reconstruction of a mock signal in order to obtain residuals of the + right correlation which are added to the given mean. + + Parameters + ---------- + mean : Field + the mean of the posterior Gaussian distribution + + Returns + ------- + sample : Field + Returns the a sample from the Gaussian of given mean and + covariance. + """ + + power = sqrt(power_analyze(self.S.diagonal())) + mock_signal = power_synthesize(power, real_signal=True) + + noise = self.N.diagonal().weight(-1) + + mock_noise = Field.from_random(random_type="normal", + domain=self.N.domain, + dtype=noise.dtype.type) + mock_noise *= sqrt(noise) + + mock_data = self.R(mock_signal) + mock_noise + + mock_j = self.R.adjoint_times(self.N.inverse_times(mock_data)) + mock_m = self.inverse_times(mock_j) + sample = mock_signal - mock_m + mean + return sample diff --git a/nifty/library/wiener_filter_energy.py b/nifty/library/wiener_filter_energy.py index 7b8ef0cac73376ff4395a4f61bbcba716797d5a8..3bf42e8046ce2087f929101422757ce23cc5283d 100644 --- a/nifty/library/wiener_filter_energy.py +++ b/nifty/library/wiener_filter_energy.py @@ -1,5 +1,4 @@ from ..minimization.energy import Energy -from ..operators.inversion_enabler import InversionEnabler from .wiener_filter_curvature import WienerFilterCurvature @@ -28,8 +27,7 @@ class WienerFilterEnergy(Energy): self.R = R self.N = N self.S = S - self._curvature = InversionEnabler(WienerFilterCurvature(R, N, S), - inverter=inverter) + self._curvature = WienerFilterCurvature(R, N, S, inverter=inverter) self._inverter = inverter if _j is None: _j = self.R.adjoint_times(self.N.inverse_times(d)) diff --git a/nifty/minimization/conjugate_gradient.py b/nifty/minimization/conjugate_gradient.py index 9de83875c32372b035b76fca2dede3a8466d1a2b..31457c17d0e7af8050d0e7763c12fc9dc896907a 100644 --- a/nifty/minimization/conjugate_gradient.py +++ b/nifty/minimization/conjugate_gradient.py @@ -32,9 +32,6 @@ class ConjugateGradient(Minimizer): ---------- controller : IterationController Object that decides when to terminate the minimization. - preconditioner : Operator *optional* - This operator can be provided which transforms the variables of the - system to improve the conditioning (default: None). References ---------- @@ -42,8 +39,7 @@ class ConjugateGradient(Minimizer): 2006, Springer-Verlag New York """ - def __init__(self, controller, preconditioner=None): - self._preconditioner = preconditioner + def __init__(self, controller): self._controller = controller def __call__(self, energy, preconditioner=None): @@ -54,6 +50,9 @@ class ConjugateGradient(Minimizer): energy : Energy object at the starting point of the iteration. Its curvature operator must be independent of position, otherwise linear conjugate gradient minimization will fail. + preconditioner : Operator *optional* + This operator can be provided which transforms the variables of the + system to improve the conditioning (default: None). Returns ------- @@ -62,9 +61,6 @@ class ConjugateGradient(Minimizer): status : integer Can be controller.CONVERGED or controller.ERROR """ - if preconditioner is None: - preconditioner = self._preconditioner - controller = self._controller status = controller.start(energy) if status != controller.CONTINUE: diff --git a/nifty/operators/inversion_enabler.py b/nifty/operators/inversion_enabler.py index fa0e3634d458089b338ef708bf60e5c9ec2f7338..719da004bc99c723be73b7f56ad0eb7c04017d7f 100644 --- a/nifty/operators/inversion_enabler.py +++ b/nifty/operators/inversion_enabler.py @@ -18,57 +18,29 @@ from ..minimization.quadratic_energy import QuadraticEnergy from ..field import Field -from .linear_operator import LinearOperator -class InversionEnabler(LinearOperator): +class InversionEnabler(object): - def __init__(self, op, inverter, preconditioner=None): - self._op = op - self._inverter = inverter - if preconditioner is None and hasattr(op, "preconditioner"): - self._preconditioner = op.preconditioner - else: - self._preconditioner = preconditioner + def __init__(self, inverter, preconditioner=None): super(InversionEnabler, self).__init__() + self._inverter = inverter + self._preconditioner = preconditioner - @property - def domain(self): - return self._op.domain - - @property - def target(self): - return self._op.target - - @property - def unitary(self): - return self._op.unitary - - @property - def op(self): - return self._op - - def _operation(self, x, o1, o2, tdom): - try: - return o1(x) - except NotImplementedError: - x0 = Field.zeros(tdom, dtype=x.dtype) - energy = QuadraticEnergy(A=o2, b=x, position=x0) - r = self._inverter(energy, preconditioner=self._preconditioner)[0] - return r.position + def _operation(self, x, op, tdom): + x0 = Field.zeros(tdom, dtype=x.dtype) + energy = QuadraticEnergy(A=op, b=x, position=x0) + r = self._inverter(energy, preconditioner=self._preconditioner)[0] + return r.position def _times(self, x): - return self._operation(x, self._op._times, - self._op.inverse_times, self.target) + return self._operation(x, self._inverse_times, self.target) def _adjoint_times(self, x): - return self._operation(x, self._op._adjoint_times, - self._op._adjoint_inverse_times, self.domain) + return self._operation(x, self._adjoint_inverse_times, self.domain) def _inverse_times(self, x): - return self._operation(x, self._op._inverse_times, - self._op._times, self.domain) + return self._operation(x, self._times, self.domain) def _adjoint_inverse_times(self, x): - return self._operation(x, self._op._adjoint_inverse_times, - self._op._adjoint_times, self.target) + return self._operation(x, self._adjoint_times, self.target) diff --git a/nifty/sugar.py b/nifty/sugar.py index 08e30f9ca83ada3705effdbb5575de2a52c07dc2..da15940964f76a90e57c45bede6b019e3a538ffa 100644 --- a/nifty/sugar.py +++ b/nifty/sugar.py @@ -27,7 +27,6 @@ __all__ = ['PS_field', 'power_synthesize_special', 'create_power_field', 'create_power_operator', - 'generate_posterior_sample', 'create_composed_fft_operator'] @@ -251,49 +250,6 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None): spaces=space) -def generate_posterior_sample(mean, covariance): - """ Generates a posterior sample from a Gaussian distribution with given - mean and covariance - - This method generates samples by setting up the observation and - reconstruction of a mock signal in order to obtain residuals of the right - correlation which are added to the given mean. - - Parameters - ---------- - mean : Field - the mean of the posterior Gaussian distribution - covariance : WienerFilterCurvature - The posterior correlation structure consisting of a - response operator, noise covariance and prior signal covariance - - Returns - ------- - sample : Field - Returns the a sample from the Gaussian of given mean and covariance. - """ - - S = covariance.op.S - R = covariance.op.R - N = covariance.op.N - - power = sqrt(power_analyze(S.diagonal())) - mock_signal = power_synthesize(power, real_signal=True) - - noise = N.diagonal().weight(-1) - - mock_noise = Field.from_random(random_type="normal", domain=N.domain, - dtype=noise.dtype.type) - mock_noise *= sqrt(noise) - - mock_data = R(mock_signal) + mock_noise - - mock_j = R.adjoint_times(N.inverse_times(mock_data)) - mock_m = covariance.inverse_times(mock_j) - sample = mock_signal - mock_m + mean - return sample - - def create_composed_fft_operator(domain, codomain=None, all_to='other'): fft_op_list = []