add sparsity overloads to AbstractSystems

docs/src/systems/AbstractSystem.md

@@ -56,6 +56,10 @@ generate_factorized_W
 generate_hessian
 ```
 
+Additionally, `jacobian_sparsity(sys)` and `hessian_sparsity(sys)`
+exist on the appropriate systems for fast generation of the sparsity
+patterns via an abstract interpretation without requiring differentiation.
+
 ## Problem Constructors
 
 At the end, the system types have `DEProblem` constructors, like `ODEProblem`,

docs/src/systems/NonlinearSystem.md

@@ -19,6 +19,7 @@ NonlinearSystem
 ```julia
 calculate_jacobian
 generate_jacobian
+jacobian_sparsity
 ```
 
 ## Problem Constructors

docs/src/systems/ODESystem.md

@@ -28,6 +28,7 @@ calculate_factorized_W
 generate_jacobian
 generate_tgrad
 generate_factorized_W
+jacobian_sparsity
 ```
 
 ## Problem Constructors

docs/src/systems/OptimizationSystem.md

@@ -21,6 +21,7 @@ calculate_gradient
 calculate_hessian
 generate_gradient
 generate_hessian
+hessian_sparsity
 ```
 
 ## Problem Constructors

docs/src/systems/SDESystem.md

@@ -24,6 +24,7 @@ calculate_factorized_W
 generate_jacobian
 generate_tgrad
 generate_factorized_W
+jacobian_sparsity
 ```
 
 ## Problem Constructors

src/systems/diffeqs/abstractodesystem.jl

@@ -88,6 +88,10 @@ function calculate_massmatrix(sys::AbstractODESystem; simplify=true)
     M == I ? I : M
 end
 
+jacobian_sparsity(sys::AbstractODESystem) =
+    jacobian_sparsity([eq.rhs for eq ∈ equations(sys)],
+                      [dv(sys.iv()) for dv in states(sys)])
+
 function DiffEqBase.ODEFunction(sys::AbstractODESystem, args...; kwargs...)
     ODEFunction{true}(sys, args...; kwargs...)
 end

src/systems/nonlinear/nonlinearsystem.jl

@@ -67,6 +67,10 @@ function generate_function(sys::NonlinearSystem, vs = states(sys), ps = paramete
                           conv = AbstractSysToExpr(sys), kwargs...)
 end
 
+jacobian_sparsity(sys::NonlinearSystem) =
+    jacobian_sparsity([eq.rhs for eq ∈ equations(sys)],
+                      [dv() for dv in states(sys)])
+
 """
 ```julia
 function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem,u0map,tspan,

src/systems/optimization/optimizationsystem.jl

@@ -73,6 +73,9 @@ end
 equations(sys::OptimizationSystem) = isempty(sys.systems) ? sys.op : sys.op + reduce(+,namespace_operation.(sys.systems))
 namespace_operation(sys::OptimizationSystem) = namespace_operation(sys.op,sys.name,nothing)
 
+hessian_sparsity(sys::OptimizationSystem) =
+    hessian_sparsity(sys.op,[dv() for dv in states(sys)])
+
 """
 ```julia
 function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem,

test/nonlinearsystem.jl

@@ -1,5 +1,5 @@
 using ModelingToolkit, StaticArrays, LinearAlgebra
-using DiffEqBase
+using DiffEqBase, SparseArrays
 using Test
 
 canonequal(a, b) = isequal(simplify(a), simplify(b))
@@ -62,4 +62,8 @@ eqs = [0 ~ σ*a,
 ns = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])
 nlsys_func = generate_function(ns, [x,y,z], [σ,ρ,β])
 jac = calculate_jacobian(ns)
+
+@test ModelingToolkit.jacobian_sparsity(ns).colptr == sparse(jac).colptr
+@test ModelingToolkit.jacobian_sparsity(ns).rowval == sparse(jac).rowval
+
 jac = generate_jacobian(ns)

test/odesystem.jl

@@ -1,6 +1,6 @@
 using ModelingToolkit, StaticArrays, LinearAlgebra
 using OrdinaryDiffEq
-using DiffEqBase
+using DiffEqBase, SparseArrays
 using Test
 
 # Define some variables
@@ -135,6 +135,9 @@ eqs = [D(x) ~ σ*a,
 de = ODESystem(eqs)
 generate_function(de, [x,y,z], [σ,ρ,β])
 jac = calculate_jacobian(de)
+@test ModelingToolkit.jacobian_sparsity(de).colptr == sparse(jac).colptr
+@test ModelingToolkit.jacobian_sparsity(de).rowval == sparse(jac).rowval
+
 f = ODEFunction(de, [x,y,z], [σ,ρ,β])
 
 @derivatives D'~t

test/optimizationsystem.jl

@@ -1,4 +1,5 @@
-using ModelingToolkit
+using ModelingToolkit, SparseArrays
+
 @variables x y
 @parameters a b
 loss = (a - x)^2 + b * (y - x^2)^2
@@ -19,3 +20,4 @@ calculate_hessian(combinedsys)
 generate_function(combinedsys)
 generate_gradient(combinedsys)
 generate_hessian(combinedsys)
+ModelingToolkit.hessian_sparsity(combinedsys)