Work on using SymEngine as libCellML's CAS

[dladd]

As I only have a few weeks left in my contract to work on libCellML, it was suggested that it might be best to put a pin in work on the Validator ( #127 ) for now. The broad strokes are in place for validation and contributors can go through and add further checks for Specification rules. There are also a few TODO: ... notes in that pull to indicate further work I was planning to get to. If we still have time at the end of my contract, hopefully I can come back and clean up the remaining TODO's in #127 and #116.

Instead, it was suggested that I should focus on getting SymEngine working as our computer algebra system (CAS) for libCellML (see #95 for previous discussion on choosing a CAS). The largest part of this will likely be getting MathML parsing/printing into SymEngine, which we can then use to manipulate our math string.

[dladd]

For anyone wanting to follow along, I'm discussing work on printing/parsing MathML with the SymEngine team here.

[dladd]

To update on this, I made a start on printing MathML in SymEngine but found there were still a few missing pieces in terms of what we'll need for libCellML. For instance, support for relational operators* is still in the works. Also, as SymEngine is primarily focused on improving efficiency, its documentation is currently a bit less approachable than SymPy's (though the dev team is very helpful).

We had a discussion at the ABI and thought it would be a good idea to first set up a few examples using SymPy that cover some of our expected use-cases. It was suggested that some of the OpenCOR tutorial examples (e.g., a simple ODE, the Lorenz attractor) would be a good starting point. Any other recommendations of relatively simple tests that highlight issues a CAS might have with supporting maths from CellML are welcome :-)

This should give us a good indication of whether SymPy/SymEngine can cover what we need in a CAS. Assuming all goes well with SymPy, we can attempt the same set of examples using SymEngine. We should then have a better idea of what we may need to contribute so that libCellML can be supported by SymEngine alone.

* Basic support for relational operators isn't necessarily a prohibitive amount of work to contribute alone but highlights that it might be worth figuring out what is needed before we jump in.

[agarny]

You should indeed start with the two tutorial examples you mentioned, then we can see from there.

[dladd]

I've gotten some simple tests going in SymPy now (code, output). In general, SymPy seems to handle manipulations of systems of algebra and ODEs quite well, is nice to work with and well-documented. A few things to note:

• Getting dsolve() to substitute initial conditions into integration constants is still in progress in SymPy (see their pull 11264). Basic functionality seems to be working in that pull now (this is what I used for the output linked above).
• Numerically solving ODEs has been scoped out of SymPy and into the mpmath library. mpmath is still a bit limited in this regard (e.g., no Runge-Kutta algorithms yet). It seems like most people go with SciPy if they are using Python and want to numerically solve ODEs (Lorenz attractor example).
• Sympy doesn't seem to explicitly flag over/underdetermined systems of equations as such but instead returns no solution / the range of possible solutions. But perhaps we can do some of our own basic constraint counting during validation once we are able to manipulate symbols.

I'll now see how far I can get with an equivalent in SymEngine. I expect initial values to be problematic since they are only now being implemented in SymPy. I will also try to define systems without relational operators for now.

[dladd]

Actually, a bigger issue could be that it doesn't appear that SymEngine has solve() or dsolve() functionality. But, on second thought, maybe it is reasonable to expect that if an analytic solution to an ODE exists, a user is aware of this before they start implementing the math in a CellML model?

[agarny]

Actually, a bigger issue could be that it doesn't appear that SymEngine has solve() or dsolve() functionality. But, on second thought, maybe it is reasonable to expect that if an analytic solution to an ODE exists, a user is aware of this before they start implementing the math in a CellML model?

Never assume anything when it comes to the end user... :)

[dladd]

Never assume anything when it comes to the end user... :)

Fair point. But is having solve()/dsolve() a hard requirement for us? It is certainly nice for a user to have but for libCellML's minimal default CAS this may not be needed?

[nickerso]

I think this depends on how we can determine the actual system to be "solved" by an application when simulating a model. In the SimPy examples above, this is achieved using solve, but is there another way to get this information out of SimPy?

[dladd]

Good point- from what I've seen so far, solve() is indeed what is used to determine if a given system is solvable in SymPy (outside of special cases like Diophantine equations).

@nickerso also pointed out that we might want solve() functionality if we had a case where a user wanted to use the same system with different input/output variables. As an example of how we might do this with SymPy, I've added a Mooney-Rivlin example (code, output), where we use solve() to get a system of equations for stress in terms of strain and vice-versa. I also attempted this with a Guccione model but, as I write this, SymPy has been attempting to get equations for strain for over an hour...

[isuruf]

Mooney-Rivlin example is a linear system and solve dispatches to linsolve and then it's a simple linear system solve. On the other hand, Guccione model probably has a solution that is too hard for SymPy to handle. (A solution probably has lots of lambertw terms). Btw, there is a E * 33 in line 176 and 200. Is that supposed to be E33?

[dladd]

Thanks @isuruf ! Good catch- that was indeed supposed to be E33. With that fix I'm now getting a maximum recursion depth error, which I assume is just because this is a very difficult system to solve for strain (I would not enjoy doing that by hand!)

[dladd]

For anyone interested, I've also added some very simple SymEngine examples here .

[isuruf]

:D. That's why you should never use import *. The error is actually a NotImplementedError in the lambertw solver. Recursion occurs when handling that exception.

[dladd]

I'd like to poll the team's preference between:

• (A) Move toward using SymPy as our CAS for more features / functionality. With the language difference, I assume this will likely necessitate one of:
  • (A.1) Wrapping libCellML to be used from Python with something like SWIG or Boost.Python
  • (A.2) Using SymPy through C++ with something like Cython or Python embedding
