Thursday, January 3, 2008

SymPy/sympycore (pure Python) up to 5x faster than Maxima (future of Sage.calculus?)

According to this test, sympycore is from 2.5x to 5x faster than Maxima. This is an absolutely fantastic result and also a perfect certificate for Python in scientific computing. Considering that we compare pure Python to LISP.

Ok, this made us excited, so we dugg deeper and ran more benchmarks. But first, let me say a few general remarks. I want a fast CAS (Computer Algebra System) in Python. General CAS, that people use, that is useful, that is easily extensible (!), that is not missing anything, that is comparable to Mathematica and Maple -- and most importantly -- I want it now and I don't care about 30 years horizons (I won't be able to do any serious programming in 30 years anyway). All right. How to do that? Well, many people tried... And failed. The only opensource CAS system, that has any chance of becoming the opensource CAS, in my own opinion, is Sage. You can read more about my impressions form Sage here. I am actually only interested in mathematical physics, so basically Sage.calculus. Currently Sage uses Maxima, because Maxima is old, proven, working system and it's reasonably fast and quite reliable, but written in LISP. Some people like LISP. I don't and I find it extremely difficult to extend Maxima. Also even though Maxima is in LISP, it uses it's own language for interacting with the user (well, that's not the way). I like python, so I want to use Python. Sage has written Python wrappers to Maxima, so Sage can do almost everything that Maxima can, plus many other things. Now. But the Sage.calculus has issues.

First, I don't know how to extend the wrappers with some new things, see my post in the sage-devel for details, it's almost 2 months old with no reaction, which shows that it's a difficult issue (or nonsense:)).

And second, it's slow. For some examples that Sage users have found out, even SymPy, as it is now, is 7x faster than Sage and sympycore 23x faster and with the recent speed improvements 40x faster than Sage.

So let's improve Sage.calculus. How? Well, no one knows for sure, but
I believe in my original idea of pure Python CAS (SymPy), possibly with some parts rewritten in C. Fortunately, quite a lot of us believe that this is the way.

What is this sympycore thing? In sympy, we wanted to have something now, instead of tomorrow, so we were adding a lot of features, not looking too much on speed. But then Pearu Peterson came and said, guys, we need speed too. So he rewrote the core (resulting in 10x to 100x speedup) and we moved to the new core. But first, the speed isn't sufficient, and second it destabilized SymPy a lot (there are still some problems with caching and assumptions half a year later). So with the next package of speed improvements, we decided to either port them to the current sympy, or wait until the new core stabilizes enough. So the new new core is called sympycore now, currently it only has the very basic arithmetics (and derivatives and simple integrals), but it's very fast. It's mainly done by Pearu. But for example the latest speed improvement using sexpressions was invented by Fredrik Johansson, another SymPy developer and the author of mpmath.

OK, let's go back to the benchmarks. First thing we realized is that Pearu was using CLISP 2.41 (2006-10-13) and compiled Maxima by hand in the above timings, but when I tried Maxima in Debian (which is compiled with GNU Common Lisp (GCL) GCL 2.6.8), I got different results, Maxima did beat sympycore.

SymPyCore:

In [5]: %time e=((x+y+z)**100).expand()
CPU times: user 0.57 s, sys: 0.00 s, total: 0.57 s
Wall time: 0.57

In [6]: %time e=((x+y+z)**20 * (y+x)**19).expand()
CPU times: user 0.25 s, sys: 0.00 s, total: 0.25 s
Wall time: 0.25

Maxima:

(%i7) t0:elapsed_real_time ()$ expand ((x+y+z)^100)$ elapsed_real_time ()-t0;
(%o9) 0.41
(%i16) t0:elapsed_real_time ()$ expand ((x + y+z)^20*(x+z)^19)$ elapsed_real_time ()-t0;
(%o18) 0.080000000000005


So when expanding, Maxima is comparable to sympycore (0.41 vs 0.57), but for general arithmetics, Maxima is 3.5x faster. We also compared GiNaC (resp. swiginac):

>>> %time e=((x+y+z)**20 * (y+x)**19).expand()
CPU times: user 0.03 s, sys: 0.00 s, total: 0.03 s
Wall time: 0.03


Then we compared just the (x+y+z)**200:

sympycore:
>>> %time e=((x+y+z)**200).expand()
CPU times: user 1.80 s, sys: 0.06 s, total: 1.86 s
Wall time: 1.92
swiginac:
>>> %time e=((x+y+z)**200).expand()
CPU times: user 0.52 s, sys: 0.02 s, total: 0.53 s
maxima:
(%i41) t0:elapsed_real_time ()$ expand ((x + y+z)^200)$ elapsed_real_time ()-t0;
(%o43) 2.220000000000027


Where GiNaC still wins, but sympycore beats Maxima, but the timings really depend on the algorithm used, sympycore uses Millers algorithm which is the most efficient.

So then we tried a fair comparison: compare expanding x * y where x and y are expanded powers (to make more terms):

