These are chat archives for symengine/symengine

27th
Oct 2015
Ralf Stephan
@rwst
Oct 27 2015 07:36

@certik: As to series in Sage there are two implementations used: GiNaC and one based on polynomials (flint). The first supports expansion of most symbolic function expressions but is buggy; the second is fast but does not support symbolic functions. I proposed some fixes to the first but met with disinterest. However, the situation with SymPy is ridiculous as well:

In [12]: %time _ = (1/(1-X)).series(n=100)
CPU times: user 3.06 s, sys: 8 ms, total: 3.06 s
Wall time: 3.06 s

So, this is what motivates me with symengine series expansion.

As to your design decisions I can't put my finger on any flaw. With series you want both lazy and indexed generation. I'm not sure how sparse the usual series really is, certainly not sparse if in one variable? The speed of series expansion of expressions would be independent from the kind of polynomial used, but the speed of manipulation would be affected. Note the above example is about expansion.

Isuru Fernando
@isuruf
Oct 27 2015 07:42
@rwst, btw the implementation @certik mentioned is not in the public API.
In [6]: from sympy.polys.ring_series import *

In [7]: %time _ = rs_series(1/(1-x), x, 100)
CPU times: user 4 ms, sys: 0 ns, total: 4 ms
Wall time: 4.13 ms
Ralf Stephan
@rwst
Oct 27 2015 07:49
@isuruf , where do you get the function rs_series? I only have ..._inversion and _from_list in SymPy git master.
Isuru Fernando
@isuruf
Oct 27 2015 07:49
sympy.polys.ring_series.rs_series
Ralf Stephan
@rwst
Oct 27 2015 07:53
Ah my fault, thanks.
wow that's fast, so no need for symengine?
Ralf Stephan
@rwst
Oct 27 2015 07:59
In [3]: from symengine import var
In [4]: var('f w x y z')
Out[4]: (f, w, x, y, z)
In [5]: %time _ = rs_series(1/(1-x), x, 100)
...
GeneratorsError: expected a string, Symbol or expression or a non-empty sequence of strings, Symbols or expressions
Isuru Fernando
@isuruf
Oct 27 2015 07:59
Yes, it's fast. We'll need to port it to C++ for use in other wrappers, C, Julia and Ruby
@shivamvats, can you look into the above error?
rs_series doesn't call sympify before using the expression and symbols
Issue similar to sympy/sympy#9963
Ralf Stephan
@rwst
Oct 27 2015 08:05
Solvable in symengne.py?
*symengine.py
Double wow, forget Sage, this is vastly superior.
Ralf Stephan
@rwst
Oct 27 2015 08:11
They haven't even started to think about Puiseux series, and there they are...
Isuru Fernando
@isuruf
Oct 27 2015 08:17
Nope. have to be done in SymPy. or call _sympy_ before passing symengine objects
Ralf Stephan
@rwst
Oct 27 2015 08:34
@isuruf , but how are rs_series expanded? The example you gave does not the same as mine:
In [5]: rs_series(1/(1-x), x, 100)
Out[5]: (1/(-x + 1))

In [6]: series(1/(1-x), x, n=100)
Out[6]:
2    3    4    5    6    7    8    9    10    11    12    13    14    15    1
1 + x + x  + x  + x  + x  + x  + x  + x  + x  + x   + x   + x   + x   + x   + x   + x
...
Isuru Fernando
@isuruf
Oct 27 2015 08:37
Yeah, it's strange. I don't know about the internals of rs_series. @certik and @shivamvats might be able to
Ralf Stephan
@rwst
Oct 27 2015 08:43
Probably sympy/sympy#9843 because muls and adds work, also some divs do (1/sin).
Still, nearly 100x faster than series with cos*sin.
So, first fix tht ticket, then translate.
Ralf Stephan
@rwst
Oct 27 2015 08:50
In [23]: from symengine import *
%time _ = rs_series(cos(x)*sin(x), x, 100)
...
AttributeError: 'symengine.lib.symengine_wrapper.Sin' object has no attribute 'has'
The whole function stuff is missing too?
Ralf Stephan
@rwst
Oct 27 2015 09:16
Or rather, symengine.Function is not an symengine.Expression, and even if so, Expression has no has().
Ralf Stephan
@rwst
Oct 27 2015 09:52
@isuruf @certik : how do you rate the need to include an Expr translation versus numeric or other ring improvements in symengine? While manipulation of series probably is possible without Expr, series expansion from symbolic expressions is not. You might also need Expr for other things like simplification of expressions involving functions.
Shivam Vats
@shivamvats
Oct 27 2015 11:20

Thank you @rswt for your interest in ring_series.

