add evaluable.Singular

This patch adds evaluable.Singular and singular_like to emit a singular matrix
warning at runtime and return nan values in accordance with numeric.inv. This
allows simplifications of inverses that are known a priori to be singular.

nutils/evaluable.py

@@ -1196,6 +1196,9 @@ def _determinant(self, axis1, axis2):
     def _inverse(self, axis1, axis2):
         if axis1 < self.ndim-1 and axis2 < self.ndim-1:
             return InsertAxis(inverse(self.func, (axis1, axis2)), self.length)
+        # either axis1 or axis2 is inserted
+        if self.length._intbounds[0] > 1: # matrix is at least 2x2
+            return singular_like(self)
 
     def _loopsum(self, index):
         return InsertAxis(loop_sum(self.func, index), self.length)
@@ -1443,11 +1446,13 @@ def __init__(self, func: Array):
         super().__init__(args=(func,), shape=func.shape, dtype=func.dtype)
 
     def _simplified(self):
+        if _equals_scalar_constant(self.func.shape[-1], 1):
+            return reciprocal(self.func)
+        if _equals_scalar_constant(self.func.shape[-1], 0):
+            return singular_like(self)
         result = self.func._inverse(self.ndim-2, self.ndim-1)
         if result is not None:
             return result
-        if _equals_scalar_constant(self.func.shape[-1], 1):
-            return reciprocal(self.func)
 
     evalf = staticmethod(numeric.inv)
 
@@ -2763,6 +2768,17 @@ def _intbounds_impl(self):
         return intbounds[self.index] if intbounds else (float('-inf'), float('inf'))
 
 
+class Singular(Array):
+
+    def __init__(self, dtype):
+        assert dtype in (float, complex), 'Singular dtype must be float or complex'
+        super().__init__(args=(), shape=(), dtype=dtype)
+
+    def evalf(self):
+        warnings.warn('singular matrix', RuntimeWarning)
+        return numpy.array(numpy.nan, self.dtype)
+
+
 class Zeros(Array):
     'zero'
 
@@ -2839,6 +2855,9 @@ def _assparse(self):
     def _intbounds_impl(self):
         return 0, 0
 
+    def _inverse(self, axis1, axis2):
+        return singular_like(self)
+
 
 class Inflate(Array):
 
@@ -4646,6 +4665,10 @@ def ones_like(arr):
     return ones(arr.shape, arr.dtype)
 
 
+def singular_like(arr):
+    return appendaxes(Singular(arr.dtype), arr.shape)
+
+
 def reciprocal(arg):
     arg = asarray(arg)
     if arg.dtype in (bool, int):