Commit ef214849 authored by Jakob Knollmueller's avatar Jakob Knollmueller

Revert "added some distribution_strategy's for mpirun, fixed some errors in WienerFilterEnergy"

This reverts commit 03cc7f7e
parent 03cc7f7e
Pipeline #14965 passed with stage
in 7 minutes and 7 seconds
...@@ -10,78 +10,73 @@ rank = comm.rank ...@@ -10,78 +10,73 @@ rank = comm.rank
np.random.seed(42) np.random.seed(42)
# class AdjointFFTResponse(LinearOperator): class AdjointFFTResponse(LinearOperator):
# def __init__(self, FFT, R, default_spaces=None): def __init__(self, FFT, R, default_spaces=None):
# super(AdjointFFTResponse, self).__init__(default_spaces) super(AdjointFFTResponse, self).__init__(default_spaces)
# self._domain = FFT.target self._domain = FFT.target
# self._target = R.target self._target = R.target
# self.R = R self.R = R
# self.FFT = FFT self.FFT = FFT
#
# def _times(self, x, spaces=None): def _times(self, x, spaces=None):
# return self.R(self.FFT.adjoint_times(x)) return self.R(self.FFT.adjoint_times(x))
#
# def _adjoint_times(self, x, spaces=None): def _adjoint_times(self, x, spaces=None):
# return self.FFT(self.R.adjoint_times(x)) return self.FFT(self.R.adjoint_times(x))
# @property @property
# def domain(self): def domain(self):
# return self._domain return self._domain
#
# @property @property
# def target(self): def target(self):
# return self._target return self._target
#
# @property @property
# def unitary(self): def unitary(self):
# return False return False
#
if __name__ == "__main__": if __name__ == "__main__":
distribution_strategy = 'equal' distribution_strategy = 'not'
# Set up position space # Set up position space
signal_space = RGSpace([256,256]) s_space = RGSpace([128,128])
# s_space = HPSpace(32) # s_space = HPSpace(32)
harmonic_space = FFTOperator.get_default_codomain(signal_space)
fft = FFTOperator(harmonic_space, target=signal_space,
domain_dtype=np.complex, target_dtype=np.float)
# Define harmonic transformation and associated harmonic space # Define harmonic transformation and associated harmonic space
fft = FFTOperator(s_space)
h_space = fft.target[0]
# Setting up power space # Setting up power space
power_space = PowerSpace(harmonic_space, p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
distribution_strategy=distribution_strategy)
# Choosing the prior correlation structure and defining correlation operator # Choosing the prior correlation structure and defining correlation operator
power_spectrum = (lambda k: (42 / (k + 1) ** 3)) p_spec = (lambda k: (42 / (k + 1) ** 3))
S = create_power_operator(harmonic_space, power_spectrum=power_spectrum, S = create_power_operator(h_space, power_spectrum=p_spec,
distribution_strategy=distribution_strategy) distribution_strategy=distribution_strategy)
# Drawing a sample sh from the prior distribution in harmonic space # Drawing a sample sh from the prior distribution in harmonic space
sp = Field(power_space, val=power_spectrum, sp = Field(p_space, val=p_spec,
distribution_strategy=distribution_strategy) distribution_strategy=distribution_strategy)
sh = sp.power_synthesize(real_signal=True) sh = sp.power_synthesize(real_signal=True)
ss = fft(sh) ss = fft.adjoint_times(sh)
# Choosing the measurement instrument # Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.05) # Instrument = SmoothingOperator(s_space, sigma=0.05)
mask = Field(signal_space, val=1, Instrument = DiagonalOperator(s_space, diagonal=1.)
distribution_strategy=distribution_strategy) # Instrument._diagonal.val[200:400, 200:400] = 0
mask.val[63:127, 63:127] = 0.
Instrument = DiagonalOperator(signal_space, diagonal=mask)
#Adding a harmonic transformation to the instrument #Adding a harmonic transformation to the instrument
R = ComposedOperator([fft, Instrument], default_spaces=[0, 0]) R = AdjointFFTResponse(fft, Instrument)
signal_to_noise = 1. signal_to_noise = 1.
N = DiagonalOperator(signal_space, diagonal=ss.var()/signal_to_noise, N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
bare=True, n = Field.from_random(domain=s_space,
distribution_strategy=distribution_strategy)
n = Field.from_random(domain=signal_space,
random_type='normal', random_type='normal',
std=ss.std()/np.sqrt(signal_to_noise), std=ss.std()/np.sqrt(signal_to_noise),
mean=0, distribution_strategy=distribution_strategy) mean=0)
# Creating the mock data # Creating the mock data
d = R(sh) + n d = R(sh) + n
...@@ -100,50 +95,34 @@ if __name__ == "__main__": ...@@ -100,50 +95,34 @@ if __name__ == "__main__":
minimizer = RelaxedNewton(convergence_tolerance=0, minimizer = RelaxedNewton(convergence_tolerance=0,
iteration_limit=1, iteration_limit=1,
callback=convergence_measure) callback=convergence_measure)
#
minimizer = VL_BFGS(convergence_tolerance=0, # minimizer = VL_BFGS(convergence_tolerance=0,
iteration_limit=500, # iteration_limit=50,
callback=convergence_measure, # callback=convergence_measure,
max_history_length=3) # max_history_length=3)
#
inverter = ConjugateGradient(convergence_level=3, inverter = ConjugateGradient(convergence_level=3,
convergence_tolerance=1e-5, convergence_tolerance=1e-5,
preconditioner=None) preconditioner=None)
# Setting starting position # Setting starting position
m0 = Field(harmonic_space, val=.0, m0 = Field(h_space, val=.0)
distribution_strategy=distribution_strategy)
# 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)
inverter=inverter)
D0 = energy.curvature D0 = energy.curvature
# Solving the problem analytically # Solving the problem analytically
m0 = D0.inverse_times(j) m0 = D0.inverse_times(j)
# solution, convergence = minimizer(energy)
# m0 = solution.position sample_variance = Field(sh.domain,val=0. + 0j)
m0_s = Field(signal_space,val=fft(m0).val.get_full_data().real) sample_mean = Field(sh.domain,val=0. + 0j)
plotter = plotting.RG2DPlotter() # sampling the uncertainty map
plotter.title = 'mock_signal.html'; n_samples = 1
plotter(ss) for i in range(n_samples):
plotter.title = 'data.html' sample = sugar.generate_posterior_sample(m0,D0)
plotter(Field(signal_space, sample_variance += sample**2
val=Instrument.adjoint_times(d).val.get_full_data()\ sample_mean += sample
.reshape(signal_space.shape))) variance = sample_variance/n_samples - (sample_mean/n_samples)
plotter.title = 'map.html'; plotter(m0_s)
#
# sample_variance = Field(sh.domain,val=0. + 0j,
# distribution_strategy=distribution_strategy)
# sample_mean = Field(sh.domain,val=0. + 0j,
# distribution_strategy=distribution_strategy)
#
# # sampling the uncertainty map
# n_samples = 1
# for i in range(n_samples):
# sample = sugar.generate_posterior_sample(m0,D0)
# sample_variance += sample**2
# sample_mean += sample
# variance = sample_variance/n_samples - (sample_mean/n_samples)**2
...@@ -29,7 +29,7 @@ class WienerFilterEnergy(Energy): ...@@ -29,7 +29,7 @@ class WienerFilterEnergy(Energy):
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=position, d=self.d, R=self.R, N=self.N, return self.__class__(position=position, d=self.d, R=self.R, N=self.N,
S=self.S, inverter=self.inverter) S=self.S, inverter=self.inverter)
...@@ -47,9 +47,8 @@ class WienerFilterEnergy(Energy): ...@@ -47,9 +47,8 @@ class WienerFilterEnergy(Energy):
@property @property
@memo @memo
def curvature(self): def curvature(self):
return WienerFilterCurvature(R=self.R, N=self.N, S=self.S, return WienerFilterCurvature(R=self.R, N=self.N, S=self.S)
inverter=self.inverter)
@property
@memo @memo
def _Dx(self): def _Dx(self):
return self.curvature(self.position) return self.curvature(self.position)
......
...@@ -71,8 +71,7 @@ class InvertibleOperatorMixin(object): ...@@ -71,8 +71,7 @@ class InvertibleOperatorMixin(object):
def _times(self, x, spaces, x0=None): def _times(self, x, spaces, x0=None):
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)
distribution_strategy=x.distribution_strategy)
(result, convergence) = self.__inverter(A=self.inverse_times, (result, convergence) = self.__inverter(A=self.inverse_times,
b=x, b=x,
...@@ -81,8 +80,7 @@ class InvertibleOperatorMixin(object): ...@@ -81,8 +80,7 @@ class InvertibleOperatorMixin(object):
def _adjoint_times(self, x, spaces, x0=None): def _adjoint_times(self, x, spaces, x0=None):
if x0 is None: if x0 is None:
x0 = Field(self.domain, val=0., dtype=x.dtype, x0 = Field(self.domain, val=0., dtype=x.dtype)
distribution_strategy=x.distribution_strategy)
(result, convergence) = self.__inverter(A=self.adjoint_inverse_times, (result, convergence) = self.__inverter(A=self.adjoint_inverse_times,
b=x, b=x,
...@@ -91,8 +89,7 @@ class InvertibleOperatorMixin(object): ...@@ -91,8 +89,7 @@ class InvertibleOperatorMixin(object):
def _inverse_times(self, x, spaces, x0=None): def _inverse_times(self, x, spaces, x0=None):
if x0 is None: if x0 is None:
x0 = Field(self.domain, val=0., dtype=x.dtype, x0 = Field(self.domain, val=0., dtype=x.dtype)
distribution_strategy=x.distribution_strategy)
(result, convergence) = self.__inverter(A=self.times, (result, convergence) = self.__inverter(A=self.times,
b=x, b=x,
...@@ -101,8 +98,7 @@ class InvertibleOperatorMixin(object): ...@@ -101,8 +98,7 @@ class InvertibleOperatorMixin(object):
def _adjoint_inverse_times(self, x, spaces, x0=None): def _adjoint_inverse_times(self, x, spaces, x0=None):
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)
distribution_strategy=x.distribution_strategy)
(result, convergence) = self.__inverter(A=self.adjoint_times, (result, convergence) = self.__inverter(A=self.adjoint_times,
b=x, b=x,
......
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