Planet SymPy

May 14, 2008

Jason Gedge

Summer Work

Well my research job started last week and I'm looking forward to some results, although they won't come till later on this summer. I will be researching two areas:
  1. Image stitching
  2. Image [feature] matching
The basic idea is to take a set of photos and generate multiple panoramas from them appropriately. After this I will be estimating camera poses and with some form of a transform allow the user to travel from one panorama to the next. This should end up being a fully automated process so that one can easily create "virtual tours". It is sort of a combination of different software (such as Microsoft PhotoSynth and Quicktime VR).

I'm hoping by the end of the summer I'll have most of the concepts straightened out and some tools/software created to generate these tours and display them. Depending on the results I may carry on with this project for my honours thesis. Such a thesis would be an extension on this work and would probably include improved algorithms along with efficiency improvements, such as utilizing the GPU for tour generation.

Whenever I start getting results, I'll make another post. Once I understand the domain better I'll probably even post some technical information (for those of us who enjoy learning)!

by Jason G (noreply@blogger.com) at May 14, 2008 06:45 PM

Saroj Adhikari

GSoC withdrawal

This is going to be weird. I came to this decision after a long long thought process which has been invading my mind for some days now.

Before applying to GSoC, I had applied to a bunch of physics REUs (Research Experience for Undergraduates); almost all rejected but one I was awaiting the decision, which eventually offered me a position. Both of these summer opportunities are great and since I do not want to miss any of these, I had a very hard time deciding what to do.

So, now, I decide to go officially with physics REU because I can still continue contributing to SymPy and work on my project without officially enrolling on GSoC, during my free time. I should confess here that for about a week, I was also thinking to do both SoC and REU; I wish I could.

Once again, I apologize for this late withdrawal.

This is what I emailed to my mentor Michael, Leslie, and Ondrej about my GSoC withdrawal. I also apologize to other SymPy applicants through this post.

by Saroj (noreply@blogger.com) at May 14, 2008 02:33 PM

May 13, 2008

Jason Gedge

What's up?

I've been quite inactive as of late. One reason being that I was truly inactive, the other being that my blog was locked temporarily because it was tagged as a "spam" blog!

Nevertheless, I'm back in action and I hope to start posting some things on one of my new interests: ray tracing. I have a ray tracer up and running and things are going pretty good so far. Current features include:
  • Standard reflection + refraction (no special laws yet like, for example, Fresnel)
  • Several geometries: triangles, planes, spheres, boxes, quads, and cylinders
  • Anti-aliasing
  • Lighting: area, point, directional
  • Diffuse reflections
And some plans for the near future:
  • Support for loading several common model files
  • A space partitioning scheme to reduce rendering time
  • Transformed geometry (i.e., geometry + transformation matrix)
  • Support for instancing (multiple instances of the same geometry)
  • Python bindings to control/define scenes and the ray tracer
  • Photon mapping
  • Radiosity rendering
But some of those will take time. The list is sorted in order of highest priority at the top, to the lowest priority at the bottom. That's about it. And just to whet your appetite, I will leave you with my most current rendering:

by Jason G (noreply@blogger.com) at May 13, 2008 01:20 PM

May 01, 2008

Saroj Adhikari

Google Summer of Code 2008

I am starting this blog to document my Summer of Code proceedings.

My application SymPy: Improving the Plotting Module and a Pure Python TeX Engine was accepted under Python Software Foundation. I would like to thank Google, PSF, all SymPy people (especially Ondrej), and my mentor for the Summer of Code Michael Droettboom.

I will do my best to complete my project and I will continue contributing to SymPy even after that.

This blog will probably go beyond summer of code and this summer. 

by Saroj (noreply@blogger.com) at May 01, 2008 04:57 PM

April 26, 2008

Official SymPy blog

SymPy 0.5.14 released

SymPy 0.5.14 has been released on April 26, 2008. It is available at

http://code.google.com/p/sympy/

About SymPy

SymPy is a Python library for symbolic mathematics. It aims to become a full-featured computer algebra system (CAS) while keeping the code as simple as possible in order to be comprehensible and easily extensible. SymPy is written entirely in Python and does not require any external libraries.



user-visible changes:
  • SymPy is now 25% faster on average compared to the previous release (see below)
  • Documentation was improved a lot (commit 1, 2, 3). See http://docs.sympy.org/
  • rsolve_poly & rsolve_hyper fixed (commit 1, 2)
  • subs and subs_dict unified to .subs() (commit 1, 2)
  • faster and more robust polynomials module (commit 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, ..., look into the hg history)
  • improved Matrix.det(), implemented Berkowitz algorithm (commit 1, 2)
  • improved isympy (interactive shell for SymPy) (commit 1, 2, 3, 4, 5, 6)
  • pretty-printing improved (commit 1, 2, 3)
  • Rel, Eq, Ne, Lt, Le, Gt, Ge implemented (commit 1)
  • Limit class represents unevaluated limits now (commit 1)
  • Bailey-Borwein-Plouffe algorithm (finds the nth hexidecimal digit of pi without calculating the previous digits) implemented (commit 1)
  • solver for transcendental equations added (commit 1)
  • .nseries() methods implemented (more robust/faster than .oseries) (commit 1, 2)
  • multivariate Lambdas implemented (commit 1)


changes that affected speed:
  • __eq__/__ne__/__nonzero__ returns True/False directly so dict lookups are not expensive anymore (commit 1, 2, 3, 4, 5, 6)
  • sum(x**i/i,i=1..400) is now 4.8x faster (commit 1, 2, 3, 4, 5)
  • isinstance(term, C.Mul) was replaced by term.is_Mul and similarly for other basic classes (commit 1, 2)



Plus a lot of smaller bugfixes, you can browse our Mercurial history for details.


This release contains patches from 15 developers, which is so far the highest number of people/release (33 people have sent patches to SymPy so far, see the list of contributors):

  • Mateusz Paprocki
  • Fredrik Johansson
  • James Aspnes
  • Friedrich Hagedorn
  • Pan Peng
  • Abderrahim Kitouni
  • Nimish Telang
  • Jurjen N.E. Bos
  • Elrond der Elbenfuerst
  • Rizgar Mella
  • Felix Kaiser
  • Roberto Nobrega
  • David Roberts
  • Ondřej Čertík
  • Kirill Smelkov



If you'd like to contribute too, you can browse our open issues here:

http://code.google.com/p/sympy/issues/list

