Variant structure and types

This document aims to describe the structure of the Algebraic Data Type (ADT) used to represent a symbolic tree, along with several utility types to allow robustly interacting with it.

SymbolicUtils uses Moshi.jl's ADT structure. The ADT is named BasicSymbolicImpl, and an alias BSImpl is available for convenience. The actual type of a variable is BSImpl.Type, aliased as BasicSymbolic. A BasicSymbolic is considered immutable. Mutating its fields is unsafe behavior.

In SymbolicUtils v3, the type T in BasicSymbolic{T} was the type represented by the symbolic variable. In other words, T was the symtype of the variable.

In SymbolicUtils v4, the symtype is not stored in the type, and is instead a field of the struct. This allows for greatly increased type-stability. The type T in BasicSymbolic{T} now represents a tag known as thw vartype. This flag determines the assumptions made about the symbolic algebra. It can take one of three values:

SymReal: The default behavior.
SafeReal: Idential to SymReal, but common factors in the numerator and denominator of a division are not cancelled.
TreeReal: Assumes nothing about the algebra, and always uses the Term variant to represent an expression.

A given expression must be pure in its vartype. In other words, no operation supports operands of different vartypes.

A short note on (im-)mutability

While ismutabletype(BasicSymbolic) returns true, symbolic types are IMMUTABLE. Any mutation is undefined behavior and can lead to very confusing and hard-to-debug issues. This includes internal mutation, such as mutating AddMul.dict. The arrays returned from TermInterface.arguments and TermInterface.sorted_arguments are read-only arrays for this reason.

Expression symtypes

The "symtype" of a symbolic variable/expression is the Julia type that the variable/expression represents. It can be queried with SymbolicUtils.symtype. Note that this query is unstable - the returned type cannot be inferred.

Expression shapes

In SymbolicUtils v4, arrays are first-class citizens. This is implemented by storing the shape of the symbolic. The shape can be queried using SymbolicUtils.shape and is one of two types.

Symbolics with known shape

The most common case is when the shape of a symbolic variable is known. For example:

@syms x[1:2] y[-3:6, 4:7] z

All of the variables created above have known shape. In this case, SymbolicUtils.shape returns a (custom) vector of UnitRange{Int} semantically equivalent to Base.axes. This does not return a Tuple since the number of dimensions cannot be inferred and thus returning a tuple would introduce type-instability. All array operations will perform validation on the shapes of their inputs (e.g. matrix multiplication) and calculates the shape of their outputs.

Scalar variables return an empty vector as their shape.

Symbolics with known `ndims`

The next most common case is when the exact shape/size of the symbolic is unknown but the number of dimensions is known. For example:

@syms x::Vector{Number} y::Matrix{Number} z::Array{Number, 3}

In this case, SymbolicUtils.shape returns a value of type SymbolicUtils.Unknown. This has a single field ndims::Int storing the number of dimensions of the symbolic. Note that a shape of SymbolicUtils.Unknown(0) does not represent a scalar. All array operations will perform as much validation as possible on their arguments. The shape of the result will be calculated on a best-effort basis.

Symbolics with unknown `ndims`

In this case, nothing is known about the symbolic except that it is an array. For example:

@syms x::Array{Number}

Symbolics.shape(x) will return SymbolicUtils.Unknown(-1). This effectively disables most shape checking for array operations.

Variants

struct Const
    const val::Any
    # ...
end

Any non-symbolic values in an expression are stored in a Const variant. This is crucial for type-stability, but it does mean that obtaining the value out of a Const is unstable and should be avoided. This value can be obtained by pattern matching using Moshi.Match.@match or using the unwrap_const utility. unwrap_const will act as an identity function for any input that is not Const, including non-symbolic inputs. Const is the only variant which does not have metadata.

SymbolicUtils.isconst can be used to check if a BasicSymbolic is the Const variant. This variant can be constructed using Const{T}(val) or BSImpl.Const{T}(val), where T is the appropriate vartype.

