From 03cc7f7e8edcc381eca505dbfc9c81feac394c1d Mon Sep 17 00:00:00 2001
From: "Knollmueller, Jakob (kjako)" <jakob@knollmueller.de>
Date: Mon, 17 Jul 2017 16:18:04 +0200
Subject: [PATCH] added some distribution_strategy's for mpirun, fixed some
errors in WienerFilterEnergy
---
demos/wiener_filter_via_hamiltonian.py | 141 ++++++++++--------
.../wiener_filter/wiener_filter_energy.py | 7 +-
.../invertible_operator_mixin.py | 12 +-
3 files changed, 93 insertions(+), 67 deletions(-)
diff --git a/demos/wiener_filter_via_hamiltonian.py b/demos/wiener_filter_via_hamiltonian.py
index 70f6a3d06..7bcf3ee21 100644
--- a/demos/wiener_filter_via_hamiltonian.py
+++ b/demos/wiener_filter_via_hamiltonian.py
@@ -10,73 +10,78 @@ rank = comm.rank
np.random.seed(42)
-class AdjointFFTResponse(LinearOperator):
- def __init__(self, FFT, R, default_spaces=None):
- super(AdjointFFTResponse, self).__init__(default_spaces)
- self._domain = FFT.target
- self._target = R.target
- self.R = R
- self.FFT = FFT
-
- def _times(self, x, spaces=None):
- return self.R(self.FFT.adjoint_times(x))
-
- def _adjoint_times(self, x, spaces=None):
- return self.FFT(self.R.adjoint_times(x))
- @property
- def domain(self):
- return self._domain
-
- @property
- def target(self):
- return self._target
-
- @property
- def unitary(self):
- return False
-
+# class AdjointFFTResponse(LinearOperator):
+# def __init__(self, FFT, R, default_spaces=None):
+# super(AdjointFFTResponse, self).__init__(default_spaces)
+# self._domain = FFT.target
+# self._target = R.target
+# self.R = R
+# self.FFT = FFT
+#
+# def _times(self, x, spaces=None):
+# return self.R(self.FFT.adjoint_times(x))
+#
+# def _adjoint_times(self, x, spaces=None):
+# return self.FFT(self.R.adjoint_times(x))
+# @property
+# def domain(self):
+# return self._domain
+#
+# @property
+# def target(self):
+# return self._target
+#
+# @property
+# def unitary(self):
+# return False
+#
if __name__ == "__main__":
- distribution_strategy = 'not'
+ distribution_strategy = 'equal'
# Set up position space
- s_space = RGSpace([128,128])
+ signal_space = RGSpace([256,256])
# 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
- fft = FFTOperator(s_space)
- h_space = fft.target[0]
# Setting up power space
- p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
+ power_space = PowerSpace(harmonic_space,
+ distribution_strategy=distribution_strategy)
# Choosing the prior correlation structure and defining correlation operator
- p_spec = (lambda k: (42 / (k + 1) ** 3))
+ power_spectrum = (lambda k: (42 / (k + 1) ** 3))
- S = create_power_operator(h_space, power_spectrum=p_spec,
+ S = create_power_operator(harmonic_space, power_spectrum=power_spectrum,
distribution_strategy=distribution_strategy)
# Drawing a sample sh from the prior distribution in harmonic space
- sp = Field(p_space, val=p_spec,
+ sp = Field(power_space, val=power_spectrum,
distribution_strategy=distribution_strategy)
sh = sp.power_synthesize(real_signal=True)
- ss = fft.adjoint_times(sh)
-
+ ss = fft(sh)
# Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.05)
- Instrument = DiagonalOperator(s_space, diagonal=1.)
-# Instrument._diagonal.val[200:400, 200:400] = 0
+ mask = Field(signal_space, val=1,
+ distribution_strategy=distribution_strategy)
+ mask.val[63:127, 63:127] = 0.
+ Instrument = DiagonalOperator(signal_space, diagonal=mask)
#Adding a harmonic transformation to the instrument
- R = AdjointFFTResponse(fft, Instrument)
+ R = ComposedOperator([fft, Instrument], default_spaces=[0, 0])
signal_to_noise = 1.
- N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
- n = Field.from_random(domain=s_space,
+ N = DiagonalOperator(signal_space, diagonal=ss.var()/signal_to_noise,
+ bare=True,
+ distribution_strategy=distribution_strategy)
+ n = Field.from_random(domain=signal_space,
random_type='normal',
std=ss.std()/np.sqrt(signal_to_noise),
- mean=0)
+ mean=0, distribution_strategy=distribution_strategy)
# Creating the mock data
d = R(sh) + n
@@ -95,34 +100,50 @@ if __name__ == "__main__":
minimizer = RelaxedNewton(convergence_tolerance=0,
iteration_limit=1,
callback=convergence_measure)
- #
- # minimizer = VL_BFGS(convergence_tolerance=0,
- # iteration_limit=50,
- # callback=convergence_measure,
- # max_history_length=3)
- #
+
+ minimizer = VL_BFGS(convergence_tolerance=0,
+ iteration_limit=500,
+ callback=convergence_measure,
+ max_history_length=3)
+
inverter = ConjugateGradient(convergence_level=3,
convergence_tolerance=1e-5,
preconditioner=None)
# Setting starting position
- m0 = Field(h_space, val=.0)
+ m0 = Field(harmonic_space, val=.0,
+ distribution_strategy=distribution_strategy)
# Initializing the Wiener Filter energy
- energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S, inverter=inverter)
+ energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S,
+ inverter=inverter)
D0 = energy.curvature
# Solving the problem analytically
m0 = D0.inverse_times(j)
-
- sample_variance = Field(sh.domain,val=0. + 0j)
- sample_mean = Field(sh.domain,val=0. + 0j)
-
- # 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)
+ # solution, convergence = minimizer(energy)
+ # m0 = solution.position
+ m0_s = Field(signal_space,val=fft(m0).val.get_full_data().real)
+
+ plotter = plotting.RG2DPlotter()
+ plotter.title = 'mock_signal.html';
+ plotter(ss)
+ plotter.title = 'data.html'
+ plotter(Field(signal_space,
+ val=Instrument.adjoint_times(d).val.get_full_data()\
+ .reshape(signal_space.shape)))
+ 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
diff --git a/nifty/library/wiener_filter/wiener_filter_energy.py b/nifty/library/wiener_filter/wiener_filter_energy.py
index 8eec64c51..4b9319374 100644
--- a/nifty/library/wiener_filter/wiener_filter_energy.py
+++ b/nifty/library/wiener_filter/wiener_filter_energy.py
@@ -29,7 +29,7 @@ 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,
S=self.S, inverter=self.inverter)
@@ -47,8 +47,9 @@ class WienerFilterEnergy(Energy):
@property
@memo
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
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 af952177e..483840cc6 100644
--- a/nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py
+++ b/nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py
@@ -71,7 +71,8 @@ class InvertibleOperatorMixin(object):
def _times(self, x, spaces, x0=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,
b=x,
@@ -80,7 +81,8 @@ class InvertibleOperatorMixin(object):
def _adjoint_times(self, x, spaces, x0=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,
b=x,
@@ -89,7 +91,8 @@ class InvertibleOperatorMixin(object):
def _inverse_times(self, x, spaces, x0=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,
b=x,
@@ -98,7 +101,8 @@ class InvertibleOperatorMixin(object):
def _adjoint_inverse_times(self, x, spaces, x0=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,
b=x,
--
GitLab