and then we suggest you to read the SymPy Patches Tutorial, that explains how to work with Mercurial effectively and how to create and send patches easily (not only) to SymPy. Any contribution is welcome, be it documentation, code, or just hanging out on our mailinglist or IRC (#sympy at freenode).

by Ondřej Čertík (noreply@blogger.com) at April 26, 2008 04:13 PM

March 25, 2008

Ondřej Čertík

SymPy accepts Google Summer of Code applications

SymPy is a pure Python library for symbolic mathematics. Last year SymPy had 5 excellent students and this year we are accepting students again.
Why should you apply? And why to SymPy?

Well, let me give you some reasons:

  • First of all, it's fun. To get some idea, read the GSoC2007 SymPy page, where you can find out what the last year students did and especially read their reports, where they describe their impressions from the summer, how they tackled problems and their overall conclusions.
  • It's not just about coding, we enjoy the social part too. There is a great community around numpy, scipy, ipython, matplotlib, Sage and similar tools and if you do scientific computing with Python, you gain a lot just being part of it, because you learn new things from the others.
  • I currently live in Prague (most people say it's a beautiful city, but I actually like Los Angeles, or the Bay Area:), if there are enough interested people, we can make a coding sprint here (plus of course some sightseeing+pubs). Anyone with a good commit history is welcome to stay at my apartment. :)
  • You earn $4500, some of which I suggest to spend on travelling to conferences/workshops, here are some tips: SciPy2008 (see also SciPy2007), EuroSciPy2008, Sage Days (you can read my impressions from SD6 and SD8), watch the numpy/scipy mailinglists for announcement of other meetings.


Read also the current status and motivation of SymPy and it's relation to Sage. If you want to apply, all the necessary information is on our wiki page.

Nevertheless, if you decide SymPy is not for you, but still you'd like to do GSoC project in a similar area, there are other good options too - one is SciPy/NumPy, the other is Sage. Unfortunately Sage was not accepted as a mentorship organization, but it has several good projects too, some of which you can do for example under the umbrella of the Python Software Foundation.

One of them is improving the Sage notebook. If you've never seen that - download Sage, start it (./sage), type "notebook()" and a nice Mathematica like notebook will popup in the browser. It allows collaborative editing ala Google Docs and many other things. If you'd like to work on it, reply to the email on sage-devel.

by Ondřej Čertík (noreply@blogger.com) at March 25, 2008 06:23 PM

March 06, 2008

Official SymPy blog

SymPy 0.5.13 released

It took us a little longer to release, so we have a longer changelog this time:

user-visible changes:
  • SymPy is now 2x faster in average compared to the previous release (see below)
  • integrate() can handle most of the basic integrals now (commit)
  • interactive experience with isympy was improved through adding support for [], () and {} to pretty-printer, and switching to it as the default ipython printer (commit 1, 2, 3)
  • new trim() function to map all non-atomic expressions, ie. functions, derivatives and more complex objects, to symbols and remove common factors from numerator and denominator. also cancel() was improved (commit 1, 2)
  • .expand() for noncommutative symbols fixed (commit)
  • bug in (x+y+sin(x)).as_independent() fixed (commit)
  • .subs_dict() improved (commit)
  • support for plotting geometry objects added (commit)
  • bug in .tangent_line() of ellipse fixed (commit 1, 2)
  • new atan2 function and assotiated fixes for .arg() and expanding rational powers (commit 1, 2, 3)
  • new .coeff() method for returning coefficient of a poly (commit 1, 2, 3)
  • pretty-printer now uses unicode by default (commit)
  • recognition of geometric sums were generalized (commit)
  • .is_positive and .is_negative now fallback to evalf() when appropriate (commit)
  • as the result oo*(pi-1) now correctly simplifies to oo (commit)
  • support for objects which provide __int__ method was added (commit)
  • we finally started SymPy User's Guide (commit 1, 2)
changes that affected speed:
  • first patches with 25% speedup (commit 1, 2, 3, 4)
  • Basic.cos et. al. removed, use C.cos instead (commit 1, 2, 3, 4)
  • sympy.core now uses direct imports (commit 1, 2)
  • sympifyit decorator (commit 1, 2, 3, 4, 5, 6, 7)
  • speedup Integers creation and arithmetic (commit 1, 2)
  • speedup unary operations for singleton numbers (commit)
  • remove silly slowdowns from fast-path of mul and div (commit 1, 2)
  • significant speedup was achieved by reusing dummy variables (commit 1, 2, 3, 4)
  • is_dummy is not an assumption anymore (commit 1, 2)
  • Symbols & Wilds are cached (commit 1, 2, 3)
  • ((2+3*I)**1000).expand() is now at least 100x faster (commit)
  • .expand() was made faster for cases where an expression is already expanded (commit)
  • rational powers of integers are now computed more efficiently (commit)
  • unknown assumptions are now cached as well as known assumptions (commit)
general cleanup:
  • BasicMeths merged into Basic (commit 1, 2, 3, 4, 5, 6, 7, 8)
  • cache subsystem was cleaned up -- now it supports only immutable objects (commit 1, 2, 3, 4, 5)


The following people have contributed to this release:
  • Kirill Smelkov
  • Mateusz Paprocki
  • Saroj Adhikari
  • Fredrik Johansson
  • Jaroslaw Tworek
  • Robert Kern
  • Pauli Virtanen
  • Ondřej Čertík

by Ondřej Čertík (noreply@blogger.com) at March 06, 2008 03:07 PM

Ondřej Čertík

Sage Days 8

Between February 29 and March 4, 2008 I attended the Sage Days 8, hosted at the Enthought headquarters in Austin, Texas. This was my 5th time in the USA and it was a marvelous experience, as with all my visits in the states.

As usual, I had some adventures in Atlanta, that interested readers can find at the end of this post. Anyway, on the Austin's airport I met Peter and his wife Crystal, Fernando, Benjamin, Jarrod, Eric and Clement. We went to have a dinner and then me and Clement were staying at Peter's house:



You can see the neighbor's cat and Peter's dog Trinity behind the window. The next day we went to Enthought, that was providing us with a breakfast and a lunch each day - and it was delicious. After the breakfast, we gathered in the room and introduced ourselves. Enthought rents 3/4 of the 21th floor in the Bank of America building, so when I looked left I saw:



When I looked behind I saw:



and in front of me, I saw all the participants (I took photos of all participants together with names). As you can see, there were really good people in there, like Travis (creator of NumPy), William (main author of Sage), Eric (CEO of Enthought), Fernando (author of IPython), Jarrod (the release manager of SciPy), Michael (the release manager of Sage) etc. See also the Fernando's welcome speech and the video of each of us introducting himself.

The views from the windows are terrific. I enjoyed working on each of the 4 sides of the skyscraper with completely different scenery, or when the sun is going down, that's also very cool.

We spent the whole Friday doing presentations, some of which you can find here. Then we went to Eric's house to have a big dinner together.

On Saturday, Sunday and Monday we were all hacking on many different things. I joined Fernando, Benjamin, Brian and Stefan on ipython1, Travis was implementing a new type (gmp integer) in NumPy, William wrote a manipulate command in Sage, Eric did the same in Traits, Gary and Michael implemented parallel testing of Sage, ...

On Tuesday we had final status reports and people left in the afternoon. In the evening we went with Clement to have a dinner and then we visited some bars on the 6th street, having a beer in each.

On Wednesday I visited John and Roy from the Computational Fluid Dynamics Lab at the University of Texas, Austin, who wrote the libMesh library, that I extensively used and also created a Debian package of. It was very influential to see the libMesh "from behind", also John and Roy are cool people (not mentioning the Debian tradition of having good relations with upstream:). Then I visited some professors at the same campus, after which I went into the Capitol and then I took the bus to the Barton Creek Square Mall to buy some ipods and jeans, so that I can say I have jeans from Texas. BTW, the ipod works excellent in Debian - I plugged it in and it just shows on my Gnome desktop. It's true that naively dragging mp3 files on it didn't make it play, but these instructions made it work.

On Thursday I fixed the remaining release blockers in SymPy and made a new release. In the evening, I am going to meet Aswin, he also uses SciPy and also is a friend of Kumar, who is now maintaining python-numpy and python-scipy Debian packages with me (Kumar also knows Prabhu, the author of Mayavi2 hosted at Enthought, so it's all connected).

Anyway, the whole workshop was an excellent experience for me. I learned a lot of new things and being able to speak with people who wrote tools that I use almost everyday is important. We also extensively discussed the future of all the projects (Sage, SciPy, NumPy, IPython, Cython, SymPy). See my summarizing email to the SymPy mailinglist.

Another thing, that I find very interesting is that Microsoft is financing the windows port of Sage, that will make basically anything that uses Python/Cython/C/Fortran very easy to install on windows (just a spkg package in sage). I find it really cool that MS is not only supporting but even financing a truly opensource project.

Finally the promised adventure in Atlanta: we took off the Prague airport on February 28th with a 2 hours delay (due to some paperwork as we were told by the captain). As I had 3 hours in Atlanta for the connection to Austin and I had to go through immigration, it was clear that I'll miss it. But I was not surprised, last time I was flying through Atlanta, they canceled my flight to LA completely. We arrived in Atlanta an hour and a half before my departure, then I was waiting for about an hour at immigration, it was incredibly slow. When I had around 20 min to departure, I had to ask people standing in front of me if they let me in, they were very nice and did. I was leaving immigration 10 min to my departure, then I was running to get my luggage and myself through customs and screening, it was 5 min to my departure when I ran down to the display with departure times. Then I was sprinting like hell to the terminal D to only see the clerk doing some final paperwork with all the people already boarded and the jetway door shut. After a little persuading he let me in too, fortunately there was still one seat left, so I made it. You can imagine my pleasant surprise in Austin when I discovered, that my luggage made it too, considering that I handed it to the Atlanta's airport personnel exactly 10 min prior the departure.

by Ondřej Čertík (noreply@blogger.com) at March 06, 2008 02:51 PM

February 24, 2008

Official SymPy blog

SymPy 0.5.12 released

Quite a lot has changed, SymPy now works with NumPy arrays, you can convert back and forth using numpy.array() and sympy.Matrix(). The integration algorithm has been improved significantly, it can handle a lot of common integrals already (some examples).

More detailed changelog:
  • SymPy works with NumPy out of the box (commit 1, 2)
  • RootOf implemented (commit)
  • Lambda support works now (commit)
  • heuristic Risch method improved (commit 1, 2, 3, 4, 5)
  • cancel function implemented (commit)
  • sqrt(x) is now equivalent to x**(1/2) (commit 1, 2)
  • Derivative is now unevaluated (commit)
  • list2numpy() implemented (commit)
  • series expansion of hyperbolic functions fixed (commit)
  • sympify('lambda x: 2*x') works, plus other fixes (commit 1, 2, 3)
  • simple maxima parser implemented (commit)
  • sin(x)[0] idiom changed to sin(x).args[0] (commit)
  • sin(x).series(x, 5) idiom changed to sin(x).series(x, 0, 5) (commit)
  • caching refactored (commit 1, 2, 3, 4, 5, 6, 7, 8, 9)
  • 2D plots now don't rotate in 3D, but translate instead (commit)
  • many bug fixes, see the Mercurial history for details


The following people have contributed to this release:

  • Kirill Smelkov
  • Mateusz Paprocki
  • Fredrik Johansson
  • Jaroslaw Tworek
  • Saroj Adhikari
  • Andrej Tokarčík
  • David Marek
  • Or Dvory
  • Bernhard R. Link
  • Ondřej Čertík

by Ondřej Čertík (noreply@blogger.com) at February 24, 2008 04:21 PM

January 10, 2008

Official SymPy blog

public wiki

We created a public wiki:

http://wiki.sympy.org/

That everyone is encouraged to participate in. The google code wiki has a huge disadvantage, that only the members of the google code sympy project are allowed to edit it. wiki.sympy.org can be edited by anyone and we hope that it will encourage SymPy users to help improve SymPy documentation and take part in the project.

We plan to move all tutorials and docs in there. Only things, that are relevant to developers, like how to do a release, etc., will stay in the google wiki.

by Ondřej Čertík (noreply@blogger.com) at January 10, 2008 06:37 AM

January 07, 2008

Official SymPy blog

SymPy 0.5.11 released

This is just a bugfix release, that fixed a problem with pyglet in the 0.5.10.

Changes:

  • ./setup.py install installs pyglet correctly now (commit)
  • var("k") fixed (commit)
  • script for automatic testing of plotting in pure environment added (commit 1, 2)

by Ondřej Čertík (noreply@blogger.com) at January 07, 2008 07:04 AM

January 04, 2008

Official SymPy blog

SymPy 0.5.10 released

Changes:

  • view renamed to preview, pngview, pdfview, dviview added (commit)
  • latex printer was rewritten, preview uses builtin pyglet instead of pygame (commit 1, 2)
  • square root denesting implemented (commit)
  • parser of simple Mathematica expressions added (commit)
  • TeXmacs interface written (commit)
  • some integration fixes (commit 1, 2, 3)
  • line width in 2D plotting can be specified (commit)
  • README was updated (commit 1, 2)
  • pyglet and mpmath were updated and moved to sympy/thirdparty (commit 1, 2, 3, 4, 5)
  • all sys.path hacks were moved to just 2 places - pyglet and examples (commit 1, 2, 3)
  • SymPy objects should work in numpy arrays now (commit 1, 2)
  • hand written sympify() parser was rewritten and simplified using Python AST (commit)
See also Changes.

by Ondřej Čertík (noreply@blogger.com) at January 04, 2008 07:06 AM

January 03, 2008

Ondřej Čertík

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.

by Ondřej Čertík (noreply@blogger.com) at January 03, 2008 11:05 AM

December 25, 2007

Ondřej Čertík

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.

by Ondřej Čertík (noreply@blogger.com) at December 25, 2007 05:05 PM

December 22, 2007

Ondřej Čertík

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.

by Ondřej Čertík (noreply@blogger.com) at December 22, 2007 09:52 AM

Official SymPy blog

Google Highly Open Participation Contest 2007-8

SymPy participates in the Google Highly Open Participation Contest 2007-8, if you are a high school student, go for it. Last chance to claim a task is January 22!

by Ondřej Čertík (noreply@blogger.com) at December 22, 2007 04:59 AM

December 21, 2007

Official SymPy blog

SymPy 0.5.9 released

Changes:

* Differential solvers were polished
* isympy now predefines f as a function
* Matrix printing improved
* Printing internals were documented

See Changes for more info.

by Ondřej Čertík (noreply@blogger.com) at December 21, 2007 04:22 PM

December 17, 2007

Official SymPy blog

Blogs

We'll use this blog for official announcements only from now on. Personal opinions will be on people's blogs, all of which are synchronised at planet.sympy.org.

We have implemented some new cool features in the Mercurial repository, so a new release will soon come out.

If you are curious, you can always find what code has been committed here:

http://hg.sympy.org/sympy

by Ondřej Čertík (noreply@blogger.com) at December 17, 2007 01:40 AM

November 16, 2007

Official SymPy blog

SAGE Days 6 wrap-up

SD6 is over, here are my final notes about the conference:

http://ondrejcertik.blogspot.com/2007/11/sage-days-6.html

SAGE is an exciting project. On SD6 I managed with the help of other SAGE developers to make SymPy play more nicely (after my patches are accepted, you will be able to freely mix SymPy and SAGE expressions and it will work).

by Ondřej Čertík (noreply@blogger.com) at November 16, 2007 08:33 AM

November 10, 2007

Official SymPy blog

SAGE Days 6

I am at SAGE Days 6, you can watch the progress on my blog.

by Ondřej Čertík (noreply@blogger.com) at November 10, 2007 03:51 PM

November 05, 2007

Official SymPy blog

Unicode printing

One very cool feature of SymPy, that Kirill recently committed, is a unicode printing. Start isympy and see for yourself:

$ bin/isympy
Python 2.4.4 console for SymPy 0.5.6-hg. These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z = symbols('xyz')
>>> k, m, n = symbols('kmn', integer=True)
>>> Basic.set_repr_level(2) # 2D output
>>> pprint_try_use_unicode() # try to setup unicode pprint


In [1]: sqrt((sqrt(x+1))+1)
Out[1]:
⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽
╱ ⎽⎽⎽⎽⎽⎽⎽
╲╱ 1 + ╲╱ 1 + x

In [2]: th, ph = Symbol("theta"), Symbol("phi")

In [3]: Integral(sin(th)/cos(ph), (th,0,pi), (ph, 0, 2*pi))
Out[3]:
2*π π
⌠ ⌠
⎮ ⎮ sin(θ)
⎮ ⎮ ────── dθ dφ
⎮ ⎮ cos(φ)
⌡ ⌡
0 0


More examples can be found on the wiki page PrettyPrinting.

Fredrik also played with the idea of ascii plotting. I am a huge fan of a terminal, because this is very versatile - I can use it over ssh, I can copy & paste, save the results to a file etc. I can easily script it. And the most important - I just prefer to write my programs in vim and execute them as a simple script. Any GUI, even the best one (or in the browser), is not giving me such a power. So I like the idea of making the terminal experience as pleasant as possible.

by Ondřej Čertík (noreply@blogger.com) at November 05, 2007 07:11 AM

October 20, 2007

Official SymPy blog

Google stopped showing sympy's webpage

This was written on October 18:

I was used to finding the sympy webpage by googling "sympy". When I did that today, the http://code.google.com/p/sympy/ isn't shown in the results (instead the sympy's freshmeat page, or sympy's Debian page is shown). I tried yahoo and msn, and sympy webpage always shows as the first or second link in there (as it used to in google).

It's interesting and it reminds me how important it is to have several search engines and do not depend on only one, so that one can switch to a different one if the first one is apparently not working correctly.

Update (October 21): Google fixed the problem and sympy's webpage is the first link again.

by Ondřej Čertík (noreply@blogger.com) at October 20, 2007 04:16 PM

October 06, 2007

Official SymPy blog

New SymPy release and beyond

Things have been a little quiet on planet.sympy.org, but the less we talk, the more we work. :) With 468 svn commits in September and 89 in October (till October 5) [1], we are working hard. We made a new release yesterday, where we simplified the internals of SymPy a lot and I also fixed the limits finally.

What remains is to fix some problematic series expansion (that also makes some limits fail, but the limit algorithm itself seems to be bugfree so far) and then I'll concentrate on making SymPy more SAGE friendly. SymPy will remain as a small Python library, doing one job and doing it well. But when used inside SAGE, it will know how to convert expressions to and from SAGE, so that people can use SymPy from SAGE ideally the same way as they can now use Maxima.

I think there are already things that SymPy is better than Maxima in terms of pure functionality, like pattern matching, limits (after I'll fix the series code, hopefully in the next release). There are however things, where Maxima is better, it still has a lot more features (more integrals, more differential equations, tensors) and it's faster. Except for the speed, all the other things are really easy to implement in SymPy. The main advantage of SymPy over Maxima is that it is written in Python and it is really easy to implement new features in it.

A huge deficiency of Maxima is that it uses lisp and it's own language for doing the calculations. SAGE fixes this problem (but it doesn't support all features of Maxima yet, for example substitutions of subexpressions, that is needed in the Gruntz limit algorithm), so SAGE is in fact like a platform, a common ground for all these mathematical libraries, allowing them to be used using the same interface and thus one can compare the real features (and speed), instead of some possibly good libraries being disqualified just because they use an unfortunate language/interface (for whatever historical reasons), which is great and I want SymPy to compete in the same ground as well, because only this will tell, if it is good or not and which things to improve. Competition is very healthy, it's always good to have a choice.

[1] http://groups.google.com/group/sympy-commits/about

by Ondřej Čertík (noreply@blogger.com) at October 06, 2007 04:49 AM

September 08, 2007

Jason Gedge

Problem Solving & TopCoder

There's no argument in the fact that problem solving is a vital part of our intellectual functions and something that helps keep the mind sharp. Even something as simple as a logic puzzle, such as Sudoku, can help influence our mental well being.

One of my favorite activities to keep my problem solving sharp is participating in competitions at TopCoder, a company dedicated towards providing online competitions in algorithms, component design, and so on. I personally participate in the algorithms competitions. I highly recommend for everyone who can program in Java/C/C++/C# to participate from time to time, even if you feel like you have no chance. If anything, it will keep your mind and your programming skills sharp. The scoring system is based on time, so you learn to gear your mind towards immediately being able to categorize problems and then developing a solution.

Currently, the TopCoder Collegiate Challenge is underway, and I have surpassed my expectations this year. Last year I only made it past the qualification rounds, but this year I have matched that and even made it past two more online rounds. If I am able to make it pass two more online rounds then I would have the opportunity to go to Orland and compete against 47 of the world's finest. I don't expect to make it this far from a realistic point-of-view, not a pessimistic one. Nevertheless, I will definitely try my best to make it even further.

Don't think it's just your physical self that needs to be kept "in shape"!

by Jason G (noreply@blogger.com) at September 08, 2007 07:44 PM

August 23, 2007

Official SymPy blog

Summer of Code summary

Summer of Code deadline has passed, now is the evaluation period, however I think I can already say, that it was a success, since all projects were finished.

Here is some statistics

ondra@fuji:~/sympy$ svn log |grep ondrej | wc -l
64
ondra@fuji:~/sympy$ svn log |grep inferno | wc -l
48
ondra@fuji:~/sympy$ svn log |grep fredrik | wc -l
39
ondra@fuji:~/sympy$ svn log |grep mattpap | wc -l
39
ondra@fuji:~/sympy$ svn log |grep lethargo | wc -l
29
ondra@fuji:~/sympy$ svn log |grep brian | wc -l
44
ondra@fuji:~/sympy$ svn log |grep pearu | wc -l
14
ondra@fuji:~/sympy$ svn log |grep Chris.Wu | wc -l
3

But one shouldn't draw too many conclusions from that, because people like me prefer to commit very often and people like Mateusz prefer to work offline and then commit when they have a bigger debugged chunk of code ready. Also the logs start with the sympy-merge branch, that was created on July 19. So everything before that (Pearu's contribution and everything that happened before the merging) is not present in the statistics. Another statistics are the Issues changes mailinglist, I did a search for all the names:

Ondrej 974
Fabian 335
Jason 227
Fredrik 151
Robert 104
Mateusz 88
Pearu 76
Brian 65
Chris 28


The above table are all issues from the beginning of the project. Fabian was the most active and then Jason is catching up recently. :)

