From 4fed5cc70a80226563f4ee548ca69cf1e0a836b3 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Thu, 8 Feb 2024 11:30:20 +0100 Subject: [PATCH 1/4] support locate for tensorial topologies This patch add support for locate on tensorial topologies, by moving the locate implementation from TransformChainsTopology to the Topology base class. This is mode possible by adding a new method _lower_args, as well as an implementation of _Mul._sample. --- nutils/topology.py | 167 +++++++++++++++++++++++------------------ tests/test_sample.py | 2 +- tests/test_topology.py | 9 +++ 3 files changed, 106 insertions(+), 72 deletions(-) diff --git a/nutils/topology.py b/nutils/topology.py index 0f2a7ad90..c893076da 100644 --- a/nutils/topology.py +++ b/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') diff --git a/tests/test_sample.py b/tests/test_sample.py index 915f1dc0d..c754577f6 100644 --- a/tests/test_sample.py +++ b/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)) diff --git a/tests/test_topology.py b/tests/test_topology.py index 2fe6aaebb..4ebcbf59e 100644 --- a/tests/test_topology.py +++ b/tests/test_topology.py @@ -966,6 +966,13 @@ 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): @@ -973,6 +980,8 @@ def test_missing_argument(self): @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)) From 128ab6c6656e00e55dd2c239554697a4326bb52f Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Thu, 8 Feb 2024 21:25:55 +0100 Subject: [PATCH 2/4] check for deprecation warnings in unit tests This patch changes deprecation warnings to errors in the testing.TestCase base class in order to catch features that are used by nutils but due to be removed by thrd parties. --- nutils/testing.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nutils/testing.py b/nutils/testing.py index c813f8068..451a97261 100644 --- a/nutils/testing.py +++ b/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) From 549fe2c6d0fff064c123bfe4eed41fbbddd82bae Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Thu, 8 Feb 2024 21:29:24 +0100 Subject: [PATCH 3/4] replace numpy.product by numpy.prod This patch changes all occurrences of numpy.product, due for removal in Numpy 2.0, by numpy.prod. --- nutils/function.py | 5 +++-- tests/test_basis.py | 4 ++-- tests/test_evaluable.py | 4 ++-- tests/test_function.py | 8 ++++---- 4 files changed, 11 insertions(+), 10 deletions(-) diff --git a/nutils/function.py b/nutils/function.py index 05e8aa073..e1482f0ad 100644 --- a/nutils/function.py +++ b/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) diff --git a/tests/test_basis.py b/tests/test_basis.py index a6675baf8..121229fc7 100644 --- a/tests/test_basis.py +++ b/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 diff --git a/tests/test_evaluable.py b/tests/test_evaluable.py index d10e4b389..a3e12649d 100644 --- a/tests/test_evaluable.py +++ b/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)) diff --git a/tests/test_function.py b/tests/test_function.py index e2f21baa5..ae69d124e 100644 --- a/tests/test_function.py +++ b/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)) From 6139944ced60aa6d9d4fc9317bd12c9ff58d38c0 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Thu, 8 Feb 2024 21:51:09 +0100 Subject: [PATCH 4/4] add SI support for numpy.prod --- nutils/SI.py | 7 +++++++ tests/test_SI.py | 5 +++++ 2 files changed, 12 insertions(+) diff --git a/nutils/SI.py b/nutils/SI.py index d0e0a627e..a86104001 100644 --- a/nutils/SI.py +++ b/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]) diff --git a/tests/test_SI.py b/tests/test_SI.py index 7c42e22e0..bb030e4dc 100644 --- a/tests/test_SI.py +++ b/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])))