merge NIFTy_5

.gitlab-ci.yml

@@ -34,18 +34,22 @@ build_docker_from_cache:
     - docker build -t $CONTAINER_TEST_IMAGE .
     - docker push $CONTAINER_TEST_IMAGE
 
-test:
+test_serial:
   stage: test
-  variables:
-    OMPI_MCA_btl_vader_single_copy_mechanism: none
   script:
-    - mpiexec -n 2 --bind-to none pytest-3 -q test
     - pytest-3 -q --cov=nifty5 test
     - >
       python3 -m coverage report --omit "*plotting*,*distributed_do*"
     - >
       python3 -m coverage report --omit "*plotting*,*distributed_do*" | grep TOTAL | awk '{ print "TOTAL: "$4; }'
 
+test_mpi:
+  stage: test
+  variables:
+    OMPI_MCA_btl_vader_single_copy_mechanism: none
+  script:
+    - mpiexec -n 2 --bind-to none pytest-3 -q test
+
 pages:
   stage: release
   before_script:

Dockerfile

@@ -9,7 +9,7 @@ RUN apt-get update && apt-get install -y \
     # Documentation build dependencies
     python3-sphinx python3-sphinx-rtd-theme python3-numpydoc \
     # Testing dependencies
-    python3-coverage python3-parameterized python3-pytest python3-pytest-cov \
+    python3-coverage python3-pytest python3-pytest-cov \
     # Optional NIFTy dependencies
     openmpi-bin libopenmpi-dev python3-mpi4py \
     # Packages needed for NIFTy

nifty5/library/los_response.py

@@ -28,14 +28,14 @@ from ..field import Field
 from ..operators.linear_operator import LinearOperator
 
 
-def _gaussian_error_function(x):
-    return 0.5/erfc(x*np.sqrt(2.))
+def _gaussian_sf(x):
+    return 0.5*erfc(x/np.sqrt(2.))
 
 
-def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
+def _comp_traverse(start, end, shp, dist, lo, mid, hi, sig, erf):
     ndim = start.shape[0]
     nlos = start.shape[1]
-    inc = np.full(len(shp), 1)
+    inc = np.full(len(shp), 1, dtype=np.int64)
     for i in range(-2, -len(shp)-1, -1):
         inc[i] = inc[i+1]*shp[i+1]
 
@@ -59,7 +59,7 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
         dmin += 1e-7
         dmax -= 1e-7
         if dmin >= dmax:  # no intersection
-            out[i] = (np.full(0, 0), np.full(0, 0.))
+            out[i] = (np.full(0, 0, dtype=np.int64), np.full(0, 0.))
             continue
         # determine coordinates of first cell crossing
         c_first = np.ceil(start[:, i]+dir*dmin)
@@ -75,7 +75,7 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
                 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))
+                add = np.append(add, np.full(len(tmp), step, dtype=np.int64))
         idx = np.argsort(cdist)
         cdist = cdist[idx]
         add = add[idx]
@@ -85,21 +85,19 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
         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)
+        wgt = apply_erf(wgt, mdist, lo[i], mid[i], hi[i], sig[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):
+def apply_erf(wgt, dist, lo, mid, hi, sig, 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))
+    mask = (dist > lo) & (dist <= hi)
+    wgt[mask] *= erf((-1/dist[mask]+1/mid)/sig)
     return wgt
 
 
@@ -119,17 +117,32 @@ class LOSResponse(LinearOperator):
         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.
+    sigmas: numpy.ndarray(float) (optional)
+        If this is not None, the inverse of the lengths of the LOSs are assumed
+        to be Gaussian distributed with these sigmas. The start point will
+        remain the same, but the endpoint is assumed to be unknown.
+        This is a typical statistical model for astrophysical parallaxes.
+        The LOS response then returns the expected integral
+        over the input given that the length of the LOS is unknown and
+        therefore the result is averaged over different endpoints.
+        default: None
+    truncation: float (optional)
+        Use only if the sigmas keyword argument is used!
+        This truncates the probability of the endpoint lying more sigmas away
+        than the truncation. Used to speed up computation and to avoid negative
+        distances. It should hold that `1./(1./length-sigma*truncation)>0`
+        for all lengths of the LOSs and all corresponding sigma of sigmas.
+        If unsure, leave blank.
+        default: 3.
 
     Notes
     -----