One should have in mind the old joke though: There are three kinds of lies: lies, damned lies, and statistics.

I'd like to thank all of you for bringing SymPy this far. Basically all major things are working, now comes the phase of polishing things and reducing the number of open Issues and improving doctests, so that SymPy is very easy to use for newcomers.

by Ondřej Čertík (noreply@blogger.com) at August 23, 2007 12:36 PM

August 21, 2007

Brian Jorgensen

Final SoC Report

I've written a final SoC report and posted it in the SymPy wiki. Click here to read it. You will also find similar documents from other SymPy SoC participants here.

by Brian Jorgensen (noreply@blogger.com) at August 21, 2007 08:34 AM

August 19, 2007

Brian Jorgensen

SymPy Plotting: Custom Colors Tutorial

One of the most-requested features for SymPy plotting has been the ability to use custom color schemes. I've now implemented this with the syntax described in this post. The upshot is that you can now use any color scheme expressible as a function of x, y, z, u, and/or v. I'll give some concrete examples to get you started, and then you can do something cool with it. We'll start by setting up a plot:

>>> from sympy import symbols, Plot
>>> x,y,z,u,v = symbols('xyzuv')
>>> p = Plot(axes='none')

Now let's plot a saddle and color it by the magnitude of its gradient:

>>> fz = x**2-y**2
>>> Fx, Fy, Fz = fz.diff(x), fz.diff(y), 0
>>> p[1] = fz, 'style=solid'
>>> p[1].color = (Fx**2 + Fy**2 + Fz**2)**(0.5)



