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.

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.

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.
  

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.

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.

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.

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)

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!

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):

• Using Visual Python (annoyingly high-level, hideous website ;-)
• Using PyOGRE (huge, probably unacceptable overhead)
• Using the DirectX binding for Python on Windows, and PyOpenGL 3.0 elsewhere (disgusting platform dependency, PyOpenGL 3.0 is quite unstable and not expected to improve anytime soon)
• Fixing PyOpenGL to work on Windows (probably a monumental task worthy of its own SoC project)
• Creating my own limited ctypes wrapper for OpenGL (jury still out)
• Hacking PyOpenGL 2.0 to work with Freeglut (possibly feasible?)
• Creating a C or C++ Python rendering extension (how Visual Python works under the hood; has hacker sex appeal)
• Writing my own event loop with PyOpenGL 2.0, but without using GLUT (after exhaustive searching I couldn't find any examples or mention of this, but I want to keep looking. You'd think this would be the most reasonable option)
  

I don't particularly like any of these, though several might work. Anything good which I missed?

Post script: I'm particularly interested in hearing about anything not listed here, though I'd still like to know if there is a gem on there which I have overlooked.

SymPy Plot: Mouse Support

Now in SVN. Left-click controls rotation, right-click controls zoom. Eventually I want to change it so that right-click controls translation and middle-scroll controls zoom. Right now this does not seem to be possible, because the scroll wheel is not supported by GLUT's mouse handler. I'm looking into it.

Although it has been a while, I never mentioned that basis vectors were added. They still need some tweaking. Also 'q' is added as a quit key.

Right now I'm at a family reunion for the 4th of July holiday. Each summer, my extended family meets for a week to act out our own miniaturized take on Burning Man. FYI, my contract delivery went well on Thursday. However, I'm not quite off the hook yet as I had speculated; it turns out they need additional work.

SymPy 0.4.1 released, with windows installer

This is the first release to include SymPy Plot.

0.4.1 source tarball
0.4.1 windows installer
SymPy homepage
Changelog

PyOpenGL is required for SymPy plotting support. While *nix users should have no problem (apt-get install python-opengl), Windows users will find it difficult, if not impossible, to use PyOpenGL with Python 2.5. I recommend setting up your Python environment as follows:

1. Install Python 2.4.4 to C:\Python24 (the default location)
   I normally like to install all of my programs in C:\Program Files, but unfortunately PyOpenGL depends on the default install path.
2. Install Numeric-23.7.win32-py2.4.exe
3. Install numarray-1.1.1.win32-py2.4.exe
4. Install PyOpenGL-2.0.2.01.py2.4-numpy23.exe
5. Install sympy-0.4.1.win32.exe

6. >>> from sympy import Symbol
   >>> from sympy.modules.graphing import Plot
   >>> x = Symbol('x')
   >>> Plot(x**2, [x, -1, 1, 10])

   

If you get a dll not found error relating to GLUT, obtain glut32.dll and put it in C:\Python24\Lib\site-packages\OpenGL.

SymPy Plot: Now available in SVN

To try it out, check out or update from SVN and look at examples\plotting.py. You'll need to download and setup PyOpenGL for your platform: on ubuntu, apt-get install python-opengl; on windows, you might have luck with the instructions here. So far, I've only tried it on Windows with Python 2.4.



Enjoy, and let me know how it goes!

SymPy Plot: First screenshots

>>> f1 = 0.1+x**2+y**2
>>> f2 = -0.1-x**2-y**2
>>> Plot( f1, f2, [x, -0.75, 0.75, 8], [y, -0.75, 0.75, 8] )

>>> Plot( x*y**3-y*x**3, [x, -1, 1, 20], [y, -1, 1, 20] )

Vision for an extensible visualization tool

In addition to providing standard 2d and 3d function graphing capabilities, the SymPy graphing module will be able to render 2d and 3d geometric entities on the plot, such as lines, vectors, and triangles. While the rendering system will be completely hidden from the casual user who simply wants to plot a few curves or surfaces, advanced users will be able to extend or create new kinds of renderable objects.

A new submodule called sympy.modules.graphing.scene makes this possible by providing a high-level rendering abstraction on top of OpenGL through the Scene class. A Scene handles the boilerplate OpenGL setup and teardown code. It also maintains a list of objects to be rendered, which are represented by objects which implement the Renderable interface. Renderable is a base class providing the render() function, which retrieves the object's primitive representation for rendering. This representation is something that OpenGL can render directly, such as a line list or a tri strip.

For example, consider the Renderable called CartesianCurve, which represents a plane curve in the form y = f(x); when CartesianCurve's render() function is called by a Scene which contains it, it returns a representative list of coordinate pairs, which are then used to render an OpenGL line list. The various other types of plots will be implemented as subclasses of Renderable in the same way. Here are some examples of possible constructor syntax for function plot Renderables:

# y = f(x)
CartesianCurve( x**3, [x, -10, 10, 20] )

# x = f(t), y = f(t)
ParametricCurve( cos(t), sin(t), [t, 0, 2*pi, 32] )

# r = f(theta)
PolarCurve( t, [t, 0, 2*pi, 32] )

# z = f(x, y)
CartesianSurface( x**2+y**2, [x, -1, 1, 5], [y, -1, 1, 5] )

# radius = f(z, theta)
CylindricalShell( 1.0, [z, 0, 1, 2], [t, 0, 2*pi, 32] )
CylindricalShell( t*sqrt(z), [z, 0, 1, 2] )

# radius = f(theta, rho)
SphericalShell( 1.0, [t], [p] )

SymPy Graphing: Starting with PyOpenGL and Feisty

