Locate missing (#775)

This PR fixes the behaviour of `SubsetTopology.locate` when
`skip_missing` is set. Formerly it would fail if points that were found
on the base topology (where `skip_missing` was honoured) but not in the
subset. With this patch the flag instead triggers a post-processing step
that prunes the initial sample.

CHANGELOG

@@ -5,6 +5,14 @@ The following overview lists user facing changes as well as newly added
 features in inverse chronological order.
 
 
+FIXED: locate points on trimmed topologies with the skip_missing flag set
+
+The `locate` method has a `skip_missing` argument that instructs the method to
+silently drop points that can not be located on the topology. This setting was
+partially ignored by trimmed topologies which could lead to a `LocateError`
+despite the flag being set. This issue is now fixed.
+
+
 CHANGED: solve, solve_withinfo arguments
 
 Solver methods newton, minimize and pseudotime have their function signature

nutils/topology.py

@@ -1494,6 +1494,8 @@ def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weigh
             coords = coords[..., _]
         if not geom.shape == coords.shape[1:] == (self.ndims,):
             raise ValueError('invalid geometry or point shape for {}D topology'.format(self.ndims))
+        if skip_missing and weights is not None:
+            raise ValueError('weights and skip_missing are mutually exclusive')
         arguments = dict(arguments or ())
         centroids = self.sample('_centroid', None).eval(geom, **arguments)
         assert len(centroids) == len(self)
@@ -1503,10 +1505,6 @@ def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weigh
         _point = evaluable.Argument('_locate_point', shape=(evaluable.constant(self.ndims),))
         egeom = geom.lower(function.LowerArgs.for_space(self.space, (self.transforms, self.opposites), _ielem, _point))
         xJ = evaluable.Tuple((egeom, evaluable.derivative(egeom, _point))).simplified
-        if skip_missing:
-            if weights is not None:
-                raise ValueError('weights and skip_missing are mutually exclusive')
-            missing = parallel.shzeros(len(coords), dtype=bool)
         with parallel.ctxrange('locating', len(coords)) as ipoints:
             for ipoint in ipoints:
                 xt = coords[ipoint]  # target
@@ -1540,14 +1538,21 @@ def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weigh
                             points[ipoint] = p
                             break
                 else:
-                    if skip_missing:
-                        missing[ipoint] = True
-                    else:
-                        raise LocateError('failed to locate point: {}'.format(xt))
-        if skip_missing:
-            ielems = ielems[~missing]
-            points = points[~missing]
-        return self._sample(ielems, points, weights)
+                    ielems[ipoint] = -1 # mark point as missing
+                    if not skip_missing:
+                        # rather than raising LocateError here, which
+                        # parallel.fork will reraise as a general Exception if
+                        # ipoint was assigned to a child process, we fast
+                        # forward through the remaining points to abandon the
+                        # loop and subsequently raise from the main process.
+                        for ipoint in ipoints:
+                            pass
+        if -1 not in ielems: # all points are found
+            return self._sample(ielems, points, weights)
+        elif skip_missing: # not all points are found and that's ok, we just leave those out
+            return self._sample(ielems[ielems != -1], points[ielems != -1])
+        else: # not all points are found and that's an error
+            raise LocateError(f'failed to locate point: {coords[ielems==-1][0]}')
 
     def _sample(self, ielems, coords, weights=None):
         index = numpy.argsort(ielems, kind='stable')
@@ -2586,16 +2591,37 @@ def basis(self, name, *args, **kwargs):
         basis = self.basetopo.basis(name, *args, **kwargs)
         return function.PrunedBasis(basis, self._indices, self.f_index, self.f_coords)
 
-    def locate(self, geom, coords, *, eps=0, **kwargs):
-        sample = self.basetopo.locate(geom, coords, eps=eps, **kwargs)
-        for isampleelem, (transforms, points) in enumerate(zip(sample.transforms[0], sample.points)):
+    def locate(self, geom, coords, *, eps=0, skip_missing=False, **kwargs):
+        sample = self.basetopo.locate(geom, coords, eps=eps, skip_missing=skip_missing, **kwargs)
+        missing = []
+        for isampleelem, (transforms, points_) in enumerate(zip(sample.transforms[0], sample.points)):
             ielem = self.basetopo.transforms.index(transforms)
             ref = self.refs[ielem]
             if ref != self.basetopo.references[ielem]:
-                for i, coord in enumerate(points.coords):
+                for i, coord in enumerate(points_.coords):
                     if not ref.inside(coord, eps):
-                        raise LocateError('failed to locate point: {}'.format(coords[sample.getindex(isampleelem)[i]]))
-        return sample
+                        if not skip_missing:
+                            raise LocateError('failed to locate point: {}'.format(coords[sample.getindex(isampleelem)[i]]))
+                        missing.append((isampleelem, i))
+        if not missing:
+            return sample
+        selection = numpy.ones(len(sample.points), dtype=bool)
+        newpoints = []
+        for isampleelem, points_ in enumerate(sample.points):
+            mymissing = [] # collect missing points for current element
+            for isampleelem_, i in missing[:points_.npoints]:
+                if isampleelem_ != isampleelem:
+                    break
+                mymissing.append(i)
+            if not mymissing: # no points are missing -> keep existing points object
+                newpoints.append(points_)
+            elif len(mymissing) < points_.npoints: # some points are missing -> create new CoordsPoints object
+                newpoints.append(points.CoordsPoints(points_.coords[~numeric.asboolean(mymissing, points_.npoints)]))
+            else: # all points are missing -> remove element from return sample
+                selection[isampleelem] = False
+            del missing[:len(mymissing)]
+        assert not missing
+        return Sample.new(sample.space, [trans[selection] for trans in sample.transforms], PointsSequence.from_iter(newpoints, sample.ndims))
 
 
 class RefinedTopology(TransformChainsTopology):

tests/test_topology.py

@@ -943,9 +943,13 @@ def test_invalidargs(self):
             self.domain.locate(self.geom, target, eps=1e-15, tol=1e-12, arguments=dict(scale=.123))
 
     def test_invalidpoint(self):
-        target = numpy.array([(.3, 1)])  # outside domain, but inside basetopo for mode==trimmed
-        with self.assertRaises(topology.LocateError):
+        target = numpy.array([(.3, 1), (.2, .3), (.1, .9), (0, 1), (.1, .3)])
+        # the first point is outside the domain, but inside basetopo for mode==trimmed
+        with self.subTest('skip_missing=False'), self.assertRaises(topology.LocateError):
             self.domain.locate(self.geom, target, eps=1e-15, tol=1e-12, arguments=dict(scale=.123))
+        with self.subTest('skip_missing=True'):
+            sample = self.domain.locate(self.geom, target, eps=1e-15, tol=1e-12, arguments=dict(scale=.123), skip_missing=True)
+            self.assertEqual(sample.npoints, 4)
 
     def test_boundary(self):
         target = numpy.array([(.2,), (.1,), (0,)])