Commit 36b815ce authored by Reimar Heinrich Leike's avatar Reimar Heinrich Leike

Merge branch 'NIFTy_4' into yango_minimizer

This merge is done in order to make the yango minimizer and the krylov sampling available in the same branch while still keeping the yango minimizer out of the master branch for NIFTy
parents e1a5bfc0 c305952a
Pipeline #27919 passed with stage
in 1 minute and 39 seconds
...@@ -3,89 +3,81 @@ import numpy as np ...@@ -3,89 +3,81 @@ import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from nifty4.sugar import create_power_operator from nifty4.sugar import create_power_operator
np.random.seed(42)
x_space = ift.RGSpace(1024) x_space = ift.RGSpace(1024)
h_space = x_space.get_default_codomain() h_space = x_space.get_default_codomain()
d_space = x_space d_space = x_space
N_hat = ift.Field(d_space,10.) N_hat = np.full(d_space.shape, 10.)
N_hat.val[400:450] = 0.0001 N_hat[400:450] = 0.0001
N = ift.DiagonalOperator(N_hat, d_space) N_hat = ift.Field.from_global_data(d_space, N_hat)
N = ift.DiagonalOperator(N_hat)
FFT = ift.HarmonicTransformOperator(h_space, target=x_space) FFT = ift.HarmonicTransformOperator(h_space, x_space)
R = ift.ScalingOperator(1., x_space) R = ift.ScalingOperator(1., x_space)
ampspec = lambda k : 1./(1.+k**2.)
S = ift.ScalingOperator(1.,h_space)
def ampspec(k): return 1. / (1. + k**2.)
S = ift.ScalingOperator(1., h_space)
A = create_power_operator(h_space, ampspec) A = create_power_operator(h_space, ampspec)
s_h = S.draw_sample() s_h = S.draw_sample()
sky = FFT*A sky = FFT * A
s_x = sky(s_h) s_x = sky(s_h)
n = N.draw_sample() n = N.draw_sample()
d = R(s_x) + n d = R(s_x) + n
R_p = R*FFT*A R_p = R * FFT * A
j = R_p.adjoint(N.inverse(d)) j = R_p.adjoint(N.inverse(d))
D_inv = R_p.adjoint*N.inverse*R_p + S.inverse D_inv = ift.SandwichOperator(R_p, N.inverse) + S.inverse
def sample(D_inv, S, j, N_samps, N_iter):
space = D_inv.domain
x = ift.Field.zeros(space)
r = j.copy()
p = r.copy()
d = p.vdot(D_inv(p))
y = []
for i in range(N_samps):
y += [S.draw_sample()]
for k in range(1,1+N_iter):
gamma = r.vdot(r)/d
if gamma == 0.:
break
x += gamma*p
print(p.vdot(D_inv(j)))
for i in range(N_samps):
y[i] -= p.vdot(D_inv(y[i]))*p/d
y[i] += np.random.randn()/np.sqrt(d)*p
#r_new = j - D_inv(x)
r_new = r - gamma*D_inv(p)
beta = r_new.vdot(r_new)/(r.vdot(r))
r = r_new
p = r + beta*p
d = p.vdot(D_inv(p))
if d == 0.:
break;
return x, y
N_samps = 200 N_samps = 200
N_iter = 10 N_iter = 100
m,samps = sample(D_inv, S, j, N_samps, N_iter) IC = ift.GradientNormController(tol_abs_gradnorm=1e-3, iteration_limit=N_iter)
m, samps = ift.library.generate_krylov_samples(D_inv, S, j, N_samps, IC)
m_x = sky(m) m_x = sky(m)
IC = ift.GradientNormController(iteration_limit = N_iter)
inverter = ift.ConjugateGradient(IC) inverter = ift.ConjugateGradient(IC)
curv = ift.library.WienerFilterCurvature(S=S,N=N,R=R_p, inverter=inverter) curv = ift.library.WienerFilterCurvature(S=S, N=N, R=R_p, inverter=inverter)
samps_old = [] samps_old = [curv.draw_sample(from_inverse=True) for i in range(N_samps)]
for i in range(N_samps):
samps_old += [curv.draw_sample(from_inverse=True)]
plt.plot(d.val,'o', label="data") plt.plot(d.to_global_data(), '+', label="data", alpha=.5)
plt.plot(s_x.val,label="original") plt.plot(s_x.to_global_data(), label="original")
plt.plot(m_x.val, label="reconstruction") plt.plot(m_x.to_global_data(), label="reconstruction")
plt.legend() plt.legend()
plt.show() plt.savefig('Krylov_reconstruction.png')
plt.close()
pltdict = {'alpha': .3, 'linewidth': .2}
for i in range(N_samps): for i in range(N_samps):
plt.plot(sky(samps_old[i]).val, color = 'b') if i == 0:
plt.plot(sky(samps[i]).val, color = 'r') plt.plot(sky(samps_old[i]).to_global_data(), color='b',
plt.plot((s_x-m_x).val,color='k') label='Traditional samples (residuals)',
**pltdict)
plt.plot(sky(samps[i]).to_global_data(), color='r',
label='Krylov samples (residuals)',
**pltdict)
else:
plt.plot(sky(samps_old[i]).to_global_data(), color='b', **pltdict)
plt.plot(sky(samps[i]).to_global_data(), color='r', **pltdict)
plt.plot((s_x - m_x).to_global_data(), color='k', label='signal - mean')
plt.legend() plt.legend()
plt.show() plt.savefig('Krylov_samples_residuals.png')
D_hat_old = ift.Field.zeros(x_space).val plt.close()
D_hat_new = ift.Field.zeros(x_space).val
for i in range(N_samps):
D_hat_old += sky(samps_old[i]).val**2
D_hat_new += sky(samps[i]).val**2
plt.plot(np.sqrt(D_hat_old/N_samps), color='b')
plt.plot(np.sqrt(D_hat_new/N_samps), color='r')
plt.plot(-np.sqrt(D_hat_old/N_samps), color='b')
plt.plot(-np.sqrt(D_hat_new/N_samps), color='r')
plt.plot((s_x-m_x).val, color='k')
plt.show()
D_hat_old = ift.Field.zeros(x_space).to_global_data()
D_hat_new = ift.Field.zeros(x_space).to_global_data()
for i in range(N_samps):
D_hat_old += sky(samps_old[i]).to_global_data()**2
D_hat_new += sky(samps[i]).to_global_data()**2
plt.plot(np.sqrt(D_hat_old / N_samps), 'r--', label='Traditional uncertainty')
plt.plot(-np.sqrt(D_hat_old / N_samps), 'r--')
plt.fill_between(range(len(D_hat_new)), -np.sqrt(D_hat_new / N_samps), np.sqrt(
D_hat_new / N_samps), facecolor='0.5', alpha=0.5,
label='Krylov uncertainty')
plt.plot((s_x - m_x).to_global_data(), color='k', label='signal - mean')
plt.legend()
plt.savefig('Krylov_uncertainty.png')
plt.close()
...@@ -393,7 +393,9 @@ def from_local_data(shape, arr, distaxis=0): ...@@ -393,7 +393,9 @@ def from_local_data(shape, arr, distaxis=0):
return data_object(shape, arr, distaxis) return data_object(shape, arr, distaxis)
def from_global_data(arr, distaxis=0): def from_global_data(arr, sum_up=False, distaxis=0):
if sum_up:
arr = np_allreduce_sum(arr)
if distaxis == -1: if distaxis == -1:
return data_object(arr.shape, arr, distaxis) return data_object(arr.shape, arr, distaxis)
lo, hi = _shareRange(arr.shape[distaxis], ntask, rank) lo, hi = _shareRange(arr.shape[distaxis], ntask, rank)
...@@ -427,7 +429,7 @@ def redistribute(arr, dist=None, nodist=None): ...@@ -427,7 +429,7 @@ def redistribute(arr, dist=None, nodist=None):
break break
if arr._distaxis == -1: # all data available, just pick the proper subset if arr._distaxis == -1: # all data available, just pick the proper subset
return from_global_data(arr._data, dist) return from_global_data(arr._data, distaxis=dist)
if dist == -1: # gather all data on all tasks if dist == -1: # gather all data on all tasks
tmp = np.moveaxis(arr._data, arr._distaxis, 0) tmp = np.moveaxis(arr._data, arr._distaxis, 0)
slabsize = np.prod(tmp.shape[1:])*tmp.itemsize slabsize = np.prod(tmp.shape[1:])*tmp.itemsize
......
...@@ -80,7 +80,7 @@ def from_local_data(shape, arr, distaxis=-1): ...@@ -80,7 +80,7 @@ def from_local_data(shape, arr, distaxis=-1):
return arr return arr
def from_global_data(arr, distaxis=-1): def from_global_data(arr, sum_up=False, distaxis=-1):
return arr return arr
......
...@@ -155,7 +155,7 @@ class Field(object): ...@@ -155,7 +155,7 @@ class Field(object):
return Field.empty(field._domain, dtype) return Field.empty(field._domain, dtype)
@staticmethod @staticmethod
def from_global_data(domain, arr): def from_global_data(domain, arr, sum_up=False):
"""Returns a Field constructed from `domain` and `arr`. """Returns a Field constructed from `domain` and `arr`.
Parameters Parameters
...@@ -165,10 +165,13 @@ class Field(object): ...@@ -165,10 +165,13 @@ class Field(object):
arr : numpy.ndarray arr : numpy.ndarray
The data content to be used for the new Field. The data content to be used for the new Field.
Its shape must match the shape of `domain`. Its shape must match the shape of `domain`.
If MPI is active, the contents of `arr` must be the same on all sum_up : bool, optional
MPI tasks. If True, the contents of `arr` are summed up over all MPI tasks
(if any), and the sum is used as data content.
If False, the contens of `arr` are used directly, and must be
identical on all MPI tasks.
""" """
return Field(domain, dobj.from_global_data(arr)) return Field(domain, dobj.from_global_data(arr, sum_up))
@staticmethod @staticmethod
def from_local_data(domain, arr): def from_local_data(domain, arr):
......
...@@ -5,3 +5,4 @@ from .nonlinear_power_energy import NonlinearPowerEnergy ...@@ -5,3 +5,4 @@ from .nonlinear_power_energy import NonlinearPowerEnergy
from .nonlinear_wiener_filter_energy import NonlinearWienerFilterEnergy from .nonlinear_wiener_filter_energy import NonlinearWienerFilterEnergy
from .poisson_energy import PoissonEnergy from .poisson_energy import PoissonEnergy
from .nonlinearities import Exponential, Linear, Tanh, PositiveTanh from .nonlinearities import Exponential, Linear, Tanh, PositiveTanh
from .krylov_sampling import generate_krylov_samples
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
import numpy as np
from ..minimization.quadratic_energy import QuadraticEnergy
def generate_krylov_samples(D_inv, S, j, N_samps, controller):
"""
Generates inverse samples from a curvature D.
This algorithm iteratively generates samples from
a curvature D by applying conjugate gradient steps
and resampling the curvature in search direction.
Parameters
----------
D_inv : EndomorphicOperator
The curvature which will be the inverse of the covarianc
of the generated samples
S : EndomorphicOperator (from which one can sample)
A prior covariance operator which is used to generate prior
samples that are then iteratively updated
j : Field, optional
A Field to which the inverse of D_inv is applied. The solution
of this matrix inversion problem is a side product of generating
the samples.
If not supplied, it is sampled from the inverse prior.
N_samps : Int
How many samples to generate.
controller : IterationController
convergence controller for the conjugate gradient iteration
Returns
-------
(solution, samples) : A tuple of a field 'solution' and a list of fields
'samples'. The first entry of the tuple is the solution x to
D_inv(x) = j
and the second entry are a list of samples from D_inv.inverse
"""
# MR FIXME: this should be synchronized with the "official" Nifty CG
j = S.draw_sample(from_inverse=True) if j is None else j
x = j*0.
energy = QuadraticEnergy(x, D_inv, j)
y = [S.draw_sample() for _ in range(N_samps)]
status = controller.start(energy)
if status != controller.CONTINUE:
return x, y
r = j.copy()
p = r.copy()
d = p.vdot(D_inv(p))
while True:
gamma = r.vdot(r)/d
if gamma == 0.:
break
x = x + gamma*p
Dip = D_inv(p)
for samp in y:
samp += (randn() * sqrt(d) - samp.vdot(Dip)) / d * p
energy = energy.at(x)
status = controller.check(energy)
if status != controller.CONTINUE:
return x, y
r_new = r - gamma * Dip
beta = r_new.vdot(r_new) / r.vdot(r)
r = r_new
p = r + beta * p
d = p.vdot(Dip)
if d == 0.:
break
return x, y
...@@ -63,6 +63,7 @@ class FFTOperator(LinearOperator): ...@@ -63,6 +63,7 @@ class FFTOperator(LinearOperator):
import pyfftw import pyfftw
pyfftw.interfaces.cache.enable() pyfftw.interfaces.cache.enable()
pyfftw.interfaces.cache.set_keepalive_time(1000.)
def apply(self, x, mode): def apply(self, x, mode):
self._check_input(x, mode) self._check_input(x, mode)
......
...@@ -74,7 +74,7 @@ class InversionEnabler(EndomorphicOperator): ...@@ -74,7 +74,7 @@ class InversionEnabler(EndomorphicOperator):
prec = self._approximation prec = self._approximation
if prec is not None: if prec is not None:
prec = prec._flip_modes(self._ilog[mode]) prec = prec._flip_modes(self._ilog[mode])
energy = QuadraticEnergy(A=invop, b=x, position=x0) energy = QuadraticEnergy(x0, invop, x)
r, stat = self._inverter(energy, preconditioner=prec) r, stat = self._inverter(energy, preconditioner=prec)
if stat != IterationController.CONVERGED: if stat != IterationController.CONVERGED:
logger.warning("Error detected during operator inversion") logger.warning("Error detected during operator inversion")
......
...@@ -48,6 +48,17 @@ class LinearOperator(NiftyMetaBase()): ...@@ -48,6 +48,17 @@ class LinearOperator(NiftyMetaBase()):
by means of a single integer number. by means of a single integer number.
""" """
# Field Operator Modes
TIMES = 1
ADJOINT_TIMES = 2
INVERSE_TIMES = 4
ADJOINT_INVERSE_TIMES = 8
INVERSE_ADJOINT_TIMES = 8
# Operator Transform Flags
ADJOINT_BIT = 1
INVERSE_BIT = 2
_ilog = (-1, 0, 1, -1, 2, -1, -1, -1, 3) _ilog = (-1, 0, 1, -1, 2, -1, -1, -1, 3)
_validMode = (False, True, True, False, True, False, False, False, True) _validMode = (False, True, True, False, True, False, False, False, True)
_modeTable = ((1, 2, 4, 8), _modeTable = ((1, 2, 4, 8),
...@@ -61,13 +72,6 @@ class LinearOperator(NiftyMetaBase()): ...@@ -61,13 +72,6 @@ class LinearOperator(NiftyMetaBase()):
_addInverse = (0, 5, 10, 15, 5, 5, 15, 15, 10, 15, 10, 15, 15, 15, 15, 15) _addInverse = (0, 5, 10, 15, 5, 5, 15, 15, 10, 15, 10, 15, 15, 15, 15, 15)
_backwards = 6 _backwards = 6
_all_ops = 15 _all_ops = 15
TIMES = 1
ADJOINT_TIMES = 2
INVERSE_TIMES = 4
ADJOINT_INVERSE_TIMES = 8
INVERSE_ADJOINT_TIMES = 8
ADJOINT_BIT = 1
INVERSE_BIT = 2
def _dom(self, mode): def _dom(self, mode):
return self.domain if (mode & 9) else self.target return self.domain if (mode & 9) else self.target
...@@ -92,9 +96,10 @@ class LinearOperator(NiftyMetaBase()): ...@@ -92,9 +96,10 @@ class LinearOperator(NiftyMetaBase()):
The domain on which the Operator's output Field lives.""" The domain on which the Operator's output Field lives."""
raise NotImplementedError raise NotImplementedError
def _flip_modes(self, mode): def _flip_modes(self, op_transform):
from .operator_adapter import OperatorAdapter from .operator_adapter import OperatorAdapter
return self if mode == 0 else OperatorAdapter(self, mode) return self if op_transform == 0 \
else OperatorAdapter(self, op_transform)
@property @property
def inverse(self): def inverse(self):
......
...@@ -23,33 +23,34 @@ import numpy as np ...@@ -23,33 +23,34 @@ import numpy as np
class OperatorAdapter(LinearOperator): class OperatorAdapter(LinearOperator):
"""Class representing the inverse and/or adjoint of another operator.""" """Class representing the inverse and/or adjoint of another operator."""
def __init__(self, op, mode): def __init__(self, op, op_transform):
super(OperatorAdapter, self).__init__() super(OperatorAdapter, self).__init__()
self._op = op self._op = op
self._mode = int(mode) self._trafo = int(op_transform)
if self._mode < 1 or self._mode > 3: if self._trafo < 1 or self._trafo > 3:
raise ValueError("invalid mode") raise ValueError("invalid mode")
@property @property
def domain(self): def domain(self):
return self._op._dom(1 << self._mode) return self._op._dom(1 << self._trafo)
@property @property
def target(self): def target(self):
return self._op._tgt(1 << self._mode) return self._op._tgt(1 << self._trafo)
@property @property
def capability(self): def capability(self):
return self._capTable[self._mode][self._op.capability] return self._capTable[self._trafo][self._op.capability]
def _flip_modes(self, mode): def _flip_modes(self, op_transform):
newmode = mode ^ self._mode newmode = op_transform ^ self._trafo
return self._op if newmode == 0 else OperatorAdapter(self._op, newmode) return self._op if newmode == 0 else OperatorAdapter(self._op, newmode)
def apply(self, x, mode): def apply(self, x, mode):
return self._op.apply(x, self._modeTable[self._mode][self._ilog[mode]]) return self._op.apply(x,
self._modeTable[self._trafo][self._ilog[mode]])
def draw_sample(self, from_inverse=False, dtype=np.float64): def draw_sample(self, from_inverse=False, dtype=np.float64):
if self._mode & self.INVERSE_BIT: if self._trafo & self.INVERSE_BIT:
return self._op.draw_sample(not from_inverse, dtype) return self._op.draw_sample(not from_inverse, dtype)
return self._op.draw_sample(from_inverse, dtype) return self._op.draw_sample(from_inverse, dtype)
...@@ -160,6 +160,14 @@ def NiftyMetaBase(): ...@@ -160,6 +160,14 @@ def NiftyMetaBase():
return with_metaclass(NiftyMeta, type('NewBase', (object,), {})) return with_metaclass(NiftyMeta, type('NewBase', (object,), {}))
def nthreads():
if nthreads._val is None:
import os
nthreads._val = int(os.getenv("OMP_NUM_THREADS", "1"))
return nthreads._val
nthreads._val = None
def hartley(a, axes=None): def hartley(a, axes=None):
# Check if the axes provided are valid given the shape # Check if the axes provided are valid given the shape
if axes is not None and \ if axes is not None and \
...@@ -169,7 +177,7 @@ def hartley(a, axes=None): ...@@ -169,7 +177,7 @@ def hartley(a, axes=None):
raise TypeError("Hartley transform requires real-valued arrays.") raise TypeError("Hartley transform requires real-valued arrays.")
from pyfftw.interfaces.numpy_fft import rfftn from pyfftw.interfaces.numpy_fft import rfftn
tmp = rfftn(a, axes=axes) tmp = rfftn(a, axes=axes, threads=nthreads())
def _fill_array(tmp, res, axes): def _fill_array(tmp, res, axes):
if axes is None: if axes is None:
...@@ -211,7 +219,7 @@ def my_fftn_r2c(a, axes=None): ...@@ -211,7 +219,7 @@ def my_fftn_r2c(a, axes=None):
raise TypeError("Transform requires real-valued input arrays.") raise TypeError("Transform requires real-valued input arrays.")
from pyfftw.interfaces.numpy_fft import rfftn from pyfftw.interfaces.numpy_fft import rfftn
tmp = rfftn(a, axes=axes) tmp = rfftn(a, axes=axes, threads=nthreads())
def _fill_complex_array(tmp, res, axes): def _fill_complex_array(tmp, res, axes):
if axes is None: if axes is None:
......
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