From ac8f1fe281df90a13b1893df4d5472da75de1b9c Mon Sep 17 00:00:00 2001 From: Theo Steininger Date: Wed, 12 Jul 2017 00:11:41 +0200 Subject: [PATCH] Refactored critical power files. --- demos/critical_filtering.py | 9 +- .../critical_filter/critical_power_energy.py | 103 +++++++++--------- .../wiener_filter/wiener_filter_energy.py | 5 +- .../invertible_operator_mixin.py | 3 +- .../smoothness_operator.py | 38 ++++--- nifty/sugar.py | 15 +-- 6 files changed, 88 insertions(+), 85 deletions(-) diff --git a/demos/critical_filtering.py b/demos/critical_filtering.py index 8ca2cf92..9a9e2103 100644 --- a/demos/critical_filtering.py +++ b/demos/critical_filtering.py @@ -123,9 +123,6 @@ if __name__ == "__main__": callback=convergence_measure, max_history_length=3) - inverter = ConjugateGradient(convergence_level=1, - convergence_tolerance=10e-4, - preconditioner=None) # Setting starting position flat_power = Field(p_space,val=10e-8) m0 = flat_power.power_synthesize(real_signal=True) @@ -144,7 +141,7 @@ if __name__ == "__main__": D0 = map_energy.curvature m0 = D0.inverse_times(j) # Initializing the power energy with updated parameters - power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=10., samples=3, inverter=inverter) + power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, smoothness_prior=10., samples=3) (power_energy, convergence) = minimizer1(power_energy) @@ -153,6 +150,8 @@ if __name__ == "__main__": t0.val = power_energy.position.val.real # Plotting current estimate - plot_parameters(m0,t0,log(sp), data_power) + print i + if i%50 == 0: + plot_parameters(m0,t0,log(sp), data_power) diff --git a/nifty/library/critical_filter/critical_power_energy.py b/nifty/library/critical_filter/critical_power_energy.py index bddac6bb..7848aee5 100644 --- a/nifty/library/critical_filter/critical_power_energy.py +++ b/nifty/library/critical_filter/critical_power_energy.py @@ -1,3 +1,5 @@ +import numpy as np + from nifty.energies.energy import Energy from nifty.operators.smoothness_operator import SmoothnessOperator from nifty.library.critical_filter import CriticalPowerCurvature @@ -10,9 +12,10 @@ from nifty import Field, exp class CriticalPowerEnergy(Energy): """The Energy of the power spectrum according to the critical filter. - It describes the energy of the logarithmic amplitudes of the power spectrum of - a Gaussian random field under reconstruction uncertainty with smoothness and - inverse gamma prior. It is used to infer the correlation structure of a correlated signal. + It describes the energy of the logarithmic amplitudes of the power spectrum + of a Gaussian random field under reconstruction uncertainty with smoothness + and inverse gamma prior. It is used to infer the correlation structure of a + correlated signal. Parameters ---------- @@ -25,20 +28,21 @@ class CriticalPowerEnergy(Energy): If not specified, the map is assumed to be no reconstruction. default : None alpha : float - The spectral prior of the inverse gamma distribution. 1.0 corresponds to - non-informative. + The spectral prior of the inverse gamma distribution. 1.0 corresponds + to non-informative. default : 1.0 q : float - The cutoff parameter of the inverse gamma distribution. 0.0 corresponds to - non-informative. + The cutoff parameter of the inverse gamma distribution. 0.0 corresponds + to non-informative. + default : 0.0 + smoothness_prior : float + Controls the strength of the smoothness prior default : 0.0 - sigma : float - The parameter of the smoothness prior. - default : ??? None? ??????? logarithmic : boolean Whether smoothness acts on linear or logarithmic scale. samples : integer - Number of samples used for the estimation of the uncertainty corrections. + Number of samples used for the estimation of the uncertainty + corrections. default : 3 w : Field The contribution from the map with or without uncertainty. It is used @@ -49,76 +53,72 @@ class CriticalPowerEnergy(Energy): default : None """ - def __init__(self, position, m, D=None, alpha =1.0, q=0., sigma=0., - logarithmic = True, samples=3, w=None, inverter=None): - super(CriticalPowerEnergy, self).__init__(position = position) + def __init__(self, position, m, D=None, alpha=1.0, q=0., + smoothness_prior=0., logarithmic=True, samples=3, w=None): + super(CriticalPowerEnergy, self).__init__(position=position) self.m = m self.D = D self.samples = samples - self.sigma = sigma + self.smoothness_prior = np.float(smoothness_prior) self.alpha = Field(self.position.domain, val=alpha) self.q = Field(self.position.domain, val=q) - self.T = SmoothnessOperator(domain=self.position.domain[0], sigma=self.sigma, + self.T = SmoothnessOperator(domain=self.position.domain[0], + strength=self.smoothness_prior, logarithmic=logarithmic) self.rho = self.position.domain[0].rho - self.inverter = inverter - self.w = w - if self.w is None: - self.w = self._calculate_w(self.m, self.D, self.samples) - self.theta = (exp(-self.position) * (self.q + self.w / 2.)) - + self._w = w if w is not None else None def at(self, position): - return self.__class__(position, self.m, D=self.D, - alpha =self.alpha, - q=self.q, - sigma=self.sigma, w=self.w, - samples=self.samples) + return self.__class__(position, self.m, D=self.D, alpha=self.alpha, + q=self.q, smoothness_prior=self.smoothness_prior, + w=self.w, samples=self.samples) @property def value(self): - energy = self._theta.vdot(Field(self.position.domain,val=1.), bare= True) + energy = self._theta.sum() energy += self.position.vdot(self._rho_prime, bare=True) energy += 0.5 * self.position.vdot(self._Tt) return energy.real @property def gradient(self): - gradient = - self._theta.weight(-1) + gradient = -self._theta.weight(-1) gradient += (self._rho_prime).weight(-1) - gradient += self._Tt + gradient += self._Tt gradient.val = gradient.val.real return gradient @property def curvature(self): - curvature = CriticalPowerCurvature(theta=self._theta.weight(-1), T=self.T, - inverter=self.inverter) + curvature = CriticalPowerCurvature(theta=self._theta.weight(-1), + T=self.T) return curvature - def _calculate_w(self, m, D, samples): - w = Field(domain=self.position.domain, val=0. , dtype=m.dtype) - if D is not None: - for i in range(samples): - posterior_sample = generate_posterior_sample(m, D, inverter = self.inverter) - projected_sample = posterior_sample.power_analyze( - logarithmic=self.position.domain[0].config["logarithmic"], - nbin= self.position.domain[0].config["nbin"]) - w += (projected_sample) * self.rho - w /= float(samples) - else: - w = m.power_analyze( - logarithmic=self.position.domain[0].config["logarithmic"], - nbin=self.position.domain[0].config["nbin"]) - w *= self.rho - - return w - + @property + def w(self): + if self._w is None: + w = Field(domain=self.position.domain, val=0., dtype=self.m.dtype) + if self.D is not None: + for i in range(self.samples): + posterior_sample = generate_posterior_sample( + self.m, self.D) + projected_sample = posterior_sample.power_analyze( + logarithmic=self.position.domain[0].config["logarithmic"], + nbin=self.position.domain[0].config["nbin"]) + w += (projected_sample) * self.rho + w /= float(self.samples) + else: + w = self.m.power_analyze( + logarithmic=self.position.domain[0].config["logarithmic"], + nbin=self.position.domain[0].config["nbin"]) + w *= self.rho + self._w = w + return self._w @property @memo def _theta(self): - return (exp(-self.position) * (self.q + self.w / 2.)) + return exp(-self.position) * (self.q + self.w / 2.) @property @memo @@ -129,4 +129,3 @@ class CriticalPowerEnergy(Energy): @memo def _Tt(self): return self.T(self.position) - diff --git a/nifty/library/wiener_filter/wiener_filter_energy.py b/nifty/library/wiener_filter/wiener_filter_energy.py index 798319e5..8eec64c5 100644 --- a/nifty/library/wiener_filter/wiener_filter_energy.py +++ b/nifty/library/wiener_filter/wiener_filter_energy.py @@ -29,7 +29,6 @@ class WienerFilterEnergy(Energy): self.R = R self.N = N self.S = S - self.inverter = inverter def at(self, position): return self.__class__(position=position, d=self.d, R=self.R, N=self.N, @@ -48,10 +47,8 @@ class WienerFilterEnergy(Energy): @property @memo def curvature(self): - return WienerFilterCurvature(R=self.R, N=self.N, S=self.S, - inverter=self.inverter) + return WienerFilterCurvature(R=self.R, N=self.N, S=self.S) - @property @memo def _Dx(self): return self.curvature(self.position) diff --git a/nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py b/nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py index e134656b..ad41ba0b 100644 --- a/nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py +++ b/nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py @@ -38,7 +38,8 @@ class InvertibleOperatorMixin(object): (default: ConjugateGradient) preconditioner : LinearOperator - Preconditions the minimizaion problem + Preconditioner that is used by ConjugateGraduent if no minimizer was + given. Attributes ---------- diff --git a/nifty/operators/smoothness_operator/smoothness_operator.py b/nifty/operators/smoothness_operator/smoothness_operator.py index ae7efd29..e9e35453 100644 --- a/nifty/operators/smoothness_operator/smoothness_operator.py +++ b/nifty/operators/smoothness_operator/smoothness_operator.py @@ -5,14 +5,17 @@ from nifty.operators.laplace_operator import LaplaceOperator class SmoothnessOperator(EndomorphicOperator): - """An operator measuring the smoothness on an irregular grid with respect to some scale. + """An operator measuring the smoothness on an irregular grid with respect + to some scale. - This operator applies the irregular LaplaceOperator and its adjoint to some Field over a - PowerSpace which corresponds to its smoothness and weights the result with a scale parameter sigma. - It is used in the smoothness prior terms of the CriticalPowerEnergy. For this purpose we - use free boundary conditions in the LaplaceOperator, having no curvature at both ends. - In addition the first entry is ignored as well, corresponding to the overall mean of the map. - The mean is therefore not considered in the smoothness prior. + This operator applies the irregular LaplaceOperator and its adjoint to some + Field over a PowerSpace which corresponds to its smoothness and weights the + result with a scale parameter sigma. It is used in the smoothness prior + terms of the CriticalPowerEnergy. For this purpose we use free boundary + conditions in the LaplaceOperator, having no curvature at both ends. In + addition the first entry is ignored as well, corresponding to the overall + mean of the map. The mean is therefore not considered in the smoothness + prior. Parameters @@ -26,7 +29,8 @@ class SmoothnessOperator(EndomorphicOperator): # ---Overwritten properties and methods--- - def __init__(self, domain, sigma, logarithmic=True, default_spaces=None): + def __init__(self, domain, strength=1., logarithmic=True, + default_spaces=None): super(SmoothnessOperator, self).__init__(default_spaces=default_spaces) @@ -37,10 +41,10 @@ class SmoothnessOperator(EndomorphicOperator): if not isinstance(self.domain[0], PowerSpace): raise TypeError("The domain must contain exactly one PowerSpace.") - if sigma <= 0: + if strength <= 0: raise ValueError("ERROR: invalid sigma.") - self._sigma = sigma + self._strength = strength self._laplace = LaplaceOperator(domain=self.domain, logarithmic=logarithmic) @@ -68,11 +72,17 @@ class SmoothnessOperator(EndomorphicOperator): return False def _times(self, x, spaces): - res = self._laplace.adjoint_times(self._laplace(x, spaces), spaces) - return (1./self.sigma)**2*res + if self._strength != 0: + result = self._laplace.adjoint_times(self._laplace(x, spaces), + spaces) + result *= self._strength**2 + else: + result = x.copy_empty() + result.val[:] = 0 + return result # ---Added properties and methods--- @property - def sigma(self): - return self._sigma + def strength(self): + return self._strength diff --git a/nifty/sugar.py b/nifty/sugar.py index 5114a949..b1835af9 100644 --- a/nifty/sugar.py +++ b/nifty/sugar.py @@ -64,11 +64,13 @@ def create_power_operator(domain, power_spectrum, dtype=None, return DiagonalOperator(domain, diagonal=f, bare=True) -def generate_posterior_sample(mean, covariance, inverter=None): - """ Generates a posterior sample from a Gaussian distribution with given mean and covariance +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. + 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 ---------- @@ -77,9 +79,6 @@ def generate_posterior_sample(mean, covariance, inverter=None): covariance : WienerFilterCurvature The posterior correlation structure consisting of a response operator, noise covariance and prior signal covariance - inverter : ConjugateGradient *optional* - the conjugate gradient used to invert the curvature for the Wiener filter. - default : None Returns ------- @@ -91,8 +90,6 @@ def generate_posterior_sample(mean, covariance, inverter=None): S = covariance.S R = covariance.R N = covariance.N - if inverter is None: - inverter = ConjugateGradient(preconditioner=S) power = S.diagonal().power_analyze()**.5 mock_signal = power.power_synthesize(real_signal=True) -- GitLab