• (B) Continue to focus on using SymEngine as our CAS . This will require more work to support the kinds of checking we want to do on systems of equations and will not (for now) have the same level of functionality or documentation. But it is a smaller dependency, written in C++, and would give us the ability to do basic symbolic manipulation.

[hsorby]

For various reasons, that I will elucidate if requested, I'm for B.

[mtcooling]

Also B.

On Mon, Aug 29, 2016 at 1:54 PM, Hugh Sorby notifications@github.com wrote:

For various reasons, that I will elucidate if requested, I'm for B.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/cellml/libcellml/issues/130#issuecomment-243015820, or mute the thread https://github.com/notifications/unsubscribe-auth/AD7gqIr9ImSxhB88uSOzNLsHrV7FlqDsks5qkju6gaJpZM4JfmO3 .

[agarny]

Definitely B.

[dladd]

Great, sounds pretty unanimous :-)

The way I now foresee the maths validation side of this working:

• 1. As part of our current math validation, we assemble a list of variables (and new bvar variables) and excise their associated units. From this we can get a variable : units map and a "clean" MathML string.
• 2. Validate the MathML using the DTD. If this returns no errors, we should be able to parse the equation(s).
• 3. To validate units consistency, for each equation we would:
  • 3.1. Parse the MathML string into a SymEngine expression
  • 3.2. Substitute all variables in the expression with their corresponding units using the map from (1)
  • 3.3. Substitute all units with base units
  • 3.4. Check that the expression now simplifies to zero.
• 4. In the case of a linear or polynomial system of equations, we can do some simple constraint counting- essentially just checking whether the number of unknowns is equivalent to the number of equations. Log any over/underdetermined system errors.

[agarny]

I would swap 3.2 and 3.3, so that a given substitution is done only once rather than several times (in case a given unit is used several times in a model). In fact, 3.3 might be done as part of 1.

Otherwise, when there isn't units consistency, I would like to be able to provide my end-users with some information as why it is the case, so that they can easily fix the problem. However, I get the feeling that with 3.4, I will only be able to tell them whether there is units consistency or not, and nothing more. Am I correct or am I missing something?

[nickerso]

My hope for 3.4 is to result in nice messages like "V/s != kg/K" for a given equation (or subpart of an equation). You would get back the result of the simplification/reduction and if its correct, that would be something like dimensionless or 0 = 0.

[dladd]

However, I get the feeling that with 3.4, I will only be able to tell them whether there is units consistency or not, and nothing more. Am I correct or am I missing something?

My hope for 3.4 is to result in nice messages like "V/s != kg/K" for a given equation (or subpart of an equation). You would get back the result of the simplification/reduction and if its correct, that would be something like dimensionless or 0 = 0.

This was what I was thinking also- we could return the string (or maybe store the SE expressions) for the equation, the units-substituted version (3.2), and the final reduced form (3.4). As a simple example, if we have

5*(x - y) = 0

and

x:units = meters y:units = seconds 5:units = dimensionless

we might log an error as:

Math in component 'component_name' contains an equation:
5*(x - y) = 0
whose unit substitution:
dimensionless*(meters - seconds) = 0
incorrectly reduces to:
meters = seconds

EDIT: Here's a SymPy example of this with our single equation ODE (code, output)

[agarny]

Thanks @dladd and @nickerso, that makes sense now. However, @dladd, when it comes to your example, I would expect/hope to be told that meters - seconds = 0, which would make it easier for an end-user to relate to. Anyway, your approach seems reasonable to me, so let's see how it all goes.

[dladd]

Thanks @agarny . Fair enough- I'd agree that whatever form this takes, it should be consistent.

I did a sneaky edit of my last comment and added a SymPy example of how this might work on our exponential decay case (code, output). This works using evaluation with SymPy's relational operators. So we'll probably want to get relationals going in SymEngine if we can. An alternative might be to look for addition or subtraction operators in the reduced expression tree, which should indicate incompatible units.

I also just realised that that example works because the exp() coefficient term simplifies to zero. As part of this, we'll also want to substitute special functions like exp(1), exp(-1), sin(1) with 1. This may also require some more thought about tracking units associated with non-variable numbers (i.e., MathML <cn> elements).

[agarny]

Indeed, things like exp(1) will have to be replaced with 1. Something similar will have to be done say a*1 which will have to be replaced with a.

Otherwise, and just FTR, in your 'incorrect' coded example, one would ideally get:

meter = (meter*second - meter) + meter

rather than:

                                        -1
meter = meter + (meter*second - meter)*e 

[isuruf]

A simpler solution would be to introduce a Dimensionless object as @nickerso mentioned. Dimensionless would be another numeric type that would have automatic simplification like D + D = D, D - D = D, D * D = D, D + 1 = D and also automatic simplification like exp(D) = D, etc. Then you would only have to do,

subMap = {
        t: ut*d,
        a: ua*d,
        b: ub*d,
        y0: uy0*d,
        y(t): uy*d
    }

and then SymEngine or SymPy would take care of the rest. You might want to expand the result before checking that it is 0. You can take a look at RealDouble class in SymEngine, for how to do this, or ping me on github or gitter.

[MichaelClerx]

What's the current status on this @hsorby @nickerso ? Are we still planning to have a CAS? I think we decided it was pretty much out of scope, no?

[hsorby]

Currently I am aiming for a limited self-implemented CAS. If we can get away with that then happy days otherwise we will look again at SymEngine.

[MichaelClerx]

Ok! Then closing this in favour of #232