diff --git a/nutils/evaluable.py b/nutils/evaluable.py
index 894bfa43f..8a3ca64e5 100644
--- a/nutils/evaluable.py
+++ b/nutils/evaluable.py
@@ -75,16 +75,51 @@ def asarray(arg):
 
 asarrays = types.tuple[asarray]
 
-def as_canonical_length(value):
-  if isarray(value):
-    return value._as_canonical_length
-  elif numeric.isint(value):
-    value = int(value) # Ensure this is an `int`, not `numpy.int64`.
-  else:
-    raise ValueError('length should be an `int` or `Array` with zero dimensions and dtype `int`, got {!r}'.format(value))
-  return value
+def asindex(arg):
+  arg = asarray(arg)
+  if arg.ndim or arg.dtype not in (int, bool): # NOTE: bool to be removed after introduction of Cast
+    raise ValueError('argument is not an index: {}'.format(arg))
+  return arg
+
+@types.apply_annotations
+def equalindex(n:asindex, m:asindex):
+  '''Compare two array indices.
+
+  Returns `True` if the two indices are certainly equal, `False` if they are
+  certainly not equal, or `None` if equality cannot be determined at compile
+  time.
+  '''
+
+  if n is m:
+    return True
+  n = n.simplified
+  m = m.simplified
+  if n is m:
+    return True
+  if n.arguments != m.arguments:
+    return False
+  if n.isconstant: # implies m.isconstant
+    return int(n) == int(m)
+
+asshape = types.tuple[asindex]
+
+@types.apply_annotations
+def equalshape(N:asshape, M:asshape):
+  '''Compare two array shapes.
+
+  Returns `True` if all indices are certainly equal, `False` if any indices are
+  certainly not equal, or `None` if equality cannot be determined at compile
+  time.
+  '''
 
-asshape = types.tuple[as_canonical_length]
+  assert len(N) == len(M)
+  retval = True
+  for eq in map(equalindex, N, M):
+    if eq == False:
+      return False
+    if eq == None:
+      retval = None
+  return retval
 
 class ExpensiveEvaluationWarning(warnings.NutilsInefficiencyWarning): pass
 
@@ -375,7 +410,7 @@ def simplified(obj):
     if isinstance(obj, Evaluable):
       retval = obj._simplified()
       if retval is not None and isinstance(obj, Array):
-        assert isinstance(retval, Array) and retval.shape == obj.shape, '{}._simplified resulted in shape change'.format(type(obj).__name__)
+        assert isinstance(retval, Array) and equalshape(retval.shape, obj.shape), '{}._simplified resulted in shape change'.format(type(obj).__name__)
       return retval
 
   def _simplified(self):
@@ -392,7 +427,7 @@ def _optimized_for_numpy1(obj: simplified.fget):
     if isinstance(obj, Evaluable):
       retval = obj._simplified() or obj._optimized_for_numpy()
       if retval is not None and isinstance(obj, Array):
-        assert isinstance(retval, Array) and retval.shape == obj.shape, '{0}._optimized_for_numpy or {0}._simplified resulted in shape change'.format(type(obj).__name__)
+        assert isinstance(retval, Array) and equalshape(retval.shape, obj.shape), '{0}._optimized_for_numpy or {0}._simplified resulted in shape change'.format(type(obj).__name__)
       return retval
 
   def _optimized_for_numpy(self):
@@ -430,7 +465,7 @@ def _combine_loop_concatenates(self, outer_exclude):
         # start, stop and the concatenation length) are formed by
         # `loop_concatenate`-ing `func.shape[-1]`. If the shape is constant,
         # this can be simplified to a `Range`.
-        data = Tuple((Tuple((lc.func, lc.start, lc.stop, asarray(lc.shape[-1]))) for lc in lcs))
+        data = Tuple((Tuple((lc.func, lc.start, lc.stop, lc.shape[-1])) for lc in lcs))
         # Combine `LoopConcatenate` instances in `data` excluding
         # `outer_exclude` and those that will be processed in a subsequent loop
         # (the remainder of `exclude`). The latter consists of loops that are
@@ -509,10 +544,10 @@ class SparseArray(Evaluable):
   'sparse array'
 
   @types.apply_annotations
-  def __init__(self, chunks:types.tuple[types.tuple[asarray]], shape:types.tuple[asarray], dtype:asdtype):
+  def __init__(self, chunks:types.tuple[asarrays], shape:asarrays, dtype:asdtype):
     self._shape = shape
     self._dtype = dtype
-    super().__init__(args=[Tuple(map(asarray, shape)), *map(Tuple, chunks)])
+    super().__init__(args=[Tuple(shape), *map(Tuple, chunks)])
 
   def evalf(self, shape, *chunks):
     length = builtins.sum(values.size for *indices, values in chunks)
@@ -697,10 +732,10 @@ def swapaxes(arg, axis1, axis2):
 class Axis(types.Immutable):
   __slots__ = 'length'
   @types.apply_annotations
-  def __init__(self, length:as_canonical_length):
+  def __init__(self, length:asindex):
     self.length = length
   def __str__(self):
-    return '?' if isarray(self.length) else str(self.length)
+    return str(int(self.length)) if self.length.isconstant else '?'
   def __repr__(self):
     return '{}{}'.format(type(self).__name__[0].lower(), self)
 
@@ -713,18 +748,18 @@ class Raveled(Axis):
   def __init__(self, shape:asshape):
     assert len(shape) == 2
     self.shape = shape
-    super().__init__(as_canonical_length(shape[0] * shape[1]))
+    super().__init__(shape[0] * shape[1])
 
 class Diagonal(Axis):
   __slots__ = ()
   @types.apply_annotations
-  def __init__(self, length:as_canonical_length, marker):
+  def __init__(self, length:asindex, marker):
     super().__init__(length)
 
 class Sparse(Axis):
   __slots__ = 'mask'
   @types.apply_annotations
-  def __init__(self, length:as_canonical_length, mask:types.arraydata=False):
+  def __init__(self, length:asindex, mask:types.arraydata=False):
     assert mask.dtype == bool
     self.mask = numpy.asarray(mask) # True for indices that are certain to be filled, used to detect dense addition
     super().__init__(length)
@@ -749,7 +784,7 @@ def _assparse(self):
         for *indices, values in chunks:
           assert len(indices) == self.ndim
           assert all(idx.dtype == int for idx in indices)
-          assert all(idx.shape == values.shape for idx in indices)
+          assert all(equalshape(idx.shape, values.shape) for idx in indices)
       elif chunks:
         assert len(chunks) == 1
         chunk, = chunks
