Automatic sparsity detection

[ChrisRackauckas]

Currently needs overload of base functionality.

using ModelingToolkit
using Test

# Define some variables
@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives Dx'~x Dy'~y Dz'~z

### Automatic conversion to symbolic

function lorenz(du,u,p,t)
 du[1] = p[1]*(u[2]-u[1])
 du[2] = u[1]*(p[2]-u[3]) - u[2]
 du[3] = u[1]*u[2] - p[3]*u[3]
end
u = [x,y,z]
du = similar(u)
p = [σ,ρ,β]
lorenz(du,u,p,t)

du

J = [Dx(du[1]) Dy(du[1]) Dz(du[1])
     Dx(du[2]) Dy(du[2]) Dz(du[2])
     Dx(du[3]) Dy(du[3]) Dz(du[3])]
J = expand_derivatives.(J)

using SparseArrays
function SparseArrays.SparseMatrixCSC(M::Matrix{T}) where {T}
    @assert !SparseArrays.has_offset_axes(M)
    I = findall(!iszero, M)
    eltypeTiI = [i[1] for i in I]
    eltypeTiJ = [i[2] for i in I]
    eltypeTvV = [M[i] for i in I]
    return SparseArrays.sparse_IJ_sorted!(eltypeTiI, eltypeTiJ, eltypeTvV, size(M)...)
end
SparseMatrixCSC(J)

[asprionj]

Once this has happened, all matrices (Jacobians etc.) generated with MT will be automatically sparse, right?

[ChrisRackauckas]

Not necessarily. We need to find the right interface. I think it should be if the user requests it, or if we detect sparsity is low enough (or structured)

[asprionj]

Agreed -- and sorry, just now found the related discussion over at DiffEqBiological.jl
. So this is just about enabling sparsity detection, if beneficial or requested.