wence/better type inference

[wence-]

Assembly of 0-forms in complex mode needs to correctly infer the type of the return variable. This is necessary so that if we're adjointing, a real-valued functional is actually returned as a real.

This is kind of a complicated patch. One other way of doing this would be to augment UFL's form preprocessing to attach a flag about the realness of the form and we then cast to real in firedrake assembly.

I think the reification of terms in the gem DAG is still good, but maybe the impero_c changes are unnecessary and the UFL-approach is TRTTD.

[miklos1]

Considering only 0-forms, it may (or may not) be simpler to do this in Firedrake assembly based on UFL type inference, but in general, I think, TSFC is the right place for this kind of inference, since TSFC has the most information naturally available: UFL knows very little about the basis functions, FIAT/FInAT knows only the basis functions and nothing about the form, TSFC has all in one place. Additionally, doing inference after optimization may naturally yield better results with no extra effort, whereas doing type inference earlier (e.g. on UFL) may give more conservative estimates in some cases, or needs more effort and special casing to retain precision.

TSFC is also in the right position to infer other things that could potentially affect the kernel interface: not just whether the result is complex-valued or not, but also for example whether the resulting matrix is diagonal or not (e.g., underintegration).

[miklos1]

With respect to complex mode and type inference (e.g. real, complex, or logical), I think the right way to go forward is a kind of typed GEM. That is, all GEM nodes should be able to tell their types, whether they are real, complex, or Boolean valued. The basis could be the existing tsfc.loopy.assign_dtypes, but I suggest caching the result on the GEM nodes so that no-brainer optimizations could casually make use of type. In other words, type could be another slot on the node with a cached_property (similar to how hash values are calculated on demand at first, then stored in a slot).

• Real, Imag, and Conj could be added as GEM nodes instead of being MathFunctions with special name strings. Real.__new__(e) and Imag.__new__(e) could immediately return e if e.type == REAL, similarly Imag.__new__(e) could return zero in this case.
• Leaf nodes should all have known types: for literals (including zero, identity matrix, Kroenecker delta) this is easy, one can just look at the numbers. Variables should get a type attribute.
• The only intermediate nodes which remain ambiguous are Sqrt and Power. They could theoretically give complex results for real arguments, yet complex values are not desired when running in real mode. These nodes should get a real attribute which, if true, signifies that the corresponding operations should rather fail than return a complex number. (Alternatively, one could make separate real and complex nodes.)

In effect, complex mode would be the internally presumed default and real mode the special case. Furthermore, the current profusion of if complex_mode could be exiled to the peripheries: translation of coefficient vectors to real or complex valued Variables, and the translation of sqrt and power in ufl2gem according to complex mode.

These complex mode branches could also be eliminated from tsfc.fem as long as they can be eliminated from tsfc.ufl_utils.preprocess_expression, which I think is fine, since ufl.algorithms.remove_complex_nodes should become redundant with typed GEM as described above.

[wence-]

That seems like a reasonable approach. One could plausibly do type assignation on construction at O(1) cost by always ensuring that nodes are typed (rather than the on-demand approach of hashing).

[miklos1]

If you look at gem.node.Node, you can see that a hash_value slot is defined, but it is only assigned when __hash__ is first called: it will do the potentially expensive calculation at the first call, save the result in hash_value, and repeated calls to hash just look up hash_value.

I think a similar approach could work for GEM node types.

[miklos1]

I know this is just a draft pull request, but this line would create a circularity in package dependency, by making gem depend on tsfc. This in undesirable: gem should not depend on anything other than the Python standard library and quasi-standard packages.

[wence-]

Yeah, that's easy to fix by doing the more general gem typing approach.

[sv2518]

Do we not need to merge this before complex?

[wence-]

No, we worked around it. To do this more generally I think miklos' comment is the right way to handle this more cleanly.