r"""
Interface to Mathematica
The Mathematica interface will only work if Mathematica is installed on your
computer with a command line interface that runs when you give the ``math``
command. The interface lets you send certain Sage objects to Mathematica,
run Mathematica functions, import certain Mathematica expressions to Sage,
or any combination of the above.
To send a Sage object ``sobj`` to Mathematica, call ``mathematica(sobj)``.
This exports the Sage object to Mathematica and returns a new Sage object
wrapping the Mathematica expression/variable, so that you can use the
Mathematica variable from within Sage. You can then call Mathematica
functions on the new object; for example::
sage: mobj = mathematica(x^2-1) # optional - mathematica
sage: mobj.Factor() # optional - mathematica
(-1 + x)*(1 + x)
In the above example the factorization is done using Mathematica's
``Factor[]`` function.
To see Mathematica's output you can simply print the Mathematica wrapper
object. However if you want to import Mathematica's output back to Sage,
call the Mathematica wrapper object's ``sage()`` method. This method returns
a native Sage object::
sage: mobj = mathematica(x^2-1) # optional - mathematica
sage: mobj2 = mobj.Factor(); mobj2 # optional - mathematica
(-1 + x)*(1 + x)
sage: mobj2.parent() # optional - mathematica
Mathematica
sage: sobj = mobj2.sage(); sobj # optional - mathematica
(x - 1)*(x + 1)
sage: sobj.parent() # optional - mathematica
Symbolic Ring
If you want to run a Mathematica function and don't already have the input
in the form of a Sage object, then it might be simpler to input a string to
``mathematica(expr)``. This string will be evaluated as if you had typed it
into Mathematica::
sage: mathematica('Factor[x^2-1]') # optional - mathematica
(-1 + x)*(1 + x)
sage: mathematica('Range[3]') # optional - mathematica
{1, 2, 3}
If you don't want Sage to go to the trouble of creating a wrapper for the
Mathematica expression, then you can call ``mathematica.eval(expr)``, which
returns the result as a Mathematica AsciiArtString formatted string. If you
want the result to be a string formatted like Mathematica's InputForm, call
``repr(mobj)`` on the wrapper object ``mobj``. If you want a string
formatted in Sage style, call ``mobj._sage_repr()``::
sage: mathematica.eval('x^2 - 1') # optional - mathematica
2
-1 + x
sage: repr(mathematica('Range[3]')) # optional - mathematica
'{1, 2, 3}'
sage: mathematica('Range[3]')._sage_repr() # optional - mathematica
'[1, 2, 3]'
Finally, if you just want to use a Mathematica command line from within
Sage, the function ``mathematica_console()`` dumps you into an interactive
command-line Mathematica session. This is an enhanced version of the usual
Mathematica command-line, in that it provides readline editing and history
(the usual one doesn't!)
Tutorial
--------
We follow some of the tutorial from
http://library.wolfram.com/conferences/devconf99/withoff/Basic1.html/.
For any of this to work you must buy and install the Mathematica
program, and it must be available as the command
``math`` in your PATH.
Syntax
~~~~~~
Now make 1 and add it to itself. The result is a Mathematica
object.
::
sage: m = mathematica
sage: a = m(1) + m(1); a # optional - mathematica
2
sage: a.parent() # optional - mathematica
Mathematica
sage: m('1+1') # optional - mathematica
2
sage: m(3)**m(50) # optional - mathematica
717897987691852588770249
The following is equivalent to ``Plus[2, 3]`` in
Mathematica::
sage: m = mathematica
sage: m(2).Plus(m(3)) # optional - mathematica
5
We can also compute `7(2+3)`.
::
sage: m(7).Times(m(2).Plus(m(3))) # optional - mathematica
35
sage: m('7(2+3)') # optional - mathematica
35
Some typical input
~~~~~~~~~~~~~~~~~~
We solve an equation and a system of two equations::
sage: eqn = mathematica('3x + 5 == 14') # optional - mathematica
sage: eqn # optional - mathematica
5 + 3*x == 14
sage: eqn.Solve('x') # optional - mathematica
{{x -> 3}}
sage: sys = mathematica('{x^2 - 3y == 3, 2x - y == 1}') # optional - mathematica
sage: print sys # optional - mathematica
2
{x - 3 y == 3, 2 x - y == 1}
sage: sys.Solve('{x, y}') # optional - mathematica
{{x -> 0, y -> -1}, {x -> 6, y -> 11}}
Assignments and definitions
~~~~~~~~~~~~~~~~~~~~~~~~~~~
If you assign the mathematica `5` to a variable `c`
in Sage, this does not affect the `c` in Mathematica.
::
sage: c = m(5) # optional - mathematica
sage: print m('b + c x') # optional - mathematica
b + c x
sage: print m('b') + c*m('x') # optional - mathematica
b + 5 x
The Sage interfaces changes Sage lists into Mathematica lists::
sage: m = mathematica
sage: eq1 = m('x^2 - 3y == 3') # optional - mathematica
sage: eq2 = m('2x - y == 1') # optional - mathematica
sage: v = m([eq1, eq2]); v # optional - mathematica
{x^2 - 3*y == 3, 2*x - y == 1}
sage: v.Solve(['x', 'y']) # optional - mathematica
{{x -> 0, y -> -1}, {x -> 6, y -> 11}}
Function definitions
~~~~~~~~~~~~~~~~~~~~
Define mathematica functions by simply sending the definition to
the interpreter.
::
sage: m = mathematica
sage: _ = mathematica('f[p_] = p^2'); # optional - mathematica
sage: m('f[9]') # optional - mathematica
81
Numerical Calculations
~~~~~~~~~~~~~~~~~~~~~~
We find the `x` such that `e^x - 3x = 0`.
::
sage: e = mathematica('Exp[x] - 3x == 0') # optional - mathematica
sage: e.FindRoot(['x', 2]) # optional - mathematica
{x -> 1.512134551657842}
Note that this agrees with what the PARI interpreter gp produces::
sage: gp('solve(x=1,2,exp(x)-3*x)')
1.512134551657842473896739678 # 32-bit
1.5121345516578424738967396780720387046 # 64-bit
Next we find the minimum of a polynomial using the two different
ways of accessing Mathematica::
sage: mathematica('FindMinimum[x^3 - 6x^2 + 11x - 5, {x,3}]') # optional - mathematica
{0.6150998205402516, {x -> 2.5773502699629733}}
sage: f = mathematica('x^3 - 6x^2 + 11x - 5') # optional - mathematica
sage: f.FindMinimum(['x', 3]) # optional - mathematica
{0.6150998205402516, {x -> 2.5773502699629733}}
Polynomial and Integer Factorization
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
We factor a polynomial of degree 200 over the integers.
::
sage: R.<x> = PolynomialRing(ZZ)
sage: f = (x**100+17*x+5)*(x**100-5*x+20)
sage: f
x^200 + 12*x^101 + 25*x^100 - 85*x^2 + 315*x + 100
sage: g = mathematica(str(f)) # optional - mathematica
sage: print g # optional - mathematica
2 100 101 200
100 + 315 x - 85 x + 25 x + 12 x + x
sage: g # optional - mathematica
100 + 315*x - 85*x^2 + 25*x^100 + 12*x^101 + x^200
sage: print g.Factor() # optional - mathematica
100 100
(20 - 5 x + x ) (5 + 17 x + x )
We can also factor a multivariate polynomial::
sage: f = mathematica('x^6 + (-y - 2)*x^5 + (y^3 + 2*y)*x^4 - y^4*x^3') # optional - mathematica
sage: print f.Factor() # optional - mathematica
3 2 3
x (x - y) (-2 x + x + y )
We factor an integer::
sage: n = mathematica(2434500) # optional - mathematica
sage: n.FactorInteger() # optional - mathematica
{{2, 2}, {3, 2}, {5, 3}, {541, 1}}
sage: n = mathematica(2434500) # optional - mathematica
sage: F = n.FactorInteger(); F # optional - mathematica
{{2, 2}, {3, 2}, {5, 3}, {541, 1}}
sage: F[1] # optional - mathematica
{2, 2}
sage: F[4] # optional - mathematica
{541, 1}
We can also load the ECM package and factoring using it::
sage: _ = mathematica.eval("<<NumberTheory`FactorIntegerECM`"); # optional - mathematica
sage: mathematica.FactorIntegerECM('932901*939321') # optional - mathematica
8396109
Long Input
----------
The Mathematica interface reads in even very long input (using
files) in a robust manner.
::
sage: t = '"%s"'%10^10000 # ten thousand character string.
sage: a = mathematica(t) # optional - mathematica
sage: a = mathematica.eval(t) # optional - mathematica
Loading and saving
------------------
Mathematica has an excellent ``InputForm`` function,
which makes saving and loading Mathematica objects possible. The
first examples test saving and loading to strings.
::
sage: x = mathematica(pi/2) # optional - mathematica
sage: print x # optional - mathematica
Pi
--
2
sage: loads(dumps(x)) == x # optional - mathematica
True
sage: n = x.N(50) # optional - mathematica
sage: print n # optional - mathematica
1.5707963267948966192313216916397514420985846996876
sage: loads(dumps(n)) == n # optional - mathematica
True
Complicated translations
------------------------
The ``mobj.sage()`` method tries to convert a Mathematica object to a Sage
object. In many cases, it will just work. In particular, it should be able to
convert expressions entirely consisting of:
- numbers, i.e. integers, floats, complex numbers;
- functions and named constants also present in Sage, where:
- Sage knows how to translate the function or constant's name from
Mathematica's, or
- the Sage name for the function or constant is trivially related to
Mathematica's;
- symbolic variables whose names don't pathologically overlap with
objects already defined in Sage.
This method will not work when Mathematica's output includes:
- strings;
- functions unknown to Sage;
- Mathematica functions with different parameters/parameter order to
the Sage equivalent.
If you want to convert more complicated Mathematica expressions, you can
instead call ``mobj._sage_()`` and supply a translation dictionary::
sage: m = mathematica('NewFn[x]') # optional - mathematica
sage: m._sage_(locals={'NewFn': sin}) # optional - mathematica
sin(x)
For more details, see the documentation for ``._sage_()``.
OTHER Examples::
sage: def math_bessel_K(nu,x):
... return mathematica(nu).BesselK(x).N(20)
...
sage: math_bessel_K(2,I) # optional - mathematica
0.180489972066962*I - 2.592886175491197 # 32-bit
-2.59288617549119697817 + 0.18048997206696202663*I # 64-bit
::
sage: slist = [[1, 2], 3., 4 + I]
sage: mlist = mathematica(slist); mlist # optional - mathematica
{{1, 2}, 3., 4 + I}
sage: slist2 = list(mlist); slist2 # optional - mathematica
[{1, 2}, 3., 4 + I]
sage: slist2[0] # optional - mathematica
{1, 2}
sage: slist2[0].parent() # optional - mathematica
Mathematica
sage: slist3 = mlist.sage(); slist3 # optional - mathematica
[[1, 2], 3.0, I + 4]
::
sage: mathematica('10.^80') # optional - mathematica
1.*^80
sage: mathematica('10.^80').sage() # optional - mathematica
1e+80
AUTHORS:
- William Stein (2005): first version
- Doug Cutrell (2006-03-01): Instructions for use under Cygwin/Windows.
- Felix Lawrence (2009-08-21): Added support for importing Mathematica lists
and floats with exponents.
"""
import os
import re
from expect import (Expect, ExpectElement, ExpectFunction,
FunctionElement, AsciiArtString)
from sage.misc.misc import graphics_filename
def clean_output(s):
if s is None:
return ''
i = s.find('Out[')
j = i + s[i:].find('=')
s = s[:i] + ' '*(j+1-i) + s[j+1:]
s = s.replace('\\\n','')
return s.strip('\n')
def _un_camel(name):
"""
Convert `CamelCase` to `camel_case`.
EXAMPLES::
sage: sage.interfaces.mathematica._un_camel('CamelCase')
'camel_case'
sage: sage.interfaces.mathematica._un_camel('EllipticE')
'elliptic_e'
sage: sage.interfaces.mathematica._un_camel('FindRoot')
'find_root'
sage: sage.interfaces.mathematica._un_camel('GCD')
'gcd'
"""
s1 = re.sub('(.)([A-Z][a-z]+)', r'\1_\2', name)
return re.sub('([a-z0-9])([A-Z])', r'\1_\2', s1).lower()
class Mathematica(Expect):
"""
Interface to the Mathematica interpreter.
"""
def __init__(self, maxread=100, script_subdirectory="", logfile=None, server=None, server_tmpdir=None):
Expect.__init__(self,
name = 'mathematica',
prompt = 'In[[0-9]+]:=',
command = "math",
maxread = maxread,
server = server,
server_tmpdir = server_tmpdir,
script_subdirectory = script_subdirectory,
verbose_start = False,
logfile=logfile,
eval_using_file_cutoff=50)
def _read_in_file_command(self, filename):
return '<<"%s"'%filename
def _keyboard_interrupt(self):
print "Interrupting %s..."%self
e = self._expect
e.sendline(chr(3))
e.expect('Interrupt> ')
e.sendline("a")
e.expect(self._prompt)
return e.before
def _install_hints(self):
"""
Hints for installing mathematica on your computer.
AUTHORS:
- William Stein and Justin Walker (2006-02-12)
"""
return """
In order to use the Mathematica interface you need to have Mathematica
installed and have a script in your PATH called "math" that runs the
command-line version of Mathematica. Alternatively, you could use a
remote connection to a server running Mathematica -- for hints, type
print mathematica._install_hints_ssh()
(1) You might have to buy Mathematica (http://www.wolfram.com/).
(2) * LINUX: The math script comes standard with your Mathematica install.
* APPLE OS X:
(a) create a file called math (in your PATH):
#!/bin/sh
/Applications/Mathematica.app/Contents/MacOS/MathKernel $@
The path in the above script must be modified if you installed
Mathematica elsewhere or installed an old version of
Mathematica that has the version in the .app name.
(b) Make the file executable.
chmod +x math
* WINDOWS:
Install Mathematica for Linux into the VMware virtual machine (sorry,
that's the only way at present).
"""
def eval(self, code, strip=True, **kwds):
s = Expect.eval(self, code, **kwds)
if strip:
return AsciiArtString(clean_output(s))
else:
return AsciiArtString(s)
def set(self, var, value):
"""
Set the variable var to the given value.
"""
cmd = '%s=%s;'%(var,value)
out = self._eval_line(cmd, allow_use_file=True)
if len(out) > 8:
raise TypeError, "Error executing code in Mathematica\nCODE:\n\t%s\nMathematica ERROR:\n\t%s"%(cmd, out)
def get(self, var, ascii_art=False):
"""
Get the value of the variable var.
AUTHORS:
- William Stein
- Kiran Kedlaya (2006-02-04): suggested using InputForm
"""
if ascii_art:
return self.eval(var, strip=True)
else:
return self.eval('InputForm[%s, NumberMarks->False]'%var, strip=True)
def _eval_line(self, line, allow_use_file=True, wait_for_prompt=True, restart_if_needed=False):
s = Expect._eval_line(self, line,
allow_use_file=allow_use_file, wait_for_prompt=wait_for_prompt)
return str(s).strip('\n')
def _function_call_string(self, function, args, kwds):
"""
Returns the string used to make function calls.
EXAMPLES::
sage: mathematica._function_call_string('Sin', ['x'], [])
'Sin[x]'
"""
return "%s[%s]"%(function, ",".join(args))
def _left_list_delim(self):
return "{"
def _right_list_delim(self):
return "}"
def _left_func_delim(self):
return "["
def _right_func_delim(self):
return "]"
def chdir(self, dir):
"""
Change Mathematica's current working directory.
EXAMPLES::
sage: mathematica.chdir('/') # optional
sage: mathematica('Directory[]') # optional
"/"
"""
self.eval('SetDirectory["%s"]'%dir)
def _true_symbol(self):
return ' True'
def _false_symbol(self):
return ' False'
def _equality_symbol(self):
return '=='
def _assign_symbol(self):
return ":="
def _exponent_symbol(self):
"""
Returns the symbol used to denote the exponent of a number in
Mathematica.
EXAMPLES::
sage: mathematica._exponent_symbol() # optional - mathematica
'*^'
::
sage: bignum = mathematica('10.^80') # optional - mathematica
sage: repr(bignum) # optional - mathematica
'1.*^80'
sage: repr(bignum).replace(mathematica._exponent_symbol(), 'e').strip() # optional - mathematica
'1.e80'
"""
return "*^"
def _object_class(self):
return MathematicaElement
def console(self, readline=True):
mathematica_console(readline=readline)
def trait_names(self):
a = self.eval('Names["*"]')
return a.replace('$','').replace('\n \n>','').replace(',','').replace('}','').replace('{','').split()
def help(self, cmd):
return self.eval('? %s'%cmd)
def __getattr__(self, attrname):
if attrname[:1] == "_":
raise AttributeError
return MathematicaFunction(self, attrname)
class MathematicaElement(ExpectElement):
def __getitem__(self, n):
return self.parent().new('%s[[%s]]'%(self._name, n))
def __getattr__(self, attrname):
self._check_valid()
if attrname[:1] == "_":
raise AttributeError
return MathematicaFunctionElement(self, attrname)
def __float__(self):
P = self.parent()
return float(P.eval('N[%s,16]'%self.name()))
def _reduce(self):
return self.parent().eval('InputForm[%s]'%self.name())
def __reduce__(self):
return reduce_load, (self._reduce(), )
def _latex_(self):
z = self.parent().eval('TeXForm[%s]'%self.name())
i = z.find('=')
return z[i+1:].strip()
def __repr__(self):
P = self._check_valid()
return P.get(self._name, ascii_art=False).strip()
def _sage_(self, locals={}):
r"""
Attempt to return a Sage version of this object.
This method works successfully when Mathematica returns a result
or list of results that consist only of:
- numbers, i.e. integers, floats, complex numbers;
- functions and named constants also present in Sage, where:
- Sage knows how to translate the function or constant's name
from Mathematica's naming scheme, or
- you provide a translation dictionary `locals`, or
- the Sage name for the function or constant is simply the
Mathematica name in lower case;
- symbolic variables whose names don't pathologically overlap with
objects already defined in Sage.
This method will not work when Mathematica's output includes:
- strings;
- functions unknown to Sage that are not specified in `locals`;
- Mathematica functions with different parameters/parameter order to
the Sage equivalent. In this case, define a function to do the
parameter conversion, and pass it in via the locals dictionary.
EXAMPLES:
Mathematica lists of numbers/constants become Sage lists of
numbers/constants::
sage: m = mathematica('{{1., 4}, Pi, 3.2e100, I}') # optional - mathematica
sage: s = m.sage(); s # optional - mathematica
[[1.0, 4], pi, 3.2*e100, I]
sage: s[1].n() # optional - mathematica
3.14159265358979
sage: s[3]^2 # optional - mathematica
-1
::
sage: m = mathematica('x^2 + 5*y') # optional - mathematica
sage: m.sage() # optional - mathematica
x^2 + 5*y
::
sage: m = mathematica('Sin[Sqrt[1-x^2]] * (1 - Cos[1/x])^2') # optional - mathematica
sage: m.sage() # optional - mathematica
(cos(1/x) - 1)^2*sin(sqrt(-x^2 + 1))
::
sage: m = mathematica('NewFn[x]') # optional - mathematica
sage: m._sage_(locals={'NewFn': sin}) # optional - mathematica
sin(x)
::
sage: var('bla') # optional - mathematica
bla
sage: m = mathematica('bla^2') # optional - mathematica
sage: bla^2 - m.sage() # optional - mathematica
0
::
sage: m = mathematica('bla^2') # optional - mathematica
sage: mb = m.sage() # optional - mathematica
sage: var('bla') # optional - mathematica
bla
sage: bla^2 - mb # optional - mathematica
0
AUTHORS:
- Felix Lawrence (2010-11-03): Major rewrite to use ._sage_repr() and
sage.calculus.calculus.symbolic_expression_from_string() for greater
compatibility, while still supporting conversion of symbolic
expressions.
"""
from sage.symbolic.pynac import symbol_table
from sage.symbolic.constants import constants_name_table as constants
from sage.calculus.calculus import symbolic_expression_from_string
from sage.calculus.calculus import _find_func as find_func
res = self._sage_repr()
if '"' in res:
raise NotImplementedError, "String conversion from Mathematica \
does not work. Mathematica's output was: %s" % res
lsymbols = symbol_table['mathematica'].copy()
lsymbols.update(locals)
autotrans = [ str.lower,
_un_camel,
lambda x: x
]
p = re.compile('(?<!\.)[a-zA-Z]\w*')
for m in p.finditer(res):
if m.group() in lsymbols:
pass
elif m.end() < len(res) and res[m.end()] == '(':
for t in autotrans:
f = find_func(t(m.group()), create_when_missing = False)
if f != None:
lsymbols[m.group()] = f
break
else:
raise NotImplementedError, "Don't know a Sage equivalent \
for Mathematica function '%s'. Please specify one \
manually using the 'locals' dictionary" % m.group()
else:
for t in autotrans:
if t(m.group()) in constants:
lsymbols[m.group()] = constants[t(m.group())]
break
try:
return symbolic_expression_from_string(res, lsymbols,
accept_sequence=True)
except:
raise NotImplementedError, "Unable to parse Mathematica \
output: %s" % res
def __str__(self):
P = self._check_valid()
return P.get(self._name, ascii_art=True)
def __len__(self):
"""
Return the object's length, evaluated by mathematica.
EXAMPLES::
sage: len(mathematica([1,1.,2])) # optional - mathematica
3
AUTHORS:
- Felix Lawrence (2009-08-21)
"""
return self.Length()
def show(self, filename=None, ImageSize=600):
r"""
Show a mathematica expression or plot in the Sage notebook.
EXAMPLES::
sage: P = mathematica('Plot[Sin[x],{x,-2Pi,4Pi}]') # optional - mathematica
sage: show(P) # optional - mathematica
sage: P.show(ImageSize=800) # optional - mathematica
sage: Q = mathematica('Sin[x Cos[y]]/Sqrt[1-x^2]') # optional - mathematica
sage: show(Q) # optional - mathematica
<html><div class="math">\frac{\sin (x \cos (y))}{\sqrt{1-x^2}}</div></html>
"""
P = self._check_valid()
if P.eval('InputForm[%s]' % self.name()).strip().startswith('Graphics['):
if filename is None:
filename = graphics_filename()
orig_dir = P.eval('Directory[]').strip()
P.chdir(os.path.abspath("."))
s = 'Export["%s", %s, ImageSize->%s]'%(filename, self.name(), ImageSize)
P.eval(s)
P.chdir(orig_dir)
else:
print '<html><div class="math">%s</div></html>' % self._latex_()
def str(self):
return str(self)
def __cmp__(self, other):
P = self.parent()
if P.eval("%s < %s"%(self.name(), other.name())).strip() == 'True':
return -1
elif P.eval("%s > %s"%(self.name(), other.name())).strip() == 'True':
return 1
elif P.eval("%s == %s"%(self.name(), other.name())).strip() == 'True':
return 0
else:
return -1
def N(self, *args):
"""
EXAMPLES::
sage: mathematica('Pi').N(10) # optional -- mathematica
3.1415926536
sage: mathematica('Pi').N(50) # optional -- mathematica
3.14159265358979323846264338327950288419716939937511
"""
return self.parent().N(self, *args)
class MathematicaFunction(ExpectFunction):
def _sage_doc_(self):
M = self._parent
return M.help(self._name)
class MathematicaFunctionElement(FunctionElement):
def _sage_doc_(self):
M = self._obj.parent()
return M.help(self._name)
mathematica = Mathematica(script_subdirectory='user')
def reduce_load(X):
return mathematica(X)
import os, sys
def mathematica_console(readline=True):
if not readline:
os.system('math')
return
f1 = os.popen('math ', 'w')
f1.flush()
try:
while True:
sys.stdout.write('')
try:
line = raw_input(' ')
f1.writelines(line+'\n')
f1.flush()
except KeyboardInterrupt:
f1.close()
break
except EOFError:
pass
sys.stdout.write('\n')