Skip to content

Instantly share code, notes, and snippets.

@shashi
Last active April 19, 2021 23:36
Show Gist options
  • Save shashi/419058d8840ac789fe0829b84157949a to your computer and use it in GitHub Desktop.
Save shashi/419058d8840ac789fe0829b84157949a to your computer and use it in GitHub Desktop.

Extensible Symbolic Programming in Julia

       Vancouver Julia Meetup

            04/20/2021
           Shashi Gowda

         PhD candidate in
Computational Science and Engineering
         Julia Lab, MIT

Outline

  • Tour of Symbolics.jl
  • Code generation
  • Rule-based rewriting of expressions
  • Interface and Interoperability

Symbolics.jl/SymbolicUtils.jl - Goals

  • Easy to use (like SymPy/Mathematica)
  • Fast!: use the fastest implementation
  • Extensible: allow packages with term-like types to use the features
  • Code-generation: turn symbolic results into executable, parallel Julia code.

Symbolics.jl - Tour

Expressions

using Symbolics
@variables t x y u(..)

x + y*x + x*y, x^2 + y

u(x,y,t)

Symbolics.jl - Tour

Expressions

julia> typeof(x+x*y)
Num # <- we will come back to this later

julia> typeof(x+x*y) <: Number
true

Symbolics.jl - Tour

Arrays of Expressions

julia> M = [t^2 + t + t^2  2t + 4t
            x + y + y + 2t  x^2 - x^2 + y^2]
2×2 Matrix{Num}:
  t + 2(t^2)   6t
 x + 2t + 2y  y^2

Symbolics.jl - Tour

Arrays of Expressions

julia> using LinearAlgebra

julia> det(M)
(t + 2(t^2))*(y^2) - (6t*(x + 2t + 2y))

julia> M \ M
2×2 Matrix{Num}:
 1  0
 0  1

julia> M * M

Symbolics.jl - Tour

Codegen!

using LinearAlgebra, SparseArrays

julia> f, f! = build_function(M, [x, y, t])

julia> f, f! = build_function(sparse(Diagonal(M)), [x, y, t])

Symbolics.jl - Tour

Differentiation

julia> Symbolics.jacobian([x + x*y, x^2 + y],
                          [x, y])
2×2 Matrix{Num}:
 1 + y  x
    2x  1

rosenbrock(X) = sum(1:length(X)-1) do i
    100 * (X[i+1] - X[i]^2)^2 + (1 - X[i])^2
end

@variables X[1:100]

rr = rosenbrock(X)

Symbolics.hessian(rr, X)

Symbolics.jl - Tour

Sparsity

using BenchmarkTools

julia> subsets = getindex.((X,), UnitRange.((1:5:100),(5:5:100)));

julia> exprs = rosenbrock.(subsets);

julia> @btime Symbolics.jacobian(exprs, X)
  1.078 s (4935567 allocations: 198.67 MiB)

julia> @btime Symbolics.sparsejacobian(exprs, X);
  23.134 ms (105622 allocations: 5.24 MiB)

julia> Symbolics.jacobian_sparsity(exprs, X)

Symbolics.jl - Tour

Sparsity

using BenchmarkTools

julia> subsets = getindex.((X,), UnitRange.((1:5:100),(5:5:100)))

julia> exprs = rosenbrock.(subsets);

julia> @btime Symbolics.jacobian(exprs, X)
  1.078 s (4935567 allocations: 198.67 MiB)

julia> @btime Symbolics.sparsejacobian(exprs, X);
  23.134 ms (105622 allocations: 5.24 MiB)

julia> Symbolics.jacobian_sparsity(exprs, X)

Symbolics.jl - Missing Parts

  • Integration (definite integral)
  • Polynomial division (easy)
  • Groebner basis
  • Polynomial equation solving
  • Simplification under assumptions (Z3)

Under the hood - What is an Expression?

julia> @variables x y

julia> sin(cos(x))

julia> typeof(sin(cos(y)))

julia> using Symbolics: value

julia> value(sin(cos(x)))

Under the hood - What is an Expression?

julia> @variables x y

julia> value(sin(cos(x)))

julia> using SymbolicUtils

julia> SymbolicUtils.print_tree(stdout,
                value(sin(cos(x))))

Under the hood - What is an Expression?

Add -- Asymptotically better for storage

using SymbolicUtils

@syms x y

x + y

z1 = 10
for i=1:1000
    z1 += sin(rand([x, y]))
end

Under the hood - What is an Expression?

Mul -- Asymptotically better for storage

z2 = 10
for i=1:1000
    z2 *= sin(rand([x, y]))
end

Rule-based rewriting

using SymbolicUtils

@syms w z α::Real β::Real

r1 = @rule sin(2(~x)) => 2sin(~x)*cos(~x)

r1(sin(2z))

r1(sin(3z))

r1(sin(2*(w-z)))

r2 = @rule sin(~x + ~y) => sin(~x)*cos(~y) + cos(~x)*sin(~y);

r2(sin+β))

How do rules look at a term?

  • istree
  • operation
  • arguments
function show_structure(x)
    @show typeof(x)
    @show istree(x)
    if istree(x)
        @show operation(x)
        @show arguments(x)
    end
end

show_structure(hypot(sin(x), cos(x)))
show_structure(z1)
show_structure(z2)

Benefits of the Interface

  • Using the term interface allows term-rewriting of any type.
  • It also allows packages like Metatheory to be an alternative term-rewriting system

Thank you!

(Thanks to Yingbo Ma for helping me prepare this talk!)

Work presented here is by:

  • me
  • Yingbo Ma
  • Chris Rackauckas
  • Mason Protter
  • many other contributors
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment