juliasymbolics/symbolics.jl

Symbolic programming for the next generation of numerical software

Symbolics.jl

Symbolics.jl is a fast and modern Computer Algebra System (CAS) for a fast and modern programming language (Julia). The goal is to have a high-performance and parallelized symbolic algebra system that is directly extendable in the same language as that of the users.

Installation

To install Symbolics.jl, use the Julia package manager:

julia> using Pkg
julia> Pkg.add("Symbolics")

Documentation

For information on using the package, see the stable documentation. Use the in-development documentation for the version of the documentation which contains the unreleased features.

Relationship to Other Packages

• SymbolicUtils.jl: This is a rule-rewriting system that is the core of Symbolics.jl. Symbolics.jl builds off of SymbolicUtils.jl to extend it to a whole symbolic algebra system, complete with support for differentiation, solving symbolic systems of equations, etc. If you're looking for the barebones to build a new CAS for specific algebras, SymbolicUtils.jl is that foundation. Otherwise, Symbolics.jl is for you.
• ModelingToolkit.jl: This is a symbolic-numeric modeling system for the SciML ecosystem. It heavily uses Symbolics.jl for its representation of symbolic equations along with tools like differentiation, and adds the representation of common modeling systems like ODEs, SDEs, and more.

Example

julia> using Symbolics

julia> @variables t x y
julia> D = Differential(t)

julia> z = t + t^2
julia> D(z) # symbolic representation of derivative(t + t^2, t)
Differential(t)(t + t^2)

julia> expand_derivatives(D(z))
1 + 2t

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

