Merge pull request #691 from evalf/array

Array

CHANGELOG.rst

@@ -9,6 +9,42 @@ features in inverse chronological order.
 New in v8.0 (in development)
 ----------------------------
 
+- Deprecated: function methods that have Numpy equivalents
+
+  The :mod:`nutils.function` methods that have direct equivalents in the
+  :mod:`numpy` module (``function.sum``, ``function.sqrt``, ``function.sin``,
+  etc) have been deprecated in favour of using Numpy's methods (``numpy.sum``,
+  ``numpy.sqrt``, ``numpy.sin``, etc) and will be removed in the next release.
+  Ultimately, only methods that relate to the variable nature of function
+  arrays and therefore have no Numpy equivalent, such as ``function.grad`` and
+  ``function.normal``, will remain in the function module.
+
+  Be aware that some functions were not 100% equivalent to their Numpy
+  counterpart. For instance, ``function.max`` is the equivalent to
+  ``numpy.maximum``, as the deprecation message helpfully points out. More
+  problematically, ``function.dot`` behaves very differently from both
+  ``numpy.dot`` and ``numpy.matmul``. Porting the code over to equivalent
+  instructions will therefore require some attention.
+
+- Deprecated: Array.dot for ndim != 1
+
+  The :func:`nutils.function.Array.dot` method is incompatible with Numpy's
+  equivalent method for arrays of ndim != 1, or when axes are specified (which
+  Numpy does not allow). Aiming for 100% compatibility, the next release cycle
+  will remove the axis argument and temporarily forbid arguments of ndim != 1.
+  The release cycle thereafter will re-enable arguments with ndim != 1, with
+  logic equal to Numpy's method. In the meantime, the advice is to rely on
+  ``numpy.dot``, ``numpy.matmul`` or the ``@`` operator instead.
+
+- Deprecated: Array.sum for ndim > 1 without axis argument
+
+  The :func:`nutils.function.Array.sum` method by default operates on the last
+  axis. This is different from Numpy's behavour, which by default sums all
+  axes. Aiming for 100% compatibility, the next release cycle will make the
+  axis argument mandatory for any array of ndim > 1. The release cycle
+  thereafter will reintroduce the default value to match Numpy's. To prepare
+  for this, relying on the current default now triggers a deprecation warning.
+
 - New: iteration count via info.niter
 
   The info struct returned by ``solve_withinfo`` newly contains the amount of

examples/adaptivity.py

@@ -38,7 +38,7 @@ def main(etype: str, btype: str, degree: int, nrefine: int):
     domain, geom = mesh.unitsquare(2, etype)
 
     x, y = geom - .5
-    exact = (x**2 + y**2)**(1/3) * function.cos(function.arctan2(y+x, y-x) * (2/3))
+    exact = (x**2 + y**2)**(1/3) * numpy.cos(numpy.arctan2(y+x, y-x) * (2/3))
     domain = domain.trim(exact-1e-15, maxrefine=0)
     linreg = util.linear_regressor()
 
@@ -71,7 +71,7 @@ def main(etype: str, btype: str, degree: int, nrefine: int):
         args = solver.solve_linear('u:v', res, constrain=cons)
 
         ndofs = len(args['u'])
-        error = function.sqrt(domain.integral(['du du dV', '∇_k(du) ∇_k(du) dV'] @ ns, degree=7)).eval(**args)
+        error = numpy.sqrt(domain.integral(['du du dV', '∇_k(du) ∇_k(du) dV'] @ ns, degree=7)).eval(**args)
         rate, offset = linreg.add(numpy.log(ndofs), numpy.log(error))
         treelog.user(f'ndofs: {ndofs}, L2 error: {error[0]:.2e} ({rate[0]:.2f}), H1 error: {error[1]:.2e} ({rate[1]:.2f})')
 

examples/cahnhilliard.py

