Commit 58378e58 authored by Jakob Knollmueller's avatar Jakob Knollmueller
Browse files

added possibility to specify inverters to energy classes

parent d28da9fe
Pipeline #14421 passed with stage
in 5 minutes and 26 seconds
...@@ -95,6 +95,8 @@ if __name__ == "__main__": ...@@ -95,6 +95,8 @@ if __name__ == "__main__":
# Creating the mock data # Creating the mock data
d = R(sh) + n d = R(sh) + n
# The information source
j = R.adjoint_times(N.inverse_times(d))
realized_power = log(sh.power_analyze(logarithmic=p_space.config["logarithmic"], realized_power = log(sh.power_analyze(logarithmic=p_space.config["logarithmic"],
nbin=p_space.config["nbin"])) nbin=p_space.config["nbin"]))
data_power = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"], data_power = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"],
...@@ -114,11 +116,15 @@ if __name__ == "__main__": ...@@ -114,11 +116,15 @@ if __name__ == "__main__":
convergence_level=2, convergence_level=2,
iteration_limit=3, iteration_limit=3,
callback=convergence_measure) callback=convergence_measure)
minimizer2 = VL_BFGS(convergence_tolerance=0, minimizer2 = VL_BFGS(convergence_tolerance=0,
iteration_limit=7, iteration_limit=7,
callback=convergence_measure, callback=convergence_measure,
max_history_length=3) max_history_length=3)
inverter = ConjugateGradient(convergence_level=1,
convergence_tolerance=10e-4,
preconditioner=None)
# Setting starting position # Setting starting position
flat_power = Field(p_space,val=10e-8) flat_power = Field(p_space,val=10e-8)
m0 = flat_power.power_synthesize(real_signal=True) m0 = flat_power.power_synthesize(real_signal=True)
...@@ -132,15 +138,12 @@ if __name__ == "__main__": ...@@ -132,15 +138,12 @@ if __name__ == "__main__":
distribution_strategy=distribution_strategy) distribution_strategy=distribution_strategy)
# Initializing the nonlinear Wiener Filter energy # Initializing the nonlinear Wiener Filter energy
map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0) map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0, inverter=inverter)
# Minimization with chosen minimizer # Solving the Wiener Filter analytically
map_energy = map_energy.analytic_solution()
# Updating parameters for correlation structure reconstruction
m0 = map_energy.position
D0 = map_energy.curvature D0 = map_energy.curvature
m0 = D0.inverse_times(j)
# Initializing the power energy with updated parameters # Initializing the power energy with updated parameters
power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=10., samples=3) power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, sigma=10., samples=3, inverter=inverter)
(power_energy, convergence) = minimizer1(power_energy) (power_energy, convergence) = minimizer1(power_energy)
......
...@@ -80,6 +80,7 @@ if __name__ == "__main__": ...@@ -80,6 +80,7 @@ if __name__ == "__main__":
# Creating the mock data # Creating the mock data
d = R(sh) + n d = R(sh) + n
j = R.adjoint_times(N.inverse_times(d))
# Choosing the minimization strategy # Choosing the minimization strategy
...@@ -101,66 +102,27 @@ if __name__ == "__main__": ...@@ -101,66 +102,27 @@ if __name__ == "__main__":
# max_history_length=3) # max_history_length=3)
# #
inverter = ConjugateGradient(convergence_level=3,
convergence_tolerance=10e-5,
preconditioner=None)
# Setting starting position # Setting starting position
m0 = Field(h_space, val=.0) m0 = Field(h_space, val=.0)
# Initializing the Wiener Filter energy # Initializing the Wiener Filter energy
energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S) energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S, inverter=inverter)
D0 = energy.curvature
# Solving the problem analytically # Solving the problem analytically
solution = energy.analytic_solution() m0 = D0.inverse_times(j)
sample_variance = Field(sh.power_analyze(logarithmic=False).domain,val=0. + 0j) sample_variance = Field(sh.power_analyze(logarithmic=False).domain,val=0. + 0j)
sample_mean = Field(sh.domain,val=0. + 0j) sample_mean = Field(sh.domain,val=0. + 0j)
n_samples = 200
# sampling the uncertainty map
n_samples = 1
for i in range(n_samples): for i in range(n_samples):
sample = sugar.generate_posterior_sample(solution.position,solution.curvature) sample = sugar.generate_posterior_sample(m0,D0)
sample_variance += sample.power_analyze(logarithmic=False) sample_variance += sample**2
sample_mean += sample sample_mean += sample
sample_variance = sample_variance/n_samples - (sample_mean/n_samples).power_analyze(logarithmic=False) variance = sample_variance/n_samples - (sample_mean/n_samples)
D = HarmonicPropagatorOperator(S=S, N=N, R=Instrument)
class DiagonalProber(DiagonalProberMixin, Prober):
pass
diagProber = DiagonalProber(domain=sh.domain)
diagProber(D)
diag = diagProber.diagonal
pr_diag = diag.power_analyze(logarithmic=False)**0.5
# Solving the problem with chosen minimization strategy
(energy, convergence) = minimizer(energy)
# Transforming fields to position space for plotting
ss = fft.adjoint_times(sh)
m = fft.adjoint_times(energy.position)
m_wf = fft.adjoint_times(solution.position)
# Plotting
d_data = d.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html')
ss_data = ss.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=ss_data)], filename='ss.html')
sh_data = sh.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=sh_data)], filename='sh.html')
m_data = m.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=m_data)], filename='map.html')
m_wf_data = m_wf.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=m_wf_data)], filename='map_wf.html')
...@@ -67,11 +67,3 @@ if __name__ == "__main__": ...@@ -67,11 +67,3 @@ if __name__ == "__main__":
D = PropagatorOperator(S=S, N=N, R=R) D = PropagatorOperator(S=S, N=N, R=R)
m = D(j) m = D(j)
d_data = d.val.get_full_data().real
m_data = m.val.get_full_data().real
ss_data = ss.val.get_full_data().real
# if rank == 0:
# pl.plot([go.Heatmap(z=d_data)], filename='data.html')
# pl.plot([go.Heatmap(z=m_data)], filename='map.html')
# pl.plot([go.Heatmap(z=ss_data)], filename='map_orig.html')
...@@ -42,10 +42,13 @@ class CriticalPowerEnergy(Energy): ...@@ -42,10 +42,13 @@ class CriticalPowerEnergy(Energy):
The contribution from the map with or without uncertainty. It is used The contribution from the map with or without uncertainty. It is used
to pass on the result of the costly sampling during the minimization. to pass on the result of the costly sampling during the minimization.
default : None default : None
inverter : ConjugateGradient
The inversion strategy to invert the curvature and to generate samples.
default : None
""" """
def __init__(self, position, m, D=None, alpha =1.0, q=0., sigma=0., def __init__(self, position, m, D=None, alpha =1.0, q=0., sigma=0.,
logarithmic = True, samples=3, w=None): logarithmic = True, samples=3, w=None, inverter=None):
super(CriticalPowerEnergy, self).__init__(position = position) super(CriticalPowerEnergy, self).__init__(position = position)
self.m = m self.m = m
self.D = D self.D = D
...@@ -56,11 +59,13 @@ class CriticalPowerEnergy(Energy): ...@@ -56,11 +59,13 @@ class CriticalPowerEnergy(Energy):
self.T = SmoothnessOperator(domain=self.position.domain[0], sigma=self.sigma, self.T = SmoothnessOperator(domain=self.position.domain[0], sigma=self.sigma,
logarithmic=logarithmic) logarithmic=logarithmic)
self.rho = self.position.domain[0].rho self.rho = self.position.domain[0].rho
self.inverter = inverter
self.w = w self.w = w
if self.w is None: if self.w is None:
self.w = self._calculate_w(self.m, self.D, self.samples) self.w = self._calculate_w(self.m, self.D, self.samples)
self.theta = (exp(-self.position) * (self.q + self.w / 2.)) self.theta = (exp(-self.position) * (self.q + self.w / 2.))
def at(self, position): def at(self, position):
return self.__class__(position, self.m, D=self.D, return self.__class__(position, self.m, D=self.D,
alpha =self.alpha, alpha =self.alpha,
...@@ -85,14 +90,14 @@ class CriticalPowerEnergy(Energy): ...@@ -85,14 +90,14 @@ class CriticalPowerEnergy(Energy):
@property @property
def curvature(self): def curvature(self):
curvature = CriticalPowerCurvature(theta=self.theta.weight(-1), T=self.T) curvature = CriticalPowerCurvature(theta=self.theta.weight(-1), T=self.T, inverter=self.inverter)
return curvature return curvature
def _calculate_w(self, m, D, samples): def _calculate_w(self, m, D, samples):
w = Field(domain=self.position.domain, val=0. , dtype=m.dtype) w = Field(domain=self.position.domain, val=0. , dtype=m.dtype)
if D is not None: if D is not None:
for i in range(samples): for i in range(samples):
posterior_sample = generate_posterior_sample(m, D) posterior_sample = generate_posterior_sample(m, D, inverter = self.inverter)
projected_sample = posterior_sample.power_analyze( projected_sample = posterior_sample.power_analyze(
logarithmic=self.position.domain[0].config["logarithmic"], logarithmic=self.position.domain[0].config["logarithmic"],
nbin= self.position.domain[0].config["nbin"]) nbin= self.position.domain[0].config["nbin"])
......
...@@ -22,12 +22,13 @@ class WienerFilterEnergy(Energy): ...@@ -22,12 +22,13 @@ class WienerFilterEnergy(Energy):
The prior signal covariance in harmonic space. The prior signal covariance in harmonic space.
""" """
def __init__(self, position, d, R, N, S): def __init__(self, position, d, R, N, S, inverter=None):
super(WienerFilterEnergy, self).__init__(position = position) super(WienerFilterEnergy, self).__init__(position = position)
self.d = d self.d = d
self.R = R self.R = R
self.N = N self.N = N
self.S = S self.S = S
self.inverter = inverter
def at(self, position): def at(self, position):
return self.__class__(position, self.d, self.R, self.N, self.S) return self.__class__(position, self.d, self.R, self.N, self.S)
...@@ -49,7 +50,7 @@ class WienerFilterEnergy(Energy): ...@@ -49,7 +50,7 @@ class WienerFilterEnergy(Energy):
@property @property
def curvature(self): def curvature(self):
curvature = WienerFilterCurvature(R=self.R, N=self.N, S=self.S) curvature = WienerFilterCurvature(R=self.R, N=self.N, S=self.S, inverter=self.inverter)
return curvature return curvature
@memo @memo
...@@ -57,10 +58,3 @@ class WienerFilterEnergy(Energy): ...@@ -57,10 +58,3 @@ class WienerFilterEnergy(Energy):
residual = self.d - self.R(self.position) residual = self.d - self.R(self.position)
return residual return residual
def analytic_solution(self):
D_inverse = self.curvature
j = self.R.adjoint_times(self.N.inverse_times(self.d))
new_position = D_inverse.inverse_times(j)
return self.at(new_position)
...@@ -72,7 +72,7 @@ class InvertibleOperatorMixin(object): ...@@ -72,7 +72,7 @@ class InvertibleOperatorMixin(object):
if x0 is None: if x0 is None:
x0 = Field(self.target, val=0., dtype=x.dtype) x0 = Field(self.target, val=0., dtype=x.dtype)
(result, convergence) = self.__inverter(A=self.inverse_times, (result, convergence) = self.__inverter(A=self.inverses_times,
b=x, b=x,
x0=x0) x0=x0)
return result return result
......
...@@ -20,7 +20,7 @@ from nifty import PowerSpace,\ ...@@ -20,7 +20,7 @@ from nifty import PowerSpace,\
Field,\ Field,\
DiagonalOperator,\ DiagonalOperator,\
sqrt sqrt
from nifty.minimization.conjugate_gradient import ConjugateGradient
__all__ = ['create_power_operator'] __all__ = ['create_power_operator']
...@@ -66,10 +66,36 @@ def create_power_operator(domain, power_spectrum, dtype=None, ...@@ -66,10 +66,36 @@ def create_power_operator(domain, power_spectrum, dtype=None,
f **= 2 f **= 2
return DiagonalOperator(domain, diagonal=f, bare=True) return DiagonalOperator(domain, diagonal=f, bare=True)
def generate_posterior_sample(mean, covariance): def generate_posterior_sample(mean, covariance, inverter = None):
""" 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
inverter : ConjugateGradient *optional*
the conjugate gradient used to invert the curvature for the Wiener filter.
default : None
Returns
-------
sample : Field
Returns the a sample from the Gaussian of given mean and covariance.
"""
S = covariance.S S = covariance.S
R = covariance.R R = covariance.R
N = covariance.N N = covariance.N
if inverter is None:
inverter = ConjugateGradient(preconditioner=S)
power = S.diagonal().power_analyze()**.5 power = S.diagonal().power_analyze()**.5
mock_signal = power.power_synthesize(real_signal=True) mock_signal = power.power_synthesize(real_signal=True)
......
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