-    `starts, `ends`, `sigmas_low`, and `sigmas_up` have to be identical on
+    `starts, `ends`, `sigmas`, and `truncation` 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):
+    def __init__(self, domain, starts, ends, sigmas=None, truncation=3.):
         self._domain = DomainTuple.make(domain)
         self._capability = self.TIMES | self.ADJOINT_TIMES
 
@@ -141,20 +154,15 @@ class LOSResponse(LinearOperator):
         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)
+        if sigmas is None:
+            sigmas = np.zeros(nlos, dtype=np.float32)
+        sigmas = np.array(sigmas)
         if starts.shape[0] != ndim:
             raise TypeError("dimension mismatch")
-        if nlos != sigmas_low.shape[0]:
+        if nlos != sigmas.shape[0]:
             raise TypeError("dimension mismatch")
         if starts.shape != ends.shape:
             raise TypeError("dimension mismatch")
-        if sigmas_low.shape != sigmas_up.shape:
-            raise TypeError("dimension mismatch")
 
         self._local_shape = dobj.local_shape(self.domain[0].shape)
         local_zero_point = (np.array(
@@ -164,7 +172,11 @@ class LOSResponse(LinearOperator):
         diffs = ends-starts
         difflen = np.linalg.norm(diffs, axis=0)
         diffs /= difflen
-        real_ends = ends + sigmas_up*diffs
+        real_distances = 1./(1./difflen - truncation*sigmas)
+        if np.any(real_distances < 0):
+            raise ValueError("parallax error truncation to high: "
+                             "getting negative distances")
+        real_ends = starts + diffs*real_distances
         lzp = local_zero_point.reshape((-1, 1))
         dist = np.array(self.domain[0].distances).reshape((-1, 1))
         localized_pixel_starts = (starts-lzp)/dist + 0.5
@@ -175,8 +187,11 @@ class LOSResponse(LinearOperator):
                              localized_pixel_ends,
                              self._local_shape,
                              np.array(self.domain[0].distances),
-                             difflen-sigmas_low, difflen, difflen+sigmas_up,
-                             _gaussian_error_function)
+                             1./(1./difflen+truncation*sigmas),
+                             difflen,
+                             1./(1./difflen-truncation*sigmas),
+                             sigmas,
+                             _gaussian_sf)
 
         boxsz = 16
         nlos = len(w_i)