The Const constructors have an additional special behavior. If given an array of symbolics (or array of array of ... symbolics), it will return a Term (see below) with hvncat as the operation. This allows standard symbolic operations (such as substitute) to work on arrays of symbolics without excessive special-case handling and improved type-stability.

struct Sym
    const name::Symbol
    const metadata::MetadataT
    const shape::ShapeT
    const type::TypeT
    # ...
end

Sym represents a symbolic quantity with a given name. This and Const are the two atomic variants. metadata is the symbolic metadata associated with this variable. type is the tag for the type of quantity represented here. shape stores the shape if the variable is an array symbolic.

metadata is either nothing or a map from DataType keys to arbitrary values. Any

interaction with metadata should be done by providing such a mapping during construction or using getmetadata, setmetadata, hasmetadata.

type is a Julia type.
shape is as described above.

These three fields are present in all subsequent variants as well.

A Sym can be constructed using Sym{T}(name::Symbol; type, shape, metadata) or BSImpl.Sym{T}(name::Symbol; type, shape, metadata).

struct Term
    const f::Any
    const args::SmallV{BasicSymbolicImpl.Type{T}}
    const metadata::MetadataT
    const shape::ShapeT
    const type::TypeT
    # ...
end

Term is the generic expression form for an operation f applied to the arguments in args. In other words, this represents f(args...). Any constant (non-symbolic) arguments (including arrays of symbolics) are converted to symbolics and wrapped in Const.

A Term can be constructed using Term{T}(f, args; type, shape, metadata) or BSImpl.Term{T}(f, args; type, shape, metadata).

struct AddMul
    const coeff::Any
    const dict::ACDict{T}
    const variant::AddMulVariant.T
    const metadata::MetadataT
    const shape::ShapeT
    const type::TypeT
    # ...
end

AddMul is a specialized representation for associative-commutative addition and multiplication. The two operations are distinguised using the AddMulVariant EnumX.jl enum. It has two variants: AddMulVariant.ADD and AddMulVariant.MUL.

For multiplication terms, coeff is a constant non-symbolic coefficient multipled with the expression. dict is a map from terms being multiplied to their exponents. For example, 2x^2 * (y + z)^3 is represented with coeff = 2 and dict = ACDict{T}(x => 2, (y + z) => 3). A valid multiplication term is subject to the following constraints:

coeff must be non-symbolic.
The values of dict must be non-symbolic.
The keys of dict must not be expressions with ^ as the operation UNLESS the exponent is symbolic. For example, x^x * y^2 is represented with dict = ACDict{T}((x^x) => 1, y => 2).
dict must not be empty.
coeff must not be zero.
If dict has only one element, coeff must not be one. Such a case should be represented as a power term (with ^ as the operation).
If dict has only one element where the key is an addition, coeff must not be negative one. Such a case should be represented by distributing the negation.

The Mul{T}(coeff, dict; type, shape, metadata) constructor validates these constraints and automatically returns the appropriate alternative form where applicable. It should be preferred. BSImpl.AddMul{T}(coeff, dict, variant; type, shape, metadata) is faster but does not validate the constraints and should be used with caution. Incorrect usage can and will lead to both invalid expressions and undefined behavior.

For addition terms, coeff is a constant non-symbolic coefficient added to the expression. dict is a map from terms being added to the constant non-symbolic coefficients they are multiplied by. For example, to represent 1 + 2x + 3y * zcoeff would be 1 and dict would be Dict(x => 2, (y * z) => 3). A valid addition term is subject to the following constraints:

coeff must be non-symbolic.
The values of dict must be non-symbolic.
The keys of dict must not be additions expressions represented with AddMul.
dict must not be empty.
If dict has only one element, coeff must not be zero. Such a case should be represented using the appropriate multiplication form.

The Add{T}(coeff, dict; type, shape, metadata) constructor validates these constraints and automatically returns the appropriate alternative form where applicable. It should be preferred. BSImpl.AddMul{T}(coeff, dict, variant; type, shape, metadata) is faster but does not validate the constraints and should be used with caution. Incorrect usage can and will lead to both invalid expressions and undefined behavior.