julia> B = simplify.([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

julia> simplify.(substitute.(B, (Dict(x => y^2),)))
2×2 Matrix{Num}:
    t + 2(t^2)   6t
 2t + y^2 + 2y  y^2

julia> substitute.(B, (Dict(x => 2.0, y => 3.0, t => 4.0),))
2×2 Matrix{Num}:
 36.0  24.0
 16.0   9.0

Citation

If you use Symbolics.jl, please cite this paper (or see the free arxiv version)

@article{10.1145/3511528.3511535,
author = {Gowda, Shashi and Ma, Yingbo and Cheli, Alessandro and Gw\'{o}\'{z}zd\'{z}, Maja and Shah, Viral B. and Edelman, Alan and Rackauckas, Christopher},
title = {High-Performance Symbolic-Numerics via Multiple Dispatch},
year = {2022},
issue_date = {September 2021},
publisher = {Association for Computing Machinery},
address = {New York, NY, USA},
volume = {55},
number = {3},
issn = {1932-2240},
url = {https://doi.org/10.1145/3511528.3511535},
doi = {10.1145/3511528.3511535},
abstract = {As mathematical computing becomes more democratized in high-level languages, high-performance symbolic-numeric systems are necessary for domain scientists and engineers to get the best performance out of their machine without deep knowledge of code optimization. Naturally, users need different term types either to have different algebraic properties for them, or to use efficient data structures. To this end, we developed Symbolics.jl, an extendable symbolic system which uses dynamic multiple dispatch to change behavior depending on the domain needs. In this work we detail an underlying abstract term interface which allows for speed without sacrificing generality. We show that by formalizing a generic API on actions independent of implementation, we can retroactively add optimized data structures to our system without changing the pre-existing term rewriters. We showcase how this can be used to optimize term construction and give a 113x acceleration on general symbolic transformations. Further, we show that such a generic API allows for complementary term-rewriting implementations. Exploiting this feature, we demonstrate the ability to swap between classical term-rewriting simplifiers and e-graph-based term-rewriting simplifiers. We illustrate how this symbolic system improves numerical computing tasks by showcasing an e-graph ruleset which minimizes the number of CPU cycles during expression evaluation, and demonstrate how it simplifies a real-world reaction-network simulation to halve the runtime. Additionally, we show a reaction-diffusion partial differential equation solver which is able to be automatically converted into symbolic expressions via multiple dispatch tracing, which is subsequently accelerated and parallelized to give a 157x simulation speedup. Together, this presents Symbolics.jl as a next-generation symbolic-numeric computing environment geared towards modeling and simulation.},
journal = {ACM Commun. Comput. Algebra},
month = {jan},
pages = {92–96},
numpages = {5}
}

symbolics.jl's Issues

Difference Operator

That's a good idea. It would be different from the derivative operator

jlperla sglyon

┆Issue is synchronized with this Trello card by Unito

Juno inline display of Variable

If the following code is evaluated in the Juno editor and the result bubble is expanded to show the parameters, the display of Variable shows the docstring, it's also noticeably slow, taking a few seconds to display. It would perhaps be better just to show the name of the variable?

using ModelingToolkit
@variables a
@parameters b
eqs = [0~a+1]
NonlinearSystem(eqs, [a], [b])

┆Issue is synchronized with this Trello card by Unito

Latexify fails for nested derivatives

The Latexify recipe fails for nested derivative operators. It seems like this case was just missed in the current implementation / tests. This postwalk will pass the Latex string of the inner derivative to Latexify.latexraw when processing the outer derivative.

The branch to fix this should be aware of SciML/ModelingToolkit.jl#502.

using ModelingToolkit, Latexify

@parameters t, x

@variables u(..)

@derivatives Dt'~t
@derivatives Dx'~x

eqs = [Dt(u(t, x)) ~ Dx(Dx(u(t, x)))]

latexify(eqs)

in Latexify.jl: 
You are trying to create latex-maths from a `String` that cannot be parsed as
an expression. 

`latexify` will, by default, try to parse any string inputs into expressions
and this parsing has just failed.

If you are passing strings that you want returned verbatim as part of your input,
try making them `LaTeXString`s first. 

If you are trying to make a table with plain text, try passing the keyword
argument `latex=false`. You should also ensure that you have chosen an output
environment that is capable of displaying not-maths objects. Try for example
`env=:table` for a latex table or `env=:mdtable` for a markdown table.


Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] latexraw(::String; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/dan/.julia/packages/Latexify/7DHiI/src/latexraw.jl:114
 [3] latexraw(::String) at /home/dan/.julia/packages/Latexify/7DHiI/src/latexraw.jl:110
 [4] (::ModelingToolkit.var"#414#430")(::Expr) at /home/dan/.julia/packages/ModelingToolkit/BHRXD/src/latexify_recipes.jl:14
 [5] walk(::Expr, ::Function, ::ModelingToolkit.var"#414#430") at /home/dan/.julia/packages/MacroTools/gME9C/src/utils.jl:112
 [6] postwalk at /home/dan/.julia/packages/MacroTools/gME9C/src/utils.jl:122 [inlined]
 [7] SciML/ModelingToolkit.jl#413 at ./none:0 [inlined]
 [8] iterate at ./generator.jl:47 [inlined]
 [9] collect(::Base.Generator{Array{Expr,1},ModelingToolkit.var"#413#429"}) at ./array.jl:686
 [10] macro expansion at /home/dan/.julia/packages/ModelingToolkit/BHRXD/src/latexify_recipes.jl:14 [inlined]
 [11] apply_recipe(::Array{Equation,1}; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/dan/.julia/packages/Latexify/7DHiI/src/recipes.jl:164
 [12] apply_recipe at /home/dan/.julia/packages/Latexify/7DHiI/src/recipes.jl:158 [inlined]
 [13] latexify(::Array{Equation,1}; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/dan/.julia/packages/Latexify/7DHiI/src/latexify_function.jl:4
 [14] latexify(::Array{Equation,1}) at /home/dan/.julia/packages/Latexify/7DHiI/src/latexify_function.jl:4
 [15] top-level scope at In[6]:12
 [16] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091

┆Issue is synchronized with this Trello card by Unito

Registering ∈ does not work completely

julia> import Base: ∈

julia> @register ∈(x, y)
in (generic function with 70 methods)

julia> x ∈ 1
in(x, 1)

julia> x ∈ 1:5
ERROR: MethodError: no method matching isinteger(::Num)

Get variables from an array of terms gives an empty set

This was unexpected:

julia> var = @variables x y
(x, y)

julia> expr1 = x + y <= 1; expr2 = x - y <= 2;

julia>  aux = [expr1, expr2]
2-element Array{Term{Bool},1}:
 (x + y) <= 1
 (x - y) <= 2

julia> get_variables(aux)
Any[]

Maybe a new method for get_variables can be added for an array of Term near
https://github.com/SciML/ModelingToolkit.jl/blob/135821e06a4b0e16adb62c0ad87410cf1d57c20e/src/utils.jl#L105.

build_function on large system, AIBECS example for testing

I would like to leverage some tools from MTK to use in AIBECS.jl, and as per discussions on slack, here is a MWE to test build_function on an AIBECS example. In the example below I used a tiny model grid, but the line defining grd (the model grid) should be changed to get larger systems for testing, so I added a bunch of comments to point this out in the code directly.

What the MWE does, essentially, is build a sparse operator A that acts on u, and tries to make A*u fast by replacing u -> A*u with an in-place function that it traces with build_function. Hopefully, this is helpful to your suite test! The MWE:

using ModelingToolkit
using AIBECS
using BenchmarkTools

# Load the model grid
# Note: Replace `Primeau_2x2x2` by `OCCA` or `OCIM2` for large systems
# sizes (of `u` for these grids):
# - `Primeau_2x2x2` -> 5
# - `OCCA` -> ≈ 80'000
# - `OCIM2` -> ≈ 200'000
# Note 2: the large systems have to be downloaded (once only, thanks to DataDeps.jl)
# They are not too big in memory though, so that should be OK.
const grd, T = Primeau_2x2x2.load() # <- change `Primeau_2x2x2` by `OCCA` or `OCIM2` here for larger systems

# Create a sparse operator (which will act on `u`)
# FYI, A is a diagonal + a far-off diagonal
A(w₀, w′) = transportoperator(grd, z -> w₀ + w′ * z, fsedremin=0.1)

# Now trace `A` with MTK to recreate a fast version of it
# Set up parameters and a dummy tracer to trace T_POP(p) * x
const nb = count(iswet(grd)) # size of `u`
@parameters w₀ w′ u[1:nb]
# trace the resulting operation of 
Au = simplify.(A(w₀, w′) * u)
Au_fun! = eval(build_function(Au, w₀, w′, u)[2])

# Now test that the fast function returns the same as the slow A * u on a random u
u = rand(nb)
du = rand(nb)
Au_fun!(du, 1.0, 2.0, u)
@show du ≈ A(1.0, 2.0) * u

# And compare the timing
@btime A(1.0, 2.0) * $u
@btime Au_fun!($du, 1.0, 2.0, $u)

Obviously, I would love to improve my code so feel free to tell me what's wrong with it if you spot anything!

┆Issue is synchronized with this Trello card by Unito

checkbounds kwarg in build_function not behaving properly

For checkbounds in build_function, the docs say:

`checkbounds`: For whether to enable bounds checking inside of the generated
  function. Defaults to false, meaning that `@inbounds` is applied

but this doesn't seem to be the case. E.g.

julia> using ModelingToolkit

julia> @variables x
(x,)

julia> build_function([x],[x], checkbounds=false)[2]
:(function (var"##out#278", var"##arg#277")
      #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:252 =#
      #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:253 =#
      let x = var"##arg#277"[1]
          #= /Users/marc/.julia/packages/ModelingToolkit/P0rIX/src/build_function.jl:253 =#
          begin
              #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:295 =#
              var"##out#278"[1] = x
              #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:297 =#
              nothing
          end
      end
  end)

there's no @inbounds When checkbounds=true

julia> build_function([x],[x], checkbounds=true)[2]
:(function (var"##out#280", var"##arg#279")
      #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:252 =#
      #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:253 =#
      let x = var"##arg#279"[1]
          #= /Users/marc/.julia/packages/ModelingToolkit/P0rIX/src/build_function.jl:253 =#
          #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:299 =# @inbounds begin
                  #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:295 =#
                  var"##out#280"[1] = x
                  #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:297 =#
                  nothing
              end
      end
  end)

there is an @inbounds. For the out-of-place function the option is ignored for both:

julia> build_function([x],[x], checkbounds=true)[1]
:(function (var"##arg#281",)
      #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:252 =#
      #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:253 =#
      let x = var"##arg#281"[1]
          #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:344 =#
          (SymbolicUtils.Code.create_array)(typeof(var"##arg#281"), nothing, Val{(1,)}(), x)
      end
  end)

julia> build_function([x],[x], checkbounds=false)[1]
:(function (var"##arg#283",)
      #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:252 =#
      #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:253 =#
      let x = var"##arg#283"[1]
          #= /Users/marc/.julia/packages/SymbolicUtils/HntGZ/src/code.jl:344 =#
          (SymbolicUtils.Code.create_array)(typeof(var"##arg#283"), nothing, Val{(1,)}(), x)
      end
  end)

I'm using ModelingToolkit 5.6.2 on julia 1.5.2.

Discussion: The problem with Expression <: Number, and a proposed fix

I wasn't sure where to put this but I figured this was a good place. Maybe discourse instead? Open to suggestions

---

@ChrisRackauckas and I discussed yesterday that julia's subtyping dispatch model makes symbolic graph languages like ModelingToolkit (MTK) difficult. For instance, in MTK we have Expression <: Number, but this is far from ideal. For instance, what if someone is more interested in symbolically representing arrays? Or strings? or some other type that doesn't have anything to do with numbers. Sure, we can get by with <: Number even if the object isn't actually a number, but it's awkward.

I think once (if?) MTK starts doing algebraic simplicifications, it's going to need to be able to deal with applying type dependant simplification rules, for example only assuming commutativity between Real and Complex numbers but not say Matrix or Quaternion.

One approach for getting around the sub-typing constraint that I've been exploring is to instead to have symbolic computations happen inside a Cassette.jl pass, inside of which we apply code transformations such that Expression{T} acts like it is a T.

It turns out that @shashi just recently opened a PR to Cassette that would make this quite easy, so I mocked up a demo of what I was thinking that I'd like to share here for discussion.

Two notes

1. The following code requires Shashi's Cassette PR JuliaLabs/Cassette.jl#157.
2. Apologies for not following MTK naming conventions in this demo but I think you'll get the idea

#--------------------------------------------------------------------------------
# Set up some symbolic types
#--------------------------------------------------------------------------------
abstract type Symbolic{T} end # Symbolic{T} will act like it is <: T

struct Sym{T} <: Symbolic{T}
    name::Symbol
end

Symbol(s::Sym) = s.name

struct SymExpr{T} <: Symbolic{T}
    op
    args::Vector{Any}
end

#--------------------------------------------------------------------------------
# Pretty printing
#--------------------------------------------------------------------------------
function Base.show(io::IO, s::Sym{T}) where {T}
    print(io, string(s.name)*"::$T")
end

expr(se::SymExpr, Ts) = Expr(:call, expr(se.op, Ts), expr.(se.args, (Ts,))...)
expr(x, Ts)           = x
expr(f::Function, Ts) = Symbol(f)
function expr(s::Sym{T}, Ts) where {T} 
    if s ∉ Ts
        push!(Ts, s)
    end
    s.name
end

function Base.show(io::IO, se::SymExpr{T}) where {T}
    sset = Set()
    ex = expr(se, sset)
    print(io, repr(ex)[2:end]*"::$T" * " where {"*repr(sset)[9:end-2]*"}")
end

#--------------------------------------------------------------------------------
# Set up the Cassette pass
#--------------------------------------------------------------------------------
using Cassette: Cassette, overdub, @context, ReflectOn
using Base.Core.Compiler: return_type
using SpecializeVarargs
@context SymContext

sym_substitute(::Type{Sym{T}})     where {T} = T
sym_substitute(::Type{SymExpr{T}}) where {T} = T
sym_substitute(::Type{T})          where {T} = T

function Cassette.overdub(ctx::SymContext, f::Function, args...)
    argsT = typeof(args)
    if any((<:).(argsT.parameters, Symbolic))
        argsT′ = [sym_substitute.(argsT.parameters)...]
        if isprimitive(f)
            SymExpr{return_type(f, Tuple{argsT′...})}(f, [args...])
        else
            overdub(ctx, ReflectOn{Tuple{typeof(f), argsT′...}}(), f, args...)
        end
    else
        f(args...)
    end
end

# If isprimitive(f) == true, then inside a pass, we won't recurse into the insides of f. 
# Primitives are stopping points for us 
for f in [:+, :-, :*, :/, :^, :exp, :log, 
          :sin, :cos, :tan, :asin, :acos, :atan,
          :sinh, :cosh, :tanh, :asinh, :acosh, :atanh, :adjoint]
    @eval isprimitive(::typeof($f)) = true
end
isprimitive(::Any) = false

# Convenience macro for wrapping any enclosed function calls in the cassette pass
using MacroTools: postwalk
macro sym(expr)
    out = postwalk(expr) do ex
        if ex isa Expr && ex.head == :call
            :(overdub(SymContext(), $(ex.args[1]), $(ex.args[2:end]...)))
        else
            ex
        end
    end
    esc(out)
end

Okay, with that setup code out of the way, lets see what this got us. Suppose we have two functions from an external library that have very strict type constraints:

f(x::Float64)     = 1 + sin(x)^(1/2)
g(x::Vector{Int}) = x'x + 2

With the above definitions, we have no problems getting inside these functions:

julia> x = Sym{Float64}(:x)
x::Float64

julia> y = Sym{Vector{Int}}(:y)
y::Array{Int64,1}

julia> @sym f(x)
(1 + sin(x) ^ 0.5)::Float64 where {x::Float64}

julia> @sym g(y)
(adjoint(y) * y + 2)::Int64 where {y::Array{Int64,1}}

---

Okay, I thought that was pretty neat, but there's still some major problems with this sort of approach. Here's a two that are forefront in my mind:

1. This approach can be too good at recursing into code. For instance, you wouldn't want to accidentally call sin(x) on a symbolic x if it wasn't registered because then you'd be exposed to the internals of sin which nobody is interested in seeing. Alternatively, we could make it so that it only recurses into registered functions instead of the other way around, but I think that might almost as problematic because it'd be incredibly annoying to have to register every function you might care about, plus it'd make symbolic differentiation more difficult.

2. If we want types associated with expressions, then we either have to build in the mapping between input and output types for every registered function, or we have to rely on Core.Compiler.return_type which brings it's own host of problems, which are not limited to the usual warnings you'll hear from complier people. For instance, return_type(+, Tuple{Real, Real}) == Any, so Syms should be concretely typed instead of abstract 😕.

Do you think this is a dead end? Do you think it's promising? Any other thoughts or comments? Questions?

Names of array of variables

Hi,
I noticed that variable with indexes have a symbol between numbers

Variable("x", 1, 2, 3)

prints

x₁̒₂̒₃

Is this on purpose?

┆Issue is synchronized with this Trello card by Unito

The alternative to `ifelse`

ifelse is a bit problematic for symbolic libraries because it's in Core and not in Base and dispatch cannot be added to it. So, while we want to represent branching, we need to do so without using ifelse. The opportunities are something like:

• ModelingToolkit.ifelse
• _ifelse (exported)
• A different name

We need to just pick one and start using it, and whatever is chosen is what the eventual tracer should transform things to.

┆Issue is synchronized with this Trello card by Unito

Register

@JuliaRegistrator register()

"ERROR: LoadError: syntax: "3" is not a valid function argument name" when including saved Julia expressions with `var"x(t)"` style names

Example doesn't work. Still working on small reproducible.

Build function efficient multiple returns

It is often the case (for instance in optimization) that one needs to calculate some functions, such as the gradient and the hessian, that have some shared calculations. Is there a way in modelingtoolkit to use symbolic computation and build_function to generate functions that act something like?

function gradhess!(grad,hess)

aux=shared calculations

if grad != nothing

  grad.=compute_grad (aux)

end

if hessian != nothing

  hess.= compute_hess(aux)

end

end

┆Issue is synchronized with this Trello card by Unito

Array-like symbolic objects for more efficient code generation

using ModelingToolkit
using BenchmarkTools
@variables a[1:1]
a_copy = copy(a)
a_test = [3.0]
function my_add!(a)
    @inbounds a[1] += 2.0
    return nothing
end
my_add!(a_copy)
my_add2! = eval(ModelingToolkit.build_function(a_copy,a)[2])
my_add2!(a) = my_add2!(a,a) #why does this not work?
my_add3!(a) = my_add2!(a,a)
@code_native my_add3!(a_test)
@code_native my_add!(a_test)
@btime my_add!($a_test) #2ns
@btime my_add2!($a_test,$a_test) #12ns
@btime my_add3!($a_test) #12ns

YingboMa could you have a look at this?

┆Issue is synchronized with this Trello card by Unito

Complex functions

It would be great to be able to extract real and imaginary parts of a complex function.
Currently this does not seem to work, e.g.

julia> @variables a b
(a, b)

julia> c = a + b * im
a + b * (im)

julia> c^2
(a + b * (im)) ^ 2

It needs to know that im^2 is -1.

Conversion error when evaluating vector-valued functions on SubArray

When using build_function() to generate Julia functions which take one or more vectors as arguments and again output a vector, one encounters a conversion error when passing a SubArray into the first argument of the generated function along the lines of: "Cannot convert object of type Array{Float64,1} to SubArray{Float64,1}".

MWE:

@variables y[1:2]
f(x) = x.^2
g = build_function(f(y), y)[1] |> eval
g(view(rand(3), 1:2))

Said conversion error does not come up if the desired function outputs a Number or a Matrix which leads me to believe that this issue is a result of trying to retain type-stability. I have not yet taken a closer look at the source code which is responsible for this error, but it seems like fixing this should be more or less straightforward.

┆Issue is synchronized with this Trello card by Unito

Simplification no able to multiply

If I have a result like this

(2.44140625*10^-4*(2048*t^4-192*t^2*(16*t^2+1)))/t^6

it won't simplify further.Whereas in Maxima, it gives

-(16*t^2+3)/(64*t^4)

[request] support for symbolic integration

something like:

@variables t
D = Differential(t)
z = t
julia> expand_integral(∫(z))
0.5(t^2) or 1//2(t^2)

also a question is handling the new constant, should you define a new parameter?

┆Issue is synchronized with this Trello card by Unito

Define intersection and union

Intersection and union should be defined to return the corresponding symbolic expressions.
Currently they fall back to a generic array method and give useless results:

julia> @variables x, y

julia> x ∩ y
Operation[]

julia> @which x ∩ y
intersect(itr, itrs...) in Base at array.jl:2566

toexpr() return function of operator in AST instead Symbol

There may be an error in the toexpr method

@parameters t
pow_f =  t^2 
expr =toexpr(pow_f)

AST contains ^ (function of type typeof(^))

julia> dump(expr) 
Expr
  head: Symbol call
  args: Array{Any}((3,))
    1: ^ (function of type typeof(^))
    2: Symbol t
    3: Int64 2

but perhaps, it should contain Symbol ^

julia> dump(:((^)(t, 2)))
Expr
 head: Symbol call
 args: Array{Any}((3,))
   1: Symbol ^
   2: Symbol t
   3: Int64 2

Derivative of conj?

@variables u[1:2]
ModelingToolkit.gradient(u'*u, u)
# ERROR: MethodError: no method matching derivative(::typeof(conj), ::Tuple{Operation}, ::Val{1})

I might be doing something naive here, but is there a way to specify that u is all real and u₁ == conj(u₁)?
For now I'm using sum(u.^2) instead of u'*u.

Variable typing: better handling of the internal datatype

We currently do Variable{T}(name) mostly because handling DataType is a bit difficult. Namely, if you have the DataType, how can you get the type?

x = typeof(rand(3,3))
eval(x)
x = getproperty(Main,Symbol(x))

Seems like you have to use eval, which could be an issue with scoping.

@YingboMa @shashi

Multiplication with numbers from Unitful in expressions

When using ModelingToolkit and Unitful expressions left multiplication of numbers with units is not defined.

using Unitful, ModelingToolkit
@variables x
x*(1u"kg") #this works
(1u"kg")*x #this does not

┆Issue is synchronized with this Trello card by Unito

Unnecessary tupling of arrays of variables

Following the instructions from the readme, this should work, but doesn't. Related to SciML/ModelingToolkit.jl#153.

julia> a,b,c = @variables X[1:3];
ERROR: BoundsError: attempt to access (Operation[X₁(), X₂(), X₃()],)
  at index [2]
Stacktrace:
 [1] indexed_iterate(::Tuple{Array{Operation,1}}, ::Int64, ::Int64) at ./tuple.jl:60
 [2] top-level scope at none:0

You can work around it with:

julia> X = @variables X[1:3]
(Operation[X₁(), X₂(), X₃()],)

julia> a,b,c = X[1]
3-element Array{Operation,1}:
 X₁()
 X₂()
 X₃()

Why does @variables X[1:3] return a tuple?

Parsing of array functions

Related to SciML/ModelingToolkit.jl#102, how should something like this be dealt with:

f(X) = sum(X .* X)

This is probably not common in the ODE context, but is very common in the optimization context.

Simplification of transposed scalars

MT fails to recognize a simplification of scalars that are transposed. In the following MWE, the transpose wrapper prevents the multiplication be 0 to be simplified.

julia> simplify(0*Ja)
ModelingToolkit.Constant(0)

julia> o
transpose(0.0)

julia> simplify(o*Ja)
transpose(0.0) * Ja

┆Issue is synchronized with this Trello card by Unito

Create dependent variable via ModelingToolkit.Variable

According to https://mtk.sciml.ai/dev/IR/ I can create a dependent variable x(t) via

t = Variable(:t)
x = Variable(:x)(t)

However, on #master I get the error message

Sym x is not callable. Use @syms x(var1, var2,...) to create it as a callable. See ?@fun for more options

Comparing the source with v3.20.1, the corresponding functor is gone:

(x::Variable)(args...) = Operation(x, collect(Expression, args))

I guess something changed and the docs are not yet up to date. I would be grateful for some hints on how to do this now. I need to programatically create dependent variables depending on a dimensional parameter d like

x = [Variable(:x, i)(t) for i in 1:d]

┆Issue is synchronized with this Trello card by Unito

Simplify exp(x/2)^2

julia> exp(x/2)^2 |> simplify
exp(0.5x)^2

It would be great if this gave exp(x).

Function generation benchmark

using Markdown
using InteractiveUtils
using ModelingToolkit, NLsolve, RuntimeGeneratedFunctions, BenchmarkTools
function buildsystem(eqs, vars, params)
    ns = NonlinearSystem(eqs, vars, params)
    fexpr = generate_function(ns, vars, params, expression=Val{true})[2]
    jexpr = generate_jacobian(ns, vars, params, expression=Val{true})[2]
    return (ns, fexpr, jexpr)
end

splitpairs(p) = (ops = first.(p), vals = tuple(last.(p)...))
function solvesystem(eqs, vars, params)
    v = splitpairs(vars)
    p = splitpairs(params)
    ns, fexpr, jexpr = buildsystem(eqs, v.ops, p.ops)
    f! = @RuntimeGeneratedFunction(fexpr)
    j! = @RuntimeGeneratedFunction(jexpr)
    nlsolve((out, x) -> f!(out, x, p.vals),
            (out, x) -> j!(out, x, p.vals),
            collect(Float64, v.vals))
end

function getsystem(eqs, vars, params)
    v = splitpairs(vars)
    p = splitpairs(params)
    ns, fexpr, jexpr = buildsystem(eqs, v.ops, p.ops)
    f! = @RuntimeGeneratedFunction(fexpr)
    j! = @RuntimeGeneratedFunction(jexpr)
    (out, x) -> f!(out, x, p.vals),
    (out, x) -> j!(out, x, p.vals),
            collect(Float64, v.vals)
end

vars = @variables f Δp′ Re
params = @parameters ε ρ v Dₕ μ
eqs = [0 ~ -2*log(ε/(3.7Dₕ) + 2.51/(Re * √f)) - 1/√f,
       0 ~ ρ * v * Dₕ / μ - Re,
       0 ~ f * ρ/2 * v^2 / Dₕ - Δp′]
v0 = collect(vars) .=> 1.0
p0 = [ε => 1e-4, ρ => 1e3, v => 20, Dₕ => 0.06, μ => 9e-4]
function f_manual!(F, x, p)
    f, Δp′, Re = x
    ε, ρ, v, Dₕ, μ = p
    F[1] = -2*log(ε/(3.7Dₕ) + 2.51/(Re * √f)) - 1/√f
    F[2] = ρ * v * Dₕ / μ - Re
    F[3] = f * ρ/2 * v^2 / Dₕ - Δp′
end
f!,j!,xx = getsystem(eqs, v0, p0)
t_MTK = @btime nlsolve(f!, j!, ones(3))
t_manual = @btime nlsolve((F, x) -> f_manual!(F, x, (1e-4, 1e3, 20, 0.05, 9e-4)), ones(3))

t_MTK = @belapsed solvesystem($eqs, $v0, $p0)

using Profile
@profile solvesystem(eqs, v0, p0)
Juno.profiler()

function getsystem2(eqs, vars, params)
    v = splitpairs(vars)
    p = splitpairs(params)
    ns, fexpr, jexpr = buildsystem(eqs, v.ops, p.ops)
    f! = @RuntimeGeneratedFunction(fexpr)
    j! = @RuntimeGeneratedFunction(jexpr)
    f!,j!
end
f!,j! = getsystem2(eqs, v0, p0)
t_manual = @btime nlsolve((F, x) -> f!(F, x, (1e-4, 1e3, 20, 0.05, 9e-4)), (F, x) -> j!(F, x, (1e-4, 1e3, 20, 0.05, 9e-4)), ones(3))

┆Issue is synchronized with this Trello card by Unito

Solving rational equation

I have the following equation which I would like to solve for the parameters kp,ki,kd,T. The variable s is the polynomial variable.

Equation(ModelingToolkit.Constant(0), (1.0 + s ^ 2 + s) * inv(0.2 * s ^ 2 + s) ^ 1 - ((kp + ki / s) + kd * s) * (1 / (s * T + 1) ^ 2))

I can't figure out how to do this without supplying a numerical value for s. It should really hold for all s. I guess I want to do something like
https://mtk.sciml.ai/dev/tutorials/nonlinear/#Solving-Nonlinear-Systems-with-NLsolve-1
but I have to manually construct the system of equations?

Below is my code

using ModelingToolkit, ControlSystems

# @parameters s
@variables kp, ki, kd, T, s

function sym_pid(; kp=0., ki=0., kd=0., time=false, series=false)
    @variables s
    if series
        return time ? kp*(1 + 1/(ki*s) + kd*s) : kp*(1 + ki/s + kd*s)
    else
        return time ? kp + 1/(ki*s) + kd*s : kp + ki/s + kd*s
    end
end

Cs = sym_pid(;kp, ki, kd) # symbolic controller
Lf = 1/(s*T + 1)^2
C = Cs*Lf #|> simplify

kp_,ki_,kd_,T_ = 1,1,1,0.1
Cn = pid(;kp=kp_, ki=ki_, kd=kd_) # Numerical controller
Lf_ = tf(1,[T_, 1])^2
C_ = Cn*Lf_
"Convert a numeric transfer function to a symbolic"
function sym_tf(G::TransferFunction)
    n = reverse(numvec(G)[])
    sn = sum(s^i*n[i+1] for i in 0:length(n)-1)
    d = reverse(denvec(G)[])
    sd = sum(s^i*d[i+1] for i in 0:length(n)-1)
    simplify(sn/sd)
end

C_ns = sym_tf(C_)

## System of equations
eqs = [0 ~ C_ns - C]
ns = NonlinearSystem(eqs, [kp,ki,kd,T],[])
nlsys_func = generate_function(ns)[2]

f = eval(nlsys_func)
du = zeros(4); u = ones(4)
params = []
f(du,u,params) # Breaks here

j_func = generate_jacobian(ns)[2] # second is in-place
j! = eval(j_func)

using NLsolve
res = nlsolve((out, x) -> f(out, x, params), (out, x) -> j!(out, x, params), ones(4))

┆Issue is synchronized with this Trello card by Unito

Easier construction of typed variables

using Unitful
@variables x(t)::SArray{3,3} y::u"m"

┆Issue is synchronized with this Trello card by Unito

Pretty print symbolic matrix

Is there a way to pretty print a symbolic matrix in the REPL? For example, is it possible to show ModelingToolkit.Constant(0) as 0 in the Matrix output?

Parse `-1 <= x^2 + y^2 <= 1`

We should allow parsing -1 < x^2 + y^2 < 1.

Implementing linear constraints

Hi,
I'm interested in applying ModelingToolkit to model simple geometrical objects such as hyperplanes and halfspaces.

As a first step, the code below is supposed to check if the given ModelingToolkit.Operation describes either a halfspace or a hyperplane:

# represents the hyperplane <a, x> = b
struct Hyperplane{N, VN<:AbstractVector{N}}
    a::VN
    b::N
end

# represents the halfspace <a, x> <= b
struct HalfSpace{N, VN<:AbstractVector{N}}
    a::VN
    b::N
end

function is_hyperplane(expr::Operation)
   return expr.op == ==
end

function is_halfspace(expr::Operation)
    if expr.op ∈ [<, >]
        return true
    elseif expr.op == |
        x, y = expr.args
        if (is_halfspace(x) && is_hyperplane(y)) || (is_halfspace(y) && is_hyperplane(x))
            return true
        end
    end
    return false
end

# all these work
vars = @variables x y

@assert is_hyperplane(2x + 3y == 5)
@assert is_hyperplane(2x == 5y)
@assert is_halfspace(2x + 3y < 5)
@assert is_halfspace(2x + 3y > 5)
@assert is_halfspace(2x + 3y ≤ 5)
@assert is_halfspace(2x + 3y ≥ 5)
@assert is_halfspace(2x ≤ 5y - 1)
@assert is_halfspace(2x ≥ 5y - 1)
@assert is_halfspace(2x <= 5y - 1)
@assert is_halfspace(2x >= 5y - 1)

Here comes the question. What is the recommended way to get the coefficient values in such linear expressions? For instance,

expr = 2x + 3y == 5
a = expr.args[1]
b = expr.args[2]
sexpr = simplify(a - b)

-5 + 2x + 3y

is fine, and one can get the value 5 with sexpr.args[1].value. I can probably recursively check the arg field in expr to see if it's a value or an operation, but i was wondering if there's a higher level way of doing this... or it is planned to be added? Taking the derivative should also do the trick, but i didn't manage to properly use such functions,

@show derivative(cos(x), 1)
@show derivative(2x, 1)

derivative(cos(x), 1) = -(sin(x))   # ok
derivative(2x, 1) = (*)(x)   # ?
(*)(x)

Thanks in advance!

┆Issue is synchronized with this Trello card by Unito

Allow symbolic indices

julia> @variables x[1:3, 1:3]
(Operation[x₁ˏ₁ x₁ˏ₂ x₁ˏ₃; x₂ˏ₁ x₂ˏ₂ x₂ˏ₃; x₃ˏ₁ x₃ˏ₂ x₃ˏ₃],)

julia> @variables i, j
(i, j)

julia> x[i, j]
ERROR: ArgumentError: invalid index: i of type Operation
Stacktrace:
 [1] to_index(::Operation) at ./indices.jl:297
 [2] to_index(::Array{Operation,2}, ::Operation) at ./indices.jl:274
 [3] to_indices at ./indices.jl:325 [inlined]
 [4] to_indices at ./indices.jl:322 [inlined]
 [5] getindex(::Array{Operation,2}, ::Operation) at ./abstractarray.jl:1043
 [6] top-level scope at REPL[43]:1

This would be useful e.g. for matrix calculus.

┆Issue is synchronized with this Trello card by Unito

Support for Distributions in the IR

I'd like for transitions in a jump system to be stochastic e.g. something like

rate₁ = α*X
affect₁ = [X ~ X + rand(Poisson(β))]

The above doesn't work, as there's no Poisson(::Operation).

┆Issue is synchronized with this Trello card by Unito

substitute and substituter replacing division with inv and adding ^ 1 powers

@parameters t
@variables A(t) B(t)
c = A/B
substitute(c, [A.op(t) => A.op(), B.op(t) => B.op()])

replaces "A(t) / B(t)" with "inv(B) ^ 1 * A" and not "A / B".

@code_native on these two cases seems to indicate these are not equivalent when compiled

g1 = eval(build_function(substitute(c, [A.op(t) => A.op(), B.op(t) => B.op()]),[A,B],Variable[],t,expression=Val{true}))
g2 = eval(build_function(A.op()/B.op(),[A,B],Variable[],t,expression=Val{true}))

@code_native g1([x,y],Float64[],t)
        .section        __TEXT,__text,regular,pure_instructions
; ┌ @ build_function.jl:125 within `#39'
; │┌ @ array.jl:809 within `getindex'
        movq    (%rdi), %rax
        movabsq $5538124600, %rcx       ## imm = 0x14A191338
        vmovsd  (%rcx), %xmm0           ## xmm0 = mem[0],zero
; │└
; │┌ @ number.jl:199 within `inv'
; ││┌ @ float.jl:407 within `/'
        vdivsd  8(%rax), %xmm0, %xmm0
; │└└
; │┌ @ float.jl:405 within `*'
        vmulsd  (%rax), %xmm0, %xmm0
; │└
        retq
        nopl    (%rax,%rax)
; └vs


@code_native g2([x,y],Float64[],t)
        .section        __TEXT,__text,regular,pure_instructions
; ┌ @ build_function.jl:125 within `#41'
; │┌ @ array.jl:809 within `getindex'
        movq    (%rdi), %rax
        vmovsd  (%rax), %xmm0           ## xmm0 = mem[0],zero
; │└
; │┌ @ float.jl:407 within `/'
        vdivsd  8(%rax), %xmm0, %xmm0
; │└
        retq
        nopl    (%rax)
; └

┆Issue is synchronized with this Trello card by Unito

Double inequalities

Is there a reason that double inequalities as a < x < b shouldn't work?

julia> @variables x
(x,)

julia> 0 < x < 1
TypeError: non-boolean (Operation) used in boolean context

Stacktrace:
 [1] top-level scope at In[4]:1
 [2] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091

as a shortcut of

julia> (0 < x) & (x < 1)
(0 < x) & (x < 1)

Support for norms

using ModelingToolkit, LinearAlgebra
@variables x
norm(x)

First version checklist

• Move code from SymbolicUtils (see JuliaSymbolics/SymbolicUtils.jl#178)
• Move tests from SymbolicUtils
• Move Num, @variables from MTK (what happens to @parameters?)
• Move CAS functions -- Differential, jacobian, hessian, *_sparsity from MTK
• Tests from MTK

Define variable as 'real'

Hi, is there a way to define a new variable as 'real' so that for example

@variables v[1:2,1:2]
v'

produces

v₁ˏ₁  v₂ˏ₁  
v₁ˏ₂  v₂ˏ₂

instead of

conj(v₁ˏ₁)  conj(v₂ˏ₁)  
conj(v₁ˏ₂)  conj(v₂ˏ₂)

?

Something very simply done in matlab which I can't figure out using this pkg in jl (trying to switch).

Simplified API functions for creating Variables

Just a proposal; from the perspective of questions like
SciML/Catalyst.jl#281
and my own the other day, the current notation for building Syms and Terms programmatically seems complicated:

t = Num(Variable{ModelingToolkit.Parameter{Real}}(:t))  # t
x = Num(Variable{ModelingToolkit.FnType{Tuple{Any},Real}}(Symbol("x"),1))(t) # x_1(t)

Would there be objections to adding a few helper functions to simplify this for common use cases? I know @ChrisRackauckas prefers to keep the API smaller, but I suspect this is going to lead to many user questions in the future as understanding the typing here requires going in and reading the SymbolicUtils docs too. (Not to mention it is much more complicated than just remembering Operation and Variable with appropriate args.) Any thoughts?

No method matching -(::Variable...)

I couldn't really figure out what went wrong inside runge_kutta_discretize, but it appears that there's simply a missing method for -(::Variable{ModelingToolkit.Parameter{Number}})

MWE

using ModelingToolkit, OrdinaryDiffEq
const MT = ModelingToolkit
@parameters t, k, k₃, Jₘ, Jₐ, c, tx
@variables θ(t) q(t) u(t)
@derivatives D'~t

Δq = q - θ
Δv = D(q) - D(θ)
eqs = [
    D(D(q)) ~ (-k * Δq - k₃ * Δq^3 - c * Δv + u) / Jₘ
    D(D(θ)) ~ (+k * Δq + k₃ * Δq^3 + c * Δv) / Jₐ
]

loss = q^2 + θ^2 + u^2

sys = ODESystem(eqs, name=:DMM) |> ode_order_lowering
sysc = ControlSystem(loss, sys.eqs, t, states(sys), [u], MT.parameters(sys), name=:CDMM)

syso = runge_kutta_discretize(sysc, 0.01, (0, 0.5))

julia> syso = runge_kutta_discretize(sysc, 0.01, (0, 0.5))
ERROR: MethodError: no method matching -(::Variable{ModelingToolkit.Parameter{Number}})
Closest candidates are:
  -(::Any, ::VectorizationBase.Static{N}) where N at /home/fredrikb/.julia/packages/VectorizationBase/so9eG/src/static.jl:153
  -(::Any, ::ChainRulesCore.DoesNotExist) at /home/fredrikb/.julia/packages/ChainRulesCore/trzfY/src/differential_arithmetic.jl:26
  -(::Any, ::ChainRulesCore.Zero) at /home/fredrikb/.julia/packages/ChainRulesCore/trzfY/src/differential_arithmetic.jl:66
  ...
Stacktrace:
 [1] (::ModelingToolkit.var"#606#607")(::Array{Operation,1}, ::Array{Operation,1}, ::Array{Variable,1}, ::Operation) at ./array.jl:0
 [2] #invokelatest#1 at ./essentials.jl:710 [inlined]
 [3] invokelatest at ./essentials.jl:709 [inlined]
 [4] (::ModelingToolkit.var"#575#597"{Int64,Float64,ModelingToolkit.var"#606#607",Array{Array{Array{Operation,1},1},1},Array{Variable,1},Operation,Array{Array{Array{Operation,1},1},1}})(::Int64) at ./none:0
 [5] iterate at ./generator.jl:47 [inlined]
 [6] collect(::Base.Generator{UnitRange{Int64},ModelingToolkit.var"#575#597"{Int64,Float64,ModelingToolkit.var"#606#607",Array{Array{Array{Operation,1},1},1},Array{Variable,1},Operation,Array{Array{Array{Operation,1},1},1}}}) at ./array.jl:686
 [7] SciML/ModelingToolkit.jl#574 at ./none:0 [inlined]
 [8] iterate at ./generator.jl:47 [inlined]
 [9] collect(::Base.Generator{UnitRange{Int64},ModelingToolkit.var"#574#596"{Float64,Int64,ModelingToolkit.var"#606#607",Array{Array{Array{Operation,1},1},1},Array{Variable,1},Operation,Array{Array{Array{Operation,1},1},1}}}) at ./array.jl:686
 [10] runge_kutta_discretize(::ControlSystem, ::Float64, ::Tuple{Int64,Float64}; tab::DiffEqBase.ImplicitRKTableau{Array{Float64,2},Array{Float64,1}}) at /home/fredrikb/.julia/dev/ModelingToolkit/src/systems/control/controlsystem.jl:179
 [11] runge_kutta_discretize(::ControlSystem, ::Float64, ::Tuple{Int64,Float64}) at /home/fredrikb/.julia/dev/ModelingToolkit/src/systems/control/controlsystem.jl:135
 [12] top-level scope at none:1

build_function() fails when processing sparse Jacobian, but not sparse Hessian

I'm not sure if this is a genuine bug or whether I'm doing something daft (I'm still fairly new to Julia).

Anyway, I'm trying to use ModelingToolkit to generate sparse Jacobians and Hessians, and then use build_function() to obtain the resulting Julia expressions.

The originating mathematical functions I'd like to process will be fairly arbitrary, and depend both on spatial variables and parameters. In my use-case, I need to calculate sparse Jacobians and Hessians with respect to the spatial variables and, sometimes, then differentiate again with respect to the parametric variables.

Calculation of the sparse Jacobian and Hessian themselves seems to work fine. Also, subsequently taking gradients, elementwise, with respect to the parametric variables works fine.

However, build_function() fails when passed the Jacobian SparseMatrix (it only processes the first element) but not when passed the SparseMatrix containing the Hessian. For the Jacobian, build_function() also seems to leave the SparseMatrix in an inconsistent state.

build_function() succeeds elementwise on the Jacobian or if the SparseMatrix is converted to a dense matrix first. So, I suspect the issue may be in processing of the SparseMatrix.

I've attached two examples from the REPL. I've tried varying syntax a bit but cannot get the first example to work.

Example generating sparse Jacobian then passing to build_function() -- fails and leaves the sparse matrix in inconsistent state

julia> @variables x[1:100] pr[1:2] y
(Num[x₁, x₂, x₃, x₄, x₅, x₆, x₇, x₈, x₉, x₁₀  …  x₉₁, x₉₂, x₉₃, x₉₄, x₉₅, x₉₆, x₉₇, x₉₈, x₉₉, x₁₀₀], Num[pr₁, pr₂], y)

julia> function g(z::Vector, pr::Vector)
           return (
               z[1]^2 * z[2]^3 * sin(z[3]) +
               5 * z[4] +
               cos(z[10]) * exp(z[11]+z[12]) +
               pr[1]^2 * z[50]^2 +
               z[80] * z[95]^11
           )
       end
g (generic function with 1 method)

julia> y = g(x, pr)
5x₄ + x₈₀*(x₉₅^11) + cos(x₁₀)*exp(x₁₁ + x₁₂) + (pr₁^2)*(x₅₀^2) + sin(x₃)*(x₁^2)*(x₂^3)

julia> J_sparse = ModelingToolkit.sparsejacobian([y], x)
1×100 SparseMatrixCSC{Num,Int64} with 10 stored entries:
  [1,  1]  =  2x₁*sin(x₃)*(x₂^3)
  [1,  2]  =  3sin(x₃)*(x₁^2)*(x₂^2)
  [1,  3]  =  cos(x₃)*(x₁^2)*(x₂^3)
  [1,  4]  =  5
  [1, 10]  =  -1sin(x₁₀)*exp(x₁₁ + x₁₂)
  [1, 11]  =  cos(x₁₀)*exp(x₁₁ + x₁₂)
  [1, 12]  =  cos(x₁₀)*exp(x₁₁ + x₁₂)
  [1, 50]  =  2x₅₀*(pr₁^2)
  [1, 80]  =  x₉₅^11
  [1, 95]  =  11x₈₀*(x₉₅^10)

julia> expr_J_sparse = build_function(J_sparse, x, pr)
(:((var"##MTKArg#253", var"##MTKArg#254")->begin
          @inbounds begin
                  let (x₁, x₂, x₃, x₄, x₅, x₆, x₇, x₈, x₉, x₁₀, x₁₁, x₁₂, x₁₃, x₁₄, x₁₅, x₁₆, x₁₇, x₁₈, x₁₉, x₂₀, x₂₁, x₂₂, x₂₃, x₂₄, x₂₅, x₂₆, x₂₇, x₂₈, x₂₉, x₃₀, x₃₁, x₃₂, x₃₃, x₃₄, x₃₅, 
x₃₆, x₃₇, x₃₈, x₃₉, x₄₀, x₄₁, x₄₂, x₄₃, x₄₄, x₄₅, x₄₆, x₄₇, x₄₈, x₄₉, x₅₀, x₅₁, x₅₂, x₅₃, x₅₄, x₅₅, x₅₆, x₅₇, x₅₈, x₅₉, x₆₀, x₆₁, x₆₂, x₆₃, x₆₄, x₆₅, x₆₆, x₆₇, x₆₈, x₆₉, x₇₀, x₇₁, x₇₂, x₇₃, x₇₄, x₇₅, x₇₆, x₇₇, x₇₈, x₇₉, x₈₀, x₈₁, x₈₂, x₈₃, x₈₄, x₈₅, x₈₆, x₈₇, x₈₈, x₈₉, x₉₀, x₉₁, x₉₂, x₉₃, x₉₄, x₉₅, x₉₆, x₉₇, x₉₈, x₉₉, x₁₀₀, pr₁, pr₂) = (var"##MTKArg#253"[1], var"##MTKArg#253"[2], var"##MTKArg#253"[3], var"##MTKArg#253"[4], var"##MTKArg#253"[5], var"##MTKArg#253"[6], var"##MTKArg#253"[7], var"##MTKArg#253"[8], var"##MTKArg#253"[9], var"##MTKArg#253"[10], var"##MTKArg#253"[11], var"##MTKArg#253"[12], var"##MTKArg#253"[13], var"##MTKArg#253"[14], var"##MTKArg#253"[15], var"##MTKArg#253"[16], var"##MTKArg#253"[17], var"##MTKArg#253"[18], var"##MTKArg#253"[19], var"##MTKArg#253"[20], var"##MTKArg#253"[21], var"##MTKArg#253"[22], var"##MTKArg#253"[23], var"##MTKArg#253"[24], var"##MTKArg#253"[25], var"##MTKArg#253"[26], var"##MTKArg#253"[27], var"##MTKArg#253"[28], var"##MTKArg#253"[29], var"##MTKArg#253"[30], var"##MTKArg#253"[31], var"##MTKArg#253"[32], var"##MTKArg#253"[33], var"##MTKArg#253"[34], var"##MTKArg#253"[35], var"##MTKArg#253"[36], var"##MTKArg#253"[37], var"##MTKArg#253"[38], var"##MTKArg#253"[39], var"##MTKArg#253"[40], var"##MTKArg#253"[41], var"##MTKArg#253"[42], var"##MTKArg#253"[43], var"##MTKArg#253"[44], var"##MTKArg#253"[45], var"##MTKArg#253"[46], var"##MTKArg#253"[47], var"##MTKArg#253"[48], var"##MTKArg#253"[49], var"##MTKArg#253"[50], var"##MTKArg#253"[51], var"##MTKArg#253"[52], var"##MTKArg#253"[53], var"##MTKArg#253"[54], var"##MTKArg#253"[55], var"##MTKArg#253"[56], var"##MTKArg#253"[57], var"##MTKArg#253"[58], var"##MTKArg#253"[59], var"##MTKArg#253"[60], var"##MTKArg#253"[61], var"##MTKArg#253"[62], var"##MTKArg#253"[63], var"##MTKArg#253"[64], var"##MTKArg#253"[65], var"##MTKArg#253"[66], var"##MTKArg#253"[67], var"##MTKArg#253"[68], var"##MTKArg#253"[69], var"##MTKArg#253"[70], var"##MTKArg#253"[71], var"##MTKArg#253"[72], var"##MTKArg#253"[73], var"##MTKArg#253"[74], var"##MTKArg#253"[75], var"##MTKArg#253"[76], var"##MTKArg#253"[77], var"##MTKArg#253"[78], var"##MTKArg#253"[79], var"##MTKArg#253"[80], var"##MTKArg#253"[81], var"##MTKArg#253"[82], var"##MTKArg#253"[83], var"##MTKArg#253"[84], var"##MTKArg#253"[85], var"##MTKArg#253"[86], var"##MTKArg#253"[87], var"##MTKArg#253"[88], var"##MTKArg#253"[89], var"##MTKArg#253"[90], var"##MTKArg#253"[91], var"##MTKArg#253"[92], var"##MTKArg#253"[93], var"##MTKArg#253"[94], var"##MTKArg#253"[95], var"##MTKArg#253"[96], var"##MTKArg#253"[97], var"##MTKArg#253"[98], var"##MTKArg#253"[99], var"##MTKArg#253"[100], var"##MTKArg#254"[1], var"##MTKArg#254"[2])
                      SparseMatrixCSC{eltype(var"##MTKArg#253"), Int}(1, 1, [1, 2, 3, 4, 5, 5, 5, 5, 5, 5  …  10, 10, 10, 10, 11, 11, 11, 11, 11, 11], [1], [(*)(2, x₁, (sin)(x₃), (^)(x₂, 3))])
                  end
              end
      end), :((var"##MTIIPVar#256", var"##MTKArg#253", var"##MTKArg#254")->begin
          @inbounds begin
                  let (x₁, x₂, x₃, x₄, x₅, x₆, x₇, x₈, x₉, x₁₀, x₁₁, x₁₂, x₁₃, x₁₄, x₁₅, x₁₆, x₁₇, x₁₈, x₁₉, x₂₀, x₂₁, x₂₂, x₂₃, x₂₄, x₂₅, x₂₆, x₂₇, x₂₈, x₂₉, x₃₀, x₃₁, x₃₂, x₃₃, x₃₄, x₃₅, 
x₃₆, x₃₇, x₃₈, x₃₉, x₄₀, x₄₁, x₄₂, x₄₃, x₄₄, x₄₅, x₄₆, x₄₇, x₄₈, x₄₉, x₅₀, x₅₁, x₅₂, x₅₃, x₅₄, x₅₅, x₅₆, x₅₇, x₅₈, x₅₉, x₆₀, x₆₁, x₆₂, x₆₃, x₆₄, x₆₅, x₆₆, x₆₇, x₆₈, x₆₉, x₇₀, x₇₁, x₇₂, x₇₃, x₇₄, x₇₅, x₇₆, x₇₇, x₇₈, x₇₉, x₈₀, x₈₁, x₈₂, x₈₃, x₈₄, x₈₅, x₈₆, x₈₇, x₈₈, x₈₉, x₉₀, x₉₁, x₉₂, x₉₃, x₉₄, x₉₅, x₉₆, x₉₇, x₉₈, x₉₉, x₁₀₀, pr₁, pr₂) = (var"##MTKArg#253"[1], var"##MTKArg#253"[2], var"##MTKArg#253"[3], var"##MTKArg#253"[4], var"##MTKArg#253"[5], var"##MTKArg#253"[6], var"##MTKArg#253"[7], var"##MTKArg#253"[8], var"##MTKArg#253"[9], var"##MTKArg#253"[10], var"##MTKArg#253"[11], var"##MTKArg#253"[12], var"##MTKArg#253"[13], var"##MTKArg#253"[14], var"##MTKArg#253"[15], var"##MTKArg#253"[16], var"##MTKArg#253"[17], var"##MTKArg#253"[18], var"##MTKArg#253"[19], var"##MTKArg#253"[20], var"##MTKArg#253"[21], var"##MTKArg#253"[22], var"##MTKArg#253"[23], var"##MTKArg#253"[24], var"##MTKArg#253"[25], var"##MTKArg#253"[26], var"##MTKArg#253"[27], var"##MTKArg#253"[28], var"##MTKArg#253"[29], var"##MTKArg#253"[30], var"##MTKArg#253"[31], var"##MTKArg#253"[32], var"##MTKArg#253"[33], var"##MTKArg#253"[34], var"##MTKArg#253"[35], var"##MTKArg#253"[36], var"##MTKArg#253"[37], var"##MTKArg#253"[38], var"##MTKArg#253"[39], var"##MTKArg#253"[40], var"##MTKArg#253"[41], var"##MTKArg#253"[42], var"##MTKArg#253"[43], var"##MTKArg#253"[44], var"##MTKArg#253"[45], var"##MTKArg#253"[46], var"##MTKArg#253"[47], var"##MTKArg#253"[48], var"##MTKArg#253"[49], var"##MTKArg#253"[50], var"##MTKArg#253"[51], var"##MTKArg#253"[52], var"##MTKArg#253"[53], var"##MTKArg#253"[54], var"##MTKArg#253"[55], var"##MTKArg#253"[56], var"##MTKArg#253"[57], var"##MTKArg#253"[58], var"##MTKArg#253"[59], var"##MTKArg#253"[60], var"##MTKArg#253"[61], var"##MTKArg#253"[62], var"##MTKArg#253"[63], var"##MTKArg#253"[64], var"##MTKArg#253"[65], var"##MTKArg#253"[66], var"##MTKArg#253"[67], var"##MTKArg#253"[68], var"##MTKArg#253"[69], var"##MTKArg#253"[70], var"##MTKArg#253"[71], var"##MTKArg#253"[72], var"##MTKArg#253"[73], var"##MTKArg#253"[74], var"##MTKArg#253"[75], var"##MTKArg#253"[76], var"##MTKArg#253"[77], var"##MTKArg#253"[78], var"##MTKArg#253"[79], var"##MTKArg#253"[80], var"##MTKArg#253"[81], var"##MTKArg#253"[82], var"##MTKArg#253"[83], var"##MTKArg#253"[84], var"##MTKArg#253"[85], var"##MTKArg#253"[86], var"##MTKArg#253"[87], var"##MTKArg#253"[88], var"##MTKArg#253"[89], var"##MTKArg#253"[90], var"##MTKArg#253"[91], var"##MTKArg#253"[92], var"##MTKArg#253"[93], var"##MTKArg#253"[94], var"##MTKArg#253"[95], var"##MTKArg#253"[96], var"##MTKArg#253"[97], var"##MTKArg#253"[98], var"##MTKArg#253"[99], var"##MTKArg#253"[100], var"##MTKArg#254"[1], var"##MTKArg#254"[2])
                      (var"##MTIIPVar#256").nzval[1] = (*)(2, x₁, (sin)(x₃), (^)(x₂, 3))
                  end
              end
          nothing
      end))

julia> J_sparse
1×100 SparseMatrixCSC{Num,Int64} with 10 stored entries:Error showing value of type SparseMatrixCSC{Num,Int64}:
ERROR: BoundsError: attempt to access 1-element Array{Int64,1} at index [1:10]
Stacktrace:
 [1] throw_boundserror(::Array{Int64,1}, ::Tuple{UnitRange{Int64}}) at .\abstractarray.jl:541
 [2] checkbounds at .\abstractarray.jl:506 [inlined]
 [3] getindex at .\array.jl:815 [inlined]
 [4] show(::IOContext{REPL.Terminals.TTYTerminal}, ::SparseMatrixCSC{Num,Int64}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\SparseArrays\src\sparsematrix.jl:230  
 [5] show(::IOContext{REPL.Terminals.TTYTerminal}, ::MIME{Symbol("text/plain")}, ::SparseMatrixCSC{Num,Int64}) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\SparseArrays\src\sparsematrix.jl:197
 [6] display(::REPL.REPLDisplay, ::MIME{Symbol("text/plain")}, ::Any) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:214
 [7] display(::REPL.REPLDisplay, ::Any) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:218
 [8] display(::Any) at .\multimedia.jl:328
 [9] #invokelatest#1 at .\essentials.jl:710 [inlined]
 [10] invokelatest at .\essentials.jl:709 [inlined]
 [11] print_response(::IO, ::Any, ::Bool, ::Bool, ::Any) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:238
 [12] print_response(::REPL.AbstractREPL, ::Any, ::Bool, ::Bool) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:223
 [13] (::REPL.var"#do_respond#54"{Bool,Bool,VSCodeServer.var"#40#41"{REPL.LineEditREPL,REPL.LineEdit.Prompt},REPL.LineEditREPL,REPL.LineEdit.Prompt})(::Any, ::Any, ::Any) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:822
 [14] #invokelatest#1 at .\essentials.jl:710 [inlined]
 [15] invokelatest at .\essentials.jl:709 [inlined]
 [16] run_interface(::REPL.Terminals.TextTerminal, ::REPL.LineEdit.ModalInterface, ::REPL.LineEdit.MIState) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\LineEdit.jl:2355
 [17] run_frontend(::REPL.LineEditREPL, ::REPL.REPLBackendRef) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:1144
 [18] (::REPL.var"#38#42"{REPL.LineEditREPL,REPL.REPLBackendRef})() at .\task.jl:356

Example generating sparse Hessian then passing to build_function() -- seems ok

julia> @variables x[1:100] pr[1:2] y
(Num[x₁, x₂, x₃, x₄, x₅, x₆, x₇, x₈, x₉, x₁₀  …  x₉₁, x₉₂, x₉₃, x₉₄, x₉₅, x₉₆, x₉₇, x₉₈, x₉₉, x₁₀₀], Num[pr₁, pr₂], y)

julia> function g(z::Vector, pr::Vector)
           return (
               z[1]^2 * z[2]^3 * sin(z[3]) +
               5 * z[4] +
               cos(z[10]) * exp(z[11]+z[12]) +
               pr[1]^2 * z[50]^2 +
               z[80] * z[95]^11
           )
       end
g (generic function with 1 method)

julia> y = g(x, pr)
5x₄ + x₈₀*(x₉₅^11) + cos(x₁₀)*exp(x₁₁ + x₁₂) + (pr₁^2)*(x₅₀^2) + sin(x₃)*(x₁^2)*(x₂^3)

julia> H_sparse = ModelingToolkit.sparsehessian(y, x)
100×100 SparseMatrixCSC{Num,Int64} with 22 stored entries:
  [1 ,  1]  =  2sin(x₃)*(x₂^3)
  [2 ,  1]  =  6x₁*sin(x₃)*(x₂^2)
  [3 ,  1]  =  2x₁*cos(x₃)*(x₂^3)
  [1 ,  2]  =  6x₁*sin(x₃)*(x₂^2)
  [2 ,  2]  =  6x₂*sin(x₃)*(x₁^2)
  [3 ,  2]  =  3cos(x₃)*(x₁^2)*(x₂^2)
  [1 ,  3]  =  2x₁*cos(x₃)*(x₂^3)
  [2 ,  3]  =  3cos(x₃)*(x₁^2)*(x₂^2)
  ⋮
  [11, 11]  =  cos(x₁₀)*exp(x₁₁ + x₁₂)
  [12, 11]  =  cos(x₁₀)*exp(x₁₁ + x₁₂)
  [10, 12]  =  -1sin(x₁₀)*exp(x₁₁ + x₁₂)
  [11, 12]  =  cos(x₁₀)*exp(x₁₁ + x₁₂)
  [12, 12]  =  cos(x₁₀)*exp(x₁₁ + x₁₂)
  [50, 50]  =  2(pr₁^2)
  [95, 80]  =  11(x₉₅^10)
  [80, 95]  =  11(x₉₅^10)
  [95, 95]  =  110x₈₀*(x₉₅^9)

julia> expr_hess = build_function(H_sparse, x, pr)
(:((var"##MTKArg#253", var"##MTKArg#254")->begin
          @inbounds begin
                  let (x₁, x₂, x₃, x₄, x₅, x₆, x₇, x₈, x₉, x₁₀, x₁₁, x₁₂, x₁₃, x₁₄, x₁₅, x₁₆, x₁₇, x₁₈, x₁₉, x₂₀, x₂₁, x₂₂, x₂₃, x₂₄, x₂₅, x₂₆, x₂₇, x₂₈, x₂₉, x₃₀, x₃₁, x₃₂, x₃₃, x₃₄, x₃₅, 
x₃₆, x₃₇, x₃₈, x₃₉, x₄₀, x₄₁, x₄₂, x₄₃, x₄₄, x₄₅, x₄₆, x₄₇, x₄₈, x₄₉, x₅₀, x₅₁, x₅₂, x₅₃, x₅₄, x₅₅, x₅₆, x₅₇, x₅₈, x₅₉, x₆₀, x₆₁, x₆₂, x₆₃, x₆₄, x₆₅, x₆₆, x₆₇, x₆₈, x₆₉, x₇₀, x₇₁, x₇₂, x₇₃, x₇₄, x₇₅, x₇₆, x₇₇, x₇₈, x₇₉, x₈₀, x₈₁, x₈₂, x₈₃, x₈₄, x₈₅, x₈₆, x₈₇, x₈₈, x₈₉, x₉₀, x₉₁, x₉₂, x₉₃, x₉₄, x₉₅, x₉₆, x₉₇, x₉₈, x₉₉, x₁₀₀, pr₁, pr₂) = (var"##MTKArg#253"[1], var"##MTKArg#253"[2], var"##MTKArg#253"[3], var"##MTKArg#253"[4], var"##MTKArg#253"[5], var"##MTKArg#253"[6], var"##MTKArg#253"[7], var"##MTKArg#253"[8], var"##MTKArg#253"[9], var"##MTKArg#253"[10], var"##MTKArg#253"[11], var"##MTKArg#253"[12], var"##MTKArg#253"[13], var"##MTKArg#253"[14], var"##MTKArg#253"[15], var"##MTKArg#253"[16], var"##MTKArg#253"[17], var"##MTKArg#253"[18], var"##MTKArg#253"[19], var"##MTKArg#253"[20], var"##MTKArg#253"[21], var"##MTKArg#253"[22], var"##MTKArg#253"[23], var"##MTKArg#253"[24], var"##MTKArg#253"[25], var"##MTKArg#253"[26], var"##MTKArg#253"[27], var"##MTKArg#253"[28], var"##MTKArg#253"[29], var"##MTKArg#253"[30], var"##MTKArg#253"[31], var"##MTKArg#253"[32], var"##MTKArg#253"[33], var"##MTKArg#253"[34], var"##MTKArg#253"[35], var"##MTKArg#253"[36], var"##MTKArg#253"[37], var"##MTKArg#253"[38], var"##MTKArg#253"[39], var"##MTKArg#253"[40], var"##MTKArg#253"[41], var"##MTKArg#253"[42], var"##MTKArg#253"[43], var"##MTKArg#253"[44], var"##MTKArg#253"[45], var"##MTKArg#253"[46], var"##MTKArg#253"[47], var"##MTKArg#253"[48], var"##MTKArg#253"[49], var"##MTKArg#253"[50], var"##MTKArg#253"[51], var"##MTKArg#253"[52], var"##MTKArg#253"[53], var"##MTKArg#253"[54], var"##MTKArg#253"[55], var"##MTKArg#253"[56], var"##MTKArg#253"[57], var"##MTKArg#253"[58], var"##MTKArg#253"[59], var"##MTKArg#253"[60], var"##MTKArg#253"[61], var"##MTKArg#253"[62], var"##MTKArg#253"[63], var"##MTKArg#253"[64], var"##MTKArg#253"[65], var"##MTKArg#253"[66], var"##MTKArg#253"[67], var"##MTKArg#253"[68], var"##MTKArg#253"[69], var"##MTKArg#253"[70], var"##MTKArg#253"[71], var"##MTKArg#253"[72], var"##MTKArg#253"[73], var"##MTKArg#253"[74], var"##MTKArg#253"[75], var"##MTKArg#253"[76], var"##MTKArg#253"[77], var"##MTKArg#253"[78], var"##MTKArg#253"[79], var"##MTKArg#253"[80], var"##MTKArg#253"[81], var"##MTKArg#253"[82], var"##MTKArg#253"[83], var"##MTKArg#253"[84], var"##MTKArg#253"[85], var"##MTKArg#253"[86], var"##MTKArg#253"[87], var"##MTKArg#253"[88], var"##MTKArg#253"[89], var"##MTKArg#253"[90], var"##MTKArg#253"[91], var"##MTKArg#253"[92], var"##MTKArg#253"[93], var"##MTKArg#253"[94], var"##MTKArg#253"[95], var"##MTKArg#253"[96], var"##MTKArg#253"[97], var"##MTKArg#253"[98], var"##MTKArg#253"[99], var"##MTKArg#253"[100], var"##MTKArg#254"[1], var"##MTKArg#254"[2])
                      SparseMatrixCSC{eltype(var"##MTKArg#253"), Int}(100, 100, [1, 4, 7, 10, 10, 10, 10, 10, 10, 10  …  21, 21, 21, 21, 23, 23, 23, 23, 23, 23], [1, 2, 3, 1, 2, 3, 1, 2, 3, 10  …  10, 11, 12, 10, 11, 12, 50, 95, 80, 95], [(*)(2, (sin)(x₃), (^)(x₂, 3)), (*)(6, x₁, (sin)(x₃), (^)(x₂, 2)), (*)(2, x₁, (cos)(x₃), (^)(x₂, 3)), (*)(6, x₁, (sin)(x₃), (^)(x₂, 2)), (*)(6, x₂, (sin)(x₃), (^)(x₁, 2)), (*)(3, (cos)(x₃), (^)(x₁, 2), (^)(x₂, 2)), (*)(2, x₁, (cos)(x₃), (^)(x₂, 3)), (*)(3, (cos)(x₃), (^)(x₁, 2), (^)(x₂, 2)), (*)(-1, (sin)(x₃), (^)(x₁, 2), (^)(x₂, 3)), (*)(-1, (cos)(x₁₀), (exp)((+)(x₁₁, x₁₂))), (*)(-1, (sin)(x₁₀), (exp)((+)(x₁₁, x₁₂))), (*)(-1, (sin)(x₁₀), (exp)((+)(x₁₁, x₁₂))), (*)(-1, (sin)(x₁₀), (exp)((+)(x₁₁, x₁₂))), (*)((cos)(x₁₀), (exp)((+)(x₁₁, x₁₂))), (*)((cos)(x₁₀), (exp)((+)(x₁₁, x₁₂))), (*)(-1, (sin)(x₁₀), (exp)((+)(x₁₁, x₁₂))), (*)((cos)(x₁₀), (exp)((+)(x₁₁, x₁₂))), (*)((cos)(x₁₀), (exp)((+)(x₁₁, x₁₂))), (*)(2, (^)(pr₁, 2)), (*)(11, (^)(x₉₅, 10)), (*)(11, (^)(x₉₅, 10)), (*)(110, x₈₀, (^)(x₉₅, 9))])
                  end
              end
      end), :((var"##MTIIPVar#256", var"##MTKArg#253", var"##MTKArg#254")->begin
          @inbounds begin
                  let (x₁, x₂, x₃, x₄, x₅, x₆, x₇, x₈, x₉, x₁₀, x₁₁, x₁₂, x₁₃, x₁₄, x₁₅, x₁₆, x₁₇, x₁₈, x₁₉, x₂₀, x₂₁, x₂₂, x₂₃, x₂₄, x₂₅, x₂₆, x₂₇, x₂₈, x₂₉, x₃₀, x₃₁, x₃₂, x₃₃, x₃₄, x₃₅, 
x₃₆, x₃₇, x₃₈, x₃₉, x₄₀, x₄₁, x₄₂, x₄₃, x₄₄, x₄₅, x₄₆, x₄₇, x₄₈, x₄₉, x₅₀, x₅₁, x₅₂, x₅₃, x₅₄, x₅₅, x₅₆, x₅₇, x₅₈, x₅₉, x₆₀, x₆₁, x₆₂, x₆₃, x₆₄, x₆₅, x₆₆, x₆₇, x₆₈, x₆₉, x₇₀, x₇₁, x₇₂, x₇₃, x₇₄, x₇₅, x₇₆, x₇₇, x₇₈, x₇₉, x₈₀, x₈₁, x₈₂, x₈₃, x₈₄, x₈₅, x₈₆, x₈₇, x₈₈, x₈₉, x₉₀, x₉₁, x₉₂, x₉₃, x₉₄, x₉₅, x₉₆, x₉₇, x₉₈, x₉₉, x₁₀₀, pr₁, pr₂) = (var"##MTKArg#253"[1], var"##MTKArg#253"[2], var"##MTKArg#253"[3], var"##MTKArg#253"[4], var"##MTKArg#253"[5], var"##MTKArg#253"[6], var"##MTKArg#253"[7], var"##MTKArg#253"[8], var"##MTKArg#253"[9], var"##MTKArg#253"[10], var"##MTKArg#253"[11], var"##MTKArg#253"[12], var"##MTKArg#253"[13], var"##MTKArg#253"[14], var"##MTKArg#253"[15], var"##MTKArg#253"[16], var"##MTKArg#253"[17], var"##MTKArg#253"[18], var"##MTKArg#253"[19], var"##MTKArg#253"[20], var"##MTKArg#253"[21], var"##MTKArg#253"[22], var"##MTKArg#253"[23], var"##MTKArg#253"[24], var"##MTKArg#253"[25], var"##MTKArg#253"[26], var"##MTKArg#253"[27], var"##MTKArg#253"[28], var"##MTKArg#253"[29], var"##MTKArg#253"[30], var"##MTKArg#253"[31], var"##MTKArg#253"[32], var"##MTKArg#253"[33], var"##MTKArg#253"[34], var"##MTKArg#253"[35], var"##MTKArg#253"[36], var"##MTKArg#253"[37], var"##MTKArg#253"[38], var"##MTKArg#253"[39], var"##MTKArg#253"[40], var"##MTKArg#253"[41], var"##MTKArg#253"[42], var"##MTKArg#253"[43], var"##MTKArg#253"[44], var"##MTKArg#253"[45], var"##MTKArg#253"[46], var"##MTKArg#253"[47], var"##MTKArg#253"[48], var"##MTKArg#253"[49], var"##MTKArg#253"[50], var"##MTKArg#253"[51], var"##MTKArg#253"[52], var"##MTKArg#253"[53], var"##MTKArg#253"[54], var"##MTKArg#253"[55], var"##MTKArg#253"[56], var"##MTKArg#253"[57], var"##MTKArg#253"[58], var"##MTKArg#253"[59], var"##MTKArg#253"[60], var"##MTKArg#253"[61], var"##MTKArg#253"[62], var"##MTKArg#253"[63], var"##MTKArg#253"[64], var"##MTKArg#253"[65], var"##MTKArg#253"[66], var"##MTKArg#253"[67], var"##MTKArg#253"[68], var"##MTKArg#253"[69], var"##MTKArg#253"[70], var"##MTKArg#253"[71], var"##MTKArg#253"[72], var"##MTKArg#253"[73], var"##MTKArg#253"[74], var"##MTKArg#253"[75], var"##MTKArg#253"[76], var"##MTKArg#253"[77], var"##MTKArg#253"[78], var"##MTKArg#253"[79], var"##MTKArg#253"[80], var"##MTKArg#253"[81], var"##MTKArg#253"[82], var"##MTKArg#253"[83], var"##MTKArg#253"[84], var"##MTKArg#253"[85], var"##MTKArg#253"[86], var"##MTKArg#253"[87], var"##MTKArg#253"[88], var"##MTKArg#253"[89], var"##MTKArg#253"[90], var"##MTKArg#253"[91], var"##MTKArg#253"[92], var"##MTKArg#253"[93], var"##MTKArg#253"[94], var"##MTKArg#253"[95], var"##MTKArg#253"[96], var"##MTKArg#253"[97], var"##MTKArg#253"[98], var"##MTKArg#253"[99], var"##MTKArg#253"[100], var"##MTKArg#254"[1], var"##MTKArg#254"[2])
                      (var"##MTIIPVar#256").nzval[1] = (*)(2, (sin)(x₃), (^)(x₂, 3))
                      (var"##MTIIPVar#256").nzval[2] = (*)(6, x₁, (sin)(x₃), (^)(x₂, 2))
                      (var"##MTIIPVar#256").nzval[3] = (*)(2, x₁, (cos)(x₃), (^)(x₂, 3))
                      (var"##MTIIPVar#256").nzval[4] = (*)(6, x₁, (sin)(x₃), (^)(x₂, 2))
                      (var"##MTIIPVar#256").nzval[5] = (*)(6, x₂, (sin)(x₃), (^)(x₁, 2))
                      (var"##MTIIPVar#256").nzval[6] = (*)(3, (cos)(x₃), (^)(x₁, 2), (^)(x₂, 2))
                      (var"##MTIIPVar#256").nzval[7] = (*)(2, x₁, (cos)(x₃), (^)(x₂, 3))
                      (var"##MTIIPVar#256").nzval[8] = (*)(3, (cos)(x₃), (^)(x₁, 2), (^)(x₂, 2))
                      (var"##MTIIPVar#256").nzval[9] = (*)(-1, (sin)(x₃), (^)(x₁, 2), (^)(x₂, 3))
                      (var"##MTIIPVar#256").nzval[10] = (*)(-1, (cos)(x₁₀), (exp)((+)(x₁₁, x₁₂)))
                      (var"##MTIIPVar#256").nzval[11] = (*)(-1, (sin)(x₁₀), (exp)((+)(x₁₁, x₁₂)))
                      (var"##MTIIPVar#256").nzval[12] = (*)(-1, (sin)(x₁₀), (exp)((+)(x₁₁, x₁₂)))
                      (var"##MTIIPVar#256").nzval[13] = (*)(-1, (sin)(x₁₀), (exp)((+)(x₁₁, x₁₂)))
                      (var"##MTIIPVar#256").nzval[14] = (*)((cos)(x₁₀), (exp)((+)(x₁₁, x₁₂)))
                      (var"##MTIIPVar#256").nzval[15] = (*)((cos)(x₁₀), (exp)((+)(x₁₁, x₁₂)))
                      (var"##MTIIPVar#256").nzval[16] = (*)(-1, (sin)(x₁₀), (exp)((+)(x₁₁, x₁₂)))
                      (var"##MTIIPVar#256").nzval[17] = (*)((cos)(x₁₀), (exp)((+)(x₁₁, x₁₂)))
                      (var"##MTIIPVar#256").nzval[18] = (*)((cos)(x₁₀), (exp)((+)(x₁₁, x₁₂)))
                      (var"##MTIIPVar#256").nzval[19] = (*)(2, (^)(pr₁, 2))
                      (var"##MTIIPVar#256").nzval[20] = (*)(11, (^)(x₉₅, 10))
                      (var"##MTIIPVar#256").nzval[21] = (*)(11, (^)(x₉₅, 10))
                      (var"##MTIIPVar#256").nzval[22] = (*)(110, x₈₀, (^)(x₉₅, 9))
                  end
              end
          nothing
      end))`

Feature Request: complex variables

Making Num <: Real (#446) lead to stuff like

using ModelingToolkit
@variables x
term = x + im

ERROR: MethodError: +(::Num, ::Complex{Bool}) is ambiguous. Candidates:
  +(x::Num, z::Complex) in ModelingToolkit at /Users/hw/.julia/packages/ModelingToolkit/FPuGy/src/ModelingToolkit.jl:70
  +(x::Real, z::Complex{Bool}) in Base at complex.jl:300
  +(a::Num, b::Number) in ModelingToolkit at /Users/hw/.julia/packages/SymbolicUtils/HntGZ/src/methods.jl

For system creation this problem can be averted by using the Syms instead of Num:

term = x.val + im

But x.val::Sym{Real} so this might lead to problems when you acctually want to put complex numbers into the generated function. Does build_function assume the variables to be <:Real?

It would be really nice to have support for complex variables/parameters.

Inversion of symbolic StaticMatrix broken for sizes larger than 4x4

Starting at 5x5 the StaticArrays package uses LU decomposition with pivoting for inverses and other operations, which is incompatible with symbolic matrices. There is already an override of lu() in operations.jl for AbstractMatrices not to use pivoting but does not seem to be activated for StaticMatrices.

  @variables M[1:5,1:5]
  
  julia> N = SMatrix{5,5}(M)
  5×5 SArray{Tuple{5,5},Operation,2,25} with indices SOneTo(5)×SOneTo(5):
  M₁ˏ₁  M₁ˏ₂  M₁ˏ₃  M₁ˏ₄  M₁ˏ₅
  M₂ˏ₁  M₂ˏ₂  M₂ˏ₃  M₂ˏ₄  M₂ˏ₅
  M₃ˏ₁  M₃ˏ₂  M₃ˏ₃  M₃ˏ₄  M₃ˏ₅
  M₄ˏ₁  M₄ˏ₂  M₄ˏ₃  M₄ˏ₄  M₄ˏ₅
  M₅ˏ₁  M₅ˏ₂  M₅ˏ₃  M₅ˏ₄  M₅ˏ₅

  julia> inv(N)
  ERROR: TypeError: non-boolean (Operation) used in boolean context
  Stacktrace:
  [1] __lu(::SArray{Tuple{5,5},Operation,2,25}, ::Val{true}) at /Users/dcabecinhas/.julia/packages/StaticArrays/mlIi1/src/lu.jl:150
  [2] macro expansion at /Users/dcabecinhas/.julia/packages/StaticArrays/mlIi1/src/lu.jl:70 [inlined]
  [3] _lu(::SArray{Tuple{5,5},Operation,2,25}, ::Val{true}, ::Bool) at /Users/dcabecinhas/.julia/packages/StaticArrays/mlIi1/src/lu.jl:68
  [4] lu(::SArray{Tuple{5,5},Operation,2,25}, ::Val{true}; check::Bool) at /Users/dcabecinhas/.julia/packages/StaticArrays/mlIi1/src/lu.jl:42
  [5] lu(::SArray{Tuple{5,5},Operation,2,25}, ::Val{true}) at /Users/dcabecinhas/.julia/packages/StaticArrays/mlIi1/src/lu.jl:42 (repeats 2 times)
  [6] macro expansion at /Users/dcabecinhas/.julia/packages/StaticArrays/mlIi1/src/inv.jl:77 [inlined]
  [7] _inv at /Users/dcabecinhas/.julia/packages/StaticArrays/mlIi1/src/inv.jl:73 [inlined]
  [8] inv(::SArray{Tuple{5,5},Operation,2,25}) at /Users/dcabecinhas/.julia/packages/StaticArrays/mlIi1/src/inv.jl:5
  [9] top-level scope at REPL[34]:1

julia> N = SMatrix{4,4}(M[1:4,1:4])
 4×4 SArray{Tuple{4,4},Operation,2,16} with indices SOneTo(4)×SOneTo(4):
 M₁ˏ₁  M₁ˏ₂  M₁ˏ₃  M₁ˏ₄
 M₂ˏ₁  M₂ˏ₂  M₂ˏ₃  M₂ˏ₄
 M₃ˏ₁  M₃ˏ₂  M₃ˏ₃  M₃ˏ₄
 M₄ˏ₁  M₄ˏ₂  M₄ˏ₃  M₄ˏ₄

 julia> inv(N)
 4×4 SArray{Tuple{4,4},Operation,2,16} with indices SOneTo(4)×SOneTo(4):
   ((M₂ˏ₂ * (M₃ˏ₃ * M₄ˏ₄ - M₄ˏ₃ * M₃ˏ₄) - M₂ˏ₃ * (M₃ˏ₂ * M₄ˏ₄ - M₄ˏ₂ * M₃ˏ₄)) [...]

Sparse vector support

There are zeros in this generated code in contrast to the matrix case:

using ModelingToolkit
using SparseArrays
@variables a,b 
c = [a, 0.0, 0.0, 0.0, b]
ModelingToolkit.build_function(sparse(c),[a,b])
d = [a 0.0 0.0 b]
ModelingToolkit.build_function(sparse(d),[a,b])

Trig identities

Simplify
cos(A) * cos(B) - sin(A) * sin(B) ⟶ cos(A + B)

and all the rest of the high-school trig identities.

┆Issue is synchronized with this Trello card by Unito

Support for variable-size arrays

I know that variables can be fixed size array, but what would be the best way to use modeling toolkit with the following function:

function purcell_single_ode(du,u,p,t)
    Δ,g,κ,γ = p
    α = u[1]
    dα = -γ/2 * α
    for i in 1:length(Δ)
        β = u[1+i]
        dα += -1im*(g[i]*β)
        du[1+i] = -1im*(Δ[i]*β + g[i]*α) - κ[i]/2*β
    end
    du[1] = dα
    return nothing
end

fun = ODEFunction{true}(purcell_single_ode);

p = (Δ=[5,-8,2], g=[0.4,0.08,0.1]*2π, κ=[2,3,1], γ=0.02)

prob = ODEProblem(fun,ComplexF64[1;zeros(length(p.Δ))],(0.0,5),p)
sol = solve(prob,DP8(),abstol=1e-9,reltol=1e-9);

Here, I can easily change the problem dimensions just by passing another set of parameters.

I am not sure how to do that in ModelingToolkit without rebuilding the whole function with different @variables u[1:N+1], Δ[1:N], ... each time.

┆Issue is synchronized with this Trello card by Unito

Better latexification of equations

A set of equations:

[
    0 ~ n1.F - D * Dk * max(n1.v - n2.v, 0)
    0 ~ n2.F - D * max(n2.v - n1.v, 0)
]

renders in LaTeX as

It would be nice to replace the primes with either arrows or dots, and indicate multiplication explicitly by \cdot (or something).

┆Issue is synchronized with this Trello card by Unito