Feedback: generalized way of symbolic integration of arbitrary pytensor graphs?

[kratsg]

Using a combination of scan and sympy's integral functionality, I wanted to understand what limitations there might be with the following approach I have in mind.

I am thinking of the concept of a function that takes as input a pytensor TensorVariable that would:

• stringify the graph
• use sympy to parse the expression
• then use sympy to perform the integral
• convert that back to a pytensor graph

The goal here is to provide a way to arbitrarily normalize pytensor computational graphs based on the domains of observables being provided at the time of evaluation (not necessarily known ahead of time). What I'm not sure about, just from the exploring I've done, is if this is a solved problem, if it's overkill, or if I should numerically integrate (e.g. gauss-legendre approach + scan) rather than trying to symbolically do it ... or even instead fallback to a numerical integratino with gauss-legendre if the symbolic approach fails.

[ricardoV94]

You may want to poke @jessegrabowski

[jessegrabowski]

@kratsg you might be able to try sympytensor, it allows you to start from a sympy expression and get out a pytensor graph. Obviously that's only the 2nd half of your plan, but depending on the broader context, it might be enough?

[kratsg]

Ahh, I solved that side of it here: https://github.com/scipp-atlas/pyhs3/blob/main/src/pyhs3/generic_parse.py#L87-L90 . I need the other direction essentially.

[kratsg]

I should probably switch to sympytensor if it plans to be well-supported and covers things I haven't ran into yet.

[kratsg]

scipp-atlas/pyhs3#182

[jessegrabowski]

I should probably switch to sympytensor if it plans to be well-supported and covers things I haven't ran into yet.

It's an open source project, so in terms of well-supported it is what it is. It should alreayd have pretty broad coverage of sympy though, I use it in several projects.

[maresb]

Do you have any concrete examples in mind @kratsg? In my experience, usually anything beyond the most trivial examples is not symbolically integrable.

[kratsg]

Everything will generally be fundamental distributions, such as gaussians, poissons, laundau, lognormal, etc. We also have asymmetric crystal ball functions which are piecewise. The structure of the code I was imagining was to provide a way known integral if possible for a given dimension of that distribution. For example: Gauss(x, mu, sigma) has three dimensions: x, mu, sigma. If one asks for the 1-dimensional integral over 'x', that would be known and can skip manually symbolically integrating it by providing a shortcircuit (somewhat similar for mu, sigma). It gets more complicated if one needs to do a two-dimensional integral over (x, mu) or (mu, sigma), or (x, sigma) -- as we normalize our distributions over the bounds of the variable(s) in question.

Evaluating a likelihood for a Gaussian whose observable 'x' is bounded over x \in [-10, 20] is normalized differently than x \in [40, 50] . So we need to be able to normalize differently depending on the domain of x in this case. Does this help...?

[jessegrabowski]

One other comment: the pytensor graph is already close to (is?) an AST representation of an expression, so you probably don't want to stringify it. Instead you can define how different pytensor ops convert to sympy, then just walk the tree applying the conversions. I'm interested in something like this (it's in scope for sympytensor), but I haven't thought carefully about how to handle dimensions. Sympy is fundamentally about scalar operations (plus some very limited support for linear algebra), while pytensor is fundamentally about tensors. So how should y = pt.expand_dims(x, 0) be parsed? We could just reject it of course. Or y = x.sum(axis=0). You can do stuff with sp.IndexBase and sp.Sum, but it gets ugly really fast.