Josh Kantor's Lorenz attractor example

Josh posted a nice example of plotting a Lorenz attractor in Sage:


Put this in a notebook cell (be careful about newlines):

Integer = int
RealNumber = float

def lorenz(t,y,params):
    return [params[0]*(y[1]-y[0]),y[0]*(params[1]-y[2])- y[1],y[0]*y[1]-params[2]*y[2]]

def lorenz_jac(t,y,params):
    return [ [-params[0],params[0],0],[(params[1]-y[2]),-1,-y[0]],[y[1],y[0],-params[2]],[0,0,0]]

T=ode_solver()
T.algorithm="bsimp"
T.function=lorenz
T.jacobian=lorenz_jac
T.ode_solve(y_0=[.5,.5,.5],t_span=[0,155],params=[10,40.5,3],num_points=10000)
l=[T.solution[i][1] for i in range(len(T.solution))]

line3d(l,thickness=0.3, viewer='tachyon', figsize=8)

and this is what you get (click to zoom):

I'm glad I chose Python for Sage -- some cool scientific computing Python projects

When brainstorming for talk ideas for a workshop, Fernando Perez came up with a bunch of random off-the-top of his head very high quality/impact scientific computing projects that involve Python. Here they are (the following was written by Fernando Perez, the author of IPython):

- I can obviously talk about ipython and related projects, but I can also give a math/technical talk off this type of work (the whole implementation is python, and uses quite a few tricks):

http://dx.doi.org/10.1016/j.acha.2007.08.001

