Sparsity detection missing values due to temporary storage?

[ChrisRackauckas]

using OrdinaryDiffEq, LinearAlgebra, Test

# Define the constants for the PDE
const α₂ = 1.0
const α₃ = 1.0
const β₁ = 1.0
const β₂ = 1.0
const β₃ = 1.0
const r₁ = 1.0
const r₂ = 1.0
const D = 100.0
const γ₁ = 0.1
const γ₂ = 0.1
const γ₃ = 0.1
const N = 32
const X = reshape([i for i in 1:N for j in 1:N],N,N)
const Y = reshape([j for i in 1:N for j in 1:N],N,N)
const α₁ = 1.0.*(X.>=4*N/5)

const Mx = Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1])
const My = copy(Mx)
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0

# Define the initial condition as normal arrays
u0 = zeros(N,N,3)

const MyA = zeros(N,N);
const AMx = zeros(N,N);
const DA = zeros(N,N);
# Define the discretized PDE as an ODE function
function f(du,u,p,t)
   A = @view  u[:,:,1]
   B = @view  u[:,:,2]
   C = @view  u[:,:,3]
  dA = @view du[:,:,1]
  dB = @view du[:,:,2]
  dC = @view du[:,:,3]
  mul!(MyA,My,A)
  mul!(AMx,A,Mx)
  @. DA = D*(MyA + AMx)
  @. dA = DA + α₁ - β₁*A - r₁*A*B + r₂*C
  @. dB = α₂ - β₂*B - r₁*A*B + r₂*C
  @. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end

using SparsityDetection, SparseArrays
input = rand(N,N,3)
output = similar(input)
sparsity_pattern = sparsity!(f,output,input,nothing,0.0)
jac_sparsity = Float64.(sparse(sparsity_pattern))

[shashi]

It seems to be an intersection of #15 and global access:

julia> const Z = Int[1,2]
2-element Array{Int64,1}:
 1
 2
# This works! --
julia> sparsity!([1,2],[2,3]) do du,u
                 Z[:] .= u[:]
                du[1] = Z[1]
              end
Explored path: SparsityDetection.Path(Bool[], 1)
2×2 SparseMatrixCSC{Bool,Int64} with 1 stored entry:
  [1, 1]  =  1

# This doesn't --
julia> sparsity!([1,2],[2,3]) do du,u
                 Z[:] .= u[:]
                du[:] .= Z[:]
              end
Explored path: SparsityDetection.Path(Bool[], 1)
2×2 SparseMatrixCSC{Bool,Int64} with 0 stored entries

[ChrisRackauckas]

Interesting. So even the mul! works, it's just @. DA = D*(MyA + AMx)?

[shashi]

Yes.

[shashi]

I fixed the above problem,

mul! didn't need to work, and it doesn't.

const N = 12

const My = Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1])
My[1,2] = 2.0
My[end,end-1] = 2.0

function f(du,u,p,t)
  mul!(du,My,u)
end

using SparsityDetection
input = rand(N,N)
output = similar(input)
sparsity_pattern = sparsity!(f,output,input,nothing,0.0)
# empty

[shashi]

Ok that was actually not a problem, I screwed up while moving between branches and found a commit in the reflog which I had missed out in the #18. :/ :/

Does this look right to you for a 16x16 problem?

[ChrisRackauckas]

It's a bit hard to do this one in my head, but that seems about right. A good test would be to calculate the full Jacobian with forward diff and then zero it out, or MTK. Of course, you can only do that on the small problems, and MTK will quickly out of memory if it's like 32x32 haha.