@@ -202,7 +217,7 @@ class LOSResponse(LinearOperator):
                 tmp += (fullidx[j]//boxsz)*fct
                 fct *= self._local_shape[j]
             tmp += cnt/float(nlos)
-            tmp += iarr[ofs:ofs+nval]/float(nlos*npix)
+            tmp += iarr[ofs:ofs+nval]/(float(nlos)*float(npix))
             pri[ofs:ofs+nval] = tmp
             ofs += nval
             cnt += 1

test/common.py

@@ -15,20 +15,12 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-from builtins import str
+import pytest
 
-import numpy as np
-from parameterized import parameterized
 
-np.seterr(all='raise', under='ignore')
+def list2fixture(lst):
+    @pytest.fixture(params=lst)
+    def myfixture(request):
+        return request.param
 
-
-def _custom_name_func(testcase_func, param_num, param):
-    return "{}_{}".format(
-        testcase_func.__name__,
-        parameterized.to_safe_name("_".join(str(x) for x in param.args)),
-    )
-
-
-def expand(*args, **kwargs):
-    return parameterized.expand(*args, func=_custom_name_func, **kwargs)
+    return myfixture

test/test_energies/test_consistency.py

-# 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-2019 Max-Planck-Society
-#
-# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-
-import unittest
-from itertools import product
-from test.common import expand
-
-import nifty5 as ift
-import numpy as np
-
-
-class Energy_Tests(unittest.TestCase):
-    def make_operator(self, **kwargs):
-        np.random.seed(kwargs['seed'])
-        S = ift.ScalingOperator(1., kwargs['space'])
-        s = S.draw_sample()
-        return ift.MultiField.from_dict({kwargs['space_key']: s})
-
-    @expand(product(
-        [ift.GLSpace(15),
-         ift.RGSpace(64, distances=.789),
-         ift.RGSpace([32, 32], distances=.789)],
-        [4, 78, 23]
-        ))
-    def testGaussian(self, space, seed):
-        op = self.make_operator(
-            space_key='s1', space=space, seed=seed)['s1']
-        energy = ift.GaussianEnergy(domain=space)
-        ift.extra.check_value_gradient_consistency(energy, op)
-
-#     @expand(product(
-#         [ift.GLSpace(15),
-#          ift.RGSpace(64, distances=.789),
-#          ift.RGSpace([32, 32], distances=.789)],
-#         [4, 78, 23]
-#         ))
-#     def testQuadratic(self, type1, space, seed):
-#         np.random.seed(seed)
-#         S = ift.ScalingOperator(1., space)
-#         s = [S.draw_sample() for _ in range(3)]
-#         energy = ift.QuadraticEnergy(s[0], ift.makeOp(s[1]), s[2])
-#         ift.extra.check_value_gradient_consistency(energy)
-
-    @expand(
-        product([
-            ift.GLSpace(15),
-            ift.RGSpace(64, distances=.789),
-            ift.RGSpace([32, 32], distances=.789)
-        ], [4, 78, 23]))
-    def testInverseGammaLikelihood(self, space, seed):
-        op = self.make_operator(space_key='s1', space=space, seed=seed)['s1']
-        op = op.exp()
-        d = np.random.normal(10, size=space.shape)**2
-        d = ift.Field.from_global_data(space, d)
-        energy = ift.InverseGammaLikelihood(d)
-        ift.extra.check_value_gradient_consistency(energy, op, tol=1e-7)
-
-    @expand(product(
-        [ift.GLSpace(15),
-         ift.RGSpace(64, distances=.789),
-         ift.RGSpace([32, 32], distances=.789)],
-        [4, 78, 23]
-        ))
-    def testPoissonian(self, space, seed):
-        op = self.make_operator(
-            space_key='s1', space=space, seed=seed)['s1']
-        op = op.exp()
-        d = np.random.poisson(120, size=space.shape)
-        d = ift.Field.from_global_data(space, d)
-        energy = ift.PoissonianEnergy(d)
-        ift.extra.check_value_gradient_consistency(energy, op, tol=1e-7)
-
-    @expand(product(
-        [ift.GLSpace(15),
-         ift.RGSpace(64, distances=.789),
-         ift.RGSpace([32, 32], distances=.789)],
-        [4, 78, 23]
-        ))
-    def testHamiltonian_and_KL(self, space, seed):
-        op = self.make_operator(
-            space_key='s1', space=space, seed=seed)['s1']
-        op = op.exp()
-        lh = ift.GaussianEnergy(domain=space)
-        hamiltonian = ift.Hamiltonian(lh)
-        ift.extra.check_value_gradient_consistency(hamiltonian, op)
-        S = ift.ScalingOperator(1., space)
-        samps = [S.draw_sample() for i in range(3)]
-        kl = ift.SampledKullbachLeiblerDivergence(hamiltonian, samps)
-        ift.extra.check_value_gradient_consistency(kl, op)
-
-    @expand(product(
-        [ift.GLSpace(15),
-         ift.RGSpace(64, distances=.789),
-         ift.RGSpace([32, 32], distances=.789)],
-        [4, 78, 23]
-        ))
-    def testBernoulli(self, space, seed):
-        op = self.make_operator(
-            space_key='s1', space=space, seed=seed)['s1']
-        op = op.sigmoid()
-        d = np.random.binomial(1, 0.1, size=space.shape)
-        d = ift.Field.from_global_data(space, d)
-        energy = ift.BernoulliEnergy(d)
-        ift.extra.check_value_gradient_consistency(energy, op, tol=1e-6)

test/test_energy_gradients.py

+# 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-2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+import numpy as np
+import pytest
+
+import nifty5 as ift
+from itertools import product
+
+# Currently it is not possible to parametrize fixtures. But this will
+# hopefully be fixed in the future.
+# https://docs.pytest.org/en/latest/proposals/parametrize_with_fixtures.html
+
+SPACES = [
+    ift.GLSpace(15),
+    ift.RGSpace(64, distances=.789),
+    ift.RGSpace([32, 32], distances=.789)
+]
+SEEDS = [4, 78, 23]
+PARAMS = product(SEEDS, SPACES)
+
+
+@pytest.fixture(params=PARAMS)
+def field(request):
+    np.random.seed(request.param[0])
+    S = ift.ScalingOperator(1., request.param[1])
+    s = S.draw_sample()
+    return ift.MultiField.from_dict({'s1': s})['s1']
+
+
+def test_gaussian(field):
+    energy = ift.GaussianEnergy(domain=field.domain)
+    ift.extra.check_value_gradient_consistency(energy, field)
+
+
+def test_inverse_gamma(field):
+    field = field.exp()
+    space = field.domain
+    d = np.random.normal(10, size=space.shape)**2
+    d = ift.Field.from_global_data(space, d)
+    energy = ift.InverseGammaLikelihood(d)
+    ift.extra.check_value_gradient_consistency(energy, field, tol=1e-7)
+
+
+def testPoissonian(field):
+    field = field.exp()
+    space = field.domain
+    d = np.random.poisson(120, size=space.shape)
+    d = ift.Field.from_global_data(space, d)
+    energy = ift.PoissonianEnergy(d)
+    ift.extra.check_value_gradient_consistency(energy, field, tol=1e-7)
+
+
+def test_hamiltonian_and_KL(field):
+    field = field.exp()
+    space = field.domain
+    lh = ift.GaussianEnergy(domain=space)
+    hamiltonian = ift.Hamiltonian(lh)
+    ift.extra.check_value_gradient_consistency(hamiltonian, field)
+    S = ift.ScalingOperator(1., space)
+    samps = [S.draw_sample() for i in range(3)]
+    kl = ift.SampledKullbachLeiblerDivergence(hamiltonian, samps)
+    ift.extra.check_value_gradient_consistency(kl, field)
+
+
+def test_bernoulli(field):
+    field = field.sigmoid()
+    space = field.domain
+    d = np.random.binomial(1, 0.1, size=space.shape)
+    d = ift.Field.from_global_data(space, d)
+    energy = ift.BernoulliEnergy(d)
+    ift.extra.check_value_gradient_consistency(energy, field, tol=1e-6)

test/test_energies/test_map.py → test/test_gaussian_metric.py

@@ -15,61 +15,67 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
-from itertools import product
-from test.common import expand
+import numpy as np
+import pytest
 
 import nifty5 as ift
-import numpy as np
 
 
 def _flat_PS(k):
     return np.ones_like(k)
 
 
-class Energy_Tests(unittest.TestCase):
-    @expand(product([ift.GLSpace(15),
-                     ift.RGSpace(64, distances=.789),
-                     ift.RGSpace([32, 32], distances=.789)],
-                    ["tanh", "exp", ""],
-                    [1, 1e-2, 1e2],
-                    [4, 78, 23]))
-    def testGaussianEnergy(self, space, nonlinearity, noise, seed):
-        np.random.seed(seed)
-        dim = len(space.shape)
-        hspace = space.get_default_codomain()
-        ht = ift.HarmonicTransformOperator(hspace, target=space)
-        binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
-        pspace = ift.PowerSpace(hspace, binbounds=binbounds)
-        Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
-        xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
-        xi0_var = ift.Linearization.make_var(xi0)
+pmp = pytest.mark.parametrize
 
-        def pspec(k): return 1 / (1 + k**2)**dim
-        pspec = ift.PS_field(pspace, pspec)
-        A = Dist(ift.sqrt(pspec))
-        N = ift.ScalingOperator(noise, space)
-        n = N.draw_sample()
-        s = ht(ift.makeOp(A)(xi0_var))
-        R = ift.ScalingOperator(10., space)
 
-        def d_model():
-            if nonlinearity == "":
-                return R(ht(ift.makeOp(A)))
-            else:
-                tmp = ht(ift.makeOp(A))
-                nonlin = getattr(tmp, nonlinearity)()
-                return R(nonlin)
+@pmp('space', [
+    ift.GLSpace(15),
+    ift.RGSpace(64, distances=.789),
+    ift.RGSpace([32, 32], distances=.789)
+])
+@pmp('nonlinearity', ["tanh", "exp", ""])
+@pmp('noise', [1, 1e-2, 1e2])
+@pmp('seed', [4, 78, 23])
+def test_gaussian_energy(space, nonlinearity, noise, seed):
+    np.random.seed(seed)
+    dim = len(space.shape)
+    hspace = space.get_default_codomain()
+    ht = ift.HarmonicTransformOperator(hspace, target=space)
+    binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
+    pspace = ift.PowerSpace(hspace, binbounds=binbounds)
+    Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
+    xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
+    # FIXME Needed?
+    xi0_var = ift.Linearization.make_var(xi0)
 
-        d = d_model()(xi0) + n
+    def pspec(k):
+        return 1/(1 + k**2)**dim
 
-        if noise == 1:
-            N = None
+    pspec = ift.PS_field(pspace, pspec)
+    A = Dist(ift.sqrt(pspec))
+    N = ift.ScalingOperator(noise, space)
+    n = N.draw_sample()
+    # FIXME Needed?
+    s = ht(ift.makeOp(A)(xi0_var))
+    R = ift.ScalingOperator(10., space)
 
-        energy = ift.GaussianEnergy(d, N)(d_model())
+    def d_model():
         if nonlinearity == "":
-            ift.extra.check_value_gradient_metric_consistency(
-                energy, xi0, ntries=10)
+            return R(ht(ift.makeOp(A)))
         else:
-            ift.extra.check_value_gradient_consistency(
-                energy, xi0, ntries=10, tol=5e-8)
+            tmp = ht(ift.makeOp(A))
+            nonlin = getattr(tmp, nonlinearity)()
+            return R(nonlin)
+
+    d = d_model()(xi0) + n
+
+    if noise == 1:
+        N = None
+
+    energy = ift.GaussianEnergy(d, N)(d_model())
+    if nonlinearity == "":
+        ift.extra.check_value_gradient_metric_consistency(
+            energy, xi0, ntries=10)
+    else:
+        ift.extra.check_value_gradient_consistency(
+            energy, xi0, ntries=10, tol=5e-8)

test/test_minimization/test_minimizers.py

-# 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-2019 Max-Planck-Society
-#
-# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-
-import unittest
-from itertools import product
-from test.common import expand
-
-import nifty5 as ift
-import numpy as np
-from unittest import SkipTest
-from numpy.testing import assert_allclose, assert_equal
-
-IC = ift.GradientNormController(tol_abs_gradnorm=1e-5, iteration_limit=1000)
-
-spaces = [ift.RGSpace([1024], distances=0.123), ift.HPSpace(32)]
-
-minimizers = ['ift.VL_BFGS(IC)',
-              'ift.NonlinearCG(IC, "Polak-Ribiere")',
-              # 'ift.NonlinearCG(IC, "Hestenes-Stiefel"),
-              'ift.NonlinearCG(IC, "Fletcher-Reeves")',
-              'ift.NonlinearCG(IC, "5.49")',
-              'ift.L_BFGS_B(ftol=1e-10, gtol=1e-5, maxiter=1000)',
-              'ift.L_BFGS(IC)',
-              'ift.NewtonCG(IC)']
-
-newton_minimizers = ['ift.RelaxedNewton(IC)']
-quadratic_only_minimizers = ['ift.ConjugateGradient(IC)',
-                             'ift.ScipyCG(tol=1e-5, maxiter=300)']
-slow_minimizers = ['ift.SteepestDescent(IC)']
-
-
-class Test_Minimizers(unittest.TestCase):
-
-    @expand(product(minimizers + newton_minimizers +
-                    quadratic_only_minimizers + slow_minimizers, spaces))
-    def test_quadratic_minimization(self, minimizer, space):
-        np.random.seed(42)