struct Div
    const num::BasicSymbolicImpl.Type{T}
    const den::BasicSymbolicImpl.Type{T}
    const simplified::Bool
    const metadata::MetadataT
    const shape::ShapeT
    const type::TypeT
    # ...
end

The Div variant represents division (where the operation is /). num is the numerator and den is the denominator expression. simplified is a boolean indicating whether this expression is in the most simplified form possible. If it is true, certain algorithms in simplify_fractions will not inspect this term. In almost all cases, it should be provided as false. A valid division term is subject to the following constraints:

Both the numerator and denominator cannot be Const variants. This should instead be represented as a Const variant wrapping the result of division.
The numerator cannot be zero. This should instead be represented as a Const wrapping the appropriate zero.
The denominator cannot be one. This should instead be represented as the numerator, possibly wrapped in a Const.
The denominator cannot be zero. This should instead be represented as a Const with some form of infinity.
The denominator cannot be negative one. This should instead be represented as the negation of the numerator.
Non-symbolic coefficients should be propagated to the numerator if it is a constant or multiplication term.

The Div{T}(num, den, simplified; type, shape, metadata) constructor can be used to build this form. If T is SymReal, the constructor will use quick_cancel to cancel trivially identifiable common factors in the numerator and denominator. It will also perform validation of the above constraints and return the appropriate alternative form where applicable. Some of the constraints can be relaxed for non-scalar algebras. The BSImpl.Div{T}(num, den, simplified; type, shape, metadata) does not perform such validation or transformation.

struct ArrayOp
    const output_idx::SmallV{Union{Int, BasicSymbolicImpl.Type{T}}}
    const expr::BasicSymbolicImpl.Type{T}
    const reduce::Any
    const term::Union{BasicSymbolicImpl.Type{T}, Nothing}
    const ranges::Dict{BasicSymbolicImpl.Type{T}, StepRange{Int, Int}}
    const metadata::MetadataT
    const shape::ShapeT
    const type::TypeT
end

ArrayOp is used to represent vectorized operations. This variant should not be created manually. Instead, the @arrayop macro constructs this using a generalized Einstein-summation notation, similar to that of Tullio.jl. Consider the following example:

ex = @arrayop (i, j) A[i, k] * B[k, j] + C[i, j]

This represents A * B + C for matrices A, B, C as a vectorized array operation. Some operations, such as broadcasts, are automatically represented as such a form internally. The following description of the fields assumes familiarity with the @arrayop macro.

When processing this macro, the indices i, j, k are converted to use a common global index variable to avoid potential name conflicts with other symbolic variables named i, j, k if ex is used in a larger expression. The output_idx field stores [i, j]. expr stores the expression A[i, k] * B[k, j] + C[i, j], with i, j, k replaced by the common global index variable. reduce is the operation used to reduce indices not present in output_idx (in this example, k). By default, it is +. term stores an expression that represents an equivalent computation to use for printing/code-generation. For example, here A * B + C would be a valid value for term. By default, term is nothing except when the expression is generated via broadcast or a similar operation. ranges is a dictionary mapping indices used in expr (converted to the common global index) to the range of indices over which they should iterate, in case such a range is explicitly provided.

Note

The common global index variable is printed as _1, _2, ... in arrayops. It is not a valid symbolic variable outside of an ArrayOp's expr.

A valid ArrayOp satisfies the following conditions:

output_idx only contains the integer 1 or variants of the common global index variable.
Any top-level indexing operations in expr use common global indices. A top-level indexing operation is a term whose operation is getindex, and which is not a descendant of any other term whose operation is getindex.
reduce must be a valid reduction operation that can be passed to Base.reduce.
If term is not nothing, it must be an expression with shape shape and type type.
The keys of ranges must be variants of the common global index variable, and must be present in expr.