sympycore:
>>> from sympy import *
>>> x,y,z=map(Symbol,'xyz')
>>> xx=((x+y+z)**20).expand()
>>> yy=((x+y+z)**21).expand()
>>> %time e=(xx*yy).expand()
CPU times: user 2.21 s, sys: 0.10 s, total: 2.32 s
Wall time: 2.31
swiginac:
>>> xx=((x+y+z)**20).expand()
>>> yy=((x+y+z)**21).expand()
>>> %time e=(xx*yy).expand()
CPU times: user 0.30 s, sys: 0.00 s, total: 0.30 s
Wall time: 0.30
maxima:
(%i44) xx:expand((x+y+z)^20)$
(%i45) yy:expand((x+y+z)^21)$
(%i46) t0:elapsed_real_time ()$ expand (xx*yy)$ elapsed_real_time ()-t0;
(%o48) 0.57999999999993

So, sympycore is 7x slower than swiginac and 3x slower than maxima. We are still using pure Python, so that's very promising.

When using sexpr functions directly then 3*(a*x+..) is 4-5x faster than Maxima in Debian/Ubuntu. So, the headline of this post is justified. :)

Conclusion

Let's build the car. Sage has the most features and it is the most complete car. It has issues, some wheels need to be improved (Sage.calculus). Let's change them then. Maybe SymPy could be the new wheel, maybe not, we'll see. SymPy is quite a reasonable car for calculus (it has plotting, it has exports to latex, nice, simple but powerfull command line with ipython and all those bells and whistles and it can also be used as a regular python library). But it also has issues, one wheel should be improved. That's the sympycore project.

All those smaller and smaller wheels show, that this is indeed the way to go, but very important thing is to put them back in the car. I.e. sympycore back to sympy and sympy back to Sage and integrate them well. While also leaving them as separate modules, so that users, that only need one particular wheel, can use them.

Tuesday, January 1, 2008

R.<a,b,c> = QQ[] what is that?

While playing with Jaap's wish, I finally got fedup and decided to study what the funny syntax "R.<a,b,c> = QQ[]" really means. So I did:

sage: R. = QQ[]
sage: preparse("R. = QQ[]")
"R = QQ['a, b, c']; (a, b, c,) = R._first_ngens(Integer(3))"

Now everything is crystal clear! Let's study what QQ is:

sage: QQ?
Type: RationalField
Base Class:
String Form: Rational Field
Namespace: Interactive
Docstring:

The class class{RationalField} represents the field Q of
rational numbers.

Cool, why not. So what is the _first_ngens method doing? Let's find out:

sage: R._first_ngens?
Type: builtin_function_or_method
Base Class:
String Form:
Namespace: Interactive

Hm, not really useful. Let's push harder:

sage: R._first_ngens??
Type: builtin_function_or_method
Base Class:
String Form:
Namespace: Interactive
Source:
def _first_ngens(self, n):
v = self.gens()
return v[:n]

So the equivalent code is this:

sage: R.gens()
(a, b, c)
sage: R.gens()[:3]
(a, b, c)

Cool, why not. So what is the R.gens() doing?

sage: R.gens?
Type: builtin_function_or_method
Base Class:
String Form:
Namespace: Interactive
Docstring:

Return the tuple of variables in self.

EXAMPLES:
sage: P. = QQ[]
sage: P.gens()
(x, y, z)

sage: P = MPolynomialRing(QQ,10,'x')
sage: P.gens()
(x0, x1, x2, x3, x4, x5, x6, x7, x8, x9)

sage: P. = MPolynomialRing(QQ,2) # weird names
sage: P.gens()
(SAGE, SINGULAR)

Ah, that's the answer we want. :)

Tuesday, December 25, 2007

Some timings of startups

I did some tests (I always start the program and then immediatelly hit ctrl-D):

$ time ipython
[...]

real 0m0.305s
user 0m0.112s
sys 0m0.024s

$ time bin/isympy
[...]

real 0m0.535s
user 0m0.200s
sys 0m0.040s

$ time ./sage
[...]

real 0m1.398s
user 0m0.916s
sys 0m0.208s

I did that repeatedly, so this is the usual time I get on my laptop. I think ipython is quite slow, because:

In [1]: %time import sympy
CPU times: user 0.08 s, sys: 0.01 s, total: 0.09 s
Wall time: 0.09

or if I import it in the python interpreter, it's immediate.

The 0.3s is acceptable, 0.5s is quite a lot for me, and 1.4s is a lot. I am maybe too demanding, but it just annoys me to wait for 1.4s on the import time. I just get some idea to try, so I fire up Sage or SymPy and try it. I don't want to wait for 1.4s.

This can be fixed, one just needs to find modules that slow things down, and late import them, or just fix them in some way. And usually, one does it several times -- I think Sage did that twice, SymPy also I think did that twice already.

Saturday, December 22, 2007

People still (and will) use windows