@@ -147,7 +147,7 @@ def main(size: Length, epsilon: Length, mobility: Time/Density, stens: Tension,
     domain, geom = mesh.unitsquare(nelems, etype)
     if circle:
         angle = (geom-.5) * (numpy.pi/2)
-        geom = .5 + function.sin(angle) * function.cos(angle)[[1, 0]] / numpy.sqrt(2)
+        geom = .5 + numpy.sin(angle) * numpy.cos(angle)[[1, 0]] / numpy.sqrt(2)
 
     bezier = domain.sample('bezier', 5)  # sample for surface plots
     grid = domain.locate(geom, numeric.simplex_grid([1, 1], 1/40), maxdist=1/nelems, skip_missing=True, tol=1e-5)  # sample for quivers

examples/coil.py

@@ -75,11 +75,11 @@ def main(nelems: int = 50, degree: int = 3, freq: float = 0., nturns: int = 1, r
     # `rcoil,zwires`.
 
     ns.zwires = (numpy.arange(nturns) - (nturns - 1) / 2) * 4 * rwire
-    ns.dwires = ns.rwire - numpy.sqrt((ns.r - ns.rcoil)**2 + functools.reduce(function.min, (ns.z - ns.zwires)**2))
+    ns.dwires = ns.rwire - numpy.sqrt((ns.r - ns.rcoil)**2 + functools.reduce(numpy.minimum, (ns.z - ns.zwires)**2))
     RZ = RZ.withsubdomain(coil=RZ[:-1, :-1].trim(ns.dwires/ns.rwire, maxrefine=4))
 
-    ns.rot = function.stack([function.scatter(function.trignormal(ns.θ), 3, [0, 1]), function.kronecker(1., 0, 3, 2)])
-    ns.eθ = function.stack(['-sin(θ)', 'cos(θ)', '0'] @ ns)
+    ns.rot = numpy.stack([function.scatter(function.trignormal(ns.θ), 3, [0, 1]), function.kronecker(1., 0, 3, 2)])
+    ns.eθ = numpy.stack(['-sin(θ)', 'cos(θ)', '0'] @ ns)
 
     X = RZ * REV
     ns.x = ns.rz @ ns.rot

examples/cylinderflow.py

@@ -19,7 +19,7 @@ def __init__(self, topo, ns, timestep, region=4., aspect=16/9, figscale=7.2, spa
         self.ns = ns
         self.figsize = aspect * figscale, figscale
         self.bbox = numpy.array([[-.5, aspect-.5], [-.5, .5]]) * region
-        self.bezier = topo.select(function.min(*(ns.x-self.bbox[:,0])*(self.bbox[:,1]-ns.x))).sample('bezier', 5)
+        self.bezier = topo.select(numpy.minimum(*(ns.x-self.bbox[:,0])*(self.bbox[:,1]-ns.x))).sample('bezier', 5)
         self.spacing = spacing
         self.timestep = timestep
         self.t = 0.

examples/drivencavity-compatible.py

@@ -46,7 +46,7 @@ def main(nelems: int, degree: int, reynolds: float):
     ns.add_field(('λ', 'μ'))
 
     ns.σ_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij'
-    ns.uwall = function.stack([domain.boundary.indicator('top'), 0])
+    ns.uwall = numpy.stack([domain.boundary.indicator('top'), 0])
     ns.N = 5 * degree * nelems  # nitsche constant based on element size = 1/nelems
     ns.nitsche_i = '(N v_i - (∇_j(v_i) + ∇_i(v_j)) n_j) / Re'
 

examples/platewithhole-nurbs.py

@@ -89,7 +89,7 @@ def main(nrefine: int, traction: float, radius: float, poisson: float):
     export.triplot('stressxx.png', x, σxx, tri=bezier.tri, hull=bezier.hull, clim=(numpy.nanmin(σxx), numpy.nanmax(σxx)))
 
     # evaluate error
-    err = function.sqrt(domain.integral(['du_k du_k dV', '∇_j(du_i) ∇_j(du_i) dV'] @ ns, degree=9)).eval(**args)
+    err = numpy.sqrt(domain.integral(['du_k du_k dV', '∇_j(du_i) ∇_j(du_i) dV'] @ ns, degree=9)).eval(**args)
     treelog.user('errors: L2={:.2e}, H1={:.2e}'.format(*err))
 
     return err, cons['u'], args['u']

examples/platewithhole.py

@@ -70,7 +70,7 @@ def main(nelems: int, etype: str, btype: str, degree: int, traction: float, maxr
     x, σxx = bezier.eval(['x_i', 'σ_00'] @ ns, **args)
     export.triplot('stressxx.png', x, σxx, tri=bezier.tri, hull=bezier.hull)
 
-    err = domain.integral(function.stack(['du_k du_k dV', '∇_j(du_i) ∇_j(du_i) dV'] @ ns), degree=max(degree, 3)*2).eval(**args)**.5
+    err = domain.integral(numpy.stack(['du_k du_k dV', '∇_j(du_i) ∇_j(du_i) dV'] @ ns), degree=max(degree, 3)*2).eval(**args)**.5
     treelog.user('errors: L2={:.2e}, H1={:.2e}'.format(*err))
 
     return err, cons['u'], args['u']

nutils/expression_v1.py

@@ -136,6 +136,7 @@
 import collections
 import functools
 import itertools
+import numpy
 import operator
 import types as builtin_types
 from typing import Any, Callable, Dict, List, Mapping, Optional, overload, Tuple, Union
@@ -1435,7 +1436,7 @@ def _eval_ast(ast, functions):
         value, = args
         return value
 
-    args = (_eval_ast(arg, functions) for arg in args)
+    args = [_eval_ast(arg, functions) for arg in args]
     if op == 'group':
         array, = args
         return array
@@ -1477,12 +1478,12 @@ def _eval_ast(ast, functions):
         return function.get(array, dim, index)
     elif op == 'trace':
         array, n1, n2 = args
-        return function.trace(array, n1, n2)
+        return numpy.trace(array, n1, n2)
     elif op == 'sum':
         array, axis = args
-        return function.sum(array, axis)
+        return numpy.sum(array, axis)
     elif op == 'concatenate':
-        return function.concatenate(args, axis=0)
+        return numpy.concatenate(args, axis=0)
     elif op == 'grad':
         array, geom = args
         return function.grad(array, geom)
@@ -1497,7 +1498,7 @@ def _eval_ast(ast, functions):
         return function.insertaxis(array, -1, length)
     elif op == 'transpose':
         array, trans = args
-        return function.transpose(array, trans)
+        return numpy.transpose(array, trans)
     elif op == 'jump':
         array, = args
         return function.jump(array)
@@ -1517,7 +1518,7 @@ def _eval_ast(ast, functions):
 def _sum_expr(arg: function.Array, *, consumes: int = 0) -> function.Array:
     if consumes == 0:
         raise ValueError('sum must consume at least one axis but got zero')
-    return function.sum(arg, range(arg.ndim-consumes, arg.ndim))
+    return numpy.sum(arg, range(arg.ndim-consumes, arg.ndim))
 
 
 def _norm2_expr(arg: function.Array, *, consumes: int = 0) -> function.Array:
@@ -1539,7 +1540,7 @@ def _J_expr(geom: function.Array, *, consumes: int = 0) -> function.Array:
 def _arctan2_expr(_a: function.Array, _b: function.Array) -> function.Array:
     a = function.Array.cast(_a)
     b = function.Array.cast(_b)
-    return function.arctan2(function._append_axes(a, b.shape), function._prepend_axes(b, a.shape))
+    return numpy.arctan2(function._append_axes(a, b.shape), function._prepend_axes(b, a.shape))
 
 
 class Namespace:
@@ -1671,13 +1672,13 @@ class Namespace:
     _re_assign = re.compile('^([a-zA-Zα-ωΑ-Ω][a-zA-Zα-ωΑ-Ω0-9]*)(_[a-z]+)?$')
 
     _default_functions = dict(
-        opposite=function.opposite, sin=function.sin, cos=function.cos,
-        tan=function.tan, sinh=function.sinh, cosh=function.cosh,
-        tanh=function.tanh, arcsin=function.arcsin, arccos=function.arccos,
-        arctan=function.arctan, arctan2=_arctan2_expr,
-        arctanh=function.arctanh, exp=function.exp, abs=function.abs,
-        ln=function.ln, log=function.ln, log2=function.log2, log10=function.log10,
-        sqrt=function.sqrt, sign=function.sign, d=function.d,
+        opposite=function.opposite, sin=numpy.sin, cos=numpy.cos,
+        tan=numpy.tan, sinh=numpy.sinh, cosh=numpy.cosh,
+        tanh=numpy.tanh, arcsin=numpy.arcsin, arccos=numpy.arccos,
+        arctan=numpy.arctan, arctan2=_arctan2_expr,
+        arctanh=numpy.arctanh, exp=numpy.exp, abs=numpy.abs,
+        ln=numpy.log, log=numpy.log, log2=numpy.log2, log10=numpy.log10,
+        sqrt=numpy.sqrt, sign=numpy.sign, d=function.d,
         surfgrad=function.surfgrad, n=function.normal, sum=_sum_expr,
         norm2=_norm2_expr, J=_J_expr,
     )