Plan for porting SymPy on top of SymEngine

[Ondřej Čertík]

Hi,

Here is a detailed step by step plan how to port SymPy on top of SymEngine:

https://github.com/symengine/symengine/wiki/SymPy-core-upgrade-to-SymEngine

The plan allows to still have a reference implementation in pure
Python if we want to (we probably do), so that SymPy stays pure
Python. The plan allows to work incrementally, making sure that
nothing gets broken in the end. It scales well to lots of developers
-- many people can help port the individual symbols, in a step by step
fashion and passing Travis tests will ensure that a given PR is
complete.

Let me know what you think, we can add more details as needed. I've
used a similar approach on another code (upgrading mesh data
structures in a big multi physics code) and the process works. We just
need to adapt it for the SymPy's case.

It's a pretty big undertaking, but it's doable and it won't result in
a broken SymPy.
And it will make SymPy very competitive in terms of speed.

I don't expect we do this plan right away, we still need to improve
SymEngine quite a bit (and also implement assumptions), but it is
important to have a solid plan in place, so that it is clear that our
approach is viable, to both SymPy and SymEngine.

Let me know what you think.

Ondrej

[Ralf Stephan]

Great. I was already puzzling about how to avoid the conversion costs

visible with the recent series implementation.

Also, I'd rather work from the bottom up than port selected functionality.

I really needed something like that.

Regards,.

--
You received this message because you are subscribed to the Google Groups "symengine" group.
To unsubscribe from this group and stop receiving emails from it, send an email to symengine+...@googlegroups.com.
To post to this group, send email to syme...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/symengine/CADDwiVC4Gget92Jm5Us4sdciAG_xDLDiecaTZnFO8WL6uX_CAw%40mail.gmail.com.
For more options, visit https://groups.google.com/d/optout.

[Ralf Stephan]

>5. Remove the conversions SymPy -> SymEngine -> SymPy

This depends on the method's arguments being available

as SymEngine.py types which might not be the case.

There might even be circular dependencies that cannot be

resolved in a single step.

[Ralf Stephan]

There might even be circular dependencies that cannot be

resolved in a single step.

Check out these visualizations:

https://github.com/rwst/sympy-coredep  

[Ondřej Čertík]

Thanks Ralf for the feedback.

Yes, this plan allows us to concentrate on SymEngine without worrying
about compatibility --- as long as there is some way to get the same
job done, i.e. some way to write a temporary compatibility wrapper.

Regarding your point:

This is a good point and we have to figure out the details. This is
already part of the step 1., i.e. figuring out which functions belong
to the core. I am thinking most of the functions that are implemented
in SymEngine (sin, cos, log, ...) should belong to the core. Some of
the special functions that we do not have in SymEngine (e.g. spherical
harmonics) would not be part of the core. As we keep improving
SymEngine, more and more functions will be part of it, but almost for
sure there will be things that are in SymPy, but not in SymEngine at
the point when we want to execute this plan. Let's take the spherical
harmonics Ylm(theta, phi) as an example.

If we have an expression, say:

x+y+Ylm(theta, phi)

then the SymPy -> SymEngine conversion converts x, y, + into SymEngine
classes, but Ylm will be the FunctionWrapper class. That already
works. FunctionWrapper already knows how to differentiate, we might
need to add a few more algorithms (perhaps series expansion).

The SymEngine -> SymPy conversion converts everything back to sympy.
We'll make sure that SymPy is always a superset of features of
SymEngine, so this conversion always works.

So then at the step 5., the outputs as well as inputs will become x,
y, + as SymEngine classes, and Ylm as a FunctionWrapper C++ class
(where the Python wrappers hook up Python implementation from SymPy).
So there is no problem.


As to circular dependencies, there are definitely circular
dependencies, e.g. I think the limit code is sometimes called in
series expansion and so on. But the idea is that SymEngine itself does
not have circular dependencies, so if we provide a series expansion
that doesn't depend on limits, then the plan works and can be
executed. In other words, I think the SymPy core can have circular
dependencies, and as long as the SymEngine core doesn't, then we
should be fine. E.g. in SymPy limits depend on series, and series
depend on limits. But if we provide series implementation using
SymEngine + a wrapper that is an equivalent replacement for SymPy's
series, and we do not use limits, then there is no problem.

