Thursday, December 13, 2007

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.
at