SymPy's `sympify` and floating point arithmetic

[lcontento]

AMICI uses SymPy for symbolic computation and mathematical formulas from SBML files are converted to strings and loaded with sympy.sympify. This function does not produce a faithful representation of the given string, but applies some transformations. Some of them are very reasonable, such as sympy.sympify('x - x') returning zero. Others are quite problematic. For example, consider the following expression for a sigmoid function:

sympy.sympify('1 / (1 + exp(800.0 - x))')

While the original formula can be evaluated numerically without any problem, sympify splits the exponential in exp(800.0) * exp(-x). SymPy has arbitrary precision arithmetic so it has no trouble with the very large number exp(800.0), but the resulting literal is interpreted as inf in floating point arithmetic (and completely rejected by libsbml MathML parser).

SymPy has an evaluate=False option for sympify to reduce the number of simplifications it applies. However, this is at the moment not propagated to functions like exp but only to the basic arithmetic operations (see SymPy issue sympy/sympy#13999). Using the evaluate context manager (marked as experimental though) is the way to go at the moment:

import sympy
from sympy.core.parameters import evaluate
with evaluate(False):
    print(sympy.sympify('1 / (1 + exp(800.0 - x))'))

Personally, I would convert all instances of sympy.sympify reading in user-written expressions to a custom one disabling evaluation. This should not be a difficult or particularly dangerous change (I think). We may lose some small opportunity for simplification, but in general I believe user-written code should be followed as accurately as possible and their intentions not second-guessed.

As a final note, I guess similar floating point arithmetic mistakes may occur when creating the code for the derivatives. I am not sure whether such errors would be possible to avoid without ditching SymPy completely, but hopefully they are not common enough to worry...

[FFroehlich]

I would argue that the anticipated result depends on context. In particular, when programmatically generating model code, we generally do want simplifications. We also want simplifications for larger models where size of generated code can become an issue. I do see the point though that simplifications are not always wanted.

We already have the infrastructure set up to perform simplification, so I don't see any issues with changing the code such that no simplification are performed by default by using with evaluate(False): as you proposed. Yet, I would keep something like sp.sympify as default for simplification. Completely avoiding simplification could be achieved by passing lambda x:x as simplification function.

How does that sound?

[JanHasenauer]

I also agree that simplifications can be very valuable. The proposed way sounds good.

As finding problems caused by simplifications can be rather difficult, not sure how to "warn" the user. Should we maybe set up a list of things which should be checked if something is going wrong?

[lcontento]

I agree that simplifications are useful, especially for generated code. But AMICI is not a model generator, so that the onus of providing nicely written formulas should be on the user (or on the software that generates the model for the user). I believe diligent users that pay attention to their formulas should not be penalized by default.

That aside, this issue concerns mainly the reading of formulas from SBML/PEtab. The default simplification function used by AMICI is powsimp which applies only very reasonable simplifications related to powers (and for example would not touch the sigmoid defined in the opening example). Instead, the changes applied by sympy.sympify are smaller and not real simplifications, but more side-effects of how SymPy constructs its expressions (a sympified string is simply converted to Python code and evaluated). In many cases these are actually "complications" (e.g., sympifying (x - 3.1)**2 results in 9.61*(0.32258064516129*x - 1)**2). This may sometimes result in actual simplification (e.g., x - x sympifies to zero), but even if left unevaluated most of those expression would be probably simplified by powsimp (this happens for x - x). On the whole, my impression is that sympify with evaluate=True has some potentially disastrous side-effects and few (if any) beneficial effects in our use case. Obviously, before changing default behaviour we could test this hypothesis a bit more.

Again not related to this issue, as @FFroehlich pointed out that sometimes the user may want deeper simplification and with the current ODE exporter they easily can. My personal experience with SymPy is that simplify will often create expressions that are not very reasonable as code for numerical evaluation (and in some cases introduce disastrous changes, like the one in the opening sigmoid example). In general creating very good expressions for floating point arithmetic is a tough problem (and depends on the range of inputs...). As @JanHasenauer suggested, some checks could be made to see if everything works correctly. On the spot I can think of:

1. Compare number of operations for simplified vs non-simplified version, and choose the shortest expression (could be done also recursively on subexpressions, probably with much better results but with a big-slowdown I am afraid)
2. Check that no invalid float literals are produced
3. Compare simplified vs non-simplified versions on some random inputs and check if the results are comparable. If something goes very wrong it should be easy to detect.

[dweindl]

Compare number of operations for simplified vs non-simplified version, and choose the shortest expression (could be done also recursively on subexpressions, probably with much better results but with a big-slowdown I am afraid)

Also afraid that, for large models, this might get extreeemly slow.

Check that no invalid float literals are produced

This would ideally check the whole (in AMICI unknown) domain. For example, the reason for the current powsimp there was mostly, that otherwise the generated expressions for Hill kinetics evaluate to NAN or Inf (forgot which) at concentrations of 0.0. That might actually be solved with evaluate=False.

[FFroehlich]

I agree that simplifications are useful, especially for generated code. But AMICI is not a model generator, so that the onus of providing nicely written formulas should be on the user (or on the software that generates the model for the user). I believe diligent users that pay attention to their formulas should not be penalized by default.

Right, but realistically that is not happening for most toolboxes. Those toolboxes also do not know about AMICI requirements for large models and the kind of simplifications. For AMICI, we want short expressions that are still numerically stable to evaluate, which is not necessarily what other toolboxes will produce.

That aside, this issue concerns mainly the reading of formulas from SBML/PEtab. The default simplification function used by AMICI is powsimp which applies only very reasonable simplifications related to powers (and for example would not touch the sigmoid defined in the opening example). Instead, the changes applied by sympy.sympify are smaller and not real simplifications, but more side-effects of how SymPy constructs its expressions (a sympified string is simply converted to Python code and evaluated). In many cases these are actually "complications" (e.g., sympifying (x - 3.1)**2 results in 9.61*(0.32258064516129*x - 1)**2). This may sometimes result in actual simplification (e.g., x - x sympifies to zero), but even if left unevaluated most of those expression would be probably simplified by powsimp (this happens for x - x). On the whole, my impression is that sympify with evaluate=True has some potentially disastrous side-effects and few (if any) beneficial effects in our use case. Obviously, before changing default behavior we could test this hypothesis a bit more.

Again not related to this issue, as @FFroehlich pointed out that sometimes the user may want deeper simplification and with the current ODE exporter they easily can. My personal experience with SymPy is that simplify will often create expressions that are not very reasonable as code for numerical evaluation (and in some cases introduce disastrous changes, like the one in the opening sigmoid example). In general creating very good expressions for floating point arithmetic is a tough problem (and depends on the range of inputs...). As @JanHasenauer suggested, some checks could be made to see if everything works correctly. On the spot I can think of:

1. Compare number of operations for simplified vs non-simplified version, and choose the shortest expression (could be done also recursively on subexpressions, probably with much better results but with a big-slowdown I am afraid)
2. Check that no invalid float literals are produced
3. Compare simplified vs non-simplified versions on some random inputs and check if the results are comparable. If something goes very wrong it should be easy to detect.

I think this makes a lot of sense. This was never really addressed in AMICI since the majority of models we have worked with so far were pure mass action models where these issues do not appear. I think having more control over the simplifications that are performed definitely makes sense, sp.simplify may indeed not be the best default simplification function. I agree with Daniel though that computational load for large models is a concern. Something that only gives a minor benefit for small models and is too complex to compute for large models, won't help much.

[lcontento]

I completely agree that computation time is going to be a big limiting factor. If a SBML model makes heavy use of assignment rules, it may advantageous to move some simplification at the stage of reading in the SBML (both simplifying the rules' expressions and the containing expressions before substitution).

I would start by comparing the results of evaluate=False/True (possibly on the PEtab benchmark collection?) so that we can have an idea of any difference of formula complexity/time. Probably not a big difference, but we can discuss about which behaviour should be the default after that. After that I would add a check for invalid string literals, which should be a rather inexpensive check. This could be a first pull request and I would consider this issue closed at that point.

After that, we could start investigating what is the best way to apply simplifications stronger than the current default powsimp (which is almost no simplification). Depending on how the things look at the end, we may even change the default. I guess this is going to be tougher; it may be a good idea to look around if there is someone else who has tried to solve similar problems and has had a modicum of success.

How does this plan sound? I don't think this is a particular pressing issue (otherwise there would be many more people other than me noticing problems) and I could start working on it after my current project becomes less busy (maybe 1/1.5 months). Obviously if there is someone else who is interested and can work on it right away (or wants to collaborate later), please feel free to do it.

[FFroehlich]

I completely agree that computation time is going to be a big limiting factor. If a SBML model makes heavy use of assignment rules, it may advantageous to move some simplification at the stage of reading in the SBML (both simplifying the rules' expressions and the containing expressions before substitution).

I would start by comparing the results of evaluate=False/True (possibly on the PEtab benchmark collection?) so that we can have an idea of any difference of formula complexity/time. Probably not a big difference, but we can discuss about which behaviour should be the default after that.

I would rather test this on Performace Test model, thats where we have extensively tested compilation/parsing speed and where changes usually have the biggest impact.

After that I would add a check for invalid string literals, which should be a rather inexpensive check. This could be a first pull request and I would consider this issue closed at that point.

After that, we could start investigating what is the best way to apply simplifications stronger than the current default powsimp (which is almost no simplification). Depending on how the things look at the end, we may even change the default. I guess this is going to be tougher; it may be a good idea to look around if there is someone else who has tried to solve similar problems and has had a modicum of success.

Agreed, sounds good!

How does this plan sound? I don't think this is a particular pressing issue (otherwise there would be many more people other than me noticing problems) and I could start working on it after my current project becomes less busy (maybe 1/1.5 months). Obviously if there is someone else who is interested and can work on it right away (or wants to collaborate later), please feel free to do it.

Sounds good!