Recognise repeated and joined contractions when sum factorising

[pbrubeck]

TLDR

Sum factorisation now plans the product tree instead of searching every
order of the contraction indices. It also keeps a contraction whole when the
expression uses it more than once, so the kernel evaluates it once.

Forms that used to fail with Too many indices for sum factorisation! now
compile. Flops fall on the cases that were doing redundant work, and stay equal
everywhere else.

Fixes #283.

What this does

Plan the product tree.
A search over index orders must multiply every factor that carries an index
before it can sum that index away. It can never sum over part of a product and
multiply the rest in later.

This PR builds the tree by dynamic programming over subsets of the factors. It
sums each index at the smallest subtree that holds every factor carrying that
index. The search costs 3^factors instead of indices!. A search over orders
stays as a fallback for very wide products.

Keep a repeated contraction whole.
Breaking a contraction up renames its indices. A coefficient evaluation used
three times therefore becomes three separate contractions, and the kernel
computes it three times. This PR finds such a contraction in the expression and
leaves it alone.

Tidy-up. One helper now holds the rule for what a product-tree walk descends
into. Two walks used to carry their own copy of it.

Interpolation on a hexahedron

Kernel flops and temporaries, as count / entries / largest.

interpolate	flops main	flops PR	temporaries main	temporaries PR	compile main	compile PR
f*f*f, CG1 -> CG1	96	96	9 / 11 / 2	9 / 13 / 2	0.011 s	0.011 s
f*f*f*f, CG1 -> CG1	104	104	9 / 11 / 2	9 / 13 / 2	0.011 s	0.011 s
f*f, CG4 -> DG3	2568	2568	4 / 51 / 25	4 / 66 / 25	0.008 s	0.009 s
dot(u, u)*u, vector CG4	23875	23875	14 / 210 / 75	14 / 270 / 75	0.017 s	0.018 s
dot(A, A), tensor CG2 -> CG3	16776	16776	81 / 240 / 4	81 / 240 / 4	0.089 s	0.163 s
inner(A, A), tensor CG2 -> CG2	4131	4131	120 / 126 / 3	120 / 180 / 3	0.126 s	0.129 s
grad(f)[0], CG4 -> CG4	14106	13981	92 / 292 / 25	91 / 375 / 25	0.066 s	0.066 s
div(u), vector CG4 -> DG3	23548	23100	110 / 391 / 25	109 / 543 / 25	0.080 s	0.077 s

dot(A, A) costs 0.078 s more to compile. Its value index joins the evaluations
into a contraction of seven indices. It is the only case here wide enough for
the planner to do real work.

Helmholtz on an extruded hexahedral mesh

inner(u, v)*dx + inner(d(u), d(v))*dx, with d = grad for CG and curl for
NCE. These forms need no repair, so this table shows what the planner costs.

	flops main	flops PR	entries main	entries PR	compile main	compile PR
CG1	14732	14732	1153	1171	0.120 s	0.126 s
CG3	611668	611668	4460	4502	0.105 s	0.103 s
CG7	37651332	37651332	90392	90482	0.101 s	0.104 s
NCE1	136618	136651	7151	7178	1.028 s	1.003 s
NCE3	4486030	4486063	22058	22109	0.754 s	0.758 s
NCE7	282933862	282933895	612438	612537	0.752 s	0.755 s

Compile time moves less than run-to-run noise. The NCE forms cost 33 more flops
at every degree. The planner finds cheaper trees for three parts of those forms,
70 operations against 90 each. It costs each connected part on its own, so it
does not see that the parts share terms, and it drops one shared subexpression.
Costing the whole expression at once is future work.

Tests

test_sum_factorise.py covers the fallback bound, a contraction that keeps a
repeated evaluation whole, and a contraction joined by a value index into more
indices than an order search can take. Each test was checked to fail when its
change is reverted.

AI assistance

Claude Code was used for implementation, benchmarking, and drafting this
section. The human contributor remains responsible for understanding,
validating, and maintaining the changes.

[pbrubeck]

This is the algorithmic realisation of splitting sum-factorization into subsets of factors. This models the splitting in dual_evaluation.

Before, we were tackling all indices at once and searching over all possible permutations = indices!.

Now we traverse a tree, by appending factors one by one. The cost is reduced to 3^factors.

[pbrubeck]

Suggested change

evaluation = gem.optimise.contraction(evaluation)

[pbrubeck]

Suggested change

evaluation = gem.optimise.contraction(evaluation)