We should add particular examples and I also plan to perhaps port one
or two things in SymPy, just to get a feel for the issues and document
it on the wiki on examples.

Let me know if you have further feedback.

Ondrej

[Isuru Fernando]

One other detail that we have to figure out is how to deal with incompatibilities of SymEngine and SymPy.

For ex.

SymEngine has Complex data type, while SymPy doesn't have any.
SymPy has singletons like half, one, while SymEngine doesn't have those.

Another approach for step 1 is fake typing. What it does is in SymEngine objects, there is a reference to its SymPy counterpart which may be initialized to None. __getattr__ is overriden in SymEngine object so that methods which are not available in SymEngine are dispatched to the SymPy object. __instancecheck__ has to be overriden in SymPy to make sure `isinstance(a, Add)` works. What this cannot fake are the singletons.

Isuru Fernando

--
You received this message because you are subscribed to the Google Groups "symengine" group.
To unsubscribe from this group and stop receiving emails from it, send an email to symengine+...@googlegroups.com.
To post to this group, send email to syme...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/symengine/CADDwiVBDfd2FYH-Orcar4dBcz%2BRUuJj9nLPGaTQuO2Lx73v2MQ%40mail.gmail.com.

[Ondřej Čertík]

On Sun, Feb 14, 2016 at 4:23 AM, Isuru Fernando <isu...@gmail.com> wrote:
> One other detail that we have to figure out is how to deal with
> incompatibilities of SymEngine and SymPy.
> For ex.
> SymEngine has Complex data type, while SymPy doesn't have any.
> SymPy has singletons like half, one, while SymEngine doesn't have those.

So I think that's a non-issue, as long as we can write wrappers that
emulate SymPy, including the singletons and such. And I think we
already do all these conversions, e.g. a complex number (1+3*I) which
is an expression in SymPy is converted to the appropriate SymEngine
class and vice versa.

Differences like these are the reason for the step 3. in the plan.
Even if we would otherwise emulate SymPy 1:1, we still have
differences like these, due to speed reasons. So for that reason, we
simply think of SymEngine having a completely incompatible API and we
write a compatibility wrapper in step 3. Of course, since SymEngine is
actually not that different, it will be quite easy to write the
wrapper in most cases.

>
> Another approach for step 1 is fake typing. What it does is in SymEngine
> objects, there is a reference to its SymPy counterpart which may be
> initialized to None. __getattr__ is overriden in SymEngine object so that
> methods which are not available in SymEngine are dispatched to the SymPy
> object. __instancecheck__ has to be overriden in SymPy to make sure
> `isinstance(a, Add)` works. What this cannot fake are the singletons.

I don't like that, since we should not be mixing SymEngine with SymPy like that.

I noticed we do something like that in the series, that we call SymPy
if SymEngine can't do it, but I really don't like it for multiple
reasons.

1) One reason is a cleanliness. I think the SymEngine.py should
simply be a thin and faithful wrapper of the SymEninge C++ library,
and functionally equivalent to SymEngine.rb or SymEngine.jl. It should
also follow such API, that makes the most sense. In other words, if
there are two equally good options (say for argument order and types
or method names), and one option is used by SymPy, we should
definitely use the one in SymPy. But if there is a good reason to use
a better API, then we should use a better API.

2) The other reason is that this introduces a circular dependency on
SymPy. So in step 3., we implement some of the wrappers to call
SymEngine series, but if that series calls SymPy series silently in
some circumstances, then even if you implement series using
SymEngine.py, you do not have guarantees, that you got rid of the
SymPy series implementation. In other words, at the end of step 3, you
must be running 100% on top of SymEngine. This will not be the case if
SymEngine.py is silently calling back to SymPy core.

So for these two reasons, we should not call any SymPy functionality
from the wrappers. The only places that we call into SymPy are:

* SymPy -> SymEngine conversion, this includes creating FunctionWrapper
* SymEngine -> SymPy conversion
* The only possible call back to SymPy from SymEngine.py besides the
previous two cases is via the FunctionWrapper. That is fine though,
since by definition, what's converted into a FunctionWrapper is not
part of the "core" as defined in step 1, so there is no circular
dependency there (i.e. if this object in Python calls into SymPy core,
it would get dispatched to use SymEngine at the end of step 3., since
we already converted all non-core code in SymPy to use
`old_core_api.py` in step 2.)



Ralf, Isuru, let me know what you think.

Ondrej

>

> https://groups.google.com/d/msgid/symengine/CA%2B01voPGJgkk9gbPFBd1fO9mBTR_qo3KWvHyZM2hq0U6ozkrrg%40mail.gmail.com.

[Ralf Stephan]

Ondrej,
what you mention with series was removed yesterday with
https://github.com/symengine/symengine.py/pull/49

To view this discussion on the web visit https://groups.google.com/d/msgid/symengine/CADDwiVCXyzKDy7n7-nLF8N5rEFAs8wPvK41LG5_5_VZGt93fqA%40mail.gmail.com.

[Ralf Stephan]

I agree with the plan, although I think the core will need
a lot of refactoring before it can even be started. Also
many things outside are called from the core, some of
which will need to be inside and vice versa.

[Aaron Meurer]

I don't see anything here about modifying SymPy's core. I think it
will be necessary, because there are a lot of things in SymPy's core
which do not really belong there. A common example is egregious use of
assumptions, but there are others as well.

The other question is how you will deal with places where SymEngine
has made different design decisions than SymPy. For example, I
remember from a discussion that SymEngine always combines
exp(x)*exp(y) -> exp(x + y) (let me know if this has changed). SymPy
does not do this. But if there is no way for SymEngine to represent
exp(x)*exp(y) without it combining to exp(x + y) then it will cause
issues when used as a core for SymPy, because there are places in
SymPy where it expects the two to be different.

Aaron Meurer

[Ondřej Čertík]

On Mon, Feb 15, 2016 at 9:56 AM, Aaron Meurer <asme...@gmail.com> wrote:
> I don't see anything here about modifying SymPy's core. I think it
> will be necessary, because there are a lot of things in SymPy's core
> which do not really belong there. A common example is egregious use of
> assumptions, but there are others as well.

So the step 1 is to figure out what the "core" is, that will be
exchanged to use SymEngine, and what is "non-core", or client of the
core.

We actually might be able to port assumptions to new assumptions using
a similar approach. But I think we already figured out how to get rid
of the old assumptions in the GSoC plan from few years ago. I'd like
to put it into a wiki.

> The other question is how you will deal with places where SymEngine
> has made different design decisions than SymPy. For example, I
> remember from a discussion that SymEngine always combines
> exp(x)*exp(y) -> exp(x + y) (let me know if this has changed). SymPy
> does not do this. But if there is no way for SymEngine to represent
> exp(x)*exp(y) without it combining to exp(x + y) then it will cause
> issues when used as a core for SymPy, because there are places in
> SymPy where it expects the two to be different.

So this is a good point. The general idea is that it doesn't matter if
we have a different design decision, as long as we can
express/represent the same computation in both. So things like how
complex numbers are represented shouldn't matter. But this
exp(x)*exp(y) -> exp(x + y) might matter. The way symengine (and
ginac/Sage) works is that it represents exp(x) as {E: x}, exp(y) as
{E: y} and exp(x)*exp(y) as {E: x+y}. It can't represent it in any
other way, since the hash table can have just one key E, so one has to
add the exponents. The exp(x + y) is also represented as {E: x+y}.
Note that exp(x)*exp(y) = exp(x + y) for all complex "x" and "y", so
mathematically there is no problem.

So I just need to look at the places in SymPy that depend or require
to distinguish between exp(x)*exp(y) and exp(x + y) and see if the
algorithm can be rewritten to not depend on this, more specifically to
simply assume the exp(x + y) form. I am hoping it is possible.

Either way, this needs to be resolved as well. Good point.

Ondrej

[Aaron Meurer]