Remember that the algorithm for coloring works like this:
  1. Evaluate the color function(s) across the curve or surface.
  2. Find the minimum and maximum value of each component.
  3. Scale each component to the color gradient.
When not specified explicitly, the default color gradient is (r,g,b) = (0.4,0.4,0.4)->(0.9,0.9,0.9). In our case, everything is gray-scale because we have applied the default color gradient uniformly for each color component. When defining a color scheme in this way, you might want to supply a color gradient as well:

>>> p[1].color = (Fx**2 + Fy**2 + Fz**2)**(0.5),
................ (0.1,0.1,0.9), (0.9,0.1,0.1)



Next, let's try a color gradient with four steps:

>>> gradient = [ 0.0, (0.1,0.1,0.9), 0.3, (0.1,0.9,0.1),
................ 0.7, (0.9,0.9,0.1), 1.0, (1.0,0.0,0.0) ]
>>> p[1].color = (Fx**2 + Fy**2 + Fz**2)**(0.5), gradient



The other way to specify a color scheme is to give a separate function for each component r, g, b. With this syntax, the default color scheme is defined:

>>> p[1].color = z,y,x, (0.4,0.4,0.4), (0.9,0.9,0.9)



This maps z->red, y->green, and x->blue. In some cases, you might prefer to use the following alternative syntax:

