Document GEM's shape and index behaviour

[ReubenHill]

See the below from @miklos1 in firedrakeproject/tsfc#237 . This needs to be in a design document to avoid the confusion that led to that PR seeming necessary.

The problem is that you get intermediate values that don't have the
expected shape, and you then have to somehow implicitly know what has
disappeared. In this case, compile_ufl is being used to evaluate an
expression at a bunch of points. The expression might also have shape
because it's vector or tensor-valued. If some of the shape/indices
disappear because of this sort of index cancellation, how do I safely
know which indices have disappeared? What if it's a 3D vector valued
expression being evaluated at 3 points?

I think we shall consider whether shapes (say, A being a 3 x 3 matrix)
and free indices (say, e being a scalar expression with free indices
i and j which both have extent 3) are equivalent and interchangeable,
or somehow different. Considering an expression in isolation, these two
ways of looking at tensor nature seem interchangeable, for they are
convertible: you can wrap up free indices into a tensor shape with
ComponentTensor and you can turn tensor shapes into scalars with free
indices by indexing (Indexed) them with ‘free’ (i.e. loop) indices.

This consideration, however, raises the question, why do we have both if
they are both expressing the same thing? Is it not redundant to have
both? Is it just an artefact of coming from UFL, which also has both?
There is more to it than that.

Shapes and free indices differ significantly when we look at them under
the aspect of substitutability (Ersetzbarkeit) in a surrounding
expression. Consider an expression such as Ax, uA, or uAx, where A is a
matrix, and u and x are vectors. Then, say, we are able to reduce the
rank of A (to a vector or a scalar), and we want to substitute this
“optimised” A into the original expression. This cannot be done
correctly without updating the surrounding vector operations, that is
whether they are an inner product, outer product, MatVec, MatMult,
multiplication by scalar etc.

This is not a problem, however, for free indices, since indices are
explicitly named rather than merely having an implicit order. Consider
u_i * A_{ij} * x_j. If manage to reduce A’s rank, we can simply
substitute the result back, and u_i * a_i * x_j, u_i * a_j * x_j, or u_i

• a * x_j are all correctly substituted expressions because the matching
  indices were explicitly named. It could be that the result of such a
  substitution gives reduced rank in the outer expression too, but that
  just means we’ve done even more optimisation.

This is why GEM is absolutely strict about shapes; a scalar [shape =
()], a vector of length 1 [shape = (1,)], and a 1 x 1 matrix [shape =
(1, 1)] are all distinct shapes, and they are in no way interchangeable
(unlike MATLAB, for example). Free indices, however, are handled in a
much more loose way. In fact, the reason GEM exists, and I could not
just use a subset of UFL was that in UFL, substituting something with
additional free indices (say, a coefficient with something that depends
on the quadrature index) may be illegal, because in UFL free indices are
really just a different way of looking at shapes.

So the convention I have usually used in API design was to define as
shape that which is absolutely non-negotiable, and define as free
indices those which are possible dependencies on loop indices.

If some of the shape/indices disappear because of this sort of index
cancellation, how do I safely know which indices have disappeared?

For shapes, this ought not happen, because that’s problematic. For free
indices: they are “named”, so you can just look at the expression and
tell which ones have disappeared.

What if it's a 3D vector valued expression being evaluated at 3
points?

It could be a reasonable approach to define the 3D vector value (or
whatever shape it has) as the shape of the expression, and communicate
the point dependency as free indices, so they can disappear as a
dependency if optimisation was successful.

Originally posted by @miklos1 in firedrakeproject/tsfc#237 (comment)