sparsity! reporting 0 stored entries

[ferrolho]

I have been trying to use SparsityDetection.jl for computing the sparsity pattern of a Julia function but I seem to be missing something. I can run the example from the README.md just fine, and the result is as expected; but I can't get the sparsity for a more complex function.

I have created a Notebook to demonstrate my issue:
https://nbviewer.jupyter.org/gist/ferrolho/5516d663963aeb67b231d4aec4d2efa1

Out [10] is where I'd expect to get the sparsity as illustrated in Out [9].

[ChrisRackauckas]

Hey, thanks for providing an example. I think we will need to really dig into it in order to find out what's going on!

[ferrolho]

Okay, here is another very simple function:

function EulerStep(dx, x)
    Δt = 0.001
    qᵢ = x[ 1: 7]
    q̇ᵢ = x[ 8:14]
    q̈ᵢ = x[15:21]

    dx[ 1: 7] = qᵢ + Δt * q̇ᵢ
    dx[ 8:14] = q̇ᵢ + Δt * q̈ᵢ

    nothing
end

After which, calling

input = rand(3*7)
output = zeros(2*7)
sparsity_pattern = sparsity!(EulerStep, output, input)

yields

Explored path: SparsityDetection.Path(Bool[], 1)
14×21 SparseArrays.SparseMatrixCSC{Bool,Int64} with 0 stored entries

Am I doing something terribly wrong?

[ferrolho]

It is weird because for another simple function it does work:

function g1(dx, x)
    dx[1] = 2x[1] + √x[3]
    dx[2] = 3x[1] + sin(x[4])
    dx[3] = 4x[1] + cos(x[5])
    nothing
end

input = rand(5)
output = zeros(3)
sparsity_pattern = sparsity!(g1, output, input)

Explored path: SparsityDetection.Path(Bool[], 1)
3×5 SparseMatrixCSC{Bool,Int64} with 6 stored entries:
  [1, 1]  =  1
  [2, 1]  =  1
  [3, 1]  =  1
  [1, 3]  =  1
  [2, 4]  =  1
  [3, 5]  =  1

[ferrolho]

Rewriting the function so as to eliminate range assignments solves the problem:

function EulerStep(dx, x)
    Δt = 0.001

    for i = 1:7
        dx[i  ] = x[i  ] + Δt * x[i+ 7]
        dx[i+7] = x[i+7] + Δt * x[i+14]
    end

    nothing
end

input = rand(3*7)
output = zeros(2*7)
sparsity_pattern = sparsity!(EulerStep, output, input)

Explored path: SparsityDetection.Path(Bool[], 1)
14×21 SparseMatrixCSC{Bool,Int64} with 28 stored entries:
  [1 ,  1]  =  1
  [2 ,  2]  =  1
  [3 ,  3]  =  1
  [4 ,  4]  =  1
  [5 ,  5]  =  1
  [6 ,  6]  =  1
  [7 ,  7]  =  1
  [1 ,  8]  =  1
  [8 ,  8]  =  1
  [2 ,  9]  =  1
  [9 ,  9]  =  1
  [3 , 10]  =  1
  ⋮
  [5 , 12]  =  1
  [12, 12]  =  1
  [6 , 13]  =  1
  [13, 13]  =  1
  [7 , 14]  =  1
  [14, 14]  =  1
  [8 , 15]  =  1
  [9 , 16]  =  1
  [10, 17]  =  1
  [11, 18]  =  1
  [12, 19]  =  1
  [13, 20]  =  1
  [14, 21]  =  1

Are range assignments not supported?

[ferrolho]

Besides range assignments, doing qᵢ = x[1:7] seems to break the exploration as well:

function EulerStep(dx, x)
    Δt = 0.001
    qᵢ = x[ 1: 7]
    q̇ᵢ = x[ 8:14]
    q̈ᵢ = x[15:21]

    for i = 1:7
        dx[i  ] = qᵢ[i] + Δt * q̇ᵢ[i]
        dx[i+7] = q̇ᵢ[i] + Δt * q̈ᵢ[i]
    end

    nothing
end

input = rand(3*7)
output = zeros(2*7)
sparsity_pattern = sparsity!(EulerStep, output, input)

Explored path: SparsityDetection.Path(Bool[], 1)
14×21 SparseMatrixCSC{Bool,Int64} with 0 stored entries

[ChrisRackauckas]

Awesome, thanks for the MWE. Yes, it looks like range indexing needs to be specialized on.

[ChrisRackauckas]

Indeed it's missing a term combination: https://github.com/JuliaDiffEq/SparsityDetection.jl/blob/master/src/terms.jl

[ferrolho]

I am trying to add this to SparsityDetection.jl/src/terms.jl but I need some help.

function Base.getindex(comb1::TermCombination, I...)
    TermCombination(comb1.terms[I...])
end

Wouldn't something along these lines suffice? How can I inspect this function during runtime? I may have gotten the signature wrong because I added a print statement to the body but it never got printed.

Edit
I do get it to run if I specifically call it with a ::TermCombination, as in:

t1 = TermCombination(Set([Dict([(1, 2), (2, 4)])]))
getindex(t1, 1:2)

But it does not seem to be called with sparsity!().

Also, I don't understand the underlying data structure of struct TermCombination and I could not find any docs. So, say I call getindex(t1, 1:2), is the correct thing to do to return a TermCombination similar to t1 but containing only the pairs whose keys are 1 and 2?

[ChrisRackauckas]

@shashi can you take a look at this?

[shashi]

TermCombination is used for Hessian sparsity. It looks like the range index is calling a unsafe_copy* function which we don't intercept... I'm playing with this now.

[jfeist]

The same happens with ArrayPartition from RecursiveArrayTools:

using SparsityDetection
using RecursiveArrayTools
function rhs(dcs,cs,ps,t=0.)
    (c1, c2) = cs.x
    (dc1, dc2) = dcs.x
    (g1,g2) = ps
    for j=1:size(dc1,1)
        dc1[j] = g1[j] * c1[j]
    end
    for j=1:size(dc2,1)
        dc2[j] = g1[j] * c2[j] + g2[j] * c1[j]
    end
    dcs
end
N = 30
cs = ArrayPartition(rand(N),rand(N))
dcs = similar(cs)
ps = rand.((N,N))
sparsity_pattern = sparsity!(rhs,dcs,cs,ps)

gives

Explored path: SparsityDetection.Path(Bool[], 1)
60×60 SparseArrays.SparseMatrixCSC{Bool,Int64} with 0 stored entries

while

function rhs(dcs,cs,ps,t=0.)
    (g1,g2) = ps
    N1 = length(dcs) ÷ 2
    for j=1:N1
        dcs[j]    = g1[j] * cs[j]
        dcs[j+N1] = g1[j] * cs[j+N1] + g2[j] * cs[j]
    end
    dcs
end
N = 30
cs = rand(2N)
dcs = similar(cs)
ps = rand.((N,N))
sparsity_pattern = sparsity!(rhs,dcs,cs,ps)

works fine:

Explored path: SparsityDetection.Path(Bool[], 1)
60×60 SparseArrays.SparseMatrixCSC{Bool,Int64} with 90 stored entries: [...]