The reason rs_series is not in public API is that the function is buggy (as you yourself figured out :) ) and not implemented fully. rs_series is basically a wrapper that calls the required functions (implemented in polys/ring_series.py) to expand a given expression. Currently it works with expressions involving sin, cos, tan and exponly (though making it work with other functions isn't difficult).

However, most of the elementary functions are fully implemented including series reversion, nth root and others. To use them, you need to create rings (a little extra work). For example, let's take your expression:

In [17]: from sympy.polys.ring_series import rs_series_inversion

In [18]: R, x = ring('x', QQ)

In [19]: %timeit rs_series_inversion(1-x, x, 100)
100 loops, best of 3: 1.7 ms per loop

In [20]: %timeit rs_series_inversion(1-x, x, 10000)
10 loops, best of 3: 115 ms per loop

In [21]: %timeit (1/(1-a)).series(a,0,100)
1 loops, best of 3: 2.38 s per loop
Shivam Vats
@shivamvats
Oct 27 2015 11:28

You might also want to read https://github.com/sympy/sympy/blob/master/doc/src/modules/polys/ringseries.rst

@certik Sympy's online doc repo isn't automatically updated. Do I need to clone sympy_doc and update it?

@isuruf Thank you for pointing out the error. It should probably give a warning in cases it can't handle.
Shivam Vats
@shivamvats
Oct 27 2015 11:35
Sorry wrong tag! I mean @rwst
Isuru Fernando
@isuruf
Oct 27 2015 12:35

@rwst, symengine.Basic and sympy.Basic are two different things. When we pass a symengine.Basic object to a sympy function, sympy first converts it into a sympy.Basic and then carries on the computation. We'd like to have functionality that doesn't need conversions.

(One way would be to have symengine objects behave as sympy objects by manipulating getattr and __instancecheck__ , but that'll be a lot of work to go through sympy codebase and changing little things that are incompatible. http://stackoverflow.com/questions/6803597/how-to-fake-type-with-python )

Ralf Stephan
@rwst
Oct 27 2015 13:44
Thanks @shivamvats @isuruf , I need to put my head around this.
Ondřej Čertík
@certik
Oct 27 2015 14:43
@rwst yes, the ring series in sympy is very fast. The reason it is fast is because it uses the polynomial module in sympy. SymPy has dense and sparse polynomials, the sparse polynomials happen to be faster. So we use those. Piranha is still quite a bit faster than the sparse (or dense) polynomials in sympy. So by using it in SymEngine to implement series expansion should make it faster than what's in sympy. However, we've done it first in SymPy, because it is easier to implement, and still useful.
@shivamvats what would it take to finish the rs_series function not to be buggy?
@rwst how can you use the series expansion based on Flint? Does it work for functions like sin(x)?
Ralf Stephan
@rwst
Oct 27 2015 14:56
@certik , Sage ring series only implement sqrt,exp,logof their own elements. GiNaC series have expansion of generic symbolic expressions (through differentiation I guess).
Ondřej Čertík
@certik
Oct 27 2015 14:57
Can you post an example how to use the ring series in Sage? Let me try the speed.
The GiNaC series works similar to sympy --- it expands inside out, e.g. if you do sin(x)*cos(x), then it first expands sin(x), cos(x) and then multiplies out the series.
Ralf Stephan
@rwst
Oct 27 2015 15:02
sage: R.<t> = PowerSeriesRing(QQ)
sage: %time _ = (t^2+1).sqrt(prec=100)
CPU times: user 2 ms, sys: 0 ns, total: 2 ms
Wall time: 1.68 ms
Wall time: 1.35 ms
Wall time: 1.37 ms
Wall time: 1.48 ms
Ralf Stephan
@rwst
Oct 27 2015 15:24
GiNaC for comparison:
sage: var('x')
x
sage: %time _ = sqrt(x^2+1).series(x,100)
CPU times: user 38 ms, sys: 0 ns, total: 38 ms
Wall time: 38.5 ms
Wall time: 37.6 ms
Wall time: 40.8 ms
Wall time: 40.6 ms
Ralf Stephan
@rwst
Oct 27 2015 15:34
@isuruf , I think the point is that only SymPy has a class Expr while I don't think e.g. Ruby has one. Also you don't want to convert SymPy Exprs to symengine and vice versa. So you need to either use SymPy Exprs from symengine (which is only useful temporarily) or you use symengine Expressions from SymPy (so you need a translation). Wrong?
Better: ...so you need Expression to have the full Expr API.
But Sage has its own expression, so no way around conversion there?
Ondřej Čertík
@certik
Oct 27 2015 15:46

@rwst here is how to do this using sympy:

In [1]: from sympy.polys import ring, QQ

In [2]: from sympy.polys.ring_series import rs_nth_root

In [3]: R, x = ring('x', QQ)

In [4]: %time _ = rs_nth_root(x**2+1, 2, x, 100)
CPU times: user 17 ms, sys: 4 ms, total: 21 ms
Wall time: 17.3 ms
Wall time: 15.1 ms
Wall time: 12.3 ms
Wall time: 15.7 ms
Wall time: 12.3 ms
Wall time: 13.3 ms
Wall time: 13.4 ms
Wall time: 13.3 ms
Wall time: 14.2 ms
Wall time: 14.9 ms

On my computer, Sage gives:

sage: R.<t> = PowerSeriesRing(QQ)
sage: %time _ = (t^2+1).sqrt(prec=100)
CPU times: user 9 ms, sys: 0 ns, total: 9 ms
Wall time: 8.49 ms
Wall time: 4.31 ms
Wall time: 4.25 ms
Wall time: 4.07 ms
Wall time: 1.98 ms
Wall time: 3.57 ms

You can check that both sympy and sage answers in _ are exactly the same. My experience with Piranha is that it can easily be 10x faster than the sympy polynomials. So from the timings above, you can see that there is a huge potential for speedup. What is stopping Sage from extending the PowerSeriesRing to work for all symbolic functions?

Ondřej Čertík
@certik
Oct 27 2015 15:53
@rwst as to the Expr thing --- I am still a bit confused about what exactly the problem is. My view is that the C++ SymEngine library doesn't care about this at all. We are just trying to implement things in some maintainable way, currently we settled on using Basic, Mul, Pow, ... hierarchy and implement most functionality using the visitor pattern or single dispatch, so Basic doesn't need many methods. We are keeping an option to perhaps do things differently if they turn out to be faster. Either way though, this shouldn't matter at all for how it is actually used from Python, Ruby or Julia. Let's talk about just Python: there the wrappers are in implemented using Cython, and they are free to introduce any kind of classes (including Expr or Sage's Expression if needed), and the point of the wrappers is to make sure that things work out of the box from SymPy and Sage. The only job of the C++ SymEngine library is to ensure that the library's C++ API is implemented in such a way so that the wrappers can be written to do what they need. And so I think your question is how to improve the wrappers to operate better with Sage and SymPy, am I right?
Ondřej Čertík
@certik
Oct 27 2015 16:00
That was the reason we split the wrappers, so now in the (pure) C++ symengine/symengine repository, we only have to worry about speed, correctness, maintainability and a usable API and we can concentrate on these things without worrying or even testing any kind of wrappers. In the wrappers (symengine/symengine.py, or .jl, .rb), we simply just use the C++ (or C) API and the only thing we care is so that the (Python) wrapper can be used from sympy/Sage (and we test that in the test suite), and that it doesn't introduce unnecessary overhead in terms of speed. Ruby or Julia wrappers then care about interoperability with other libraries in those languages.
Shivam Vats
@shivamvats
Oct 27 2015 17:00
@certik Right now there are two open issues: sympy/sympy#9843 and sympy/sympy#10050. I'll send a PR for the 2nd one soon. However, the more people use it, the more confident we can be about it. And we also need to make rs_series work with rest of the ring_series functions.
Ralf Stephan
@rwst
Oct 27 2015 17:07

@certik , as to why Sage series are not being improved, maybe the fact that people would then have to integrate symbolic expressions into the PowerSeriesRingwhich means working on symbolics (most shy away from that also because reviewers are lacking).

My problem with Expr is just that it doesn't exist in symengine, to use for a translation of rs_series, and I worry about how much of SymPy I need to translate to be able to translate rs_series.

Ondřej Čertík
@certik
Oct 27 2015 18:06
@shivamvats perfect, thanks! Before we start implementing it in SymEngine, we should finish up the SymPy implementation and try to make it available using the default series() function.
Ondřej Čertík
@certik
Oct 27 2015 18:13
@rwst thanks for the background regarding PowerSeriesRing. As to Expr, we can easily introduce it into the wrappers, by simply introducing the Expr class and make all the other classes subclass it.
As to translation of rs_series, are you talking about translation into SymEngine?
For that, my plan is to first get Piranha working so that we can call it from SymEngine, I think that's mostly done. Then we need to translate the rs_* functions from SymPy into C++, they will call Piranha underneath. That will give us an immediate benchmark that we can run and see what the speed is compared to Sage or other libraries, and if this approach is worth pursuing.
Finally, we then need to create something like rs_series that calls the primitive functions underneath, and hook it all up with the rest of SymEngine, and also expose it in Python wrappers.
@shivamvats so I guess we can actually start translating some of the primitive functions into C++ now. Those are well defined and it should work.
Ralf Stephan
@rwst
Oct 27 2015 19:21
@certik , exactly, that is what I'm interested in, and what's sorely missing. I'm glad you gave hints on what could be the next steps.
Ondřej Čertík
@certik
Oct 27 2015 19:54
@rwst cool. We can work on it together. I'll try to put up a minimal benchmark using Piranha, for example the sin(x)*cos(x) or something like that, and we can compare it against what's in Sage, or what could be in Sage in terms of Flint. If you have other ideas or other designs, let us know.
Shivam Vats
@shivamvats
Oct 27 2015 20:42
@certik Sure! Which poly module will we be using?
Ondřej Čertík
@certik
Oct 27 2015 22:26
@shivamvats For symengine? I would use Piranha as the underlying polynomial library, and then port your rs_* functions form sympy into symengine.
Let's port rs_sin, rs_cos and rs_mul and then run the sin(x)*cos(x) benchmark against other software and see how we are doing.
Then we'll go from there.