Blog‎ > ‎

Converting Mathematica to Python

posted Sep 19, 2017, 3:46 AM by Juan Jose Garcia-Ripoll   [ updated Sep 19, 2017, 3:51 AM ]
I like Mathematica for working with symbolic expressions. I find that there are no free software alternatives that achieve its power in the tasks of solving complex integrals, working with limits, and helping me make sense of sophisticated differential equations. Until that is solved, I will find the need to export Mathematica expressions to other software that is more efficient on the numerical computing side. Or which simply does good plots in a less verbose way.

I have created a simple Mathematica notebook which you can download here (see attachment below). The notebook uses Mathematica's pattern matching to build a sophisticated parsing routine that converts Mathematica expressions to Python code.

See for instance the following examples which detect addition, division and exponentiation, and convert those to strings, either directly, or using wrappers that automate the case of binary operators and unitary expressions:

ToPython2[a_ - b_] := ToPythonOp["-", a, b]
ToPython2[a_/b_] := ToPythonOp["/", a, b]
ToPython2[a_^2] := "(" <> PythonWrap[a] <> "**2)"
ToPython2[a_^b_] := PythonWrap[ToPythonOp["**", a, b]]

The actual implementation makes use of ugly hacks, such as defining its own square root function (which accepts negative real arguments), or parenthesising more than needed, but it is fully functional and works nicely for my own purposes. For instance, the following complex Mathematica expression


is converted to 13 lines of Python code:

    π=math.pi;

    def mysqrt(x): return np.sqrt((1.+0j)*x)

    y=π*γ;
    x=mysqrt((-1.+(y**2)));
    z=1.+((2.*(x*y))+(-2.*(y**2)));
    w=-1.+((2.*(x*y))+(2.*(y**2)));
    u=mysqrt((z/w));
    aux0=((Δ**5.)*(np.log(((z**-0.5)/Ω))))+((x-y)*(y*((mysqrt(z))*((Ω**5.)*(np.log((w*(Ω**2))))))));
    aux1=((g**2)*(((2.*(π*(u*(Δ**5.))))+(w*(Ω**3.)))*(np.log(((-z**-0.5)/Ω)))))+(2.*aux0);
    aux2=(((Δ**4.)+((2.*((-1.+(2.*(y**2)))*((Δ**2)*(Ω**2))))+(Ω**4.)))**-2.)*aux1;
    output=np.exp((((-((-1.+(y**2))**-0.5)*((z**-0.5)*(Ω*aux2))))/π));

ċ
Export to Python.nb
(40k)
Juan Jose Garcia-Ripoll,
Sep 19, 2017, 3:46 AM