@@ -808,7 +843,7 @@ class Array(Evaluable, metaclass=_ArrayMeta):
   '''
 
   __slots__ = '_axes', 'dtype', '__index'
-  __cache__ = 'blocks', 'assparse', '_assparse', '_as_canonical_length'
+  __cache__ = 'blocks', 'assparse', '_assparse'
 
   __array_priority__ = 1. # http://stackoverflow.com/questions/7042496/numpy-coercion-problem-for-left-sided-binary-operator/7057530#7057530
 
@@ -1005,17 +1040,6 @@ def as_evaluable_array(self):
 
     return self
 
-  @property
-  def _as_canonical_length(self):
-    'convert to an :class:`int` if possible'
-
-    if self.ndim != 0 or self.dtype != int:
-      raise ValueError('expected an `Array` with dimension zero and dtype `int` but got {!r}'.format(self))
-    value = self.simplified
-    if value.isconstant:
-      value = int(value.eval()) # Ensure this is an `int`, not `numpy.int64`.
-    return value
-
 class NPoints(Array):
   'The length of the points axis.'
 
@@ -1057,13 +1081,13 @@ class Normal(Array):
 
   @types.apply_annotations
   def __init__(self, lgrad:asarray):
-    assert lgrad.ndim >= 2 and lgrad.shape[-2] == lgrad.shape[-1]
+    assert lgrad.ndim >= 2 and equalindex(lgrad.shape[-2], lgrad.shape[-1])
     self.lgrad = lgrad
     super().__init__(args=[lgrad], shape=lgrad.shape[:-1], dtype=float)
 
   def evalf(self, lgrad):
     n = lgrad[...,-1]
-    if n.shape[-1] == 1: # geom is 1D
+    if equalindex(n.shape[-1], 1): # geom is 1D
       return numpy.sign(n)
     # orthonormalize n to G
     G = lgrad[...,:-1]
@@ -1074,7 +1098,7 @@ def evalf(self, lgrad):
     return numeric.normalize(n - v3)
 
   def _derivative(self, var, seen):
-    if self.shape[-1] == 1:
+    if equalindex(self.shape[-1], 1):
       return zeros(self.shape + var.shape)
     G = self.lgrad[...,:-1]
     m, n = G.shape[-2:]
@@ -1182,9 +1206,7 @@ class InsertAxis(Array):
   __slots__ = 'func', 'length'
 
   @types.apply_annotations
-  def __init__(self, func:asarray, length:asarray):
-    if length.ndim != 0 or length.dtype != int:
-      raise Exception('invalid length argument')
+  def __init__(self, func:asarray, length:asindex):
     self.func = func
     self.length = length
     super().__init__(args=[func, length], shape=func._axes+(Inserted(length),), dtype=func.dtype)
@@ -1441,7 +1463,7 @@ def __init__(self, func:asarray):
     super().__init__(args=[func], shape=func.shape[:-1], dtype=int if func.dtype == bool else func.dtype)
 
   def _simplified(self):
-    if self.func.shape[-1] == 1:
+    if equalindex(self.func.shape[-1], 1):
       return get(self.func, self.ndim, 0)
     return self.func._product()
 
@@ -1484,16 +1506,17 @@ def _derivative(self, var, seen):
 
 class LinearFrom(Array):
 
-  __slots__ = ()
+  __slots__ = 'todims', 'fromdims'
 
   @types.apply_annotations
   def __init__(self, trans:types.strict[TransformChain], todims:types.strictint, fromdims:types.strictint):
+    self.todims = todims
+    self.fromdims = fromdims
     super().__init__(args=[trans], shape=(todims, fromdims), dtype=float)
 
   def evalf(self, chain):
-    todims, fromdims = self.shape
-    assert not chain or chain[0].todims == todims
-    return transform.linearfrom(chain, fromdims)
+    assert not chain or chain[0].todims == self.todims
+    return transform.linearfrom(chain, self.fromdims)
 
 class Inverse(Array):
   '''
@@ -1505,7 +1528,7 @@ class Inverse(Array):
 
   @types.apply_annotations
   def __init__(self, func:asarray):
-    assert func.ndim >= 2 and func.shape[-1] == func.shape[-2]
+    assert func.ndim >= 2 and equalindex(func.shape[-1], func.shape[-2])
     self.func = func
     super().__init__(args=[func], shape=func.shape, dtype=float)
 
@@ -1570,7 +1593,7 @@ class Determinant(Array):
 
   @types.apply_annotations
   def __init__(self, func:asarray):
-    assert isarray(func) and func.ndim >= 2 and func.shape[-1] == func.shape[-2]
+    assert isarray(func) and func.ndim >= 2 and equalindex(func.shape[-1], func.shape[-2])
     self.func = func
     super().__init__(args=[func], shape=func.shape[:-2], dtype=_jointdtype(func.dtype, float))
 
@@ -1602,7 +1625,7 @@ class Multiply(Array):
   def __init__(self, funcs:types.frozenmultiset[asarray]):
     self.funcs = funcs
     func1, func2 = funcs