>>> p[1].color = z,(0.4,0.9), y,(0.4,0.9), x,(0.4,0.9)

You can still use multi-step gradients with three-function color schemes. When somebody uses this to visualize something useful like curvature, I'd really like to hear about it.

by Brian Jorgensen (noreply@blogger.com) at August 19, 2007 06:29 PM

August 17, 2007

Jason Gedge

WADS 2007: Day #1

I figure I should write a blog entry now before I actually forget what happened yesterday!

Session 1 (Invited Talk)
The day started with a presentation whose major title was "Finding small holes" given by Jeff Erickson. It was a small talk on computational topology. It was surprisingly interesting and I was very intrigued overall by the area. The presentation placed an emphasis on a practical application of itself, which was surface reconstruction from potentially noisy data. The processes described showed how surface reconstruction could be done much better by taking advantage of computational topology. No actual algorithms were presented, but the ideas to establish the algorithms were.

For all of you with low G.P.A.s and thinking you wouldn't really be able to go beyond an undergraduate degree because of it, think again. Jeff Erickson cracked a joke about him being potentially one of the lower G.P.A.s amongst most professors with a 2.4/4.0. Doing well in undergraduate studies doesn't necessarily mean you're prepared to do research and be a professor, and similarly for not doing so well you could be perfect through graduate studies. If it's something you would really like to do, get out there and do it!

