Basic overview
An interface between Julia and the SymPy library of Python requires a connection between the two languages. The PythonCall package provides a means to call into an underlying sympy library in Python. (As well as PyCall.) For example
julia> import PythonCalljulia> _sympy_ = PythonCall.pyimport("sympy");
The _sympy_ object holds references to the underlying sympy library. As an example, the following creates a symbolic variable, x, and calls the sin function on it:
julia> x = _sympy_.symbols("x")Python: xjulia> _sympy_.sin(x)Python: sin(x)julia> x.is_commutativePython: Truejulia> x.conjugate()Python: conjugate(x)julia> typeof(x)PythonCall.Core.Py
The PythonCall package provides some basic operations for Py objects, such as basic math operations:
julia> x + xPython: 2*xjulia> x * xPython: x**2julia> x ^ xPython: x**x
SymPyPythonCall, which uses SymPyCore, wraps the Py objects in its Sym class to provide a means to dispatch familiar Julia generics:
julia> using SymPyPythonCall
julia> x = symbols("x") # or @syms xxjulia> simplify(2sin(x)*cos(x))sin(2⋅x)
By wrapping the Python object in a Julia struct, there are many advantages, such as:
- The type can be used for dispatch, allowing
Juliageneric functions to have methods specialized for symbolic objects. - The
getpropertymethod can be specialized. This allows object methods, likex.conjugateto have arguments translated fromSymobjects toPythonobjects when being called. - The
showmethod can be used to adjust printing of objects
The package also provides methods for some sympy methods, such as simplify above. To make this work, there needs to be a means to take Sym objects to their Py counterparts and a means to take Py objects to a symbolic type. As these conversions may be type dependent two operators (↓ and ↑) are used internally to allow the definition along these lines:
simplify(x::Sym, args...; kwargs...) = ↑(sympy.simplify(↓(x), ↓(args)...; ↓(kwargs)...))(Though there is some overhead introduced, it does not seem to be significant compared to computational cost of most symbolic computations.)
The expand_log function is not wrapped as such, but can still be called from the sympy object exported by SymPyPythonCall:
julia> @syms x::real(x,)julia> simplify(log(2x))log(2⋅x)julia> sympy.expand_log(log(2x))log(x) + log(2)
Methods of sympy are also called using the conversion operators above.
Some functions of SymPy are Julia methods
A few select functions of SymPy are available as Julia methods where typically the dispatch is based on the first element being symbolic. Other methods must by qualified, as in sympy.expand_log. Historically, SymPy was pretty expansive with what it chose to export, but that is no longer the case. Currently, only a few foundational function are exported. E.g., expand but not expand_log; simplify but not trigsimp. This should reduce name collisions with other packages. The drawback, is we can add more Julian interfaces when a method is defined. Open an issue, should you think this advantageous.
Using other SymPy modules
We now show how to access one of the numerous external modules of sympy beyond those exposed immediately by SymPy. In this case, the stats module.
julia> _stats_ = PythonCall.pyimport("sympy.stats");
The stats module holds several probability functions, similar to the Distributions package of Julia. This set of commands creates a normally distributed random variable, X, with symbolic parameters:
julia> 𝑋,mu = _sympy_.symbols("X, mu")Python: (X, mu)julia> sigma = _sympy_.symbols("sigma", positive=true)Python: sigmajulia> X = _stats_.Normal(𝑋, mu, sigma)Python: Xjulia> _stats_.E(X)Python: mujulia> _stats_.E(X^2)Python: mu**2 + sigma**2julia> _stats_.variance(X)Python: sigma**2
The one thing to note is the method calls return Py objects, as there is no intercepting of the method calls done the way there is for the sympy module. Wrapping _stats_ in Sym uses the getproperty specialization:
julia> stats = Sym(_stats_);julia> @syms 𝑋, μ, σ::positive(𝑋, μ, σ)julia> X = stats.Normal(𝑋, μ, σ)𝑋julia> stats.variance(X)2 σ
Next statements like $P(X > \mu)$ can be answered by specifying the inequality using Gt in the following:
julia> stats.P(Gt(X, μ))1/2
The unicode ≧ operator (\geqq[tab]) is an infix alternative to Gt.
A typical calculation for the normal distribution is the area one or more standard deviations larger than the mean:
julia> stats.P(X ≧ μ + 1 * σ)
sqrt(2)*(-sqrt(2)*pi*exp(1/2)*erf(sqrt(2)/2)/2 + sqrt(2)*pi*exp(1/2)/2)*exp(-1/2)/(2*pi)The familiar answer could be found by calling N or evalf.
One more distribution is illustrated, the uniform distribution over a symbolic interval $[a,b]$:
julia> @syms 𝑈 a::real b::real(𝑈, a, b)julia> U = stats.Uniform(𝑈, a, b)𝑈julia> stats.E(U)a b ─ + ─ 2 2julia> stats.variance(U) |> factor2 (a - b) ──────── 12
The mpmath library
According to its website: mpmath is a free (BSD licensed) Python library for real and complex floating-point arithmetic with arbitrary precision. It has been developed by Fredrik Johansson since 2007, with help from many contributors. For Sympy, versions prior to 1.0 included mpmath, but SymPy now depends on itmpmath as an external dependency.
The mpmath library provides numeric, not symbolic routines. To access these directly, the mpmath library can be loaded, in the manner just described.
julia> _mpmath_ = PythonCall.pyimport("mpmath")Python: <module 'mpmath' from '/home/runner/work/SymPyCore.jl/SymPyCore.jl/docs/.CondaPkg/env/lib/python3.12/site-packages/mpmath/__init__.py'>julia> mpmath = Sym(_mpmath_)<module 'mpmath' from '/home/runner/work/SymPyCore.jl/SymPyCore.jl/docs/.Conda ↪ ↪ Pkg/env/lib/python3.12/site-packages/mpmath/__init__.py'>
To compare values, consider the besselj function which has an implementation in Julia's SpecialFunctions package, SymPy itself, and mpmath. SymPyCore provides an overloaded method for a symbolic second argument which is available once SpecialFunctions is loaded.
julia> using SpecialFunctionsjulia> @syms x(x,)julia> nu, c = 1//2, pi/3(1//2, 1.0471975511965976)julia> (mpmath.besselj(nu,c), besselj(nu, x)(c).evalf(), besselj(nu, c))(0.67523723711783, 0.675237237117830, 0.6752372371178303)
Different output types
SymPyPythonCall provides a few conversions into containers of symbolic objects, like for lists, tuples, finite sets, and matrices . Not all outputs are so simple to incorporate and are simply wrapped in the Sym type.
Conversion to a workable Julia structure can require some massaging. This example shows how to get at the pieces of the Piecewise type.
The output of many integration problems is a piecewise function:
julia> @syms n::integer x::real(n, x)julia> u = integrate(x^n, x)⎧ n + 1 ⎪x ⎪────── for n ≠ -1 ⎨n + 1 ⎪ ⎪log(x) otherwise ⎩
The u object is of type Sym, but there are no methods for working with it. The .args call will break this into its argument, which again will by symbolic. The Introspection.args function will perform the same thing.
julia> as = Introspection.args(u)((x^(n + 1)/(n + 1), Ne(n, -1)), (log(x), True))
The as object is a tuple of Sym objects. Consider the first one:
julia> c1 = first(as)⎛ n + 1 ⎞ ⎜x ⎟ ⎜──────, n ≠ -1⎟ ⎝n + 1 ⎠
The value of c1 prints as a tuple, but is of type Sym and sympy type CondExprPair. Though the underlying python type is iterable or indexable, the wrapped type is not. It can be made iterable in a manner of ways: by calling ↓(c1); by finding the underlying Python object through the attribute .o, as in c1.o; or by calling the Py method of PythonCall, as in PythonCall.Py(c1). More generically, PythonCall provides a PyIterable wrapper. As there is no defined lastindex, the latter is a bit more cumbersome to use. This pattern below, expands the conditions into a dictionary:
julia> [Sym(↓(a)[1]) => Sym(↓(a)[0]) for a ∈ as]2-element Vector{Pair{SymPyCore.Sym{PythonCall.Core.Py}, SymPyCore.Sym{PythonCall.Core.Py}}}: Ne(n, -1) => x^(n + 1)/(n + 1) True => log(x)
The Python object, ↓(a) is indexed, so 0-based indexing is used above. These pieces are then converted to Sym objects for familiarity.