-    assert func1.shape == func2.shape
+    assert equalshape(func1.shape, func2.shape)
     axes = [axis1 if axis1 == axis2
        else axis1 if isinstance(axis1, Sparse)
        else axis2 if isinstance(axis2, Sparse)
@@ -1658,7 +1681,7 @@ def evalf(self, arr1, arr2):
 
   def _sum(self, axis):
     func1, func2 = self.funcs
-    if self.shape[axis] == 1:
+    if equalindex(self.shape[axis], 1):
       return multiply(get(func1, axis, 0), get(func2, axis, 0))
     if isinstance(func1._axes[axis], Inserted):
       return multiply(func1._uninsert(axis), func2.sum(axis))
@@ -1678,7 +1701,7 @@ def _add(self, other):
   def _determinant(self, axis1, axis2):
     func1, func2 = self.funcs
     axis1, axis2 = sorted([axis1, axis2])
-    if self.shape[axis1] == self.shape[axis2] == 1:
+    if equalindex(self.shape[axis1], 1) and equalindex(self.shape[axis2], 1):
       return Multiply([determinant(func1, (axis1, axis2)), determinant(func2, (axis1, axis2))])
     if all(isinstance(func1._axes[axis], Inserted) for axis in (axis1, axis2)):
       return Multiply([func1._uninsert(axis2)._uninsert(axis1)**self.shape[axis1], determinant(func2, (axis1, axis2))])
@@ -1763,7 +1786,7 @@ class Add(Array):
   def __init__(self, funcs:types.frozenmultiset[asarray]):
     self.funcs = funcs
     func1, func2 = funcs
-    assert func1.shape == func2.shape
+    assert equalshape(func1.shape, func2.shape)
     axes = [axis1 if axis1 == axis2 else Axis(axis1.length) for axis1, axis2 in zip(func1._axes, func2._axes)]
     sparse = [i for i, (axis1, axis2) in enumerate(zip(func1._axes, func2._axes)) if isinstance(axis1, Sparse) and isinstance(axis2, Sparse)]
     if len(sparse) > 1: # an addition of multiple sparse axes is always sparse
@@ -1829,7 +1852,7 @@ class Einsum(Array):
   __slots__ = 'args', 'out_idx', 'args_idx', '_einsumfmt', '_has_summed_axes'
 
   @types.apply_annotations
-  def __init__(self, args:types.tuple[asarray], args_idx:types.tuple[types.tuple[types.strictint]], out_idx:types.tuple[types.strictint]):
+  def __init__(self, args:asarrays, args_idx:types.tuple[types.tuple[types.strictint]], out_idx:types.tuple[types.strictint]):
     if len(args_idx) != len(args):
       raise ValueError('Expected one list of indices for every argument, but got {} and {}, respectively.'.format(len(args_idx), len(args)))
     for iarg, (idx, arg) in enumerate(zip(args_idx, args), 1):
@@ -1843,7 +1866,7 @@ def __init__(self, args:types.tuple[asarray], args_idx:types.tuple[types.tuple[t
       for i, length in zip(idx, arg.shape):
         if i not in lengths:
           lengths[i] = length
-        elif lengths[i] != length:
+        elif not equalindex(lengths[i], length):
           raise ValueError('Axes with index {} have different lengths.'.format(i))
     try:
       shape = [lengths[i] for i in out_idx]
@@ -1942,14 +1965,14 @@ class TakeDiag(Array):
   def __init__(self, func:asarray):
     if func.ndim < 2:
       raise Exception('takediag requires an argument of dimension >= 2')
-    if func.shape[-1] != func.shape[-2]:
+    if not equalindex(func.shape[-1], func.shape[-2]):
       raise Exception('takediag axes do not match')
     self.func = func
     shape = func._axes[:-2]+(Axis(func.shape[-1]),)
     super().__init__(args=[func], shape=shape, dtype=func.dtype)
 
   def _simplified(self):
-    if self.shape[-1] == 1:
+    if equalindex(self.shape[-1], 1):
       return Take(self.func, 0)
     return self.func._takediag(self.ndim-1, self.ndim)
 
@@ -2052,7 +2075,7 @@ class Power(Array):
 
   @types.apply_annotations
   def __init__(self, func:asarray, power:asarray):
-    assert func.shape == power.shape
+    assert equalshape(func.shape, power.shape)
     self.func = func
     self.power = power
     dtype = float if func.dtype == power.dtype == int else _jointdtype(func.dtype, power.dtype)
@@ -2123,9 +2146,10 @@ class Pointwise(Array):
   @types.apply_annotations
   def __init__(self, *args:asarrays):
     retval = self.evalf(*[numpy.ones((), dtype=arg.dtype) for arg in args])
-    shapes = set(arg.shape for arg in args)
-    assert len(shapes) == 1, 'pointwise arguments have inconsistent shapes'
+    shape0 = args[0].shape
+    assert all(equalshape(arg.shape, shape0) for arg in args[1:]), 'pointwise arguments have inconsistent shapes'
     shape = tuple(axes[0] if len(set(axes)) == 1 and not isinstance(axes[0], Sparse) else Axis(axes[0].length) for axes in zip(*(arg._axes for arg in args)))
+    # TODO check shape definition
     self.args = args
     super().__init__(args=args, shape=shape, dtype=retval.dtype)
 
@@ -2364,7 +2388,7 @@ class Eig(Evaluable):
 
   @types.apply_annotations
   def __init__(self, func:asarray, symmetric:bool=False):
-    assert func.ndim >= 2 and func.shape[-1] == func.shape[-2]
+    assert func.ndim >= 2 and equalindex(func.shape[-1], func.shape[-2])
     self.symmetric = symmetric
     self.func = func
     self._w_dtype = float if symmetric else complex
@@ -2418,7 +2442,7 @@ class Zeros(Array):
 
   @types.apply_annotations
   def __init__(self, shape:asshape, dtype:asdtype):
-    super().__init__(args=[asarray(sh) for sh in shape], shape=map(Sparse, shape), dtype=dtype)
+    super().__init__(args=shape, shape=map(Sparse, shape), dtype=dtype)
 
   def evalf(self, *shape):
     return numpy.zeros(shape, dtype=self.dtype)
@@ -2474,7 +2498,7 @@ def _determinant(self, axis1, axis2):
     shape = list(self.shape)
     assert axis1 != axis2
     length, = set(map(shape.pop, sorted((axis1, axis2), reverse=True)))
-    if length == 0:
+    if iszero(length):
       return ones(shape, self.dtype)
     else:
       return Zeros(shape, self.dtype)
@@ -2489,13 +2513,13 @@ class Inflate(Array):
   __cache__ = '_assparse'
 
   @types.apply_annotations
-  def __init__(self, func:asarray, dofmap:asarray, length:asarray):
-    if func.shape[func.ndim-dofmap.ndim:] != dofmap.shape:
+  def __init__(self, func:asarray, dofmap:asarray, length:asindex):
+    if not equalshape(func.shape[func.ndim-dofmap.ndim:], dofmap.shape):
       raise Exception('invalid dofmap')
     self.func = func
     self.dofmap = dofmap
     self.length = length
-    mask = not any(isinstance(ax, Sparse) for ax in func._axes[func.ndim-dofmap.ndim:]) and dofmap.isconstant and length.isconstant and numeric.asboolean(dofmap.eval().ravel(), length.eval(), ordered=False)
+    mask = not any(isinstance(ax, Sparse) for ax in func._axes[func.ndim-dofmap.ndim:]) and dofmap.isconstant and length.isconstant and numeric.asboolean(dofmap.eval().ravel(), int(length), ordered=False)
     shape = func._axes[:-1] + (Sparse(length, mask),)
     self.warn = not dofmap.isconstant
     super().__init__(args=[func,dofmap,length], shape=func._axes[:func.ndim-dofmap.ndim]+(Sparse(length, mask),), dtype=func.dtype)
@@ -2504,7 +2528,7 @@ def _simplified(self):
     if self.dofmap == Range(self.length):
       return self.func
     for axis in range(self.dofmap.ndim):
-      if self.dofmap.shape[axis] == 1:
+      if equalindex(self.dofmap.shape[axis], 1):
         return Inflate(_take(self.func, 0, self.func.ndim-self.dofmap.ndim+axis), _take(self.dofmap, 0, axis), self.length)
       if isinstance(self.func._axes[self.func.ndim-self.dofmap.ndim+axis], Sparse):
         items = [Inflate(f, _take(self.dofmap, ind, axis), self.shape[-1]) for ind, f in self.func._desparsify(self.func.ndim-self.dofmap.ndim+axis)]
@@ -2643,7 +2667,7 @@ def evalf(self, inflateidx, takeidx):
         for j in [subinflate[k]] if uniqueinflate else numpy.equal(inflateidx.ravel(), n).nonzero()[0]:
           newinflate.append(i)
           newtake.append(j)
-    return numpy.array(newtake, dtype=int), numpy.array(newinflate, dtype=int), numpy.array(len(newtake))
+    return numpy.array(newtake, dtype=int), numpy.array(newinflate, dtype=int), numpy.array(len(newtake), dtype=int)
 
 class Diagonalize(Array):
 
@@ -2865,7 +2889,7 @@ def evalf(self, evalargs):
       raise ValueError('argument {!r} missing'.format(self._name))
     else:
       value = numpy.asarray(value)
-      assert value.shape == self.shape, 'invalid argument shape: got {}, expected {}'.format(value.shape, self.shape)
+      assert equalshape(value.shape, self.shape)
       value = value.astype(self.dtype, casting='safe', copy=False)
       return value
 
@@ -2919,9 +2943,9 @@ def __init__(self, func:asarray):
     super().__init__(args=[func], shape=func._axes[:-2]+(axisprop,), dtype=func.dtype)
 
   def _simplified(self):
-    if self.func.shape[-2] == 1:
+    if equalindex(self.func.shape[-2], 1):
       return get(self.func, -2, 0)
-    if self.func.shape[-1] == 1:
+    if equalindex(self.func.shape[-1], 1):
       return get(self.func, -1, 0)
     if isinstance(self._axes[-1], Sparse):
       return self._resparsify(self.ndim-1)
@@ -2933,12 +2957,12 @@ def evalf(self, f):
     return f.reshape(f.shape[:-2] + (f.shape[-2]*f.shape[-1],))
 
   def _multiply(self, other):
-    if isinstance(other, Ravel) and other.func.shape[-2:] == self.func.shape[-2:]:
+    if isinstance(other, Ravel) and equalshape(other.func.shape[-2:], self.func.shape[-2:]):
       return Ravel(Multiply([self.func, other.func]))
     return Ravel(Multiply([self.func, Unravel(other, *self.func.shape[-2:])]))
 
   def _add(self, other):
-    if isinstance(other, Ravel) and other.func.shape[-2:] == self.func.shape[-2:]:
+    if isinstance(other, Ravel) and equalshape(other.func.shape[-2:], self.func.shape[-2:]):
       return Ravel(Add([self.func, other.func]))
 
   def _sum(self, axis):
@@ -2961,7 +2985,7 @@ def _take(self, index, axis):
   def _unravel(self, axis, shape):
     if axis != self.ndim-1:
       return Ravel(unravel(self.func, axis, shape))
-    elif shape == self.func.shape[-2:]:
+    elif equalshape(shape, self.func.shape[-2:]):
       return self.func
 
   def _inflate(self, dofmap, length, axis):
@@ -3007,18 +3031,18 @@ class Unravel(Array):
   __slots__ = 'func'
 
   @types.apply_annotations
-  def __init__(self, func:asarray, sh1:as_canonical_length, sh2:as_canonical_length):
+  def __init__(self, func:asarray, sh1:asindex, sh2:asindex):
     if func.ndim == 0:
       raise Exception('cannot unravel scalar function')
-    if func.shape[-1] != as_canonical_length(sh1 * sh2):
+    if not equalindex(func.shape[-1], sh1 * sh2):
       raise Exception('new shape does not match axis length')
     self.func = func
-    super().__init__(args=[func, asarray(sh1), asarray(sh2)], shape=func._axes[:-1]+(sh1, sh2), dtype=func.dtype)
+    super().__init__(args=[func, sh1, sh2], shape=func._axes[:-1]+(sh1, sh2), dtype=func.dtype)
 
   def _simplified(self):
-    if self.shape[-2] == 1:
+    if equalindex(self.shape[-2], 1):
       return insertaxis(self.func, self.ndim-2, 1)
-    if self.shape[-1] == 1:
+    if equalindex(self.shape[-1], 1):
       return insertaxis(self.func, self.ndim-1, 1)
     return self.func._unravel(self.ndim-2, self.shape[-2:])
 
@@ -3054,9 +3078,7 @@ class Range(Array):
   __slots__ = 'length', 'offset'
 
   @types.apply_annotations
-  def __init__(self, length:asarray, offset:asarray=Zeros((), int)):
-    assert length.ndim == 0 and length.dtype == int
-    assert offset.ndim == 0 and offset.dtype == int
+  def __init__(self, length:asindex, offset:asindex=Zeros((), int)):
     self.length = length
     self.offset = offset
     super().__init__(args=[length, offset], shape=[length], dtype=int)
@@ -3112,9 +3134,9 @@ class Polyval(Array):
   def __init__(self, coeffs:asarray, points:asarray, ngrad:types.strictint=0):
     if points.ndim < 1:
       raise ValueError('argument `points` should have at least one axis')
-    self.points_ndim = points.shape[-1]
-    if not numeric.isint(self.points_ndim):
+    if not points.shape[-1].isconstant:
       raise ValueError('the last axis of argument `points` should be a constant integer')
+    self.points_ndim = int(points.shape[-1])
     ndim = coeffs.ndim - self.points_ndim
     if ndim < 0:
       raise ValueError('argument `coeffs` should have at least one axis per spatial dimension')
@@ -3165,7 +3187,11 @@ def _simplified(self):
 class PolyOuterProduct(Array):
 
   def __init__(self, left, right):
-    shape = (left.shape[0] * right.shape[0],) + (max(left.shape[1:], default=0) + max(right.shape[1:], default=0) - 1,) * (left.ndim + right.ndim - 2)
+    nleft = left.shape[1]
+    assert all(n == nleft for n in left.shape[2:])
+    nright = right.shape[1]
+    assert all(n == nright for n in right.shape[2:])
+    shape = (left.shape[0] * right.shape[0],) + (nleft + nright - 1,) * (left.ndim + right.ndim - 2)
     super().__init__(args=[left, right], shape=shape, dtype=float)
 
   def evalf(self, left, right):
@@ -3212,12 +3238,12 @@ class Choose(Array):
   '''Function equivalent of :func:`numpy.choose`.'''
 
   @types.apply_annotations
-  def __init__(self, index:asarray, choices:types.tuple[asarray]):
+  def __init__(self, index:asarray, choices:asarrays):
     if index.dtype != int:
       raise Exception('index must be integer valued')
     dtype = _jointdtype(*[choice.dtype for choice in choices])
     shape = index.shape
-    if not all(choice.shape == shape for choice in choices):
+    if not all(equalshape(choice.shape, shape) for choice in choices):
       raise Exception('shapes vary')
     self.index = index
     self.choices = choices
@@ -3284,7 +3310,7 @@ class LoopSum(Array):
   def __init__(self, func: asarray, index:types.strict[Argument], length: asarray):
     if index.dtype != int or index.ndim != 0:
       raise ValueError('expected an index with dtype int and dimension zero but got {}'.format(index))
-    if any(index in asarray(n).arguments for n in func.shape):
+    if any(index in n.arguments for n in func.shape):
       raise ValueError('the shape of the function must not depend on the index')
     if index in length.arguments:
       raise ValueError('the length of the loop must not depend on the index')
@@ -3293,7 +3319,7 @@ def __init__(self, func: asarray, index:types.strict[Argument], length: asarray)
     self.index = index
     self.length = length
 
-    invariants = [*map(asarray, func.shape), length]
+    invariants = [*func.shape, length]
     dependencies = []
     _populate_dependencies_sans_invariants(func, index, invariants, dependencies, set())
     indices = {d: i for i, d in enumerate(itertools.chain(invariants, [index], dependencies))}
@@ -3334,7 +3360,7 @@ def _node(self, cache, subgraph, times):
     loopgraph = Subgraph('Loop', subgraph)
     subcache[self.index] = RegularNode('LoopIndex', (), dict(length=self.length._node(cache, subgraph, times)), (type(self).__name__, _Stats()), loopgraph)
     subtimes = times.get(self, collections.defaultdict(_Stats))
-    sum_kwargs = {'shape[{}]'.format(i): asarray(n)._node(cache, subgraph, times) for i, n in enumerate(self.shape)}
+    sum_kwargs = {'shape[{}]'.format(i): n._node(cache, subgraph, times) for i, n in enumerate(self.shape)}
     sum_kwargs['func'] = self.func._node(subcache, loopgraph, subtimes)
     cache[self] = node = RegularNode('LoopSum', (), sum_kwargs, (type(self).__name__, subtimes['sum']), loopgraph)
     return node
@@ -3373,7 +3399,7 @@ def _assparse(self):
           *elem_indices, elem_values = (Transpose.to_end(arr, *variable) for arr in (*elem_indices, elem_values))
           for i in variable[:-1]:
             *elem_indices, elem_values = map(Ravel, (*elem_indices, elem_values))
-          assert all(numeric.isint(n) for n in elem_values.shape[:-1])
+          assert all(n.isconstant for n in elem_values.shape[:-1])
         chunks.append(tuple(loop_concatenate(arr, self.index, self.length) for arr in (*elem_indices, elem_values)))
     return tuple(chunks)
 
@@ -3395,12 +3421,12 @@ def _simplified(self):
 class LoopConcatenate(Array):
 
   @types.apply_annotations
-  def __init__(self, func:asarray, start:asarray, stop:asarray, cc_length:asarray, index:types.strict[Argument], length:asarray):
+  def __init__(self, func:asarray, start:asindex, stop:asindex, cc_length:asindex, index:types.strict[Argument], length:asindex):
     if index.dtype != int or index.ndim != 0:
       raise ValueError('expected an index with dtype int and dimension zero but got {}'.format(index))
     if not func.ndim:
       raise ValueError('expected an array with at least one axis')
-    if any(index in asarray(n).arguments for n in func.shape[:-1]):
+    if any(index in n.arguments for n in func.shape[:-1]):
       raise ValueError('the shape of the function must not depend on the index')
     if index in length.arguments:
       raise ValueError('the length of the loop must not depend on the index')
@@ -3439,9 +3465,9 @@ def _simplified(self):
     elif self.index not in self.func.arguments:
       return Ravel(Transpose.from_end(InsertAxis(self.func, self.length), -2))
     for iaxis, axis in enumerate(self.func._axes[:-1]):
-      if isinstance(axis, Inserted) and self.index not in asarray(axis.length).arguments:
+      if isinstance(axis, Inserted) and self.index not in axis.length.arguments:
         return insertaxis(loop_concatenate(self.func._uninsert(iaxis), self.index, self.length), iaxis, axis.length)
-    if isinstance(self.func._axes[-1], Inserted) and self.index not in asarray(self.func.shape[-1]).arguments and self.func.shape[-1] != 1:
+    if isinstance(self.func._axes[-1], Inserted) and self.index not in self.func.shape[-1].arguments and not equalindex(self.func.shape[-1], 1):
       return Ravel(InsertAxis(loop_concatenate(InsertAxis(self.func._uninsert(self.func.ndim-1), 1), self.index, self.length), self.func.shape[-1]))
 
   def _takediag(self, axis1, axis2):
@@ -3471,12 +3497,12 @@ def _loop_concatenate_deps(self):
 class LoopConcatenateCombined(Evaluable):
 
   @types.apply_annotations
-  def __init__(self, funcdata:types.tuple[types.tuple[asarray]], index:types.strict[Argument], length:asarray):
+  def __init__(self, funcdata:types.tuple[asarrays], index:types.strict[Argument], length:asindex):
     if index.dtype != int or index.ndim != 0:
       raise ValueError('expected an index with dtype int and dimension zero but got {}'.format(index))
     if any(not func.ndim for func, start, stop, cc_length in funcdata):
       raise ValueError('expected an array with at least one axis')
-    if any(index in n.arguments for func, start, stop, cc_length in funcdata for n in (*map(asarray, func.shape[:-1]), cc_length)):
+    if any(index in n.arguments for func, start, stop, cc_length in funcdata for n in (*func.shape[:-1], cc_length)):
       raise ValueError('the shape of the function must not depend on the index')
     if index in length.arguments:
       raise ValueError('the length of the loop must not depend on the index')
@@ -3487,7 +3513,7 @@ def __init__(self, funcdata:types.tuple[types.tuple[asarray]], index:types.stric
 
     invariants = []
     for func, start, stop, cc_length in funcdata:
-      invariants.extend(map(asarray, func.shape[:-1]))
+      invariants.extend(func.shape[:-1])
       invariants.append(cc_length)
     invariants.append(length)
 
@@ -3550,7 +3576,7 @@ def _node_tuple(self, cache, subgraph, times):
     subtimes = times.get(self, collections.defaultdict(_Stats))
     concats = []
     for func, start, stop, cc_length in zip(self._funcs, self._starts, self._stops, self._cc_lengths):
-      concat_kwargs = {'shape[{}]'.format(i): asarray(n)._node(cache, subgraph, times) for i, n in enumerate((*func.shape[:-1], cc_length))}
+      concat_kwargs = {'shape[{}]'.format(i): n._node(cache, subgraph, times) for i, n in enumerate((*func.shape[:-1], cc_length))}
       concat_kwargs['start'] = start._node(subcache, loopgraph, subtimes)
       concat_kwargs['stop'] = stop._node(subcache, loopgraph, subtimes)
       concat_kwargs['func'] = func._node(subcache, loopgraph, subtimes)
@@ -3580,18 +3606,36 @@ def _gathersparsechunks(chunks):
 
 def _numpy_align(*arrays):
   '''reshape arrays according to Numpy's broadcast conventions'''
+
+  # NOTE: this routine largely repeats the logic of equalshape to reduce the
+  # number of comparisons. Ultimately _numpy_align will be removed from
+  # evaluable so this situation is temporary.
+
   arrays = [asarray(array) for array in arrays]
   if len(arrays) > 1:
     ndim = max([array.ndim for array in arrays])
     for idim in range(ndim):
-      lengths = [array.shape[idim] for array in arrays if array.ndim == ndim and array.shape[idim] != 1]
-      length = lengths[0] if lengths else 1
-      assert all(l == length for l in lengths), 'incompatible shapes: {}'.format(' != '.join(str(l) for l in lengths))
-      for i, a in enumerate(arrays):
-        if a.ndim < ndim:
-          arrays[i] = insertaxis(a, idim, length)
-        elif a.shape[idim] != length:
-          arrays[i] = repeat(a, length, idim)
+      lengths = {}
+      insert = []
+      for iarray, array in enumerate(arrays):
+        if array.ndim < ndim:
+          insert.append(iarray) # mark array for axis insertion
+        else:
+          lengths.setdefault(array.shape[idim], []).append(iarray)
+      if len(lengths) > 1: # shapes differ -> simplify to canonical form
+        for length in list(lengths):
+          simple = length.simplified
+          if length is not simple:
+            lengths.setdefault(simple, []).extend(lengths.pop(length))
+      expand = () if len(lengths) == 1 \
+          else lengths.pop(Constant(1), ()) # NOTE: this uses the fact that Constant(1) is simple
+      assert len(lengths) == 1 or all(l.isconstant for l in lengths) and len(set(map(int, lengths))) == 1, \
+        'incompatible shapes: {}'.format(' != '.join(str(l) for l in lengths))
+      length, _ = lengths.popitem()
+      for i in insert:
+        arrays[i] = insertaxis(arrays[i], idim, length)
+      for i in expand:
+        arrays[i] = repeat(arrays[i], length, idim)
   return arrays
 
 def _inflate_scalar(arg, shape):
@@ -3763,7 +3807,7 @@ def stack(args, axis=0):
 
 def repeat(arg, length, axis):
   arg = asarray(arg)
-  assert arg.shape[axis] == 1
+  assert equalindex(arg.shape[axis], 1)
   return insertaxis(get(arg, axis, 0), axis, length)
 
 def get(arg, iax, item):
@@ -3783,8 +3827,8 @@ def jacobian(geom, ndims):
 
   assert geom.ndim >= 1
   J = localgradient(geom, ndims)
-  cndims = geom.shape[-1]
-  assert J.shape[-2:] == (cndims,ndims), 'wrong jacobian shape: got {}, expected {}'.format(J.shape[-2:], (cndims, ndims))
+  cndims = int(geom.shape[-1])
+  assert int(J.shape[-2]) == cndims and int(J.shape[-1]) == ndims, 'wrong jacobian shape: got {}, expected {}'.format((int(J.shape[-2]), int(J.shape[-1])), (cndims, ndims))
   assert cndims >= ndims, 'geometry dimension < topology dimension'
   detJ = abs(determinant(J)) if cndims == ndims \
     else ones(J.shape[:-2]) if ndims == 0 \
@@ -3819,7 +3863,7 @@ def derivative(func, var, seen=None):
   else:
     result = func._derivative(var, seen)
     seen[func] = result
-  assert result.shape == func.shape+var.shape, 'bug in {}._derivative'.format(type(func).__name__)
+  assert equalshape(result.shape, func.shape+var.shape), 'bug in {}._derivative'.format(type(func).__name__)
   return result
 
 def localgradient(arg, ndims):
@@ -3851,7 +3895,7 @@ def _takeslice(arg:asarray, s:types.strict[slice], axis:types.strictint):
     if start == 0 and stop == n:
       return arg
     index = Range(stop-start, start)
-  elif numeric.isint(n):
+  elif n.isconstant:
     index = Constant(numpy.arange(*s.indices(arg.shape[axis])))
   else:
     raise Exception('a non-unit slice requires a constant-length axis')
@@ -3862,22 +3906,22 @@ def take(arg:asarray, index:asarray, axis:types.strictint):
   assert index.ndim == 1
   length = arg.shape[axis]
   if index.dtype == bool:
-    assert index.shape[0] == length
+    assert equalindex(index.shape[0], length)
     index = Find(index)
   elif index.isconstant:
     index_ = index.eval()
     ineg = numpy.less(index_, 0)
-    if not numeric.isint(length):
+    if not length.isconstant:
       if ineg.any():
         raise IndexError('negative indices only allowed for constant-length axes')
     elif ineg.any():
-      if numpy.less(index_, -length).any():
-        raise IndexError('indices out of bounds: {} < {}'.format(index_, -length))
-      return _take(arg, Constant(index_ + ineg * length), axis)
-    elif numpy.greater_equal(index_, length).any():
-      raise IndexError('indices out of bounds: {} >= {}'.format(index_, length))
+      if numpy.less(index_, -int(length)).any():
+        raise IndexError('indices out of bounds: {} < {}'.format(index_, -int(length)))
+      return _take(arg, Constant(index_ + ineg * int(length)), axis)
+    elif numpy.greater_equal(index_, int(length)).any():
+      raise IndexError('indices out of bounds: {} >= {}'.format(index_, int(length)))
     elif numpy.greater(numpy.diff(index_), 0).all():
-      return mask(arg, numeric.asboolean(index_, length), axis)
+      return mask(arg, numeric.asboolean(index_, int(length)), axis)
   return _take(arg, index, axis)
 
 @types.apply_annotations
@@ -3886,9 +3930,9 @@ def _take(arg:asarray, index:asarray, axis:types.strictint):
   return Transpose.from_end(Take(Transpose.to_end(arg, axis), index), *range(axis, axis+index.ndim))
 
 @types.apply_annotations
-def _inflate(arg:asarray, dofmap:asarray, length:asarray, axis:types.strictint):
+def _inflate(arg:asarray, dofmap:asarray, length:asindex, axis:types.strictint):
   axis = numeric.normdim(arg.ndim+1-dofmap.ndim, axis)
-  assert dofmap.shape == arg.shape[axis:axis+dofmap.ndim]
+  assert equalshape(dofmap.shape, arg.shape[axis:axis+dofmap.ndim])
   return Transpose.from_end(Inflate(Transpose.to_end(arg, *range(axis, axis+dofmap.ndim)), dofmap, length), axis)
 
 def mask(arg, mask, axis=0):
@@ -3931,7 +3975,7 @@ def appendaxes(func, shape):
 
 def _loop_concatenate_data(func, index, length):
   func = asarray(func)
-  chunk_size = asarray(func.shape[-1])
+  chunk_size = func.shape[-1]
   if chunk_size.isconstant:
     chunk_sizes = InsertAxis(chunk_size, length)
   else:
@@ -3939,16 +3983,16 @@ def _loop_concatenate_data(func, index, length):
   offsets = _SizesToOffsets(chunk_sizes)
   start = _take(offsets, index, 0)
   stop = _take(offsets, index+1, 0)
-  cc_length = asarray(as_canonical_length(_take(offsets, length, 0)))
+  cc_length = _take(offsets, length, 0)
   return func, start, stop, cc_length
 
 def loop_concatenate(func, index, length):
-  length = asarray(length)
+  length = asindex(length)
   func, start, stop, cc_length = _loop_concatenate_data(func, index, length)
   return LoopConcatenate(func, start, stop, cc_length, index, length)
 
 def loop_concatenate_combined(funcs, index, length):
-  length = asarray(length)
+  length = asindex(length)
   unique_funcs = []
   unique_funcs.extend(func for func in funcs if func not in unique_funcs)
   unique_func_data = tuple(_loop_concatenate_data(func, index, length) for func in unique_funcs)
@@ -3977,7 +4021,7 @@ def replace_arguments(value, arguments):
   '''
   if isinstance(value, Argument) and value._name in arguments:
     v = asarray(arguments[value._name])
-    assert value.shape == v.shape
+    assert equalshape(value.shape, v.shape)
     return v
 
 if __name__ == '__main__':
diff --git a/nutils/function.py b/nutils/function.py
index f0263add6..bcf1a457d 100644
--- a/nutils/function.py
+++ b/nutils/function.py
@@ -52,9 +52,12 @@ def lower(self, *, transform_chains: Tuple[evaluable.TransformChain] = (), coord
   def _lower(self, **kwargs):
     result = self._ArrayMeta__lower(**kwargs)
     assert isinstance(result, evaluable.Array)
-    offset = kwargs['coordinates'][0].ndim-1 if kwargs.get('coordinates', ()) else 0
+    kwargs_nocoords = kwargs.copy()
+    coordinates = kwargs_nocoords.pop('coordinates', ())
+    offset = kwargs['coordinates'][0].ndim-1 if coordinates and not type(self).__name__ == '_WithoutPoints' else 0
     assert result.ndim == self.ndim + offset
-    assert all(m == n for m, n in zip(result.shape[offset:], self.shape) if isinstance(n, int)), 'shape mismatch'
+    self_shape = tuple(asarray(n).lower(**kwargs_nocoords) for n in self.shape)
+    assert evaluable.equalshape(result.shape[offset:], self_shape) != False, 'shape mismatch'
     return result
 
   class _ArrayMeta(_ArrayMeta):
@@ -400,9 +403,9 @@ def derivative(self, __var: Union[str, 'Argument']) -> 'Array':
     if isinstance(__var, str):
       for arg in self.as_evaluable_array().arguments:
         if isinstance(arg, evaluable.Argument) and arg._name == __var:
-          if not all(isinstance(n, int) for n in arg.shape):
+          if not all(n.isconstant for n in arg.shape):
             raise ValueError('arguments with variable shapes are not supported')
-          __var = Argument(__var, arg.shape, dtype=arg.dtype)
+          __var = Argument(__var, tuple(map(int, arg.shape)), dtype=arg.dtype)
           break
       else:
         raise ValueError('no such argument: {}'.format(__var))
@@ -426,10 +429,10 @@ def argshapes(self) -> Mapping[str, Tuple[int, ...]]:
         if arg._name in shapes:
           if shapes.get(arg._name, arg.shape) != arg.shape:
             raise Exception('non-matching arguments shapes encountered')
-        elif not all(isinstance(n, int) for n in arg.shape):
+        elif not all(n.isconstant for n in arg.shape):
           raise ValueError('arguments with variable shapes are not supported')
         else:
-          shapes[arg._name] = arg.shape
+          shapes[arg._name] = tuple(map(int, arg.shape))
     return shapes
 
 def _prepend_points(__arg: evaluable.Array, *, coordinates: Sequence[evaluable.Array] = (), **kwargs: Any) -> evaluable.Array:
@@ -633,7 +636,7 @@ def __init__(self, geom: Array) -> None:
 
   def lower(self, *, coordinates: Tuple[evaluable.Array] = (), **kwargs: Any) -> evaluable.Array:
     assert coordinates
-    ndims = coordinates[0].shape[-1]
+    ndims = int(coordinates[0].shape[-1])
     return evaluable.jacobian(self._geom.lower(coordinates=coordinates, **kwargs), ndims)
 
 class _Elemwise(Array):
@@ -2607,7 +2610,7 @@ def __init__(self, ndofs: int, nelems: int, index: Array, coords: Array) -> None
     self._arg_dofs, self._arg_coeffs = [f.optimized_for_numpy for f in self.f_dofs_coeffs(_index)]
     assert self._arg_dofs.ndim == 1
     assert self._arg_coeffs.ndim == 1 + coords.shape[0]
-    assert self._arg_dofs.shape[0] == self._arg_coeffs.shape[0]
+    assert evaluable.equalindex(self._arg_dofs.shape[0], self._arg_coeffs.shape[0])
     self._arg_ndofs = evaluable.asarray(self._arg_dofs.shape[0])
 
   def lower(self, **kwargs: Any) -> evaluable.Array:
diff --git a/nutils/solver.py b/nutils/solver.py
index 315501b75..fa7de4684 100644
--- a/nutils/solver.py
+++ b/nutils/solver.py
@@ -875,7 +875,7 @@ def _parse_lhs_cons(lhs0, constrain, targets, argobjs, arguments):
   for target in targets:
     if target not in argobjs:
       raise SolverError('target does not occur in functional: {!r}'.format(target))
-    shape = argobjs[target].shape
+    shape = tuple(map(int, argobjs[target].shape))
     if target not in arguments:
       arguments[target] = numpy.zeros(shape)
     elif arguments[target].shape != shape:
diff --git a/nutils/topology.py b/nutils/topology.py
index 7a3d66e3d..6264131e1 100644
--- a/nutils/topology.py
+++ b/nutils/topology.py
@@ -197,7 +197,7 @@ def integrate_elementwise(self, funcs, *, degree, asfunction=False, ischeme='gau
     retvals = [sparse.toarray(retval) for retval in self.sample(ischeme, degree).integrate_sparse(
       [function.kronecker(func, pos=self.f_index, length=len(self), axis=0) for func in funcs], arguments=arguments)]
     if asfunction:
-      return [function.Elemwise(retval, self.f_index, dtype=float) for retval in retvals]
+      return [function.get(retval, 0, self.f_index) for retval in retvals]
     else:
       return retvals
 
diff --git a/tests/test_evaluable.py b/tests/test_evaluable.py
index 0bab45695..db1e20b12 100644
--- a/tests/test_evaluable.py
+++ b/tests/test_evaluable.py
@@ -19,6 +19,9 @@ def setUp(self):
     self.pairs = [(i, j) for i in range(self.op_args.ndim-1) for j in range(i+1, self.op_args.ndim) if self.op_args.shape[i] == self.op_args.shape[j]]
     _builtin_warnings.simplefilter('ignore', evaluable.ExpensiveEvaluationWarning)
 
+  def assertShapes(self):
+    self.assertEqual(self.n_op_argsfun.shape, tuple(map(int, self.op_args.shape)))
+
   def assertArrayAlmostEqual(self, actual, desired, decimal):
     if actual.shape != desired.shape:
       self.fail('shapes of actual {} and desired {} are incompatible.'.format(actual.shape, desired.shape))
@@ -83,7 +86,7 @@ def test_eval_withtimes(self):
 
   def test_getitem(self):
     for idim in range(self.op_args.ndim):
-      for item in range(self.op_args.shape[idim]):
+      for item in range(self.n_op_argsfun.shape[idim]):
         s = (Ellipsis,) + (slice(None),)*idim + (item,) + (slice(None),)*(self.op_args.ndim-idim-1)
         self.assertFunctionAlmostEqual(decimal=15,
           desired=self.n_op_argsfun[s],
@@ -132,14 +135,14 @@ def test_determinant(self):
 
   def test_take(self):
     indices = [0,-1]
-    for iax, sh in enumerate(self.op_args.shape):
+    for iax, sh in enumerate(self.n_op_argsfun.shape):
       if sh >= 2:
         self.assertFunctionAlmostEqual(decimal=15,
           desired=numpy.take(self.n_op_argsfun, indices, axis=iax),
           actual=evaluable.take(self.op_args, indices, axis=iax))
 
   def test_take_block(self):
-    for iax, sh in enumerate(self.op_args.shape):
+    for iax, sh in enumerate(self.n_op_argsfun.shape):
       if sh >= 2:
         indices = [[0,sh-1],[sh-1,0]]
         self.assertFunctionAlmostEqual(decimal=15,
@@ -147,7 +150,7 @@ def test_take_block(self):
           actual=evaluable._take(self.op_args, indices, axis=iax))
 
   def test_take_nomask(self):
-    for iax, sh in enumerate(self.op_args.shape):
+    for iax, sh in enumerate(self.n_op_argsfun.shape):
       if sh >= 2:
         indices = [0,sh-1]
         self.assertFunctionAlmostEqual(decimal=15,
@@ -156,14 +159,14 @@ def test_take_nomask(self):
 
   def test_take_reversed(self):
     indices = [-1,0]
-    for iax, sh in enumerate(self.op_args.shape):
+    for iax, sh in enumerate(self.n_op_argsfun.shape):
       if sh >= 2:
         self.assertFunctionAlmostEqual(decimal=15,
           desired=numpy.take(self.n_op_argsfun, indices, axis=iax),
           actual=evaluable.take(self.op_args, indices, axis=iax))
 
   def test_take_duplicate_indices(self):
-    for iax, sh in enumerate(self.op_args.shape):
+    for iax, sh in enumerate(self.n_op_argsfun.shape):
       if sh >= 2:
         indices = [0,sh-1,0,0]
         self.assertFunctionAlmostEqual(decimal=15,
@@ -171,16 +174,16 @@ def test_take_duplicate_indices(self):
           actual=evaluable.take(self.op_args, evaluable.Guard(evaluable.asarray(indices)), axis=iax))
 
   def test_inflate(self):
-    for iax, sh in enumerate(self.op_args.shape):
-      dofmap = evaluable.Constant(numpy.arange(sh) * 2)
-      desired = numpy.zeros(self.n_op_argsfun.shape[:iax] + (sh*2-1,) + self.n_op_argsfun.shape[iax+1:], dtype=self.n_op_argsfun.dtype)
+    for iax, sh in enumerate(self.n_op_argsfun.shape):
+      dofmap = evaluable.Constant(numpy.arange(int(sh)) * 2)
+      desired = numpy.zeros(self.n_op_argsfun.shape[:iax] + (int(sh)*2-1,) + self.n_op_argsfun.shape[iax+1:], dtype=self.n_op_argsfun.dtype)
       desired[(slice(None),)*iax+(slice(None,None,2),)] = self.n_op_argsfun
       self.assertFunctionAlmostEqual(decimal=15,
         desired=desired,
         actual=evaluable._inflate(self.op_args, dofmap=dofmap, length=sh*2-1, axis=iax))
 
   def test_inflate_duplicate_indices(self):
-    for iax, sh in enumerate(self.op_args.shape):
+    for iax, sh in enumerate(self.n_op_argsfun.shape):
       dofmap = numpy.arange(sh) % 2
       desired = numpy.zeros(self.n_op_argsfun.shape[:iax] + (2,) + self.n_op_argsfun.shape[iax+1:], dtype=self.n_op_argsfun.dtype)
       numpy.add.at(desired, (slice(None),)*iax+(dofmap,), self.n_op_argsfun)
@@ -203,9 +206,9 @@ def test_product(self):
 
   def test_getslice(self):
     for idim in range(self.op_args.ndim):
-      if self.op_args.shape[idim] == 1:
+      if self.n_op_argsfun.shape[idim] == 1:
         continue
-      s = (Ellipsis,) + (slice(None),)*idim + (slice(0,self.op_args.shape[idim]-1),) + (slice(None),)*(self.op_args.ndim-idim-1)
+      s = (Ellipsis,) + (slice(None),)*idim + (slice(0,int(self.n_op_argsfun.shape[idim])-1),) + (slice(None),)*(self.op_args.ndim-idim-1)
       self.assertFunctionAlmostEqual(decimal=15,
         desired=self.n_op_argsfun[s],
         actual=self.op_args[s])
@@ -243,7 +246,7 @@ def test_power(self):
       actual=(self.op_args**3))
 
   def test_power0(self):
-    power = (numpy.arange(self.op_args.size) % 2).reshape(self.op_args.shape)
+    power = (numpy.arange(self.n_op_argsfun.size) % 2).reshape(self.n_op_argsfun.shape)
     self.assertFunctionAlmostEqual(decimal=13,
       desired=self.n_op_argsfun**power,
       actual=self.op_args**power)
@@ -256,11 +259,11 @@ def test_sign(self):
 
   def test_mask(self):
     for idim in range(self.op_args.ndim):
-      if self.op_args.shape[idim] <= 1:
+      if self.n_op_argsfun.shape[idim] <= 1:
         continue
-      mask = numpy.ones(self.op_args.shape[idim], dtype=bool)
+      mask = numpy.ones(self.n_op_argsfun.shape[idim], dtype=bool)
       mask[0] = False
-      if self.op_args.shape[idim] > 2:
+      if self.n_op_argsfun.shape[idim] > 2:
         mask[-1] = False
       self.assertFunctionAlmostEqual(decimal=15,
         desired=self.n_op_argsfun[(slice(None,),)*idim+(mask,)],
@@ -306,7 +309,7 @@ def test_desparsify(self):
     args = []
     for arg in self.args:
       for i in range(arg.ndim):
-        arg = evaluable._inflate(arg, evaluable.Guard(numpy.arange(arg.shape[i])), arg.shape[i], i)
+        arg = evaluable._inflate(arg, evaluable.Guard(numpy.arange(int(arg.shape[i]))), arg.shape[i], i)
       args.append(arg)
     op_args = self.op(*args).simplified
     evalargs = dict(zip(self.arg_names, self.arg_values))
@@ -444,7 +447,7 @@ def _check(name, op, n_op, shapes, hasgrad=True, zerograd=False, ndim=2, low=-1,
 _check('loopconcatenate2', lambda: evaluable.loop_concatenate(evaluable.Elemwise([numpy.arange(48).reshape(4,4,3)[:,:,a:b] for a, b in util.pairwise([0,2,3])], evaluable.Argument('index', (), int), int), evaluable.Argument('index', (), int), 2), lambda: numpy.arange(48).reshape(4,4,3), [])
 _check('loopconcatenatecombined', lambda a: evaluable.loop_concatenate_combined([a+evaluable.prependaxes(evaluable.Argument('index', (), int), a.shape)], evaluable.Argument('index', (), int), 3)[0], lambda a: a+numpy.arange(3)[None], [(2,1)], hasgrad=False)
 
-_polyval_mask = lambda shape, ndim: 1 if ndim == 0 else numpy.array([sum(i[-ndim:]) < shape[-1] for i in numpy.ndindex(shape)], dtype=int).reshape(shape)
+_polyval_mask = lambda shape, ndim: 1 if ndim == 0 else numpy.array([sum(i[-ndim:]) < int(shape[-1]) for i in numpy.ndindex(shape)], dtype=int).reshape(shape)
 _polyval_desired = lambda c, x: sum(c[(...,*i)]*(x[(slice(None),*[None]*(c.ndim-x.shape[1]))]**i).prod(-1) for i in itertools.product(*[range(c.shape[-1])]*x.shape[1]) if sum(i) < c.shape[-1])
 _check('polyval_1d_p0', lambda c, x: evaluable.Polyval(c*_polyval_mask(c.shape,x.shape[1]), x), _polyval_desired, [(1,),(4,1)], ndim=1)
 _check('polyval_1d_p1', lambda c, x: evaluable.Polyval(c*_polyval_mask(c.shape,x.shape[1]), x), _polyval_desired, [(2,),(4,1)], ndim=1)
@@ -594,7 +597,7 @@ def test_derivative(self):
     self.assertTrue(evaluable.iszero(evaluable.localgradient(self.func, self.domain.ndims).simplified))
 
   def test_shape_derivative(self):
-    self.assertEqual(evaluable.localgradient(self.func, self.domain.ndims).shape, self.func.shape+(self.domain.ndims,))
+    self.assertEqual(evaluable.localgradient(self.func, self.domain.ndims).shape, self.func.shape+(evaluable.Constant(self.domain.ndims),))
 
 class jacobian(TestCase):
 
@@ -723,21 +726,23 @@ def test_loop_concatenate(self):
                      '└ B = Loop\n'
                      'NODES\n'
                      '%B0 = LoopConcatenate\n'
-                     '├ shape[0] = 2\n'
-                     '├ start = %B1 = Take; i:\n'
+                     '├ shape[0] = %A1 = Take; i:\n'
                      '│ ├ %A2 = _SizesToOffsets; i:a3\n'
                      '│ │ └ %A3 = InsertAxis; i:i2\n'
                      '│ │   ├ 1\n'
                      '│ │   └ 2\n'
-                     '│ └ %B4 = LoopIndex\n'
+                     '│ └ 2\n'
+                     '├ start = %B4 = Take; i:\n'
+                     '│ ├ %A2\n'
+                     '│ └ %B5 = LoopIndex\n'
                      '│   └ length = 2\n'
-                     '├ stop = %B5 = Take; i:\n'
+                     '├ stop = %B6 = Take; i:\n'
                      '│ ├ %A2\n'
-                     '│ └ %B6 = Add; i:\n'
-                     '│   ├ %B4\n'
+                     '│ └ %B7 = Add; i:\n'
+                     '│   ├ %B5\n'
                      '│   └ 1\n'
-                     '└ func = %B7 = InsertAxis; i:i1\n'
-                     '  ├ %B4\n'
+                     '└ func = %B8 = InsertAxis; i:i1\n'
+                     '  ├ %B5\n'
                      '  └ 1\n')
 
   @unittest.skipIf(sys.version_info < (3, 6), 'test requires dicts maintaining insertion order')
@@ -750,21 +755,23 @@ def test_loop_concatenatecombined(self):
                      '└ B = Loop\n'
                      'NODES\n'
                      '%B0 = LoopConcatenate\n'
-                     '├ shape[0] = 2\n'
-                     '├ start = %B1 = Take; i:\n'
+                     '├ shape[0] = %A1 = Take; i:\n'
                      '│ ├ %A2 = _SizesToOffsets; i:a3\n'
                      '│ │ └ %A3 = InsertAxis; i:i2\n'
                      '│ │   ├ 1\n'
                      '│ │   └ 2\n'
-                     '│ └ %B4 = LoopIndex\n'
+                     '│ └ 2\n'
+                     '├ start = %B4 = Take; i:\n'
+                     '│ ├ %A2\n'
+                     '│ └ %B5 = LoopIndex\n'
                      '│   └ length = 2\n'
-                     '├ stop = %B5 = Take; i:\n'
+                     '├ stop = %B6 = Take; i:\n'
                      '│ ├ %A2\n'
-                     '│ └ %B6 = Add; i:\n'
-                     '│   ├ %B4\n'
+                     '│ └ %B7 = Add; i:\n'
+                     '│   ├ %B5\n'
                      '│   └ 1\n'
-                     '└ func = %B7 = InsertAxis; i:i1\n'
-                     '  ├ %B4\n'
+                     '└ func = %B8 = InsertAxis; i:i1\n'
+                     '  ├ %B5\n'
                      '  └ 1\n')
 
 class simplify(TestCase):