SimpleExpressions.jl#

SimpleExpressions

A very lightweight means to create callable functions using expressions. This uses CallableExpressions as a backend. See also DynamicExpressions for a performant package with similar abilities.

The @symbolic macro, the lone export, can create a symbolic variable and optional symbolic parameter. When expressions are created with these variables, evaluation is deferred until the expression is called like a function. The expressions subtype Function so are intended to be useful with Julia's higher-order functions.

The expressions can be evaluated as a univariate function, u(x), a univariate function with parameter, u(x, p), or as a bivariate function, u(x,y) (with y being a parameter). These are all typical calling patterns when a function is passed to a numeric routine. For expressions without a symbolic value (as can happen through substitution) u() will evaluate the value.

To substitute in for either the variable or the parameter, leaving a symbolic expression, we have the calling patterns u(:,p), u(x,:) to substitute in for the parameter and variable respectively. The colon can also be nothing or missing.

When using positional arguments in a call, as above, all symbolic variables are treated identically, as are all symbolic parameters.

There are also methods for replace that allow more complicated substitutions. For replace, symbolic objects are returned. For replace, variables are distinct and identified by their symbol. Pairs may be specified to the call notation as a convenience for replace.

There are no performance claims, this package is all about convenience. Similar convenience is available in some form with SymPy, Giac, SymEngine, Symbolics, etc. As well, placeholder syntax is available in Underscores.jl, Chain.jl, DataPipes.jl etc., This package only has value in that it is very lightweight and, hopefully, intuitively simple.

Performance is good though, as CallableExpressions is performant. A benchmark case of finding a zero of a function runs without allocations in 931 ns with 0 allocations, with a symbolic expression in 4.65 μs with 286 allocations, SymEngine is two orders of magnitude slower (76.329 μs with 18911 allocations), and SymPy much slower (32.318 ms with 79796) allocations.

Extensions are provided for SpecialFunctions, AbstractTrees, Latexify, and RecipesBase.

Example

using SimpleExpressions
@symbolic x       # (x,)
u = sin(x) - (x - x^3/6)
u(0.5)  # 0.000258...
u = u - x^5/120
u(0.5) # -1.544...e-6

map(x^2, (1, 2))  # (1, 4)

using Plots
@symbolic x p     # (x, p)
u = x^5 - x - p   # (x ^ 5) + (-1 * x) + (-1 * p)
plot(u(:, 1), 0, 1.5)
plot!(u(:, 2))    # or plot(u.(:, 1:2), 0, 1.5)
eq = cos(x) ~ 2x
plot(eq, 0, pi/2) # like plot([eq...], 0, pi/2)

a ~ b

Create a SymbolicEquation.

The equation has a left and right-hand side, which can be found by tuple destructing; calling first and last; by index; or field access using .lhs and .rhs.

Symbolic equations can be evaluated, in which case the value of a-b is returned.

When a symbolic equation is passed as an argument to a symbolic expression, the pair a => b is passed to replace.

The diff function differentiates both sides.

The solve function tries to move x terms to the left-hand side; and non-x terms to the right-hand side.

isnumeric(x::AbstractSymbolic)

Check if symbolic value contains a variable or a parameter. Puns on isnumeric(::Char). Use

match(pattern::Expr, subject::AbstractSymbolic)::Union{MatchDict, Nothing}

Match subject against a pattern given as a Julia expression containing wildcards.

Return a dictionary mapping wildcard names (as symbols) to the matched values for the first match found, or nothing if there is none. Use eachmatch to get all identified matches.

Examples

julia> @symbolic x p
(x, p)

julia> match(:(~x * cos(~y)), p*cos(x))
Base.ImmutableDict{Symbol, SimpleExpressions.AbstractSymbolic} with 2 entries:
  :y => x
  :x => p

Extended help

The algorithm