I sometimes look at the download statistics of SymPy and the funny thing is that people download more the Windows installer, rather than the multiplatform source tarball. Part of this can be due to the fact, that SymPy is in Debian, Ubuntu and Sage, but still.

So here is how to make a windows installer in Debian:

$ ./setup.py bdist_wininst
[...]

and a file dist/sympy-0.5.9-hg.win32.exe is created. I only have regular python tools installed. So we just upload the exe to our site and that's it -- kind of scary, since I don't have means to test it.

But one time I tried that on one of brother's computers with windows and not only it worked, but also the 3D plotting using pyglet worked out of the box! If you are curious how it looks like:



So it has it's advantages to develop pure Python programs - they really run everywhere.

I think installing things like Python and playing with SymPy on windows must be pain (maybe people just download SymPy for windows and then run away with disgust). But just the fact that they try means, that people doing science do use windows a lot. And I don't think this is going to change any time soon. I also fully agree with Michael Abshoff's recent post about this issue.

Wednesday, December 5, 2007

Mérida - wrap up

From November 28 till December 2, 2007, I attended a Debian QA and release teams work session in Extremadura, which is an autonomous community of western Spain, that managed to install Debian on 90000 computers in every school (technically a Debian based distribution called gnuLinEx) and Junta de Extremadura also sponsored this meeting.

I took photos of all participants, see my first, second and third posts. When I arrived at the Madrid airport where we first met, I set myself a goal to remember all names and faces, so I used my blog to help me and I think I succeeded in the end. :)

See also our wiki page that we used prior and during the meeting.

So what did we do besides throwing candies?


Lucas Nussbaum will send a summary email soon about the meeting, so I'll just speak for myself:


I worked with Gonéri on svnbuildstat, that is a service for building packages and show statistics about lintian/linda/piuparts checks. It for example contains all packages of the Debian Python Modules Team (that I am a member of) and many others. We discussed and started to work on how to create robust buildbots, that can be installed as a regular Debian package with zero (if possible) configuration, so that many people can just install them without pain, thus providing a huge scalability to the project.

We wrote a preliminary patch to pbuilder for killing the build if it exceeds given memory/disk usage. I had to learn the internals of pbuilder and I lost quite some time squashing some stupid bug I caused while writing the patch.

I spent most of my time with svnbuildstat, mostly learning and discussing things. This will be important for the future, but to also have some real results, I also fixed some packages I comaintain:

Together with Kumar Appaiah we fixed the python-numpy package and I had it uploaded, then I learned how to work with quilt instead of dpatch to handle patches in Debian packages, thanks to Holger's webpage, that contains a nice tutorial. Then I switched from dpatch to quilt in python-scipy and backported a patch from upstream svn to fix a segfault bug and had the package uploaded.

Then I finished the Cython package and had it uploaded. Cython is a marvelous package to speed up Python programs and interface C/C++ programs. I greatly recommend to try. If you don't like it, you can try some of at least 10 other ways to wrap C code in Python. I also used quilt in there to backport a patch from the upstream Mercurial repository to implement parsing @classmethods. Quilt is really a pleasure to work with.

Impressions from the meeting

I've been using Debian since 2001 as my only operating system on all of my computers, so I am not a complete beginner. But it never occured to me I could get involved in Debian more than a user and an occasional bug reporter. What a mistake.

I started packaging new things and fixing packages that I need for my work and that didn't work. This got me involved quite a bit in Debian. But in Mérida it was the first time I could dring a beer (well, especially wine) with Debian Developers and I found out they are really cool people. They are all very skilled. Also something, that I love about Debian, is that the people involved in it share two common features, that are very important for them - respect to democracy and personal freedom. When I think about it, those are probably the first two items on my presonal list of values.

Of course, everytime there is a group of 1000+ people, there are good and bad people, more and less skilled, but important is the overall atmosphere - and that is as I described. I think Debian is truly unique. There is Gentoo, that has maybe 40 (?) active developers. There is Ubuntu, that has maybe 100 (?) developers, but it's basically a comercial distribution and there is not so many interesting work for non employees of Canonical. There is opensuse and fedora, where I am not sure about the numbers. The atmosphere in Debian can change in the future, one never knows, but as of the end of 2007, I think it's very cool to get involved.

Maybe it's not for everyone, but it's the right place for me.

Tuesday, December 4, 2007

Mérida - remaining photo

This post is dedicated to dato (blog), whose picture I forgot to take (thanks Cyril Brulebois for taking this one).

From left to right: Gonéri Le Bouder, me, Lars Wirzenius, Holger Levsen and Adeodato "dato" Simó!

Saturday, December 1, 2007

Merida

We spent the whole Friday hacking, Lucas will send a summary email soon. More people joined.

Ana Guerrero:


Amaya Rodrigo:


Mark Purcell:


Kilian Krause and César Gómez Martín


and Holger Levsen with red hair:


On Saturday we were again working and in the late afternoon we took a walk in Mérida, visited the famous ancient Roman monuments. Now it's 3:30am and we will soon go to Madrid and back home.