Commit 1025b2e7 authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'NIFTy_4' into yango_minimizer

parents 9eb848ac 484ddb59
Pipeline #28498 passed with stages
in 13 minutes and 6 seconds
image: parras/nifty:latest
image: $CONTAINER_TEST_IMAGE
variables:
CONTAINER_TEST_IMAGE: gitlab-registry.mpcdf.mpg.de/ift/nifty:$CI_BUILD_REF_NAME
stages:
- build_docker
- test
- release
variables:
DOCKER_DRIVER: overlay
build_docker_from_scratch:
only:
- schedules
image: docker:stable
stage: build_docker
script:
- docker login -u gitlab-ci-token -p $CI_BUILD_TOKEN gitlab-registry.mpcdf.mpg.de
- docker build -t $CONTAINER_TEST_IMAGE --no-cache .
- docker push $CONTAINER_TEST_IMAGE
build_docker_from_cache:
except:
- schedules
image: docker:stable
stage: build_docker
script:
- docker login -u gitlab-ci-token -p $CI_BUILD_TOKEN gitlab-registry.mpcdf.mpg.de
- docker build -t $CONTAINER_TEST_IMAGE .
- docker push $CONTAINER_TEST_IMAGE
test_python2_scalar:
stage: test
......
......@@ -2,14 +2,28 @@ FROM debian:testing-slim
RUN apt-get update
# Debian package installations
RUN apt-get install -y git libfftw3-dev openmpi-bin libopenmpi-dev python python-pip python-dev python-nose python-numpy python-matplotlib python-future python-mpi4py python-scipy
RUN apt-get install -y python3 python3-pip python3-dev python3-nose python3-numpy python3-matplotlib python3-future python3-mpi4py python3-scipy
# Needed for gitlab tests
RUN apt-get install -y git
# Packages needed for NIFTy
RUN apt-get install -y libfftw3-dev
RUN apt-get install -y python python-pip python-dev python-future python-scipy
RUN apt-get install -y python3 python3-pip python3-dev python3-future python3-scipy
RUN pip install pyfftw
RUN pip3 install pyfftw
# Optional NIFTy dependencies
RUN apt-get install -y openmpi-bin libopenmpi-dev python-mpi4py python3-mpi4py
RUN pip install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
RUN pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
# Documentation build dependencies
RUN apt-get install -y python-sphinx python-sphinx-rtd-theme python-numpydoc
# Python module installations
RUN pip install parameterized coverage pyfftw git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
RUN pip3 install parameterized coverage pyfftw git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
# Testing dependencies
RUN apt-get install -y python-nose python-parameterized
RUN apt-get install -y python3-nose python3-parameterized
RUN pip install coverage
# Create user (openmpi does not like to be run as root)
RUN useradd -ms /bin/bash testinguser
......
......@@ -38,7 +38,7 @@ Installation
### Requirements
- [Python](https://www.python.org/) (v2.7.x or 3.5.x)
- [NumPy](https://www.numpy.org/)
- [SciPy](https://www.scipy.org/)
- [pyFFTW](https://pypi.python.org/pypi/pyFFTW)
Optional dependencies:
......@@ -46,7 +46,6 @@ Optional dependencies:
transforms involving domains on the sphere)
- [mpi4py](https://mpi4py.scipy.org) (for MPI-parallel execution)
- [matplotlib](https://matplotlib.org/) (for field plotting)
- [SciPy](https://www.scipy.org/) (for additional minimization algorithms)
### Sources
......@@ -83,10 +82,6 @@ MPI support is added via:
sudo apt-get install openmpi-bin libopenmpi-dev
pip install --user mpi4py
Scipy-based minimizers are enabled via:
pip install --user scipy
### Installation for Python 3
If you want to run NIFTy with Python 3, you need to make the following changes
......
import nifty4 as ift
import numpy as np
if __name__ == "__main__":
np.random.seed(13)
s_space = ift.RGSpace([512, 512])
h_space = s_space.get_default_codomain()
ht = ift.HarmonicTransformOperator(h_space)
pow_spec = (lambda k: 42 / (k+1)**4)
S = ift.create_power_operator(h_space, power_spectrum=pow_spec)
sh = S.draw_sample()
ss = ht(sh)
def make_random_los(n):
starts = np.random.random((2, n))
ends = np.random.random((2, n))
return starts, ends
nlos = 1000
starts, ends = make_random_los(nlos)
R = ift.library.LOSResponse(
s_space, starts=starts, ends=ends,
sigmas_up=np.full(nlos, 0.1), sigmas_low=np.full(nlos, 0.1))
Rh = R*ht
signal_to_noise = 10
noise_amplitude = Rh(sh).val.std()/signal_to_noise
N = ift.ScalingOperator(noise_amplitude**2, Rh.target)
n = N.draw_sample()
d = Rh(sh) + n
j = Rh.adjoint_times(N.inverse_times(d))
ctrl = ift.GradientNormController(name="Iter", tol_abs_gradnorm=1e-10,
iteration_limit=300)
inverter = ift.ConjugateGradient(controller=ctrl)
Di = ift.library.WienerFilterCurvature(S=S, R=Rh, N=N, inverter=inverter)
mh = Di.inverse_times(j)
m = ht(mh)
ift.plot(m, name="reconstruction.png")
ift.plot(ss, name="signal.png")
ift.plot(ht(j), name="j.png")
......@@ -28,10 +28,6 @@ MPI support is added via::
sudo apt-get install openmpi-bin libopenmpi-dev
pip install --user mpi4py
Scipy-based minimizers are enabled via::
pip install --user scipy
Installation for Python 3
-------------------------
......@@ -40,4 +36,3 @@ to the instructions above:
- in all `apt-get` commands, replace `python-*` by `python3-*`
- in all `pip` commands, replace `pip` by `pip3`
......@@ -24,6 +24,9 @@ from .utilities import memo
from .logger import logger
from .multi import *
__all__ = ["__version__", "dobj", "DomainTuple"] + \
domains.__all__ + operators.__all__ + minimization.__all__ + \
["DomainTuple", "Field", "sqrt", "exp", "log"]
["DomainTuple", "Field", "sqrt", "exp", "log"] + \
multi.__all__
......@@ -19,6 +19,7 @@
import numpy as np
from .random import Random
from mpi4py import MPI
import sys
_comm = MPI.COMM_WORLD
ntask = _comm.Get_size()
......@@ -185,80 +186,11 @@ class data_object(object):
else:
return data_object(self._shape, tval, self._distaxis)
def __add__(self, other):
return self._binary_helper(other, op='__add__')
def __radd__(self, other):
return self._binary_helper(other, op='__radd__')
def __iadd__(self, other):
return self._binary_helper(other, op='__iadd__')
def __sub__(self, other):
return self._binary_helper(other, op='__sub__')
def __rsub__(self, other):
return self._binary_helper(other, op='__rsub__')
def __isub__(self, other):
return self._binary_helper(other, op='__isub__')
def __mul__(self, other):
return self._binary_helper(other, op='__mul__')
def __rmul__(self, other):
return self._binary_helper(other, op='__rmul__')
def __imul__(self, other):
return self._binary_helper(other, op='__imul__')
def __div__(self, other):
return self._binary_helper(other, op='__div__')
def __rdiv__(self, other):
return self._binary_helper(other, op='__rdiv__')
def __idiv__(self, other):
return self._binary_helper(other, op='__idiv__')
def __truediv__(self, other):
return self._binary_helper(other, op='__truediv__')
def __rtruediv__(self, other):
return self._binary_helper(other, op='__rtruediv__')
def __pow__(self, other):
return self._binary_helper(other, op='__pow__')
def __rpow__(self, other):
return self._binary_helper(other, op='__rpow__')
def __ipow__(self, other):
return self._binary_helper(other, op='__ipow__')
def __lt__(self, other):
return self._binary_helper(other, op='__lt__')
def __le__(self, other):
return self._binary_helper(other, op='__le__')
def __ne__(self, other):
return self._binary_helper(other, op='__ne__')
def __eq__(self, other):
return self._binary_helper(other, op='__eq__')
def __ge__(self, other):
return self._binary_helper(other, op='__ge__')
def __gt__(self, other):
return self._binary_helper(other, op='__gt__')
def __neg__(self):
return data_object(self._shape, -self._data, self._distaxis)
def __abs__(self):
return data_object(self._shape, np.abs(self._data), self._distaxis)
return data_object(self._shape, abs(self._data), self._distaxis)
def all(self):
return self.sum() == self.size
......@@ -269,6 +201,20 @@ class data_object(object):
def fill(self, value):
self._data.fill(value)
for op in ["__add__", "__radd__", "__iadd__",
"__sub__", "__rsub__", "__isub__",
"__mul__", "__rmul__", "__imul__",
"__div__", "__rdiv__", "__idiv__",
"__truediv__", "__rtruediv__", "__itruediv__",
"__floordiv__", "__rfloordiv__", "__ifloordiv__",
"__pow__", "__rpow__", "__ipow__",
"__lt__", "__le__", "__gt__", "__ge__", "__eq__", "__ne__"]:
def func(op):
def func2(self, other):
return self._binary_helper(other, op=op)
return func2
setattr(data_object, op, func(op))
def full(shape, fill_value, dtype=None, distaxis=0):
return data_object(shape, np.full(local_shape(shape, distaxis),
......@@ -302,6 +248,7 @@ def vdot(a, b):
def _math_helper(x, function, out):
function = getattr(np, function)
if out is not None:
function(x._data, out=out._data)
return out
......@@ -309,24 +256,14 @@ def _math_helper(x, function, out):
return data_object(x.shape, function(x._data), x._distaxis)
def abs(a, out=None):
return _math_helper(a, np.abs, out)
def exp(a, out=None):
return _math_helper(a, np.exp, out)
def log(a, out=None):
return _math_helper(a, np.log, out)
def tanh(a, out=None):
return _math_helper(a, np.tanh, out)
_current_module = sys.modules[__name__]
def sqrt(a, out=None):
return _math_helper(a, np.sqrt, out)
for f in ["sqrt", "exp", "log", "tanh", "conjugate"]:
def func(f):
def func2(x, out=None):
return _math_helper(x, f, out)
return func2
setattr(_current_module, f, func(f))
def from_object(object, dtype, copy, set_locked):
......
......@@ -20,7 +20,7 @@
import numpy as np
from numpy import ndarray as data_object
from numpy import full, empty, empty_like, sqrt, ones, zeros, vdot, abs, \
from numpy import full, empty, empty_like, sqrt, ones, zeros, vdot, \
exp, log, tanh
from .random import Random
......
......@@ -66,6 +66,8 @@ class DomainTuple(object):
"""
if isinstance(domain, DomainTuple):
return domain
if isinstance(domain, dict):
return domain
domain = DomainTuple._parse_domain(domain)
obj = DomainTuple._tupleCache.get(domain)
if obj is not None:
......
......@@ -182,10 +182,9 @@ class RGSpace(StructuredDomain):
"not be the same.")
# Check if the distances match, i.e. dist' = 1 / (num * dist)
if not np.all(
np.absolute(np.array(self.shape) *
np.array(self.distances) *
np.array(codomain.distances)-1) < 1e-7):
if not np.all(abs(np.array(self.shape) *
np.array(self.distances) *
np.array(codomain.distances)-1) < 1e-7):
raise AttributeError("The grid-distances of domain and codomain "
"do not match.")
......
......@@ -23,6 +23,7 @@ from . import utilities
from .domain_tuple import DomainTuple
from functools import reduce
from . import dobj
import sys
__all__ = ["Field", "sqrt", "exp", "log", "conjugate"]
......@@ -518,7 +519,7 @@ class Field(object):
float
The L2-norm of the field values.
"""
return np.sqrt(np.abs(self.vdot(x=self)))
return np.sqrt(abs(self.vdot(x=self)))
def conjugate(self):
""" Returns the complex conjugate of the field.
......@@ -539,7 +540,7 @@ class Field(object):
return Field(self._domain, -self.val)
def __abs__(self):
return Field(self._domain, dobj.abs(self.val))
return Field(self._domain, abs(self.val))
def _contraction_helper(self, op, spaces):
if spaces is None:
......@@ -745,89 +746,6 @@ class Field(object):
raise ValueError("domains are incompatible.")
self.local_data[()] = other.local_data[()]
def _binary_helper(self, other, op):
# if other is a field, make sure that the domains match
if isinstance(other, Field):
if other._domain != self._domain:
raise ValueError("domains are incompatible.")
tval = getattr(self.val, op)(other.val)
return self if tval is self.val else Field(self._domain, tval)
if np.isscalar(other) or isinstance(other, dobj.data_object):
tval = getattr(self.val, op)(other)
return self if tval is self.val else Field(self._domain, tval)
return NotImplemented
def __add__(self, other):
return self._binary_helper(other, op='__add__')
def __radd__(self, other):
return self._binary_helper(other, op='__radd__')
def __iadd__(self, other):
return self._binary_helper(other, op='__iadd__')
def __sub__(self, other):
return self._binary_helper(other, op='__sub__')
def __rsub__(self, other):
return self._binary_helper(other, op='__rsub__')
def __isub__(self, other):
return self._binary_helper(other, op='__isub__')
def __mul__(self, other):
return self._binary_helper(other, op='__mul__')
def __rmul__(self, other):
return self._binary_helper(other, op='__rmul__')
def __imul__(self, other):
return self._binary_helper(other, op='__imul__')
def __div__(self, other):
return self._binary_helper(other, op='__div__')
def __truediv__(self, other):
return self._binary_helper(other, op='__truediv__')
def __rdiv__(self, other):
return self._binary_helper(other, op='__rdiv__')
def __rtruediv__(self, other):
return self._binary_helper(other, op='__rtruediv__')
def __idiv__(self, other):
return self._binary_helper(other, op='__idiv__')
def __pow__(self, other):
return self._binary_helper(other, op='__pow__')
def __rpow__(self, other):
return self._binary_helper(other, op='__rpow__')
def __ipow__(self, other):
return self._binary_helper(other, op='__ipow__')
def __lt__(self, other):
return self._binary_helper(other, op='__lt__')
def __le__(self, other):
return self._binary_helper(other, op='__le__')
def __ne__(self, other):
return self._binary_helper(other, op='__ne__')
def __eq__(self, other):
return self._binary_helper(other, op='__eq__')
def __ge__(self, other):
return self._binary_helper(other, op='__ge__')
def __gt__(self, other):
return self._binary_helper(other, op='__gt__')
def __repr__(self):
return "<nifty4.Field>"
......@@ -836,36 +754,48 @@ class Field(object):
self._domain.__str__() + \
"\n- val = " + repr(self.val)
# Arithmetic functions working on Fields
def _math_helper(x, function, out):
if not isinstance(x, Field):
raise TypeError("This function only accepts Field objects.")
if out is not None:
if not isinstance(out, Field) or x._domain != out._domain:
raise ValueError("Bad 'out' argument")
function(x.val, out=out.val)
return out
else:
return Field(domain=x._domain, val=function(x.val))
for op in ["__add__", "__radd__", "__iadd__",
"__sub__", "__rsub__", "__isub__",
"__mul__", "__rmul__", "__imul__",
"__div__", "__rdiv__", "__idiv__",
"__truediv__", "__rtruediv__", "__itruediv__",
"__floordiv__", "__rfloordiv__", "__ifloordiv__",
"__pow__", "__rpow__", "__ipow__",
"__lt__", "__le__", "__gt__", "__ge__", "__eq__", "__ne__"]:
def func(op):
def func2(self, other):
# if other is a field, make sure that the domains match
if isinstance(other, Field):
if other._domain != self._domain:
raise ValueError("domains are incompatible.")
tval = getattr(self.val, op)(other.val)
return self if tval is self.val else Field(self._domain, tval)
if np.isscalar(other) or isinstance(other, dobj.data_object):
tval = getattr(self.val, op)(other)
return self if tval is self.val else Field(self._domain, tval)
return NotImplemented
return func2
setattr(Field, op, func(op))
def sqrt(x, out=None):
return _math_helper(x, dobj.sqrt, out)
def exp(x, out=None):
return _math_helper(x, dobj.exp, out)
def log(x, out=None):
return _math_helper(x, dobj.log, out)
def tanh(x, out=None):
return _math_helper(x, dobj.tanh, out)
# Arithmetic functions working on Fields
def conjugate(x, out=None):
return _math_helper(x, dobj.conjugate, out)
_current_module = sys.modules[__name__]
for f in ["sqrt", "exp", "log", "tanh", "conjugate"]:
def func(f):
def func2(x, out=None):
fu = getattr(dobj, f)
if not isinstance(x, Field):
raise TypeError("This function only accepts Field objects.")
if out is not None:
if not isinstance(out, Field) or x._domain != out._domain:
raise ValueError("Bad 'out' argument")
fu(x.val, out=out.val)
return out
else:
return Field(domain=x._domain, val=fu(x.val))
return func2
setattr(_current_module, f, func(f))
......@@ -6,3 +6,4 @@ from .nonlinear_wiener_filter_energy import NonlinearWienerFilterEnergy
from .poisson_energy import PoissonEnergy
from .nonlinearities import Exponential, Linear, Tanh, PositiveTanh
from .krylov_sampling import generate_krylov_samples
from .los_response import LOSResponse
......@@ -80,8 +80,8 @@ def generate_krylov_samples(D_inv, S, j, N_samps, controller):
logger.error("Error: ConjugateGradient: alpha<0.")
return energy.position, y
for samp in y:
samp += (np.random.randn()*np.sqrt(ddotq) - samp.vdot(q))/ddotq * d
for i in range(len(y)):
y[i] += (np.random.randn()*np.sqrt(ddotq) - y[i].vdot(q))/ddotq * d
q *= -alpha
r = r + q
......
import numpy as np
from scipy.special import erfc
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import aslinearoperator
from ..operators.linear_operator import LinearOperator
from ..domain_tuple import DomainTuple
from ..domains.rg_space import RGSpace
from ..domains.unstructured_domain import UnstructuredDomain
from ..field import Field
from .. import dobj
def _gaussian_error_function(x):
return 0.5*erfc(x*np.sqrt(2.))
def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
ndim = start.shape[0]
nlos = start.shape[1]
inc = np.full(len(shp), 1)
for i in range(-2, -len(shp)-1, -1):
inc[i] = inc[i+1]*shp[i+1]
pmax = np.array(shp)
out = [None]*nlos
for i in range(nlos):
dir = end[:, i]-start[:, i]
dirx = np.where(dir == 0., 1e-12, dir)
d0 = np.where(dir == 0., ((start[:, i] > 0)-0.5)*1e12,
-start[:, i]/dirx)
d1 = np.where(dir == 0., ((start[:, i] < pmax)-0.5)*-1e12,
(pmax-start[:, i])/dirx)
(dmin, dmax) = (np.minimum(d0, d1), np.maximum(d0, d1))
dmin = dmin.max()
dmax = dmax.min()
dmin = np.maximum(0., dmin)
dmax = np.minimum(1., dmax)
dmax = np.maximum(dmin, dmax)
# hack: move away from potential grid crossings
dmin += 1e-7
dmax -= 1e-7
if dmin > dmax: # no intersection
out[i] = (np.full(0, 0), np.full(0, 0.))
continue
# determine coordinates of first cell crossing
c_first = np.ceil(start[:, i]+dir*dmin)
c_first = np.where(dir > 0., c_first, c_first-1.)
c_first = (c_first-start[:, i])/dirx
pos1 = np.asarray((start[:, i]+dmin*dir), dtype=np.int)
pos1 = np.sum(pos1*inc)
cdist = np.empty(0, dtype=np.float64)
add = np.empty(0, dtype=np.int)
for j in range(ndim):
if dir[j] != 0:
step = inc[j] if dir[j] > 0 else -inc[j]
tmp = np.arange(start=c_first[j], stop=dmax,
step=abs(1./dir[j]))
cdist = np.append(cdist, tmp)
add = np.append(add, np.full(len(tmp), step))
idx = np.argsort(cdist)
cdist = cdist[idx]
add = add[idx]
cdist = np.append(np.full(1, dmin), cdist)
cdist = np.append(cdist, np.full(1, dmax))
corfac = np.linalg.norm(dir*dist)
cdist *= corfac
wgt = np.diff(cdist)
mdist = 0.5*(cdist[:-1]+cdist[1:])
wgt = apply_erf(wgt, mdist, lo[i], mid[i], hi[i], erf)
add = np.append(pos1, add)
add = np.cumsum(add)
out[i] = (add, wgt)
return out
def apply_erf(wgt, dist, lo, mid, hi, erf):
wgt = wgt.copy()
mask = dist > hi
wgt[mask] = 0.
mask = (dist > mid) & (dist <= hi)
wgt[mask] *= erf((dist[mask]-mid)/(hi-mid))
mask = (dist <= mid) & (dist > lo)
wgt[mask] *= erf((dist[mask]-mid)/(mid-lo))
return wgt
class LOSResponse(LinearOperator):
"""Line-of-sight response operator
This operator transforms from a single RGSpace to an unstructured domain
with as many entries as there were lines of sight passed to the
constructor. Adjoint application is also provided.
Parameters
----------
domain : RGSpace or DomainTuple
The operator's input domain. This must be a single RGSpace.
starts, ends : numpy.ndarray(float) with two dimensions
Arrays containing the start and end points of the individual lines
of sight. The first dimension must have as many entries as `domain`
has dimensions. The second dimensions must be identical for both arrays
and indicated the total number of lines of sight.
sigmas_low, sigmas_up : numpy.ndarray(float) (optional)
For expert use. If unsure, leave blank.
Notes
-----
`starts, `ends`, `sigmas_low`, and `sigmas_up` have to be identical on
every calling MPI task (i.e. the full LOS information has to be provided on
every task).
"""
def __init__(self, domain, starts, ends, sigmas_low=None, sigmas_up=None):
super(LOSResponse, self).__init__()
self._domain = DomainTuple.make(domain)
if ((not isinstance(self.domain[0], RGSpace)) or
(len(self._domain) != 1)):
raise TypeError("The domain must be exactly one RGSpace instance.")
ndim = len(self.domain[0].shape)
starts = np.array(starts)
nlos = starts.shape[1]
ends = np.array(ends)
if sigmas_low is None:
sigmas_low = np.zeros(nlos, dtype=np.float32)
if sigmas_up is None:
sigmas_up = np.zeros(nlos, dtype=np.float32)
sigmas_low = np.array(sigmas_low)
sigmas_up = np.array(sigmas_up)
assert starts.shape[0] == ndim, "dimension mismatch"
assert nlos == sigmas_low.shape[0], "dimension mismatch"
assert starts.shape == ends.shape, "dimension mismatch"
assert sigmas_low.shape == sigmas_up.shape, "dimension mismatch"
self._local_shape = dobj.local_shape(self