Session 2B
Nothing particularly interesting here besides for the last talk of the session. At first I thought it was going to be ridiculous because they themed their talk around Donald Duck and Uncle Scrooge, but it was actually quite interesting. The paper was called "The Stackelberg minimum spanning tree game." Their problem was to take a graph with some known and unknown edge weights, and fill in the edges with unknown weights such that the MST maximizes the sum of the values of the weights given to the unknown edges (sorry if I lost you there). One of the characters which presented this talk was Erik Demaine, a guy who came to Dalhousie at age 14 (I think that's right). He's 26 now and is a professor at MIT. The man is brilliant and quite comical also!

Session 3A
I probably should have went to session B here because this session was completely about graph drawing, which isn't a real particular interest of mine. It was cute, but that's about it.

Session 4A
A couple of interesting computational geometry presentations in this session: steiner trees, art gallery problems, and Delaunay triangulations. The last one was particularly interesting because of the reference to one of its practical applications, which is triangulating GIS data for 3D representation. I enjoyed these talks and I think it gives computational geometry a +1 on my list of things I would potentially do for my thesis/grad studies.

All in all, a great day!

by Jason G (noreply@blogger.com) at August 17, 2007 11:08 PM

WADS 2007: Days 2 & 3

D A Y 2

Session 5 (Invited Talk)
The day started off with an invited talk from Michael Langston from the University of Tennesse. The talk was on challenges in data analysis, and how sometimes you can have beautiful data to work with (good correlation, easy to work with, e.t.c.) and sometimes you can have tough data to work with. Although not particularly intriguing to myself (I'm still leaning towards a more theoretical area of research), the things he described are exactly the things I have experienced this summer with my cross-disciplinary work. All and all this was a good talk.

Session 6B
I missed the very first presentation in this session, but the second one was about Steiner trees again (everyone loves a Steiner tree). Not quite computing, because that is NP-Complete The session ended off with computing minimum-depth planar graph embeddings in quartic time. I enjoy graphs so I also enjoyed both of these talks!

Session 7A/7b
This session really wasn't doing it for me so I skipped it. I probably should have looked into the computational geometry presentations taking place in 7A though because it is starting to interest me somewhat.

Session 8B
I checked out some parameterized algorithms and kernelization algorithms. Nothing much to report besides the fact that it wasn't overly exciting! For those who are unfamiliar with kernelization, it's simply taking a problem instance and breaking it down into a smaller problem instance. From here we solve this smaller problem instance and from it we can extract the solution to the original problem. This is a common technique in creating algorithms that have good bounds for running times when analyzed parametrically.

Conference Dinner
Lobster! Nothing to say about the dinner besides for the fact that it was delicious. I joined a few of the PHD/Masters students - Aaron Lee from Carleton University, Cora Borradaile from Brown University and another guy from Calgary whose name escapes me at the moment - at the pub next to the Alexander Keith's brewery after for a couple of drinks. Good times!


D A Y 3

Session 9B
I sort of slept in through these...

Session 10B
Approximation of graph distances and shortest paths. Nothing particularly interesting to point out. I just enjoyed going to all of these graph theory talks because I'm not completely out of the blue when they mention crazy terms and the like.

Session 11B
To end off the conference were two sessions with four presentations. The first two involved suffix arrays. One of the grad students which I was fairly familiar with over the course of the past couple of days, Orgad Keller from Israel, done an excellent job at presenting his work (clarity, leaving out the overly complicated details, e.t.c.). The third talk was interesting; it dealt with streaming algorithms, those which receive input in a stream and use much less space (a fuzzy term, I know) than the actual amount of data. They proposed optimal (or maybe ti was near optimal) ways of solving the straggler identification problem. The last presentation I left for, because it was on TCP acknowledgment.

It was all over and Jason was saddened because he enjoyed his time at Dalhousie. But alas, all good things must come to and end. Back to the ol' grind now and in a couple of weeks school starts again. That's kind of exciting too. . . for now.

by Jason G (noreply@blogger.com) at August 17, 2007 11:02 PM

Chris Wu

Homestretch

Well, it seems we're in the homestretch of the SoC. I've been working on the FFLU decomp. Strangely I ran into some problems with modules not being installed properly but that finally regulated itself. Now I'm in the debugging phase and hopefully we're have a nice new fraction-free LU function.

I'm thinking that maybe I can squeeze in another Fraction-free inversion algorithm before next week...

*edit* Did I say inversion? I meant QR. So LU seems to be working but my QR one is still buggy. Oh well.

by cwu (noreply@blogger.com) at August 17, 2007 10:55 AM

August 15, 2007

Robert Schwarz

Faster arithmetic

A quick proof-of-concept finite field arithmetic and factorization was done. And the factorization over the finite field was reasonably fast, too, but the big prime reconstruction of the integer factors took some more time. Also, the resulting complete algorithm spent a lot of time with square-free factorization and other preparations. So, it was quite disappointing.

I still started over again, rewriting everything from scratch, now with tests. First the modular integers, with a custom class for each modulus. Then a generic light-weight polynomial class, for one variable, using Python dictionaries. Then some algorithms for polynomials with finite prime field coefficients, with complete factorization, which is also quite fast. Today, I've started to do some more integer polynomial arithmetic, with a gcd using a small primes modular approach. Should be one of the fastest, asymptotically, but probably needs some more tweaking in my implementations. I should now check in which cases these modular methods really are faster than working with rational numbers. All this won't be finished for the Summer of Code deadline, but can be further advanced, while I learn from my new bought book, "Modern Computer Algebra" (Gathen, Gerhard).

by Robert (noreply@blogger.com) at August 15, 2007 10:20 PM

Jason Gedge

GSoC + WADS 2007

The GSoC is nearing an end and there's not much to do really. I think once I get back home on Saturday I'm going to completely look over my Geometry module and see what little things I can do to fix it up (documentation, refactoring, etc).

The WADS 2007 conference reception was last night, and I met some interesting people from all over the world. The talks seem like they'll be fairly interesting overall, with some neat algorithms coming forth. For now I'm lacking something to blog about, but after the day is over and I have taken in some of these presentations, hopefully I'll have something interesting to write about! Note that I'm really liking this Dalhousie CS building over the CS department at MUN (i.e., our corner in the engineering building)

by Jason G (noreply@blogger.com) at August 15, 2007 06:10 AM

August 06, 2007

Robert Schwarz

Finite field factorization

I've finished updating the documentation and have started to move to next point, faster univariate factorization. This is done modular, with factorization over a finite field first. Fortunately, this point is a major part of Garthen and Gerhard's "Modern Computer Algebra" and thus well explained. Unfortunately, I've also read there:
Computaionally, (fast) polynomial factorization over a finite field is a much more advanced task than, for example, multiplication or even gcd computation. Before implementing a particular algorithm, one has to implement carefully many other routines for basic polynomial arithmetic, and designing a polynomial factorization package usually requires several woman-years.
I still hope to get something done this week, doesn't have to be optimal, then. Let's see if it is a lot if work to integrate modular (number) arithmetic, or if this is better done "inline".

by Robert (noreply@blogger.com) at August 06, 2007 10:01 AM

Nice literature

I forgot to mention some interesting books I use.


The most important one is Cox, Little, O'Shea: Ideals, Varieties, and Algorithms, which is an introduction to algebraic geometry (and commutative algebra) using a relatively constructive approach. That is, for most of the objects, there are not only proofs of existence (and uniqueness) but also algorithms (or at least references). I like it very much and have learned a lot so far; it was already a plan for the summer to revise some algebraic geometry, this fits very well.


Yesterday, I found two other books in the university, more directly concerned with (general) computer algebra:


Davenport, Siret, Tournier: Computer Algebra is a gentle introduction to different features of CAS, problems of data representation and some advanced algorithms. It could be a bit outdated on current implementation details and running times, but gives a quick overview and offers confrontation with problems otherwise ignored. (Already half-way through, not looking at all details)


Gathen, Gerhard: Modern Computer Algebra. Have not looked at it in detail, but looks really nice. One of the authors is an employee of MapleSoft and all methods are presented with a huge list of applications and examples. The content is presented with some information on historic development and stuffed with nice quotations (also in the original language, way to go!)

by Robert (noreply@blogger.com) at August 06, 2007 10:01 AM

Gröbner basis

Hey, my first blog entry so far.


My name is Robert Schwarz, I am a student of mathematics and computer science at the university of Heidelberg, Germany. For the next two weeks, I'll still be staying in Paris, where I passed a year with the Erasmus student exchange program.


The GSoC project is mostly about Gröbner bases and application. The euclidean algorithm for finding the gcd of 2 integers is one of the most popular examples in programming. The same thing also works for polynomials, but unfortunately only in one variable. The Gröbner basis is a set of polynomials that means to replace the role of the gcd.


When looking for the solutions of some system of polynomial equations f(x,y)=0, g(x,y)=0 it is clear that we could add more equations like f+g or f*h for some polynomial h without affecting the set of solutions. The set of such equations/polynomials is called the ideal I = , generated by f and g. If I want to know if some polynomial h(x,y) lies in I, I have to check if it can be written as h = a*f + b*g, that is, I try to divide h by f and g and look at the remainder, similar to the division with integers. Unfortunately, the choice of f,g and even their order can influence the remainder here, and it is the special characterization of the Gröbner bases, that it works with them. Hope, it is a bit clearer for those who have never heard of them.


The natural application lies of course in the solving of polynomial equations, with simplification of the equations and systematic elimination of variables. I think I've read somewhere that Gröbner bases are also needed for some other feature in sympy, I'm very interested in that, please tell me, so that I can start to help.


Yesterday, I finished the implementation of Buchberger's algorithm that computes Gröbner bases, so in theory, it is already available, even if it's still a bit slow.


Items on my to-do list: Speed improvements in groebner(), real-life examples for polynomial equations, elimination, factorization of polynomials, solving for (real) roots.

by Robert (noreply@blogger.com) at August 06, 2007 10:01 AM

August 04, 2007

Brian Jorgensen

New Plotting Examples

Try running examples/plotting.py in an interactive shell:



Other updates:
  • New experimental option 'use_lambda', which speeds up calculation by an order of magnitude (used in most of the examples, but doesn't work in every case yet).

  • Custom coloring can be applied in real-time:
    >>> p[1].color = 0.3,0.3,0.9 # any solid color
    >>> p[1].color = 'z,y,x' # lambda-based (fast)
    >>> p[1].color = x,y,z # sympy expr of x,y,z
    >>> p[1].color = u,v,u*v, (u,v) # or of parameters

  • Style property:
    >>> p[1].style = 'wireframe'
    >>> p[1].style = 'solid'
    >>> p[1].style = 'both'

  • Natural right-click translation in the plane of the screen.
  • Middle-click drag zoom.
  • Progress percentage displayed in window caption for long calculations.
  • z,c and numpad 1,3 rotate about the z axis regardless of camera angle.
  • x and numpad 5 reset view.

by Brian Jorgensen (noreply@blogger.com) at August 04, 2007 02:32 PM

August 03, 2007

Chris Wu

Typo

So for whatever reason I kept saying block LU, when I really meant fraction-free LU. Those are, of course, significantly different. I've been looking at this paper suggested by someone on the mailing list. It seems reasonably straight-forward (even the proof, which is surprising). I hope to finish the LU one soon and hopefully the QR one before the project finishes up.

by cwu (noreply@blogger.com) at August 03, 2007 11:48 AM

August 01, 2007

Robert Schwarz

Refactoring

Now, that the planned features are roughly present and we have moved to a new core with SymPy, there were some suggestions concerning the structure of the polynomials class. Most of them about the wrapper functions and code duplication, or the internals of the Polynomials class.
The Polynomials class is now immutable, as are all Basic-inherited classes, and tries hard to pretend to be a Basic object itself. We'll see how that will turn out. The wrapper classes essentially only extract the underlying sympy expression of the resulting Polynomial.

The plan for the next days includes:
Documentation (proper docstrings everywhere, algorithms explained)
Testing (more complete tests, with comments)

Then there are some more important features, like fast univariate integer factorization with finite fields. We'll have to look at some number theoretic python package, probably nzmath. Also, the equation system solver could tell us, when we are missing some solutions, that can't be computed but still exist.

by Robert (noreply@blogger.com) at August 01, 2007 04:17 PM

Official SymPy blog

Progress

The midterm evaluation is over and noone has failed, which is great. We started to work more closely on the #sympy at freenode IRC channel and SymPy has made a tremendous progress in past couple of weeks.

Our team was joined by Fredrik and Pearu, so now there are 8 regular contributors.

Fredrik wrote fast floating point arithmetic (in pure Python), faster than the python Decimal module. For example 10 000 decimal digits of Pi can be computed in 0.3s (and one million in less than 20 minutes) with SymPy - just run the pidigits.py in the examples directory. He among other things wrote a fast algorithm for calculating the gamma function, that in some cases is faster than in Mathematica.

We moved to the new core written by Pearu, that has made everything 10x faster or more. Basically, Pearu very cleverly designed everything so that no operation is made in vain. For example my little hobby relativity.py for calculating the Schwarzschild solution now calculates it in 1s (compared to more than 10s before)!


ondra@pc232:~/sympy/examples$ time python relativity.py
----------------------------------------
Christoffel symbols:
(1/2)*(FD(1,)(nu))(r)
(1/2)*(FD(1,)(nu))(r)

(1/2)*exp(-lam(r) + nu(r))*(FD(1,)(nu))(r)
(1/2)*(FD(1,)(lam))(r)
-r/exp(lam(r))
-r*sin(\theta)**2/exp(lam(r))

1/r
1/r
-sin(\theta)*cos(\theta)

1/sin(\theta)*cos(\theta)
1/sin(\theta)*cos(\theta)
1/r
1/r
----------------------------------------
Ricci tensor:
(1/2)*exp(-lam(r) + nu(r))*(FD(1, 1)(nu))(r) - 1/4*exp(-lam(r) + nu(r))*(FD(1,)(nu))(r)**2 + 1/r*exp(-lam(r) + nu(r))*(FD(1,)(nu))(r) + (1/2)*exp(-lam(r) + nu(r))*(-(FD(1,)(lam))(r) + (FD(1,)(nu))(r))*(FD(1,)(nu))(r) + (1/4)*exp(-lam(r) + nu(r))*(FD(1,)(lam))(r)*(FD(1,)(nu))(r)
-1/2*(FD(1, 1)(nu))(r) - 1/4*(FD(1,)(nu))(r)**2 + 1/r*(FD(1,)(lam))(r) + (1/4)*(FD(1,)(nu))(r)*(FD(1,)(lam))(r)
1 - 1/exp(lam(r)) + (1/2)*r/exp(lam(r))*(FD(1,)(lam))(r) - 1/2*r/exp(lam(r))*(FD(1,)(nu))(r)
sin(\theta)**2 - sin(\theta)**2/exp(lam(r)) + (1/2)*r*sin(\theta)**2/exp(lam(r))*(FD(1,)(lam))(r) - 1/2*r*sin(\theta)**2/exp(lam(r))*(FD(1,)(nu))(r)
----------------------------------------
solve the Einstein's equations:
lam(r) = -log(C1 + C2/r)
metric:
-C1 - C2/r 0 0 0
0 1/(C1 + C2/r) 0 0
0 0 r**2 0
0 0 0 r**2*sin(\theta)**2

real 0m1.092s
user 0m1.024s
sys 0m0.068s


Now we are in the phase of polishing things so that everything that used to work before works in the new core the same way (or better). Things that remain are printing and some limits are failing, everything else was already fixed.

One conclusion that can be made is, that Python on today computers is already fast enough to do something useful. Concerning the Pi, it's amusing to read some articles about people trying to calculate many digits, for example in one article, it took them a day to calculate one million digits on Pentium 90 (SymPy can do that in 20 minutes).

Also I found some page, that I cannot rediscover at the moment, that stated it took like 15 minutes to calculate the Schwarzschild with Maple on some older computer, so I am glad SymPy can do that in 1s (and actually the code is very slow, I am doing some calculations again and again, but I didn't find time yet to optimize it).

Anyway, the overall progress on SymPy looks very promising, I have no doubts that it will become a useful symbolic python library.

by Ondřej Čertík (noreply@blogger.com) at August 01, 2007 12:12 PM

July 31, 2007

Jason Gedge

End of July Update

Yet again, not much to update on. We're just beating away at the newly merged trunk and trying to fix up some issues before the 0.5 release. After that it'll just be more issue fixing. Considering this summer I've been taking a course on ordinary differential equations I may consider working on the facilities in SymPy that solve ODEs. Only one more month left to the summer but a lot has gotten done. I'm looking forward to it though, to see what state SymPy is in.

I'm also looking forward to school starting again. Ignoring work, this summer has been far from what I hoped it would be, but I'm probably partly to blame for that one. As for school, I'll be taking these courses:
  • Integration and Metric Spaces
  • Abstract Algebra
  • Combinatorial Analysis
  • Computer Graphics
  • Computational Complexity
I'm not looking forward to the first one, but the rest of them should be alright. Although it might be a bit of work, I expect Computer Graphics to be somewhat enjoyable, or as enjoyable as a school-related item can be! This year will be my last year of regular courses. I'm hoping to have another research award next Summer and then the Fall following it will be my honors thesis, and I have no idea what to do for that at the moment. So at the end of '08 I'll officially have my undergraduate degree from MUN: a joint major in Computer Science and Pure Mathematics.

by Jason G (noreply@blogger.com) at July 31, 2007 08:20 PM

July 29, 2007

Brian Jorgensen

Mathematicians are a Hungry Bunch

After yesterday's donut, today it's the Ding Dong Surface. Thanks to Alex who originally pointed me to the Kiss Surface (also named after junk food).

>>> p = Plot()
>>> p[1] = sqrt(1-y)*y, [x,0,2*pi,60], [y,-1,4,100], 'mode=cylindrical'




By the way, the torus shown yesterday can be plotted with:

>>> a,b = 1, 0.5
>>> p[2] = (a+b*cos(x))*cos(y), (a+b*cos(x))*sin(y), b*sin(x), [x,0,2*pi,20], [y,0,2*pi,20]


However, at writing there is a hang bug with trig functions in the new SymPy core, so this may or may not work when you try it.

by Brian Jorgensen (noreply@blogger.com) at July 29, 2007 05:05 AM

July 28, 2007

Brian Jorgensen

Just in time for the Simpson's Movie

After weeks of working mostly on the UI, I've spent three or four days reworking the calculation and display code. One result of that effort is support for parametric surfaces (x,y,z) = f(u,v).



Other things:
  • Now using GL display lists for faster rendering.
  • Support for wireframe, solid, and wireframe superimposed on solid rendering modes.
  • Smaller and more readable code base, since all surface and curve modes are now defined in terms of base parametric modes.
  • Preliminary support for custom coloring.
The new code is not hooked up to the Plot interface yet, so that is my goal for today. Hopefully you will see this in SVN by the end of the weekend, possibly tonight.

by Brian Jorgensen (noreply@blogger.com) at July 28, 2007 01:50 PM

July 27, 2007

Chris Wu

new core

So the new core was added by the SymPy crew and it seems much faster than before - at least according to my test cases. Mateusz tells me my module worked out-of-the-box, which is always nice to hear.

I managed to put together a nice (in my opinion) tutorial for my linear algebra pack. Though I will say the inability to have a table of contents of html anchors make it very hard to read and navigate the page. This is usually something quite easily done in a wiki.

The tutorial brought to light a couple of stupid mistakes and inconsistencies I had made so it was helpful in the debug process too.

I'm thinking block LU decomp next week, kind of an advanced thing as we come to the end of the project....

by cwu (noreply@blogger.com) at July 27, 2007 02:18 PM

July 24, 2007

Brian Jorgensen

SymPy Plot: First Screenshot of Labeled Axes



Good thing a picture is worth 10**3 words; I'm burnt out. Among other things, it took me way too long to figure out how to automatically position the axis labels after camera rotations. It still isn't perfect. I also need to add in grid lines. I need a break, I've been coding for a nearly unbroken week.

by Brian Jorgensen (noreply@blogger.com) at July 24, 2007 01:42 AM

July 21, 2007

Brian Jorgensen

Some thoughts on Plot coloring (with pics)

I want to support two modes for custom coloring schemes in SymPy Plot. The first mode will allow you to supply a single function and a color gradient (pseudo-code):

>>> p[1].color_function = z, (0.0,0.0,1.0), (1.0,0.0,0.0)



This calculates the min and max z-value across the surface, then interpolates from (min, blue) to (max, red). A variation on this syntax allows you to specify more nodes on the color gradient (the fourth component is the linear position [0,1] on the gradient):

>>> p[1].color_function = z, (0.0,0.0,1.0,0.0), (0.0,1.0,0.0,0.5), (1.0,0.0,0.0,1.0)



As you can see, this fades from (min, blue) to (mid, green) to (max, red).

The second coloring mode will allow you to map each RGB component to separate functions. This can be used to colormap three-space vectors, such as surface normals. The simplest usage of this would be:

>>> p[1].color_function = x, y, z

However, this would look pretty bad. In fact, all of the above images used a slightly different color range than was written in the accompanying pseudocode, but they were simplified for clarity. For example, the top one was really more like:

>>> p[1].color_function = z, (0.3,0.3,0.9), (0.9,0.3,0.3)



So in RGB-parametric mode, there will be an enhanced syntax for playing with the dynamic range independently of whatever functions you are using. The current default coloring could be represented with something like:

>>> p[1].color_function = x, (0.5,0.9), y, (0.5,0.9), z, (0.5,0.9)

by Brian Jorgensen (noreply@blogger.com) at July 21, 2007 10:59 PM

July 19, 2007

Jason Gedge

GSoC Update

It's been quite some time since I've had an update, so I felt like it is only due time to write a little something. There's been a lot of "off to the side" work. Preparations are currently under way to merge the research branch of SymPy with the actual trunk. This is a fair amount of work, and hence we're trying to coordinate everything so that things run smoothly. I have rewritten unit tests to be supported by py.test, but there is also the moving/updating of tests from the trunk and ensuring modules still work along with their tests. This whole merge will be my major task for the next little while.

I will be heading to WADS 2007, the Workshop on Algorithms and Data Structures, during mid August. It is a conference that [sort of] deals with my current research area so I was advised to go. I'm doing so for not only that reason, but also because I would like to see if the academia lifestyle is the one for me. The conference looks to have some interesting papers, and I bet there will be some interesting people to meet.

by Jason G (noreply@blogger.com) at July 19, 2007 08:51 PM

July 17, 2007

Brian Jorgensen

SymPy Plotting: Improved Usability with Pyglet

One of my personal goals for the SymPy plotting module is to introduce non-programmers to Python. This means it needs to have a familiar interface, and that it should "just work" straight out-of-the-box. With the latest code in SVN, you can now use Plot in Python 2.5 without installing any additional dependencies. In Python 2.4 and 2.3, the only external dependency is ctypes.

PyOpenGL is no longer used at all (no offense to Mike Fletcher and his loyal followers). After it was suggested last week, I'm now including a (stripped-down) version of Pyglet, which provides a ctypes-only OpenGL wrapper, in addition to excellent windowing, keyboard, and mouse support. Pyglet's event model is also compatible with the multi-threaded interface I described last week. That means you can now do this:
>>> from sympy import Symbol
>>> from sympy.modules.plotting import Plot
>>> x,y = Symbol('x'),Symbol('y')
>>> p = Plot(width=300, height=250, bbox=True)

>>> p[1] = x**2-y**2, [x,-1,1], [y,-1,1]
>>> 

>>> p[1] =  x**2+y**2, [x,-1.0,1.0], [y,-1.0,1.0]
>>> p[2] = -x**2-y**2, [x,-1.0,1.0], [y,-1.0,1.0]
>>> print p
[1]: x**2+y**2, [x,-1.0,1.0], [y,-1.0,1.0]
[2]: -x**2-y**2, [x,-1.0,1.0], [y,-1.0,1.0]
Also new:
  • Rotation mechanism based on a virtual trackball.
  • Scroll-wheel zoom.
  • Holding the shift key allows more precise movement.
  • Smooth time-based keyboard movement.
  • Added arrow and page up/down key control (ASDWRF is still there).
  • X,Y, and Z-blended coloring scheme.
In progress:
  • Configurable coloring.
  • Labeled coordinate axes.
  • Translation in the plane of the screen.
  • Automatic intervals.
  • Intuitive 2d support, including dynamic intervals.
  • Display lists (maybe faster for higher-resolution rendering, we'll see).
  • Anti-aliasing.
Let me know when you get a chance to try it out!

by Brian Jorgensen (noreply@blogger.com) at July 17, 2007 11:57 PM

Chris Wu

Back to bugs

Had to write some extra code for reducing to row echelon form but now it seems like my eigenvalues and eigenvector functions work properly (with help from Robert's nice polynomial work).

So it's back to that pesky expression bug from before the midterm evaluation. To recall, for whatever reason my Gram-Schmidt function isn't properly evaluating expressions. Back to the grind...

by cwu (noreply@blogger.com) at July 17, 2007 01:36 PM

Robert Schwarz

Finally, multivariate factorization

Ok, I gave up with Wang's algorithm, for now. It is not that hard, in general, but the article describing it left quite some questions open for me. So instead of procrastinating the implementation any further, I did another algorithm, just today. It is a lot simpler, also due to Kronecker (20th century!), and not too efficient, of course.

The last feature I want to add now, is solving polynomial equations. It is already possible now to compute a Gröbner base in a convenient order, like lexicographic, to get eliminiation working. Then, for each variable, you have to solve the univariate polynomial in that variable with solve(). You could insert each of these roots in the rest of the equation and proceed with the next variable. The problem with that is, that this approach only works for equations with a finite number of solution, which has to be decided first.

Then there is still some time left for polishing: Documentation and efficiency allow for improvement!

by Robert (noreply@blogger.com) at July 17, 2007 02:12 AM

July 11, 2007

Brian Jorgensen

PyOpenGL is Dead to Me (Where Are All the 3d Graphics Python Enthusiasts?)

Right now, displaying a plot permanently blocks the main thread; in other words, you can't do:
p = Plot( x**2, [x, -5, 5, 10] ) #never returns
p.append(x**3) #execution never gets here!
p.rotate(pi/2)
p.save("graph.jpg")
p.close()
I've tried desperately to get this to work using threading, but GLUT is not meant to be used in this way. GLUT basically forces the user to put all application logic into event callback functions which are called by the non-reentrant glutMainLoop. This won't work for either of my use cases, which are graphing from interactively from the console and embedding plotting functionality in larger applications/scripts.

The bottom line is that GLUT is not meant for real-world OpenGL applications, it's meant for learning OpenGL and and tangentially for writing simple demos and games. TOGL, a Tk Widget for OpenGL rendering which is also supported in PyOpenGL (though I couldn't get it to work right in less than an hour), suffers from essentially the same problem (Tk.mainloop is non-reentrant). As far as I can tell, these are really the only two interfaces which are available out-of-the-box in PyOpenGL.

So imagine my hope when I learned that another implementation of GLUT called Freeglut has a workaround for the inflexible event model problem. Tantalizingly, PyOpenGL 3.0 (a rewrite which uses ctypes instead of SWIG) supports Freeglut, as well as Python 2.5 (not currently supported by any functional version of PyOpenGL). I spent most of yesterday trying to get PyOpenGL 3.0 to work on Windows. Unfortunately, it just doesn't work yet, and PyOpenGL development is painfully slow.

As far as I can tell, PyOpenGL 3.0 has only one core developer, and he doesn't seem to have much time to work on it (search for 'OpenGL' on his blog, most of the results are from 2004 and 2005). I don't blame him personally, I'm just utterly shocked that Python doesn't have an actively maintained OpenGL interface library (next year's SoC?). A lag time of more than two years in producing a build for Python 2.5 on Windows is unacceptable in my view, given the size of the Python community and the traditional strength of its library support. [/Rant]

Today, I investigated a variety of new approaches (vaguely in reverse order of my preference):