The basic algorithm comes from that of rule2.jl from SymbolicIntegration and [Krebber](https://arxiv.org/pdf/1705.00907).

The pattern and the subject are walked together, top down.

• A literal in the pattern (a number, a symbol, a constant such as ℯ) matches only an equal value in the subject.

• A wildcard is bound to the part of the subject it is compared with. A wildcard that appears more than once must be bound to equal values each time, so :(~a + ~a) matches x + x but not x + y.

• A call in the pattern, such as cos(~y), matches a call in the subject with the same operation, whose arguments are then matched in turn. The pattern's operation can itself be a wildcard, see (~f)(~x) below.

• For a + or * pattern (which are commutative) the arguments of the subject may be matched in any order, so :(~x + 2) matches 2 + y. Different assignments can give different matches, which is why there can be more than one; match returns the first and eachmatch returns them all.

• Powers are matched up to their representation: for example sqrt(x) and x^(1//2) can match each other, and a pattern ~a / ~b matches x * (1 / y).

• Matching proceeds through a list of candidate bindings. A binding that conflicts with one already made, or that fails a predicate, is discarded; if no candidates remain the match fails.

Wildcards

A wildcard is written with a leading ~. A bare variable in the pattern is not a wildcard: it must match literally.

Slots and predicates

~x matches one subexpression. A predicate restricts the match; it is any function (or expression naming one) that returns a Bool, called on the value to be bound:

julia> match(:(~a::iseven * ~b), 2x)    # :a => 2, :b => x
julia> match(:(~a::iseven * ~b), 3x)    # nothing

An error thrown by a predicate counts as false.

Default slots

~!x is a slot that may be absent. In a sum its default is 0, in a product, a power, or a division it is 1. Thus :(~!a * ~b) matches x with a => 1 and b => x, and :((~b)^(~!n)) matches x with n => 1. If the term is present it is bound as usual. When the same default slot appears several times in a pattern, all occurrences must agree. Within a replacement, a default slot is written ~a once bound.

Segments

~~x and ~~~x stand for several arguments of a call, and are bound to a tuple of those arguments. ~~x allows none, ~~~x requires at least one:

julia> match(:(~x + ~~~y), x + y + z)   # :x => x, :y => (y, z)
julia> match(:(~x + ~~~y), x)           # nothing

If a segment is the only argument of the call, it is bound to all the arguments, so :(*(~~a)) matches (x + y) * z with a => (x + y, z). When there are several segments in one call the remaining arguments are divided among them, each ~~~ segment receiving at least one; several divisions may be possible, but not all are enumerated.

Operation wildcards

(~f)(~x) matches any call with one argument, binding f to the operation and x to the argument:

julia> match(:((~f)(~x)), sin(y))       # :f => sin, x => y

Caveats

• Matching is syntactic up to commutativity and the power equivalences above. It does not use other algebraic identities, so :(~x * ~x) does not match x^2.

• Matching does not consider associativity. For example, :(~x + ~y) does not match a + b + c, even though it could be argued to match (a+b) + c or a + (b + c). That matching is to expensive. The AssociativeCommutativePatternMatching implements an algorithm that does this matching.

• Pattern constants are compared by value after unwrapping, so 2 matches the symbolic number 2 (and 2.0).

replace(ex::SymbolicExpression, args::Pair...)

Replace parts of the expression with something else.

Returns a symbolic object.

The pattern/replacement is specified using variable => value; these are processed left to right.

There are different methods depending on the type of key in the the key => value pairs specified:

• A symbolic variable is replaced by the right-hand side, like ex(val,:), though the latter is more performant

• A symbolic parameter is replaced by the right-hand side, like ex(:,val)

• A function is replaced by the corresponding specified function, as the head of the sub-expression

• A sub-expression is replaced by the new expression.

• sub-expression containing a wildcard is replaced by the new expression. The specification of the replacements can use the wildcards. Wildcards expressions may be specified within SimpleExpressions or as Expr objects.

Extended help

The first two styles are straightforward.

julia> using SimpleExpressions

julia> @symbolic x p
(x, p)

julia> ex = cos(x) - x*p
cos(x) + (-1 * x * p)

julia> replace(ex, x => 2) == ex(2, :)
true

julia> replace(ex, p => 2) == ex(:, 2)
true

The third, is illustrated by:

julia> replace(sin(x + sin(x + sin(x))), sin => cos)
cos(x + cos(x + cos(x)))

The fourth is similar to the third, only an entire expression (not just its head) is replaced

julia> ex = cos(x)^2 + cos(x) + 1
(cos(x) ^ 2) + cos(x) + 1

julia> @symbolic u
(u,)

julia> replace(ex, cos(x) => u)
(u ^ 2) + u + 1

Replacements occur only if an entire node in the expression tree is matched:

julia> u = 1 + x
1 + x

julia> replace(u + exp(-u), u => x^2)
1 + x + exp(-1 * (x ^ 2))

(As this addition has three terms, 1+x is not a subtree in the expression tree.)

The fifth needs more explanation, as there can be wildcards in the expression. Wildcards can be specified as SimpleExpression expressions or as Expr objects.

First, we describe the use of symbolic wildcards.

Wildcards have a naming convention using trailing underscores:

• one matches a single subexpression or term;

• two matches zero, one, or more subexpressions.

• three matches one or more subexpressions.

In addition, the special symbol ⋯ (entered with \cdots[tab] is wild.

julia> @symbolic x p; @symbolic x_
(x_,)

julia> replace(cos(pi + x^2), cos(pi + x_) => -cos(x_))
-1 * cos(x ^ 2)

julia> ex = log(sin(x)) + tan(sin(x^2))
log(sin(x)) + tan(sin(x ^ 2))

julia> replace(ex, sin(x_) => tan((x_) / 2))
log(tan(x / 2)) + tan(tan((x ^ 2) / 2))

julia> replace(ex, sin(x_) => x_)
log(x) + tan(x ^ 2)

julia> replace(x*p, (x_) * x => x_)
p

Pattern and replacements can also be specified with Julia expressions. The details are in the docstring for match. This method is more expressive as patterns can use default slots and slots can use guards.

julia> ex = log(sin(x)) + tan(sin(x^2))
log(sin(x)) + tan(sin(x ^ 2))

julia> replace(ex, :(sin(~x)) => :(tan(~x/2)))
log(tan(x / 2)) + tan(tan((x ^ 2) / 2))

julia> replace(ex, :(sin((~x)^2)) => :(tan(~x)))
log(sin(x)) + tan(tan(x))

julia> replace(ex, :(sin((~x)^(~n::iseven))) => :(tan(~x)))

Unlike symbolic wildcards, Expr objects can have default slot (specified as ~!x) and predicates (specified after :: to test.

julia> pat, replacement = :(~!a * (~x)^(~n::(!=(-1)))), :(~a * (~x)^(~n+1)/(~n + 1))
(:(~(!a) * (~x) ^ ~(n::(!=)(-1))), :((~a * (~x) ^ (~n + 1)) / (~n + 1)))

julia> replace(x^2, pat => replacement)
(x ^ 3) / 3

julia> replace(x^(-1), pat => replacement)
1 / x

Picture

The AbstractTrees package can print this tree-representation of the expression ex = sin(x + x*log(x) + cos(x + p + x^2)):

julia> print_tree(ex;maxdepth=10)
sin
└─ +
   ├─ x
   ├─ *
   │  ├─ x
   │  └─ log
   │     └─ x
   └─ cos              <--
      └─ +             ...
         ├─ x          <--
         ├─ p          ...
         └─ ^          ...
            ├─ x       ...
            └─ 2       ...

The command wildcard expression cos(x + ...) looks at the part of the tree that has cos as a node, and the lone child is an expression with node + and child x. The ⋯ then matches p + x^2.

solve(eq::SymbolicEquation, x; conditions::Bool=false, interval=nothing)

Very simple symbolic equations can be solved with the unexported solve method.

• eq::SymbolicEquation:. equations are specified with ~, as in x^2 ~ p

• x: variable to solve for

• conditions: set to true to reveal conditions on parameters for which the identified solution(s) hold.

• interval: If specified, through a tuple (a,b) all solutions in this interval are found.

Return value:

• an equation x ~ r when there is a single solution;

• a tuple of equations, e.g. (x ~ r₁, x ~ r₂), when there are more (quadratics, even powers);

• () when numeric values show there are no real solutions;

• nothing when x does not appear in the equation (or cancels out);

• an equation whose left-hand side is not x when x could not be isolated.

Examples

import SimpleExpressions: solve
SimpleExpressions.@symbolic_variables A a B b # (A, a, B, b)
solve(sin(A)/a ~ sin(B)/b, b) # b ~ (a * sin(B)) / sin(A)
solve(sin(A)/a ~ sin(B)/b, A) # A ~ asin((sin(B) * a) / b)

SimpleExpressions.@symbolic_variables w p h  # (w, p, h)
constraint = p ~ 2w + 2h;
A = w * h;
u = solve(constraint, h)      # h ~ (p / 2) + (-1 * w)
A = A(u)                      # w * ((p / 2) + (-1 * w)), use equation in replacement
v = solve(diff(A, w) ~ 0, w)  # w ~ p / 4

Extended help

The equation is rearranged by moving terms free of x to the right-hand side and applying inverse functions (e.g. exp(u) ~ r becomes u ~ log(r), abs(u) ~ r becomes u ~ ±r). A product equal to 0 is solved factor by factor. Polynomials in x are solved when they are linear or quadratic, after factoring out powers of x, when they are a polynomial in xᵏ (e.g. x^4 - 5x^2 + 4), or when their numeric rational coefficients allow roots to be found by the rational root theorem.

Answers are kept exact where possible: rational values stay rational and other values are left as symbolic numbers, e.g. x^2 ~ 8 gives x ~ ±2⋅sqrt(2), exp(x) ~ 2 gives x ~ log(2), 2^x ~ 8 gives x ~ 3 and sin(x) ~ 1//2 gives x ~ π/6 (standard values of asin, acos and atan are recognized). Calling a symbolic number, e.g. solve(x^2 ~ 2, x)[2].rhs(), gives its floating point value. When the equation has floating point values, so does the answer.

Some further rewriting is done:

• an injective function on both sides is removed, e.g. log(u) ~ log(v) becomes u ~ v; exponentials on both sides are compared through their logarithms, e.g. 2^x ~ 3^(x+1);

• sums of logarithms are combined, log(u) + log(v) ~ r becomes log(u⋅v) ~ r;

• rational functions, e.g. 1/x + 1/(x+1) ~ 1, are solved through their numerator, discarding roots of the denominator;

• a polynomial in some g(x), e.g. exp(2x) - 3exp(x) + 2 ~ 0 or sqrt(x) + x ~ 6, is solved for g(x) first;

• logarithms in any base (log, log10, log2, log(b, u)) are combined, Σ kᵢ⋅log_b(uᵢ) ~ r becomes Π uᵢ^kᵢ ~ b^r;

• a square root is isolated and both sides squared, e.g. 2x ~ sqrt(x+3);

• abs(u) is split into the cases u ≥ 0 and u ≤ 0;

• a polynomial factor common to all terms is factored out, e.g. 2x - x*exp(6x-1) ~ 0;

• a⋅cos(u) + b⋅sin(u) ~ 0 becomes tan(u) ~ -a/b;

• a sum of exponentials of one sign equal to 0 has no solutions;

• equations mixing x with an exponential or logarithm of a linear function of x, e.g. x*exp(x) ~ 1, 2^x ~ x + 1, x + log(x) ~ 2 or x^x ~ 4, are solved with the Lambert W function: lambertw(z) and, for -1/e ≤ z < 0, the second real branch lambertw(z, -1) (from LambertW.jl). Simple values are recognized, e.g. x^x ~ 4 gives x ~ 2.

When the equation has no symbolic parameters, numeric solutions are checked in the original equation and extraneous ones are discarded.

With symbolic parameters (or other variables), each solution is valid only under some conditions: logarithms of positive values, square roots of nonnegative ones, nonzero denominators, and conditions from squaring or from splitting abs. These are evaluated at the solution: a solution failing a numeric condition is discarded, as is one which fails at every sampled value of the parameters. Conditions holding at every sampled value are dropped (this is a heuristic, not a proof). Pass conditions=true to have each solution returned as a pair eq => conds of the equation and a tuple of the remaining conditions, e.g.

solve(sqrt(x) ~ p, x; conditions=true)     # x ~ p^2 => (p >= 0,)
solve(abs(x - p) ~ 2x, x; conditions=true) # (x ~ -p => (p <= 0,), x ~ p/3 => (p >= 0,))

The check uses Julia's definitions, so e.g. x^(1//3) is not defined for negative x ((-8.0)^(1//3) is a DomainError); use cbrt(x) for the real cube root.

Pass interval=(a, b) to get every real solution with a ≤ x ≤ b (a closed interval). Periodic functions then give all their branches, e.g. sin(u) ~ r gives u = asin(r) + 2πn and u = π - asin(r) + 2πn, and those in the interval are listed:

solve(4sin(3x) ~ 2, x; interval=(0, 4π/3))     # (π/18, 5π/18, 13π/18, 17π/18)
solve(x^2 - 5x + 6 ~ 0, x; interval=(0, 2.5))  # x ~ 2; other equations are filtered

The equation must not have parameters when interval is given. Without interval, inverse trigonometric functions give principal values.

isconstant(x)::Bool

Check if x has a symbolic variable. (Symbolic parameters are treated as constant for this method.)

isvariable(x)::Bool

Check if x is a symbolic variable or parameter

map_matched(ex, is_match, f)

Traverse expression. If is_match is true, apply f to that part of expression tree and reassemble.

Basically CallableExpressions.expression_map_matched.

Not exported.

solve_equation(l, r, x; conditions=false)

Solve l ~ r for x. This is the generic solver behind solve; see its docstring for what is solved and what is returned. It walks expressions with TermInterface, builds them with ordinary operators, and uses the hooks of solve_interface.jl for everything else.

sympify(ex::String)

Turn a string into an expression, converting symbols to symbolic variables.

Name comes from SymPy.

@symbolic x [p]

Create a symbolic variable and optional symbolic parameter.

Expressions created using these variables subclass Function so may be used where functions are expected.

The ~ infix operator can be used to create equations, which, by default, are treated as lhs - rhs when called as functions.

Extended help

Calling or substituting into expressions

To call a symbolic expression regular call notation with positional arguments are used. The first argument maps to any symbolic variable; the second – when given – to any symbolic parameter. It is an error to call an expression with a parameter using just a single argument; for that substitution is needed.

Example

using SimpleExpressions
@symbolic x p
u = x^5 - x - 1
u(2) # 29 call is u(x)

u.((0,1,2)) # (-1, -1, 29)

u = 2x + p
u(1)    # errors!
u(1, 2) # 2(1)+2 or 4

u = sum(x .* p)
u(2, [1,2]) # 6  call is u(x, p)

Calling with nothing, missing, or : in a slot substitutes in the specified value leaving a symbolic expression, possibly with no variable or parameter.

@symbolic x p
u = cos(x) - p*x
u(nothing, 2)  # cos(x) - 2 * x
u(:, 2)        #  cos(x) - 2 * x, alternate calling form
u(pi, nothing) # -1.0 - p * π
v = u(1,:)(:,2)    # (cos(1)-(2*1)),

The latter can be evaluated using a zero-argument call, e.g. v().

With substitution in this manner, any symbolic variable and any symbolic parameters will receive the same substituted value.

The replace generic for symbolic objects takes pairs of values and replaces the left one with the right one working from left to right, leaving a symbolic expression. The replace method treats symbolic variables and symbolic parameters with different symbols as unique.

A symbolic equation, defined through ~, may also be used to specify a left- and right-hand value.

The main use is as an easier-to-type replacement for anonymous functions, though with differences:

1 |> sin(x) |> x^2  # 0.708… from sin(1)^2
u = cos(x) - p*x
2 |> u(:, 3) # -6.4161…, a alternative to u(2,3)

map(x^2, (1, 2)) # (1,4)

Symbolic expressions an be used with other packages, to simplify some function calls at the expense of being non-idiomatic:

using Roots
@symbolic x p
find_zero(x^5 - x - 1, 1)       # 1.167…
find_zero(x^5 - x ~ p, 1; p=4)  # 1.401…

using ForwardDiff
Base.adjoint(𝑓::Function) = x -> ForwardDiff.derivative(𝑓, x)
u = x^5 - x - 1
find_zero((u,u'), 1, Roots.Newton()) # 1.167…

Or

using Plots
plot(x^5 - x - 1, 0, 1.5)

Or using both positions, so that we call as a bivariate function:

@symbolic x y
xs = ys = range(-5, 5, length=100)
contour(xs, ys, x^2 - y^2 + 2x*y)

Symbolic derivatives can be taken with respect to the symbolic value, symbolic parameters are treated as constant. We use diff as an interface

@symbolic x p
u = x^5 - p*x - 1
diff(u, x)           # (5 * (x ^ 4)) - p
u = u(:, 1)    # set parameter
a, b = 1, 2
find_zeros(diff(u,x) ~ (u(b)-u(a)) / (b-a), (a,b)) # [1.577…]

Idiosyncrasies

Using this is a convenience for simple cases. It is easy to run into idiosyncrasies.

Expressions are not functions in terms of scope

Unlike functions, expressions are defined with variables at the time of definition, not when called. For example, with a clean environment:

@symbolic x
u = m*x + b    # errors, `m` not defined
f(x) = m*x + b # ok
m, b = 1, 2
u = m*x + b    # defined using `m` amd `b` at time of assignment
u(3)           # 1 * 3 + 2
f(3)           # 1 * 3 + 2 values of `m` and `b` when called
m, b = 3, 4
u(3)           # still computing 1 * 3 + 2
f(3)           # computing 3 * 3 + 4, using values of `m` and `b` when called

Symbolic values are really singletons when calling by position

Though one can make different symbolic variables, the basic call notation by position treats them as the same:

@symbolic x
@symbolic y    # both x, y are `SymbolicVariable` type
u = x + 2y
u(3)           # 9 coming from 3 + 2*(3)

However, this is only to simplify the call interface. Using keyword arguments allows evaluation with different values:

u(;x=3, y=2)   # 7

Using replace, we have:

u(x=>3, y=>2)  # 3 + (2 * 2); evaluate with u(x=>3, y=>2)()

The underlying CallableExpressions object is directly called in the above manner; that package does not have the narrowed design of this package.

Containers

The variables may be used as placeholders for containers, e.g.

u = sum(xi*pi for (xi, pi) in zip(x,p))
u((1,2),(3,4))  # 11

Broadcasting as a function

Broadcasting a function call works as expected

@symbolic x
u = x^2
u.((1,2)) # (1, 4)

Symbolic expressions can also be constructed that will broadcast the call

u = x.^2 .+ sin.(p)
u((1,2),3)

u = @. x^2 + sin(p)
u((1,2),(3,4))

@symbolic_variables w x[1:3] y() z=>"𝑧" Ω::isinteger

Define multiple symbolic variables or symbolic functions. Guards are currently ignored.

Not exported.

julia> using SimpleExpressions: @symbolic_variables

julia> @symbolic_variables  w x[1:3] y() z=>"𝑧" Ω::isinteger
(w, SimpleExpressions.SymbolicVariable[x₁, x₂, x₃], y, 𝑧, Ω)

julia> x
3-element Vector{SimpleExpressions.SymbolicVariable}:
 x₁
 x₂
 x₃