extend support for nested _mul.tri (#797)

For a multiplied sample A * B, the .tri implementation added trivial
implementations for the situation that either A or B consisted of only
one point, or A is one-dimensional. This PR extends this logic to any
multiplication chain A0 * A1 * ... with support for all left and right
multiplications of single-point or one-dimensional samples. In
particular this extends support to 3D structured boundaries.

nutils/function.py

@@ -4336,7 +4336,16 @@ def take(array: IntoArray, indices: IntoArray, axis: Optional[int] = None) -> Ar
             axis = 0
         else:
             axis = numeric.normdim(array.ndim, axis)
-        indices = _Wrapper.broadcasted_arrays(evaluable.NormDim, array.shape[axis], Array.cast(indices, dtype=int))
+        length = array.shape[axis]
+        indices = util.deep_reduce(numpy.stack, indices)
+        if isinstance(indices, Array):
+            indices = _Wrapper.broadcasted_arrays(evaluable.NormDim, length, indices)
+        else:
+            indices = numpy.array(indices)
+            indices[indices < 0] += length
+            if (indices < 0).any() or (indices >= length).any():
+                raise ValueError('indices out of bounds')
+            indices = _Constant(indices)
         transposed = _Transpose.to_end(array, axis)
         taken = _Wrapper(evaluable.Take, transposed, _WithoutPoints(indices), shape=(*transposed.shape[:-1], *indices.shape), dtype=array.dtype)
         return _Transpose.from_end(taken, *range(axis, axis+indices.ndim))

nutils/points.py

@@ -91,6 +91,7 @@ def hull(self):
         defines a simplex by mapping vertices into the list of points.
         '''
 
+        assert self.ndims > 0, 'hull is not defined for a zero-dimensional point set'
         edge_vertices = numpy.arange(self.ndims+1).repeat(self.ndims).reshape(self.ndims, self.ndims+1).T  # ndims+1 x ndims
         edge_simplices = numpy.sort(self.tri, axis=1)[:, edge_vertices]  # nelems x ndims+1 x ndims
         elems, edges = divmod(numpy.lexsort(edge_simplices.reshape(-1, self.ndims).T), self.ndims+1)

nutils/sample.py

@@ -15,7 +15,7 @@
 set.
 '''
 
-from . import types, points, _util as util, function, evaluable, parallel, numeric, matrix, sparse, warnings
+from . import types, points, _util as util, function, evaluable, parallel, matrix, sparse, warnings
 from .pointsseq import PointsSequence
 from .transformseq import Transforms
 from ._backports import cached_property
@@ -595,6 +595,90 @@ def __call__(self, func: function.IntoArray) -> function.Array:
         return numpy.concatenate([self._sample1(func), self._sample2(func)])
 
 
+def _simplex_strip(strip):
+    # Helper function that creates simplices for an extruded simplex, with
+    # vertices arranged in a [2,n] shape (prepended with an arbitrary number of
+    # axes). The Strategy is to create the first simplex from the first vertex
+    # in layer 1 and all vertices from layer 2, the second from all but the
+    # first of layer 2 and the first two of layer 1, and so on until the last
+    # simplex consists of all vertices of layer 1 and the last of layer 2.
+
+    assert strip.dtype == int
+    *shape, m, n = strip.shape
+    assert m == 2
+    flat = strip.reshape((*shape, 2*n)) # ravel last two axes, reallocates if necessary
+    flat = numpy.ascontiguousarray(flat) # required for use as buffer
+    *strides, s = flat.strides
+    return numpy.ndarray(buffer=flat, dtype=int, shape=(*shape, n, n+1), strides=(*strides, s, s))
+
+
+def _mul_tri(tri1, tri2):
+    # Helper function that computes the tri1 x tri2 product. The indices should
+    # be pre-multiplied with the appropriate strides.
+
+    if tri2 is None: # multiplication with 'empty' right hand side
+        return tri1
+
+    if tri1.shape[1] > tri2.shape[1]: # swap to reduce cases below
+        tri1, tri2 = tri2, tri1
+
+    ndims1 = tri1.shape[1] - 1
+    ndims2 = tri2.shape[1] - 1
+
+    tri_outer = tri1[:,None,:,None] + tri2[None,:,None,:]
+
+    if ndims1 == 0: # Left multiplication by a 0D sample
+        # Multiply the left 0D tri by the right tri and maintain the latter's
+        # triangulation.
+        tri = tri_outer.reshape(-1, ndims2+1)
+    elif ndims1 == 1: # Left multiplication by a 1D sample
+        # Multiply the left 1D tri by the right tri and triangulate using a
+        # simplex strip.
+        tri = _simplex_strip(tri_outer).reshape(-1, ndims2+2)
+    else:
+        raise NotImplementedError(f'tri not supported for {ndims1}D x {ndims2}D multiplication')
+
+    assert tri.shape[1] == ndims1 + ndims2 + 1
+    return tri
+
+
+def _mul_hull(tri1, tri2, hull1, hull2):
+    # Helper function that computes the hull1 x hull2 product. The indices
+    # should be pre-multiplied with the appropriate strides. If either tri1 or
+    # tri2 represents a 0D triangulation (i.e. a point) then the corresponding
+    # hull value will be ignored.
+
+    if tri2 is None: # multiplication with 'empty' right hand side
+        return hull1
+
+    if tri1.shape[1] > tri2.shape[1]: # swap to reduce cases below
+        tri1, tri2 = tri2, tri1
+        hull1, hull2 = hull2, hull1
+
+    ndims1 = tri1.shape[1] - 1
+    ndims2 = tri2.shape[1] - 1
+
+    if ndims1 == 0: # Left multiplication by a 0D sample
+        # Multiply the left 0D tri by the right hull and maintain the latter's
+        # triangulation.
+        hull_outer = tri1[:,None,:,None] + hull2[None,:,None,:] # ...,1,ndims2
+        hull = hull_outer.reshape(-1, ndims2)
+    elif ndims1 == 1: # Left multiplication by a 1D sample
+        # 1. Multiply the left 1D tri by the right hull and triangulate using a
+        # simplex strip.
+        hull_outer = tri1[:,None,:,None] + hull2[None,:,None,:] # ...,2,ndims2
+        hull = _simplex_strip(hull_outer).reshape(-1, ndims2+1)
+        # 2. Multiply the left 0D hull by the right tri and maintain the
+        # latter's triangulation.
+        hull_outer = hull1[:,None,:,None] + tri2[None,:,None,:] # ...,1,ndims2+1
+        hull = numpy.concatenate([hull_outer.reshape(-1, ndims2+1), hull])
+    else:
+        raise NotImplementedError(f'hull not supported for {ndims1}D x {ndims2}D multiplication')
+
+    assert hull.shape[1] == ndims1 + ndims2
+    return hull
+
+
 class _Mul(_TensorialSample):
 
     def __init__(self, sample1: Sample, sample2: Sample) -> None:
@@ -625,49 +709,84 @@ def get_lower_args(self, __ielem: evaluable.Array) -> function.LowerArgs:
         ielem1, ielem2 = evaluable.divmod(__ielem, self._sample2.nelems)
         return self._sample1.get_lower_args(ielem1) | self._sample2.get_lower_args(ielem2)
 
+    @property
+    def _reversed_factors(self):
+        # Helper method that generates the factors of arbitrarily nested
+        # multiplications in reverse order.
+
+        for s in self._sample2, self._sample1:
+            if isinstance(s, _Mul):
+                yield from s._reversed_factors
+            else:
+                yield s
+
+    def _get_element_tri_hull(self, ielem: int, with_hull: bool) -> numpy.ndarray:
+        # Helper method that returns the element_tri and element_hull for a
+        # given element index, used by get_element_tri and get_element_hull.
+        #
+        # To save work in case only the element_tri is required, a None value
+        # is returned for the latter if the with_hull flag is set to False. The
+        # converse (returning only the hull) is not possible as construction of
+        # the hull implies construction of the tri.
+
+        if not 0 <= ielem < self.nelems:
+            raise IndexError('element number is out of bounds')
+
+        # We loop from the final factor back to the first because of the order
+        # in which both the element index and the element vertices are raveled.
+
+        tri = hull = None
+        stride = 1
+        for sample in self._reversed_factors:
+            ielem, i = divmod(ielem, sample.nelems) # i is the unraveled element index in sample
+            nverts = len(sample.getindex(i))
+            sample_tri = sample.get_element_tri(i) * stride
+            if with_hull:
+                sample_hull = sample.ndims and sample.get_element_hull(i) * stride
+                hull = _mul_hull(sample_tri, tri, sample_hull, hull)
+            tri = _mul_tri(sample_tri, tri)
+            stride *= nverts # update stride to include the element's vertex count
+        assert ielem == 0
+        return tri, hull
+
     def get_element_tri(self, ielem: int) -> numpy.ndarray:
-        if self._sample1.ndims == 1:
-            ielem1, ielem2 = divmod(ielem, self._sample2.nelems)
-            npoints2 = len(self._sample2.getindex(ielem2))
-            tri12 = self._sample1.get_element_tri(ielem1)[:, None, :, None] * npoints2 + self._sample2.get_element_tri(ielem2)[None, :, None, :]  # ntri1 x ntri2 x 2 x ndims
-            return numeric.overlapping(tri12.reshape(-1, 2*self.ndims), n=self.ndims+1).reshape(-1, self.ndims+1)
-        elif self._sample1.npoints == 1:
-            return self._sample2.get_element_tri(ielem)
-        elif self._sample2.npoints == 1:
-            return self._sample1.get_element_tri(ielem)
-        else:
-            return super().get_element_tri(ielem)
+        return self._get_element_tri_hull(ielem, with_hull=False)[0]
 
     def get_element_hull(self, ielem: int) -> numpy.ndarray:
-        if self._sample1.ndims == 1:
-            ielem1, ielem2 = divmod(ielem, self._sample2.nelems)
-            npoints2 = len(self._sample2.getindex(ielem2))
-            hull1 = self._sample1.get_element_hull(ielem1)[:, None, :, None] * npoints2 + self._sample2.get_element_tri(ielem2)[None, :, None, :]  # 2 x ntri2 x 1 x ndims
-            hull2 = self._sample1.get_element_tri(ielem1)[:, None, :, None] * npoints2 + self._sample2.get_element_hull(ielem2)[None, :, None, :]  # ntri1 x nhull2 x 2 x ndims-1
-            return numpy.concatenate([hull1.reshape(-1, self.ndims), numeric.overlapping(hull2.reshape(-1, 2*(self.ndims-1)), n=self.ndims).reshape(-1, self.ndims)])
-        elif self._sample1.npoints == 1:
-            return self._sample2.get_element_hull(ielem)
-        elif self._sample2.npoints == 1:
-            return self._sample1.get_element_hull(ielem)
-        else:
-            return super().get_element_hull(ielem)
+        return self._get_element_tri_hull(ielem, with_hull=True)[1]
+
+    def _tri_hull(self, with_hull) -> numpy.ndarray:
+        # Helper method that returns the tri and hull of the sample, used by
+        # the tri and hull. These properties replace the default implementation
+        # via get_element_tri and get_element_hull by a faster algorithm that
+        # applies the product structure directly on the level of the sample.
+        #
+        # To save work in case only tri is required, a None value is returned
+        # for hull if the with_hull flag is set to False. The converse
+        # (returning only hull) is not possible as construction of hull implies
+        # construction of tri.
+        #
+        # We loop from the final factor back to the first because of the order
+        # in which the sample points are raveled.
+
+        tri = hull = None
+        stride = 1
+        for sample in self._reversed_factors:
+            sample_tri = sample.tri * stride
+            if with_hull:
+                sample_hull = sample.ndims and sample.hull * stride
+                hull = _mul_hull(sample_tri, tri, sample_hull, hull)
+            tri = _mul_tri(sample_tri, tri)
+            stride *= sample.npoints # update stride to include the sample's point count
+        return tri, hull
 
     @property
     def tri(self) -> numpy.ndarray:
-        if self._sample1.ndims == 1:
-            tri12 = self._sample1.tri[:, None, :, None] * self._sample2.npoints + self._sample2.tri[None, :, None, :]  # ntri1 x ntri2 x 2 x ndims
-            return numeric.overlapping(tri12.reshape(-1, 2*self.ndims), n=self.ndims+1).reshape(-1, self.ndims+1)
-        else:
-            return super().tri
+        return self._tri_hull(with_hull=False)[0]
 
     @property
     def hull(self) -> numpy.ndarray:
-        if self._sample1.ndims == 1:
-            hull1 = self._sample1.hull[:, None, :, None] * self._sample2.npoints + self._sample2.tri[None, :, None, :]  # 2 x ntri2 x 1 x ndims
-            hull2 = self._sample1.tri[:, None, :, None] * self._sample2.npoints + self._sample2.hull[None, :, None, :]  # ntri1 x nhull2 x 2 x ndims-1
-            return numpy.concatenate([hull1.reshape(-1, self.ndims), numeric.overlapping(hull2.reshape(-1, 2*(self.ndims-1)), n=self.ndims).reshape(-1, self.ndims)])
-        else:
-            return super().hull
+        return self._tri_hull(with_hull=True)[1]
 
     def integral(self, func: function.IntoArray) -> function.Array:
         return self._sample1.integral(self._sample2.integral(func))

tests/test_sample.py

@@ -441,51 +441,90 @@ class Dummy(Sample):
             Dummy(('a', 'b'), 1, 1, 1) * Dummy(('b', 'c'), 1, 1, 1)
 
 
+@parametrize
 class rectilinear(TestCase):
 
+    _nsimplex = 1, 1, 2, 6, 24 # number of simplices required to cover n-cube
+
     def setUp(self):
         super().setUp()
-        self.domain, self.geom = mesh.rectilinear([2, 1])
-        self.bezier2 = self.domain.sample('bezier', 2)
-        self.bezier3 = self.domain.sample('bezier', 3)
-        self.gauss2 = self.domain.sample('gauss', 2)
+        self.ndims = len(self.shape)
+        self.nelems = numpy.prod(self.shape, dtype=int)
+        self.nbelems = 2 * sum(self.nelems // n for n in self.shape)
+        self.topo, self.geom = mesh.rectilinear(self.shape)
 
     def test_integrate(self):
-        area = self.gauss2.integrate(1)
-        self.assertLess(abs(area-2), 1e-15)
+        area = self.topo.integrate(1, degree=1)
+        self.assertAlmostEqual(area, self.nelems, places=15)
 
     def test_integral(self):
-        area = self.gauss2.integral(function.asarray(1)).eval()
-        self.assertLess(abs(area-2), 1e-15)
+        area = self.topo.integral(function.asarray(1), degree=1).eval()
+        self.assertAlmostEqual(area, self.nelems, places=15)
 
     def test_eval(self):
-        x = self.bezier3.eval(self.geom)
-        self.assertEqual(x.shape, (self.bezier3.npoints,)+self.geom.shape)
+        for n in 1, 2:
+            bezier = self.topo.sample('bezier', n+1)
+            x = bezier.eval(self.geom)
+            self.assertEqual(x.shape, (bezier.npoints, *self.geom.shape))
+
+    def test_bezier(self):
+        for n in 1, 2:
+            bezier = self.topo.sample('bezier', n+1)
+            self.assertEqual(bezier.npoints, self.nelems * (n+1)**self.ndims)
 
     def test_tri(self):
-        self.assertEqual(len(self.bezier2.tri), 4)
-        self.assertEqual(len(self.bezier3.tri), 16)
+        for n in 1, 2:
+            bezier = self.topo.sample('bezier', n+1)
+            tri = bezier.tri
+            self.assertEqual(len(tri), self.nelems * self._nsimplex[self.ndims] * n**self.ndims)
+            self.assertAllEqual(numpy.unique(tri), numpy.arange(bezier.npoints))
+
+    def test_bnd_tri(self):
+        for n in 1, 2:
+            bezier = self.topo.boundary.sample('bezier', n+1)
+            tri = bezier.tri
+            self.assertEqual(len(tri), self.nbelems * self._nsimplex[self.ndims-1] * n**(self.ndims-1))
+            self.assertAllEqual(numpy.unique(tri), numpy.arange(bezier.npoints))
 
     def test_hull(self):
-        self.assertEqual(len(self.bezier2.hull), 8)
-        self.assertEqual(len(self.bezier3.hull), 16)
+        for n in 1, 2:
+            bezier = self.topo.sample('bezier', n+1)
+            hull = bezier.hull
+            self.assertEqual(len(hull), self.nelems * self._nsimplex[self.ndims-1] * 2 * self.ndims * n**(self.ndims-1))
+            if n == 1:
+                self.assertAllEqual(numpy.unique(hull), numpy.arange(bezier.npoints))
+
+    @parametrize.enable_if(lambda shape: len(shape) >= 3)
+    def test_bnd_hull(self):
+        for n in 1, 2:
+            bezier = self.topo.boundary.sample('bezier', n+1)
+            hull = bezier.hull
+            self.assertEqual(len(bezier.hull), self.nbelems * self._nsimplex[self.ndims-2] * n**(self.ndims-2) * 2 * (self.ndims-1))
+            if n == 1:
+                self.assertAllEqual(numpy.unique(hull), numpy.arange(bezier.npoints))
 
     def test_subset(self):
-        subset1 = self.bezier2.subset(numpy.eye(8)[0])
-        subset2 = self.bezier2.subset(numpy.eye(8)[1])
-        self.assertEqual(subset1.npoints, 4)
-        self.assertEqual(subset2.npoints, 4)
-        self.assertEqual(subset1, subset2)
+        bezier = self.topo.sample('bezier', 2)
+        subset1, subset2 = [bezier.subset(mask) for mask in numpy.eye(bezier.npoints, dtype=bool)[:2]]
+        self.assertEqual(subset1.npoints, 2**self.ndims)
+        self.assertEqual(subset2, subset1)
 
     def test_asfunction(self):
-        func = self.geom[0]**2 - self.geom[1]**2
-        values = self.gauss2.eval(func)
-        sampled = self.gauss2.asfunction(values)
+        func = sum(self.geom**2)
+        gauss = self.topo.sample('gauss', 2)
+        values = gauss.eval(func)
+        sampled = gauss.asfunction(values)
+        bezier = self.topo.sample('bezier', 2)
         with self.assertRaises(ValueError):
-            self.bezier2.eval(sampled)
-        self.assertAllEqual(self.gauss2.eval(sampled), values)
+            bezier.eval(sampled)
+        self.assertAllEqual(gauss.eval(sampled), values)
         arg = function.Argument('dofs', [2, 3])
-        self.assertTrue(evaluable.iszero(evaluable.asarray(self.gauss2(function.derivative(sampled, arg)))))
+        self.assertTrue(evaluable.iszero(evaluable.asarray(gauss(function.derivative(sampled, arg)))))
+
+rectilinear(shape=(4,))
+rectilinear(shape=(4,3))
+rectilinear(shape=(4,3,2))
+rectilinear(shape=(4,3,2,1))
 
 
 class integral(TestCase):