On Tue, Feb 16, 2016 at 12:49 PM, Ondřej Čertík <ondrej...@gmail.com> wrote:
> On Mon, Feb 15, 2016 at 9:56 AM, Aaron Meurer <asme...@gmail.com> wrote:
>> I don't see anything here about modifying SymPy's core. I think it
>> will be necessary, because there are a lot of things in SymPy's core
>> which do not really belong there. A common example is egregious use of
>> assumptions, but there are others as well.
>
> So the step 1 is to figure out what the "core" is, that will be
> exchanged to use SymEngine, and what is "non-core", or client of the
> core.
>
> We actually might be able to port assumptions to new assumptions using
> a similar approach. But I think we already figured out how to get rid
> of the old assumptions in the GSoC plan from few years ago. I'd like
> to put it into a wiki.

I prefer my plan at
https://groups.google.com/forum/#!msg/sympy/PHR136kdxc4/C2qs5obPGoMJ,
which was started by Sudhanshu in his GSoC last year (you can also
look at his proposal
https://github.com/sympy/sympy/wiki/GSoC-2015-Application-Sudhanshu-Mishra:-Assumptions.

The key difference is that I no longer want to remove the syntax for
the old assumptions. Instead, we should merge the two assumptions
together, so that there is only one assumptions system. Some work on
this is already done, namely, ask(Q.positive(Symbol('x',
positive=True))) works.

Additionally, the use of assumptions should be more or less banned
from the core, since they impact performance too much, and greatly
increase its complexity (making it very hard for something like
symengine to build a hot swappable replacement). This is not simple,
because it means reducing the amount of automatic simplification that
happens based on assumptions, and removing automatic simplification is
always difficult.

>
>> The other question is how you will deal with places where SymEngine
>> has made different design decisions than SymPy. For example, I
>> remember from a discussion that SymEngine always combines
>> exp(x)*exp(y) -> exp(x + y) (let me know if this has changed). SymPy
>> does not do this. But if there is no way for SymEngine to represent
>> exp(x)*exp(y) without it combining to exp(x + y) then it will cause
>> issues when used as a core for SymPy, because there are places in
>> SymPy where it expects the two to be different.
>
> So this is a good point. The general idea is that it doesn't matter if
> we have a different design decision, as long as we can
> express/represent the same computation in both. So things like how
> complex numbers are represented shouldn't matter. But this
> exp(x)*exp(y) -> exp(x + y) might matter. The way symengine (and
> ginac/Sage) works is that it represents exp(x) as {E: x}, exp(y) as
> {E: y} and exp(x)*exp(y) as {E: x+y}. It can't represent it in any
> other way, since the hash table can have just one key E, so one has to
> add the exponents. The exp(x + y) is also represented as {E: x+y}.
> Note that exp(x)*exp(y) = exp(x + y) for all complex "x" and "y", so
> mathematically there is no problem.

Well it used to be the case, that SymPy did this transformation
automatically as well (you can check old versions). But then I put in
the work to remove it, because it was getting in my way, and what I
found was that things were implicitly depending on it happening. It
was a bit of a pain to change it, because I had to go into all the
test failures and figure out in the code (generally code I didn't
understand, like the limit algorithm) where it was implicitly
depending on this.

That's always the case with automatic simplification. If you put it in
there, then things will inevitably depend on it, and it makes it
difficult to remove. There's a reason that 2*(x + y) still
automatically distributes in SymPy, even though it's a huge gotcha for
tons of users.

>
> So I just need to look at the places in SymPy that depend or require
> to distinguish between exp(x)*exp(y) and exp(x + y) and see if the
> algorithm can be rewritten to not depend on this, more specifically to
> simply assume the exp(x + y) form. I am hoping it is possible.

But how will you regression test this?

Another way would be to implement another representation in SymEngine
that can represent exp(x)*exp(y) and use that when the user enters
such an expression. Assumedly it would be slower, so one would need to
switch to the better one automatically when such an expression isn't
there.

It's not a trivial thing, though, because then you have to make sure
that all your algorithms can handle both forms (basically the same
problem that I had with SymPy). It would mean symengine automatically
makes tradeoffs on its internal representation based on its form (a
similar example would be representing a polynomial expression
efficiently as a polynomial, but falling back to a general expression
tree if a non-polynomial term is added to it).

Aaron Meurer

[Ondřej Čertík]

Would you mind adding the plan to the SymPy's wiki somewhere? Even if
we don't have time to work on it, at least let's have a plan.

This would be helpful also so that we can design the assumptions in
SymEngine in a compatible way.

So in the above, we are talking about the internal data structures for
the Mul class. E.g. exp(x)*exp(y)*a^z is represented as {E: x+y, a: z}
in SymEngine. This data structure doesn't allow to treat exp(x)*exp(y)
as two separate terms. You are right that other data structures, like
the one in SymPy, allow to treat exp(x)*exp(y) separately, but it is
slower.

The polynomials are exactly the same issue. In fact, there are
multiple ways to represent polynomials themselves, with various speed
trade-offs.

What we decided is not to do automatic conversions from a polynomial
to Basic, not even to do automatic conversions from say polynomials
with integer coefficients to polynomials with rational coefficients.
Rather, the user decides what kind of representation to use at the
beginning, and that will be used through out, and if a conversion is
needed, the user does it explicitly.

The algorithms also depend on the data structure. For example,
multiplying polynomials with integer coefficients and integer
exponents allow you to pack the exponents into one big integers and
then use an <integer, integer> hash table. For other polynomials you
have to use a different algorithm.

Another example is a series expansion. There are multiple ways to
represent a series internally and then one writes specialized
algorithms for each given representation.

We definitely want to have all this.

So currently we only have MulHashTable (i.e. a Mul implementation
based on a hashtable, which is fast, but doesn't allow to represent
exp(x)*exp(y)). We can add other representations, say MulList (a Mul
implementation that represents the terms as a list, as SymPy does it).

Using the above philosophy, if the user specifies MulList, it will
*not* be automatically converted to MulHashTable. It will simply be
treated as a special function, essentially. The internal algorithms,
when they encounter such a term will simply call a specialized
algorithm for MulList. E.g. if you do <MulList> * exp(z) then it keeps
it as MulList.


Aaron, do you remember what applications prefer MulList to
MulHashTable? When is it useful to have this flexibility of keeping
the powers separate?

Ondrej

[Aaron Meurer]

On Tue, Feb 16, 2016 at 2:11 PM, Ondřej Čertík <ondrej...@gmail.com> wrote:
> On Tue, Feb 16, 2016 at 11:13 AM, Aaron Meurer <asme...@gmail.com> wrote:
>> On Tue, Feb 16, 2016 at 12:49 PM, Ondřej Čertík <ondrej...@gmail.com> wrote:
>>> On Mon, Feb 15, 2016 at 9:56 AM, Aaron Meurer <asme...@gmail.com> wrote:
>>>> I don't see anything here about modifying SymPy's core. I think it
>>>> will be necessary, because there are a lot of things in SymPy's core
>>>> which do not really belong there. A common example is egregious use of
>>>> assumptions, but there are others as well.
>>>
>>> So the step 1 is to figure out what the "core" is, that will be
>>> exchanged to use SymEngine, and what is "non-core", or client of the
>>> core.
>>>
>>> We actually might be able to port assumptions to new assumptions using
>>> a similar approach. But I think we already figured out how to get rid
>>> of the old assumptions in the GSoC plan from few years ago. I'd like
>>> to put it into a wiki.

I'll put it on my todo list, which is quite large so don't hold your
breath. But likely if a good student this year is interested in
working on assumptions for GSoC I will write stuff down for his/her
benefit (which I would likely do anyway in a mailing list thread, but
I'll try to remember to put stuff on the wiki if I write a lot).

For now, that mailing list thread from last year is the last time I
did a brain dump on the assumptions. It's a little out of date thanks
to Sudhanshu's project, but the longer term stuff still needs to be
done.

I agree that doing it explicitly is better. But how would it work for
a core for SymPy? Currently this is explicit in the sense that if you
want, say, polynomials in SymPy, you have to use Poly or ring() or
whatever. The default mode is to use expressions with Symbol and so
on, which are very generic, but slow (although some algorithms convert
to and from Poly automatically, because that is the only place they
are implemented, e.g., factor()).

>
> The algorithms also depend on the data structure. For example,
> multiplying polynomials with integer coefficients and integer
> exponents allow you to pack the exponents into one big integers and
> then use an <integer, integer> hash table. For other polynomials you
> have to use a different algorithm.
>
> Another example is a series expansion. There are multiple ways to
> represent a series internally and then one writes specialized
> algorithms for each given representation.
>
> We definitely want to have all this.
>
> So currently we only have MulHashTable (i.e. a Mul implementation
> based on a hashtable, which is fast, but doesn't allow to represent
> exp(x)*exp(y)). We can add other representations, say MulList (a Mul
> implementation that represents the terms as a list, as SymPy does it).
>
> Using the above philosophy, if the user specifies MulList, it will
> *not* be automatically converted to MulHashTable. It will simply be
> treated as a special function, essentially. The internal algorithms,
> when they encounter such a term will simply call a specialized
> algorithm for MulList. E.g. if you do <MulList> * exp(z) then it keeps
> it as MulList.
>
>
> Aaron, do you remember what applications prefer MulList to
> MulHashTable? When is it useful to have this flexibility of keeping
> the powers separate?

To some degree it's just aesthetics. You want expressions to look like
the way you write them down. Of course, not all aesthetics are valid
reasons to do this sort of thing (for instance, there's no way in
SymPy to distinguish between 1/exp(x) and exp(-x)). For exponentials,
my original motivation was to be able to split expressions into
products for dsolve().

For exp(x)*exp(y) -> exp(x + y) it might not actually be a big deal.
You would have to try changing it in sympy.Mul and seeing what tests
break.

Aaron Meurer

>
> Ondrej
>
> --
> You received this message because you are subscribed to the Google Groups "symengine" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to symengine+...@googlegroups.com.
> To post to this group, send email to syme...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/symengine/CADDwiVCz_DVjfFfFoviL8ao-5N81n-EGWjefuY5-kcMy%3DMpo-g%40mail.gmail.com.

[Ondřej Čertík]

Thanks. I'll try to understand this myself, so that I can push this. I
think that summer students can definitely do this, but this is the
kind of project that requires me or you to be completely on top of it
and know exactly what needs to be done and how.

SymPy already has two kinds of polynomials: dense and sparse. The user
just constructs them with different functions/constructors. It doesn't
matter what the default mode is, as long as you can explicitly
construct a different data structure and use it in expressions.

Looks like those are not enough reasons to keep MulList around. It
complicates things a lot to have multiple data structures in the
general symbolic parts, as we know in SymPy with regards to non
commutative mul and other things like handling of Order. What I want
to have for the general symbolics is just one data structure + related
algorithms, the one that seems to be the fastest. Which is what we
currently have in SymEngine. Then for polynomials or series expansion,
we have more specialized (and thus even faster) data structures.

Mathematica also doesn't seem to be able to represent Exp[x]*Exp[y],
it automatically converts to exp(x+y). You can hold it via
Hold[Exp[x]*Exp[y]], but my understanding is that's essentially
equivalent to a Python string "exp(x)*exp(y)", before calling eval()
on it. I.e. you can pass it around, but the minute you release the
hold, it becomes exp(x+y).

As to splitting products for dsolve, you can still split it by going
into the exponents and split Adds there.

Ondrej

[Ralf Stephan]

On Tuesday, February 16, 2016 at 8:29:50 PM UTC+1, Aaron Meurer wrote:

I agree that doing it explicitly is better. But how would it work for
a core for SymPy? Currently this is explicit in the sense that if you
want, say, polynomials in SymPy, you have to use Poly or ring() or
whatever. The default mode is to use expressions with Symbol and so
on, which are very generic, but slow (although some algorithms convert
to and from Poly automatically, because that is the only place they
are implemented, e.g., factor()). 

I don't think there will be a different data structure than an expression tree

as the default for a symbolic-heavy CAS like SymPy. They may be slow but

I can relate here that most of the time spent in symbolic expression computation

with Sage is spent in Python, not in Pynac which usually spends between 0.1 and 10%.

OTOH, if we can talk about removal of bits from the core, maybe things can

be added that need speed like ntheory/factor_.py. AFAI can see everything

apart from the divisors functions is already in symengine.