support locate on tensorial topologies (#851)

nutils/SI.py

@@ -391,6 +391,13 @@ def __stack_like(op, *args, **kwargs):
             raise TypeError(f'incompatible arguments for {op.__name__}: ' + ', '.join(dim.__name__ for dim in dims))
         return dims[0].wrap(op(arg0, *args[1:], **kwargs))
 
+    @register('prod', 'product')
+    def __prod(op, a, axis=None, *args, **kwargs):
+        (dim, arg), = Quantity.__unpack(a)
+        axes = range(arg.ndim) if axis is None else axis if isinstance(axis, tuple) else (axis,)
+        n = functools.reduce(operator.mul, [arg.shape[axis] for axis in axes], 1)
+        return (dim**n).wrap(op(arg, axis, *args, **kwargs))
+
     @register('curvature')
     def __evaluate(op, *args, **kwargs):
         (dim0, arg0), = Quantity.__unpack(args[0])

nutils/__init__.py

@@ -1,4 +1,4 @@
 'Numerical Utilities for Finite Element Analysis'
 
-__version__ = version = '9a14'
+__version__ = version = '9a15'
 version_name = 'jook-sing'

nutils/function.py

@@ -353,7 +353,7 @@ def prod(self, __axis: int) -> 'Array':
         :class:`Array`
         '''
 
-        return numpy.product(self, __axis)
+        return numpy.prod(self, __axis)
 
     def dot(self, __other: IntoArray, axes: Optional[Union[int, Sequence[int]]] = None) -> 'Array':
         '''Return the inner product of the arguments over the given axes, elementwise over the remanining axes.
@@ -3113,7 +3113,8 @@ def sum(arg: IntoArray, axis: Optional[Union[int, Sequence[int]]] = None) -> Arr
         return summed
 
     @implements(numpy.product)
-    def product(arg: IntoArray, axis: int) -> Array:
+    @implements(numpy.prod)
+    def prod(arg: IntoArray, axis: int) -> Array:
         arg = Array.cast(arg)
         if arg.dtype == bool:
             arg = arg.astype(int)

nutils/testing.py

@@ -184,6 +184,7 @@ def setUp(self):
         self.enter_context(treelog.set(treelog.LoggingLog('nutils')))
         self.enter_context(_builtin_warnings.catch_warnings())
         _builtin_warnings.simplefilter('error', warnings.NutilsWarning)
+        _builtin_warnings.simplefilter('error', DeprecationWarning)
 
     def assertAllEqual(self, actual, desired):
         actual = numpy.asarray(actual)

nutils/topology.py

@@ -785,6 +785,78 @@ def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weigh
         located : :class:`nutils.sample.Sample`
         '''
 
+        if max(tol, eps) <= 0:
+            raise ValueError('locate requires either tol or eps to be strictly positive')
+        coords = numpy.asarray(coords, dtype=float)
+        if geom.ndim == 0:
+            geom = geom[_]
+            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)
+        ielems = parallel.shempty(len(coords), dtype=int)
+        points = parallel.shempty((len(coords), len(geom)), dtype=float)
+        _ielem = evaluable.Argument('_locate_ielem', shape=(), dtype=int)
+        _point = evaluable.Argument('_locate_point', shape=(evaluable.constant(self.ndims),))
+        egeom = geom.lower(self._lower_args(_ielem, _point))
+        xJ = evaluable.Tuple((egeom, evaluable.derivative(egeom, _point))).simplified
+        with parallel.ctxrange('locating', len(coords)) as ipoints:
+            for ipoint in ipoints:
+                xt = coords[ipoint]  # target
+                dist = numpy.linalg.norm(centroids - xt, axis=1)
+                for ielem in numpy.argsort(dist) if maxdist is None \
+                        else sorted((dist < maxdist).nonzero()[0], key=dist.__getitem__):
+                    ref = self.references[ielem]
+                    arguments['_locate_ielem'] = ielem
+                    arguments['_locate_point'] = p = numpy.array(ref.centroid)
+                    ex = ep = numpy.inf
+                    iiter = 0
+                    while ex > tol and ep > eps:  # newton loop
+                        if iiter > maxiter > 0:
+                            break  # maximum number of iterations reached
+                        iiter += 1
+                        xp, Jp = xJ.eval(**arguments)
+                        dx = xt - xp
+                        ex0 = ex
+                        ex = numpy.linalg.norm(dx)
+                        if ex >= ex0:
+                            break  # newton is diverging
+                        try:
+                            dp = numpy.linalg.solve(Jp, dx)
+                        except numpy.linalg.LinAlgError:
+                            break  # jacobian is singular
+                        ep = numpy.linalg.norm(dp)
+                        p += dp  # NOTE: modifies arguments['_locate_point'] in place
+                    else:
+                        if ref.inside(p, max(eps, ep)):
+                            ielems[ipoint] = ielem
+                            points[ipoint] = p
+                            break
+                else:
+                    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 _lower_args(self, ielem, point):
+        raise NotImplementedError
+
+    def _sample(self, ielems, coords, weights=None):
         raise NotImplementedError
 
     @cached_property
@@ -985,8 +1057,11 @@ def withgroups(self, vgroups: Mapping[str, Union[str, Topology]] = {}, bgroups:
         def indicator(self, subtopo: Union[str, Topology]) -> Topology:
             raise SkipTest('`{}` does not implement `Topology.indicator`'.format(type(self).__qualname__))
 
-        def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, ischeme=None, scale=None, skip_missing=False) -> Sample:
-            raise SkipTest('`{}` does not implement `Topology.locate`'.format(type(self).__qualname__))
+        def _lower_args(self, ielem, point):
+            raise SkipTest('`{}` does not implement `Topology._lower_args`'.format(type(self).__qualname__))
+
+        def _sample(self, ielems, coords, weights=None):
+            raise SkipTest('`{}` does not implement `Topology._sample`'.format(type(self).__qualname__))
 
 else:
     _TensorialTopology = Topology
@@ -1241,6 +1316,23 @@ def basis(self, name: str, degree: Union[int, Sequence[int]], **kwargs) -> funct
     def sample(self, ischeme: str, degree: int) -> Sample:
         return self.topo1.sample(ischeme, degree) * self.topo2.sample(ischeme, degree)
 
+    def _lower_args(self, ielem, point):
+        ielem1, ielem2 = evaluable.divmod(ielem, len(self.topo2))
+        largs1 = self.topo1._lower_args(ielem1, point[:self.topo1.ndims])
+        largs2 = self.topo2._lower_args(ielem2, point[self.topo1.ndims:])
+        return largs1 | largs2
+
+    def _sample(self, ielems, coords, weights=None):
+        ielems1, ielems2 = divmod(ielems, len(self.topo2))
+        sample1 = self.topo1._sample(ielems1, coords[...,:self.topo1.ndims], weights)
+        sample2 = self.topo2._sample(ielems2, coords[...,self.topo1.ndims:])
+        # NOTE: zip creates a sample with ndims set equal to that of the first
+        # argument, assuming (without verifying) that all samples are of the
+        # same dimension. Here this assumption is incorrect, as ndims should be
+        # the sum of the two samples, but since zipped samples do not support
+        # tri and hull the mistake is without consequences.
+        return Sample.zip(sample1, sample2)
+
 
 class _Take(_TensorialTopology):
 
@@ -1534,75 +1626,8 @@ def indicator(self, subtopo):
         values[numpy.fromiter(map(self.transforms.index, subtopo.transforms), dtype=int)] = 1
         return function.get(values, 0, self.f_index)
 
-    @log.withcontext
-    def locate(self, geom, coords, *, tol=0, eps=0, maxiter=0, arguments=None, weights=None, maxdist=None, skip_missing=False):
-        if max(tol, eps) <= 0:
-            raise ValueError('locate requires either tol or eps to be strictly positive')
-        coords = numpy.asarray(coords, dtype=float)
-        if geom.ndim == 0:
-            geom = geom[_]
-            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)
-        ielems = parallel.shempty(len(coords), dtype=int)
-        points = parallel.shempty((len(coords), len(geom)), dtype=float)
-        _ielem = evaluable.Argument('_locate_ielem', shape=(), dtype=int)
-        _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
-        with parallel.ctxrange('locating', len(coords)) as ipoints:
-            for ipoint in ipoints:
-                xt = coords[ipoint]  # target
-                dist = numpy.linalg.norm(centroids - xt, axis=1)
-                for ielem in numpy.argsort(dist) if maxdist is None \
-                        else sorted((dist < maxdist).nonzero()[0], key=dist.__getitem__):
-                    ref = self.references[ielem]
-                    arguments['_locate_ielem'] = ielem
-                    arguments['_locate_point'] = p = numpy.array(ref.centroid)
-                    ex = ep = numpy.inf
-                    iiter = 0
-                    while ex > tol and ep > eps:  # newton loop
-                        if iiter > maxiter > 0:
-                            break  # maximum number of iterations reached
-                        iiter += 1
-                        xp, Jp = xJ.eval(**arguments)
-                        dx = xt - xp
-                        ex0 = ex
-                        ex = numpy.linalg.norm(dx)
-                        if ex >= ex0:
-                            break  # newton is diverging
-                        try:
-                            dp = numpy.linalg.solve(Jp, dx)
-                        except numpy.linalg.LinAlgError:
-                            break  # jacobian is singular
-                        ep = numpy.linalg.norm(dp)
-                        p += dp  # NOTE: modifies arguments['_locate_point'] in place
-                    else:
-                        if ref.inside(p, max(eps, ep)):
-                            ielems[ipoint] = ielem
-                            points[ipoint] = p
-                            break
-                else:
-                    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 _lower_args(self, ielem, point):
+        return function.LowerArgs.for_space(self.space, (self.transforms, self.opposites), ielem, point)
 
     def _sample(self, ielems, coords, weights=None):
         index = numpy.argsort(ielems, kind='stable')

tests/test_SI.py

@@ -145,6 +145,11 @@ def test_sum(self):
         self.assertTrue(numpy.all(numpy.sum(SI.units.kg * numpy.arange(6).reshape(2, 3), 0) == SI.units.kg * numpy.array([3, 5, 7])))
         self.assertTrue(numpy.all(numpy.sum(SI.units.kg * numpy.arange(6).reshape(2, 3), 1) == SI.units.kg * numpy.array([3, 12])))
 
+    def test_prod(self):
+        self.assertTrue(numpy.all(numpy.prod(SI.units.m * numpy.arange(1,7).reshape(2, 3), axis=0) == numpy.array([4, 10, 18]) * SI.units.m**2))
+        self.assertTrue(numpy.all(numpy.prod(SI.units.m * numpy.arange(1,7).reshape(2, 3), axis=1) == numpy.array([6, 120]) * SI.units.m**3))
+        self.assertTrue(numpy.all(numpy.prod(SI.units.m * numpy.arange(1,7).reshape(2, 3)) == 720 * SI.units.m**6))
+
     def test_mean(self):
         self.assertTrue(numpy.all(numpy.mean(SI.units.kg * numpy.arange(6).reshape(2, 3), 0) == SI.units.kg * numpy.array([1.5, 2.5, 3.5])))
         self.assertTrue(numpy.all(numpy.mean(SI.units.kg * numpy.arange(6).reshape(2, 3), 1) == SI.units.kg * numpy.array([1, 4])))

tests/test_basis.py

@@ -275,7 +275,7 @@ def setUp(self):
         elif self.variant == 'tensor':
             structured, geom = mesh.rectilinear([numpy.linspace(0, 1, 5-i) for i in range(self.ndims)])
             domain = topology.ConnectedTopology(structured.space, structured.references, structured.transforms, structured.opposites, structured.connectivity)
-            nverts = numpy.product([5-i for i in range(self.ndims)])
+            nverts = numpy.prod([5-i for i in range(self.ndims)])
         elif self.variant == 'simplex':
             numpy.random.seed(0)
             nverts = 20
@@ -299,7 +299,7 @@ def test_ndofs(self):
         elif self.btype == 'discont':
             ndofs_by_ref = {
                 element.getsimplex(1)**self.ndims: (self.degree+1)**self.ndims,
-                element.getsimplex(self.ndims): numpy.product(self.degree+numpy.arange(self.ndims)+1) // numpy.product(numpy.arange(self.ndims)+1)}
+                element.getsimplex(self.ndims): numpy.prod(self.degree+numpy.arange(self.ndims)+1) // numpy.prod(numpy.arange(self.ndims)+1)}
             ndofs = sum(ndofs_by_ref[reference] for reference in self.domain.references)
         elif self.degree == 1:
             ndofs = self.nverts

tests/test_evaluable.py

@@ -512,8 +512,8 @@ def _check(name, op, n_op, *arg_values, hasgrad=True, zerograd=False, ndim=2):
 _check('arctan-complex', evaluable.arctan, numpy.arctan, ANC(4, 4))
 _check('ln', evaluable.ln, numpy.log, POS(4, 4))
 _check('ln-complex', evaluable.ln, numpy.log, NZC(4, 4))
-_check('product', lambda a: evaluable.product(a, 2), lambda a: numpy.product(a, 2), ANY(4, 3, 4))
-_check('product-complex', lambda a: evaluable.product(a, 2), lambda a: numpy.product(a, 2), ANC(4, 3, 4))
+_check('product', lambda a: evaluable.product(a, 2), lambda a: numpy.prod(a, 2), ANY(4, 3, 4))
+_check('product-complex', lambda a: evaluable.product(a, 2), lambda a: numpy.prod(a, 2), ANC(4, 3, 4))
 _check('sum', lambda a: evaluable.sum(a, 2), lambda a: a.sum(2), ANY(4, 3, 4))
 _check('sum-complex', lambda a: evaluable.sum(a, 2), lambda a: a.sum(2), ANC(4, 3, 4))
 _check('transpose1', lambda a: evaluable.transpose(a, [0, 1, 3, 2]), lambda a: a.transpose([0, 1, 3, 2]), ANY(2, 3, 4, 5))

tests/test_function.py

@@ -342,10 +342,10 @@ def _check(name, op, n_op, *args):
 _check('sum-bool', lambda a: numpy.sum(function.Array.cast(a > 0), 2), lambda a: (a > 0).sum(2), ANY(4, 3, 4))
 _check('sum-complex', lambda a: numpy.sum(function.Array.cast(a), 2), lambda a: a.sum(2), ANC(4, 3, 4))
 _check('Array_sum', lambda a: function.Array.cast(a).sum(2), lambda a: a.sum(2), ANY(4, 3, 4))
-_check('product', lambda a: numpy.product(function.Array.cast(a), 2), lambda a: numpy.product(a, 2), ANY(4, 3, 4))
-_check('product-bool', lambda a: numpy.product(function.Array.cast(a > 0), 2), lambda a: numpy.product((a > 0), 2), ANY(4, 3, 4))
-_check('product-complex', lambda a: numpy.product(function.Array.cast(a), 2), lambda a: numpy.product(a, 2), ANC(4, 3, 4))
-_check('Array_prod', lambda a: function.Array.cast(a).prod(2), lambda a: numpy.product(a, 2), ANY(4, 3, 4))
+_check('product', lambda a: numpy.prod(function.Array.cast(a), 2), lambda a: numpy.prod(a, 2), ANY(4, 3, 4))
+_check('product-bool', lambda a: numpy.prod(function.Array.cast(a > 0), 2), lambda a: numpy.prod((a > 0), 2), ANY(4, 3, 4))
+_check('product-complex', lambda a: numpy.prod(function.Array.cast(a), 2), lambda a: numpy.prod(a, 2), ANC(4, 3, 4))
+_check('Array_prod', lambda a: function.Array.cast(a).prod(2), lambda a: numpy.prod(a, 2), ANY(4, 3, 4))
 
 _check('dot', lambda a, b: numpy.dot(a, function.Array.cast(b)), numpy.dot, ANY(1, 2, 5), ANY(3, 5, 4))
 _check('dot-complex', lambda a, b: numpy.dot(a, function.Array.cast(b)), numpy.dot, ANC(1, 2, 5), ANC(3, 5, 4))

tests/test_sample.py

@@ -280,7 +280,7 @@ def test_integrate(self):
         self.assertAlmostEqual(self.stitched.integrate(function.J(self.geomX)), 5/9)  # NOTE: != norm(slope)
 
     def test_nested(self):
-        with self.assertRaisesRegex(ValueError, 'Nested integrals or samples in the same space: X, Y.'):
+        with self.assertRaisesRegex(ValueError, 'Nested integrals or samples in the same space: X.*, Y.'):
             self.stitched.integral(self.stitched.integral(1)).eval()
         topoZ, geomZ = mesh.line(2, space='Z')
         inner = self.stitched.integral((geomZ - self.geomX) * function.J(self.geomY))

tests/test_topology.py

@@ -966,13 +966,22 @@ def test_boundary_scalar(self):
         located = sample.eval(self.geom[1], scale=.123)
         self.assertAllAlmostEqual(located, target)
 
+    def test_integrate(self):
+        target = numpy.array([(.2, .3), (.1, .9), (0, 1), (.1, .3)])
+        weights = numpy.array([.1, .2, .3, .4])
+        sample = self.domain.locate(self.geom, target, eps=1e-15, tol=1e-12, arguments=dict(scale=.123), weights=weights)
+        integrated = sample.integrate(self.geom, scale=.123)
+        self.assertAllAlmostEqual(integrated, weights @ target)
+
     def test_missing_argument(self):
         target = numpy.array([(.2, .3)])
         with self.assertRaises(Exception):
             self.domain.locate(self.geom, target, eps=1e-15, tol=1e-12)
 
     @parametrize.enable_if(lambda etype, mode, **kwargs: etype == 'square' and mode != 'nonlinear')
     def test_detect_linear(self):
+        if os.environ.get('NUTILS_TENSORIAL'):
+            self.skipTest('linear detection not supported yet in tensorial mode')
         target = numpy.array([(.2, .3)])
         with self.assertLogs('nutils', level='DEBUG') as cm:
             self.domain.locate(self.geom, target, eps=1e-15, tol=1e-12, arguments=dict(scale=.123))