The @arrayop macro should be heavily preferred for creating ArrayOps. In case this is not possible (such as in recursive functions like substitute) the ArrayOp constructor should be preferred. This does not allow specifying the type and shape, since these values are tied to the fields of the variant and are thus determined. The BSImpl.ArrayOp constructor should be used with extreme caution, since it does not validate input.

Array arithmetic

SymbolicUtils implements a simple array algebra in addition to the default scalar algebra. Similar to how SymbolicUtils.promote_symtype given a function and symtypes of its arguments returns the symtype of the result, SymbolicUtils.promote_shape does the same for the shapes of the arguments. Implementing both methods is cruicial for correctly using custom functions in symbolic expressions. Without promote_shape, SymbolicUtils will use Unknown(-1) as the shape.

The array algebra implemented aims to mimic that of base Arrays as closely as possible. For example, a symbolic adjoint(::Vector) * (::Vector) will return a symbolic scalar instead of a one-element symbolic vector. promote_shape implementations will propagate the shape information on a best-effort basis. Invalid shapes (such as attempting to multiply a 3-dimensional array) will error. Following are notable exceptions to Base-like behavior:

map and mapreduce require that all input arrays have the same shape
promote_symtype and promote_shape is not implemented for map and mapreduce, since doing so requires the function(s) passed to map and mapreduce instead of their types or shapes.
Since ndims information is not present in the type, eachindex, iterate, size, axes, ndims, collect are type-unstable. safe_eachindex is useful as a type-stable iteration alternative.
ifelse requires that both the true and false cases have identical shape.
Symbolic arrays only support cartesian indexing. For example, given @syms x[1:3, 1:3] accessing x[4] is invalid and x[1, 2] should be used. Valid indices are Int, Colon, AbstractRange{Int} and symbolic expressions with integer symtype. A single CartesianIndex of appropriate dimension can also be used to symbolically index arrays.

Symbolic array operations are also supported on arrays of symbolics. However, at least one of the arguments to the function must be a symbolic (instead of an array of symbolics) to allow the dispatches defined in SymbolicUtils to be targeted instead of those in Base. To aid in constructing arrays of symbolics, the BS utility is provided. Similar to the T[...] syntax for constructing an array of etype T, BS[...] will construct an array of BasicSymbolics. At least one value in the array must be a symbolic value to infer T in Array{BasicSymbolic{T}, N}. To explicitly specify the vartype, use BS{T}[...].

Symbolic functions and dependent variables

SymbolicUtils defines FnType{A, R, T} for symbolic functions and dependent variables. Here, A is a Tuple{...} of the symtypes of arguments and R is the type returned by the symbolic function. T is the type that the function itself subtypes, or Nothing.

The syntax

@syms f(::T1, ::T2)::R

creates f with a symtype of FnType{Tuple{T1, T2}, R, Nothing}. This is a symbolic function taking arguments of type T1 and T2, and returning R. Nothing is a sentinel indicating that the supertype of the function is unknown. By contrast,

@syms f(..)::R

creates f with a symtype of FnType{Tuple, R, Nothing}. SymbolicUtils considers this case to be a dependent variable with as-yet unspecified independent variables. In other words,

@syms x f1(::Real)::Real f2(..)::Real

Here, f1(x) is considered a symbolic function f1 called with the argument x and f2(x) is considered a dependent variable that depends on x. The utilities SymbolicUtils.is_function_symbolic, SymbolicUtils.is_function_symtype, symbolicUtils.is_called_function_symbolic can be used to differentiate between these cases.

API

Basics

Missing docstring.

Missing docstring for SymbolicUtils.BasicSymbolic. Check Documenter's build log for details.

SymbolicUtils.@syms — Macro

@syms <lhs_expr>[::T1] <lhs_expr>[::T2]...

For instance:

@syms foo::Real bar baz(x, y::Real)::Complex

Create one or more variables. <lhs_expr> can be just a symbol in which case it will be the name of the variable, or a function call in which case a function-like variable which has the same name as the function being called. The Sym type, or in the case of a function-like Sym, the output type of calling the function can be set using the ::T syntax.