- Brian Granger (from Tech-X http://txcorp.org) has a NASA grant to develop distributed arrays for Python, using IPython and numpy. That would make for an excellent talk, I think (matlab, interactivesupercomputing.com, and all the 'big boys' are after distributed arrays).

- Trilinos (http://trilinos.sandia.gov), a large set of parallel solvers developed at Sandia National Lab, has a great set of Python bindings (even usable interactively via ipython).

- MPI4Py is an excellent set of MPI bindings for Python, and its author Lisandro Dalcin is also the developer of Petsc4py:
* http://mpi4py.scipy.org/
* http://code.google.com/p/petsc4py/
If Lisandro can't come, I can contact one of the Petsc guys who's a great python developer and see if he's coming or can make it, he's an excellent speaker (Matt Knepley from Argonne Nat. Lab http://www-unix.mcs.anl.gov/~knepley/).

- The Hubble space telescope people are all pyhton based, and have done enormous amounts of work on both their internal image processing pipeline and contrbuting to Matplotlib (they have currently a developer working full time on matplotlib).

- The NetworkX (Sage uses this) guys from Los Alamos will probably be coming: https://networkx.lanl.gov/wiki.

- The Scripps institute has an extremely strong Python team doing molecular visualization and visual programming. Their work is very impressive, and they're already in San Diego, so it's a no brainer to attend. Their presentations are always of very high quality: http://mgltools.scripps.edu.

- JPL uses python extensively (they've contracted out work for matplotlib to be extended to suit their specific needs).

- My new job at UC Berkeley is on neuroscience, and we could present the work that's being developed there for fMRI analysis (all Python based, fully open source, NIH funded).

- Andrew Straw's work (http://www.its.caltech.edu/~astraw/) on real-time 3d tracking of fruit flies is very, very impressive. All python based, hardware control, real-time parallel computing.

- The CACR group at Caltech has the contract for DANSE (http://wiki.cacr.caltech.edu/danse/index.php/Main_Page), the Spallation Neutron Source's data analysis framework, all python. This is currently the largest experiment being funded in the USA.

AMS Meeting Day 2: The Competition

Today I violated the Mathematica license agreement in front of Eric Weisstein (a famous Mathematica developer), I talked with the developers of Wiris (rhymes with virus) which is a commercial competitor to the Sage notebook, discussed Mupad with a longtime Mupad developer, and gave away $1000 of tutorials, DVD's, and other goodies at the Sage AMS exhibit booth.

There is a project called "Wiris", which I had never heard of until today. So Tom Boothby and I had a very interesting talk with the people at the Wiris Booth (http://www.wiris.com/). Wiris is a closed-source commercial math software company in Barcelona that makes a web-based interface to their own custom mathematical software (interestingly, one of their main developers took Calculus from Jordi Quer, who wrote a lot of modular-forms related code for Sage recently). Their primary audience and market right now is European Government Agencies who use their software for high school and beginning college education. Their software is much different than the Sage notebook, since it is written entirely in Java instead of being an AJAX javascript application. They knew about Sage and asked if it used OpenMath or MathML, and I explained that it didn't use either, that it shouldn't and that for our purposes (i.e., interfacing math software) those technologies do not solve the problems we have -- in fact, they are worse than useless. They said that us not taking the OpenMath route was disappointing. They ended the discussion by telling us that their web-based interface is much better than ours :-). I guess they were trying to intimidate us.

Eric Weisstein -- who told us that he is now an official Mathematica developer (doing graph theory among other things) in addition to his "Math World" came over to the Sage booth and asked a lot of pointed questions about Sage, mostly "how does this get funded?" We ended up talking for about an hour. There are numerous people at Mathematica who are well aware of Sage, and he claims he doesn't see Sage as a threat to Mathematica as a company. He said Mathematica has had a recent explosion in hiring as a result of greatly increased sales because of the new "Demonstrations" feature in Mathematica 6. We also talked a lot about graph theory and how Sage has a complete implementation of graph isomorphism testing, etc., which greatly impressed him (thanks Robert Miller!). Jason Grout asked about the Mathematica end use license agreement and Eric pleaded IANAL and I demonstrated using Mathematica via the Sage notebook -- all locally over localhost on a machine with a valid Mathematica license -- that this violates the Mathematica license. He got annoyed when I did a Mathematica graph through the Sage notebook... Later Eric talked about how he "hoped" Sage would continue to have momentum and not just die like other free projects. He gave Maxima as an example of a dead project, and seemed quite shocked when I mentioned that they are very much alive and have regular releases, etc. He then said that there is no such thing as free.

I met somebody from MuPad who has worked on that closed source project for a decade. They used to be a German national government funded project, then had an academic research branch until one year ago, and now are 100% commercial and private. I asked about their vision for the next 5 years, and he said "MuPad will survive", and said they were mainly happy to be stable. They also mentioned wanting to do more numerical and applied functionality.

I talked with a Unix guy who works on some Scientific Workplace. There is only a Windows version, and he is working on doing an OS X/Linux port, which they will finish "in the future"!

And that's just a little of what happened today...

AMS Meeting Part 1

I am completely exhausted right now, having put in a huge amount of time with little sleep during the last few days writing 3d graphics code for Sage (with Robert Bradshaw), making the Sage-2.9.2 release (with Michael Abshoff), printing fliers, tutorials, and preparing DVDs (with my brother) for the massive joint AMS meeting where there will be a Sage exhibit booth.

Our booth with have a huge banner that says "Creating a viable free open source alternative to Magma, Maple, Mathematica, and Matlab."
An Austrian named Harald Schilly did most of the work creating this poster.

I have no clue what to expect during the next few days. I've never run an exhibit booth before; I'm a mathematician not a "vendor", and my product is free. This should be very interesting. I've been to numerous AMS meetings before, and every booth I remember is commercial... TI, book publishers, math games, Maple, Mathematica, etc.

Sage 3.0?

The main goals for 3.0, which we should aim to complete in Feb 2008 are:

(1) SOLARIS: Full Solaris 10 support on Opteron and Sparc,

(2) 3D: Excellent interactive 3d graphics from the notebook and command line by working with the Jmol project.

(3) DEMOS: A demos server with interactive web-based use of Sage for specific applications, along with a catalogue of 500 examples that illustrate Sage and drive traffic to sagemath.org.

(4) DOCTESTS: Raise the doctest coverage of the Sage library to 50%. (Right now it is at 35.5%).

The Growth of Open Source Mathematics Software

With help from Martin Albrecht, I've been tracking the traffic to sagemath.org.

Some of the above plots compare Sagemath.org to Maplesoft.com. Before the UW press release (inspired by Sage winning Trophees du Libre), sagemath.org traffic was insignificant when compared to Maple. On Dec 3 when the press release went out until Dec 11, there was a huge spike in traffic, so that sagemath.org went from being a speck to over twice the amount of traffic as Maple gets. NOTE -- this increase in traffic started two days before the slashdotting. On Dec 5 during the slashdotting, the number of downloads of sage itself went up a huge amount (thousands per day). This hump died down after about a week on December 11. However, the number of Sage users has gone up by quite a bit, and perhaps a nontrivial fraction of the people who tried Sage a week ago are now Sage users. Traffic to sagemath.org has gone from basically nothing as measured by Alexa to nontrivial in comparison to the other big Ma*'s.

Perhaps soon trac #1000 can be closed and replaced by a ticket to have 100,000 users?

Sage, Mathematica and Hiking to Vultee Arch

Here is an example of implementing a special function in Sage using Mathematica to do the real work:

def math_bessel_K(nu,x):
    return mathematica(nu).BesselK(x).N(20).sage()

sage: math_bessel_K(2,I)
0.180489972066962*I - 2.592886175491197

When I showed this to a Sage user they exclaimed:

> This is incredibly cool! I can't believe I missed this in the Sage
> reference manual, but there it is. In wonder how it works?

This blog post is about how the above works.

I figured out a clean way to have this sort of notation work in
Mathematica, etc., from Sage when I was hiking with my wife in Sedona to Vultee Arch a couple years ago. Basically when you do say

s = mathematica(nu)

Sage converts nu into a Mathematica-readable string by calling nu._mathematica_init_():

sage: nu = sqrt(2) + I*sin(3)
sage: nu._mathematica_init_()
'(Sqrt[2]) + ((I) * (Sin[3]))'

Note that _mathematica_init_ on the nu in the above example calls other _mathematica_init_'s, e.g., sqrt(2) has a mathematica_init_ method that calls that _mathematica_init_ method for 2, etc. Mathematica itself is not used at all yet. Incidentally, just like you can latex expressions, e.g.,

sage: latex(nu)
{i \cdot \sin \left( 3 \right)} + \sqrt{ 2 }

for inclusion in a paper or something, you can also mathematica them as above, to share with mathematica users...

Anyway Sage sends this string to mathematica:

'sage0 := (Sqrt[2]) + ((I) * (Sin[3]))'

Mathematica evaluates it and Sage creates a Python wrapper object, which knows the variable name 'sage0':

sage: s = mathematica(nu)
sage: s.name( )
'sage0'
sage: type(s)

sage: s
Sqrt[2] + I*Sin[3]

Now there is a running instance of mathematica, and it has an object defined in it
called sage0, which is equal to our nu:

sage: mathematica.eval('sage0')
Sqrt[2] + I Sin[3]


You can even play around with that copy of Mathematica more directly from the command
line (or with %mathematica in a cell in the Sage notebook):

sage: %mathematica

--> Switching to Mathematica <--

''
mathematica: sage0
Sqrt[2] + I*Sin[3]
mathematica: Sqrt[sage0]
Sqrt[Sqrt[2] + I*Sin[3]]
mathematica: N[sage0]
1.4142135623730951 + 0.1411200080598672*I

[[Note this is like the Mathematica command line but better since it has readline support, etc.]]

mathematica: %sage
--> Exiting back to SAGE <--

Finally, if you want to call a function and give the version of nu that is in Mathematica as an argument, you just call the function on it using Pythonic notation:

sage: s.Sqrt()
Sqrt[Sqrt[2] + I*Sin[3]]
sage: s.N()
1.4142135623730951 + 0.1411200080598672*I

What happens here is that say s.Sqrt creates a Python class that wraps the Mathematica function Sqrt with argument s not yet evaluated, so even this works:

sage: s.Sqrt?
Docstring:
Sqrt[z] or Sqrt[z] gives the square root of z.

When we did s.Sqrt? Sage queries Mathematica for how Sqrt works.

Everything above applies equally to the Maple, Matlab, Maxima, PARI, Singular, etc., interfaces. There is *much* that can be done to make the Mathematica interface in particular even more cool. Ideas appreciated!

By the way, it snowed very nicely on the hike, and my wife and I had a great time.