Merge remote-tracking branch 'origin/NIFTy_5' into mr_docs

.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:

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/test_operators/test_adjoint.py

@@ -38,9 +38,10 @@ class Consistency_Tests(unittest.TestCase):
     def testLOSResponse(self, sp, dtype):
         starts = np.random.randn(len(sp.shape), 10)
         ends = np.random.randn(len(sp.shape), 10)
-        sigma_low = 1e-4*np.random.randn(10)
-        sigma_ups = 1e-5*np.random.randn(10)
-        op = ift.LOSResponse(sp, starts, ends, sigma_low, sigma_ups)
+        sigmas = 1e-4*np.random.randn(10)**2
+        while np.any(1./(1./np.linalg.norm(ends-starts,axis=0)-2.*sigmas)<0):
+            sigmas /= 2.
+        op = ift.LOSResponse(sp, starts, ends, sigmas, 2.)
         ift.extra.consistency_check(op, dtype, dtype)
 
     @expand(product(_h_spaces + _p_spaces + _pow_spaces,