Separating symbols and data

[mlange05]

The current design of Devito function symbols is centered around the SymbolicData objects that always have an associated .data property. The original rationale for this is that we want control over data allocation, so that we may implement data layout changes (padding, vector folding, etc.). On the other hand, various use cases require us to be able to assign different data arrays to a particular symbol on operator invocation (Operator/StencilKernel.apply()). As an example we might want to do:

shape  = (...)
m0 = DenseData(name='m0', shape=shape)
m1 = DenseData(name='m1', shape=shape)

eqn = some_equation_involving_m
op = Operator(eqn, subs=...)

op.apply(m=m0)
op.apply(m=m1)

In the current and previous seismic examples (gradient_test.py) this problem is avoided by creating the m symbol during the operator setup and copying the data into m.data, without any operator re-use. In order to avoid the potentially costly copy operation this implies, we need to adjust our API to allow multiple data objects to be associated with symbols at runtime.

[mlange05]

Some possible suggestions (although this is not an exhaustive list):

a) We could turn DenseData/TimeData into "pure symobl" objects that still carry data layout and shape information (required for the compiler to adhere to things like folding or padding). But instead of having a single .data property these would provide an .allocate() function that returns TensorData objects to encapsulate actual data. This would put the job of managing and associating data into the users hands, but allow us to do enough sanity checking that the data provided at .apply() agrees with all assumptions made by the operator during codegen.

b) Somehow over-write our function symbol caching so that multiple explicit invocations of DenseData(name='m', ...) don't alias. This would allow us to natively inject different m objects into the operator with op.apply(m=m0) and op.apply(m=m1). I'm not sure how to trick the caching to allow that safely though.

c) Another variation that somewhat combines a) and b), where we allow aliasing at operator invocation time, ie. we allow:

m0 = DenseData(name='m0')
m1 = DenseData(name='m1')
...
op.apply(m0=m0)
op.apply(m0=m1)

Food for thought, other suggestions or opinions are very welcome.

[navjotk]

I agree with the general rationale behind this issue and this is indeed a problem.

However, I would suggest that c) above can be done without a), i.e. the operator could be created with symbol m, and then at run-time, it is a matter of our apply method checking whether m1 is compatible with the assumptions made about m at code-generation time and then just passing a pointer to m1.data instead of m.data while calling the C function. Syntactically, this would be equivalent to the code snippet shown in c), however, DenseData/TimeData can still maintain their data properties. This is my favourite possibility because it improves generality without completely changing everything about the API.

About b), I am not a fan of having things like m=DenseData(name='u'), i.e. not using the same symbol name as the python variable name. I think this might create hard-to-debug issues at some point in the future. In fact, in my ideal API, the name parameter would be dropped and the symbol's internal name decided by reflection or something.

If we want to go for something like a), I would suggest going further and using pure symbols (no shape or layout information). This is because, although I have said in presentations that the equation doesn't need to change for a change from 2D to 3D models, I realise now that the equation is made of symbols that have shape information already attached to them. To truly reuse a symbolic equation between 2D and 3D, inference of the shape needs to be deferred at least. For handling potential issues pointed out in a) above, there is the possibility of generating multiple versions of the code and deciding the correct one to call at run-time (just like a "real" compiler). This way the symbol can indeed be a pure symbol without shape or layout information.

Another related issue is that the order of discretisation is "hard-coded" at the creation time of the symbolic data object. Therefore, the application developer is required to create new symbols and rewrite the equation every time they change the discretisation. If the symbol can be created without choosing the discretisation, u.dx for example could behave as a symbol and not as a function call like it is now. This way, one can create multiple Operators with the same symbol, yet different discretisations. The FD derivative will only be expanded at code generation time this way so the discretisation information can be provided to the Operator and not the symbol.

[mlange05]

@navjotk Thanks, that is indeed a quite insightful analysis. I agree with you that c) looks the most pain-free and general avenue, since a) will either create a property-checking mess or prevent us from doing data layout transformations, and b) looks dangerous because we interfere with SymPy caching.

Your final point though (although slightly off-topic) I would slightly disagree with: 1) Derivative order is already parameterisable, so there is no need to "re-write" an equation. 2) The fact that we store a single derivative order per symbol is sufficient for current purposes, but can easily be over-written with things like u.deriv(x, order=4). This would allow us, in theory at least, to even combine multiple stencil order in one equation, eg.: Eq(u.deriv(x, order=4) + some_other_term + u.deriv(x, order=8)). In this case only the symbolic layer needs a small adjustment (adding .deriv()), but everything else stands. 3) What you seem to be arguing for sounds like a full-blown static form compiler (much like UFL), where the interpretation of a high-level expression is delayed and performed in one swoop at the end. This contradicts our current dynamic approach to utilising SymPy slightly, and the original premise of this project was to not do this until we have a relevant use case that forces us to do so. Do you have any concrete examples that require delayed analysis?

[navjotk]

Just like I don't understand how delaying the evaluation of derivatives requires a full form compiler here. Completely orthogonal.

[mlange05]

This is now addressed in #223, so closing.