Examples:

@syms foo bar::Real baz::Int will create

variable foo of symtype Number (the default), bar of symtype Real and baz of symtype Int

@syms f(x) g(y::Real, x)::Int h(a::Int, f(b)) creates 1-arg f 2-arg g

and 2 arg h. The second argument to h must be a one argument function-like variable. So, h(1, g) will fail and h(1, f) will work.

Formal syntax

Following is a semi-formal CFG of the syntax accepted by this macro:

# any variable accepted by this macro must be a `var`.
# `var` can represent a quantity (`value`) or a function `(fn)`.
var = value | fn
# A `value` is represented as a name followed by a suffix
value = name suffix
# A `name` can be a valid Julia identifier
name = ident |
# Or it can be an interpolated variable, in which case `ident` is assumed to refer to
# a variable in the current scope of type `Symbol` containing the name of this variable.
# Note that in this case the created symbolic variable will be bound to a randomized
# Julia identifier.
       "$" ident |
# Or it can be of the form `Foo.Bar.baz` referencing a value accessible as `Foo.Bar.baz`
# in the current scope.
       getproperty_literal
getproperty_literal = ident "." getproperty_literal | ident "." ident
# The `suffix` can be empty (no suffix) which defaults the type to `Number`
suffix = "" |
# or it can be a type annotation (setting the type of the prefix). The shape of the result
# is inferred from the type as best it can be. In particular, `Array{T, N}` is inferred
# to have shape `Unknown(N)`.
         "::" type |
# or it can be a shape annotation, which sets the shape to the one specified by `ranges`.
# The type defaults to `Array{Number, length(ranges)}`
         "[" ranges "]" |
# lastly, it can be a combined shape and type annotation. Here, the type annotation
# sets the `eltype` of the symbolic array.
         "[" ranges "]::" type
# `ranges` is either a single `range` or a single range followed by one or more `ranges`.
ranges = range | range "," ranges
# A `range` is simply two bounds separated by a colon, as standard Julia ranges work.
# The range must be non-empty. Each bound can be a literal integer or an identifier
# representing an integer in the current scope.
range = (int | ident) ":" (int | ident) |
# Alternatively, a range can be a Julia expression that evaluates to a range. All identifiers
# used in `expr` are assumed to exist in the current scope.
        expr |
# Alternatively, a range can be a Julia expression evaluating to an iterable of ranges,
# followed by the splat operator.
        expr "..."
# A function is represented by a function-call syntax `fncall` followed by the `suffix`
# above. The type and shape from `suffix` represent the type and shape of the value
# returned by the symbolic function.
fn = fncall suffix
# a function call is a call `head` followed by a parenthesized list of arguments.
fncall = head "(" args ")"
# A function call head can be a name, representing the name of the symbolic function.
head = ident |
# Alternatively, it can be a parenthesized type-annotated name, where the type annotation
# represents the intended supertype of the function. In other words, if this symbolic
# function were to be replaced by an "actual" function, the type-annotation constrains the
# type of the "actual" function.
       "(" ident "::" type ")"
# Arguments to a function is a list of one or more arguments
args = arg | arg "," args
# An argument can take the syntax of a variable (which means we can represent functions of
# functions of functions of...). The type of the variable constrains the type of the
# corresponding argument of the function. The name and shape information is discarded.
arg = var |
# Or an argument can be an unnamed type-annotation, which constrains the type without
# requiring a name.
      "::" type |
# Or an argument can be the identifier `..`, which is used as a stand-in for `Vararg{Any}`
      ".." |
# Or an argument can be a type-annotated `..`, representing `Vararg{type}`. Note that this
# and the previous version of `arg` can only be the last element in `args` due to Julia's
# `Tuple` semantics.
      "(..)::" type |
# Or an argument can be a Julia expression followed by a splat operator. This assumes the
# expression evaluates to an iterable of symbolic variables whose `symtype` should be used
# as the argument types. Note that `expr` may be evaluated multiple times in the macro
# expansion.
      expr "..."