Just a brief update. Progress has been a little slower this week because I've begun working on the new graphing implementation using PyOpenGL and the sample function I wrote about last week. After about a day of tinkering, I had no luck with PyOpenGL on Windows using Python 2.5 (not officially supported yet, apparently). As a result, I've been working in Feisty for the last few days. Unfortunately, I'm having a hard time convincing Feisty to connect to my wpa2-protected wifi, and this is putting a damper on my work flow. Once I get that and a few other nags sorted out, I should be on track.

New direction for SymPy Graphing

In addition to gaining some new perspective on my project this weekend, I've gotten a lot done.



On Friday, I said I wanted to stop supporting matplotlib and to begin focusing on a fresh OpenGL plotting solution. It turns out that matplotlib has better 3d support than I suggested, though it is buggy and still doesn't have everything I need. Fortunately, it turns out we can have it both ways.

I've moved matplotlib 2d plotting support out of modules/graphing.py and into examples/mplot2d.py. I've also added 3d plotting support in examples/mplot3d.py. In accomplishing this, I created a generic sampling function in modules/graphing.py which generates plotting data from a 2d or 3d function over an interval or rectangle:


>>> from sympy.modules.graphing import sample
>>> from sympy import Symbol
>>> x = Symbol('x')
>>> sample(x, (x, 0, 2, 2))
[[0,1,2],[0,1,2]]
>>> sample(x**2, (x, 1, 3, 2))
[[1,2,3], [1,4,9]]


This approach will allow me to support multiple plotting solutions. You should now be able to graph many varieties of functions using matplotlib with the implementations in examples/mplot2d.py and examples/mplot3d.py.


>>> mplot3d(x**2-y**2, (x, -10.0, 10.0, 20), (y, -10.0, 10.0, 20))


The Case for Dropping Matplotlib

My pet motivation for SymPy Graphing is to aid in the creation of interactive calculus visualizations. I want to be able to write an example program which visualizes the normal and tangential components of acceleration along a curve. I want to be able to graph a 3d surface subject to a constraint, perhaps to see only the portion of a surface which falls within a cylindrical boundary. Ondrej mentioned that he would like to be able to plot curvature as the vertex color across a 3d surface. All of these things should be doable. At the same time, SymPy Graphing should be simple and intuitive enough for a non-programming Calc III student to use in Python's interactive mode as a graphing calculator.



In my work on 2d graphing this week, I've utilized matplotlib's pylab interface. Using matplotlib to plot functions from SymPy is pretty straightforward, though I find it clunky and have had many frustrations with it. But as the name implies, matplotlib isn't meant to do much beyond plotting, and has only limited 3d support. I want to do a lot more with SymPy Graphing.

As mentioned above, I'd like to support plotting an arbitrary function as the vertex color across a surface; matplotlib has no mechanism for this. Jason and I have discussed allowing geometric figures and line segments from the geometry module to be drawn through the graphing interface; again, not possible. I'd also like to support interactive rotation, zooming, translation, etc., and possibly someday (outside the scope of GSoC) point-picking; to my knowledge, you cannot do these things in matplotlib.

As I said in my proposal, 3d graphing will be supported through PyOpenGL. In general, this is not a very difficult task (there is a good summary of what it entails at the plot3d website) and this approach will provide the flexibility I want.

The upshot of all this is that I am considering a shift away from matplotlib altogether and toward supporting all graphing, 2d and 3d, through a unified OpenGL interface. I only want to support a fraction of the features of matplotlib. (I'm not trying to banish matplotlib from sympy. In any event, one can easily use matplotlib directly. See the current revision of modules/graphing.py in SVN for an example.) If my project goal was simply to write a SymPy wrapper for matplotlib, I'd be done before the official coding period even begins. To do this right, I'm going to need to get my hands a bit dirtier.

Update:
The 3d support in Matplotlib is better than I'd thought (it can do interactive rotation and scaling). There is a good example here. It also looks like it might be able to do colormaps. I'm looking into it in more depth right now.

Sympy Graphing: Usable Progress, Early Frustrations

Over the weekend, I did a bit of refactoring and added a couple of improvements to sympy/modules/graphing.py (where Fabian had laid a good foundation already). First, I've extended plot() to accept a list of functions in addition to a single function.

from sympy import Symbol
from sympy.modules.graphing import plot
x = Symbol('x')
plot([sqrt(x), log(x), sin(x)], xlabel="x label")


I also added in some code so that functions with discontinuities, such as sqrt where x < 0, are plotted correctly.

f_list = [sqrt(x), sqrt(-x), -sqrt(x), -sqrt(-x)]
plot(f_list, [x, -10.0, 10.0], title="Square Roots")


Unfortunately, this code has hit a snag getting through automated testing. At this point, my best hypothesis is that there is some incompatibility between py.test/doctest and the Tk backend for pylab. I'm going to try a different backend next. Ondrej has opened an issue here.

Google Summer of Code 2007 with SymPy

My project proposal Python as a Symbolic Graphing Calculator with SymPy, PIL, and PyOpenGL has been accepted for Google Summer of Code 2007 under the patronage of the Python Software Foundation!

I'll be working with Ondrej Certik, who started the SymPy project, a CAS for Python. I wanted to contribute to the SymPy project for several reasons. Most importantly, I love Python. Even when I'm coding C#, C++, or Java I usually end up metaprogramming in Python much of the time. Secondly, I am a math student. SymPy is useful to me already, and it's only in version 0.3. I don't think I'm the only one who is unsatisfied with other CAS offerings. I think SymPy has a lot of potential, and I'm glad to be part of the team.

Once I get school over with, this project will have my full attention. I'll be updating the blog with my progress and anything interesting I find along the way. I'll also be writing about other summer projects, including my summer reading. More on this soon.