Small numeric error introduced in sparse jacobians

[TorkelE]

Something that popped up. The actual error is really small, so probably not practically serious. Still might be worth investigating what is going one.

using Catalyst, Test

rn = @reaction_network begin
    (p,d), 0 <--> (X,Y,Z)
    k1, X + Y --> XY
    k2, X + 2Z --> XZ2
    k3, Y3 +X2 --> Y3Z2
    k4, X + Y + Z --> XYZ
    k5, XZ2 + Y3Z2 --> XY3Z4
    k6, XYZ + XYZ --> X2Y2Z2
    d, (XY3Z4, X2Y2Z2) --> 0
    X + Y, V --> 0
    k7/(1 + t), 2V --> V2
    Z, V2 --> 0
end

# Create problem
u0 = [sp => rand() for sp in species(rn)]
ps = [p => rand() for p in parameters(rn)]
prob = ODEProblem(rn, u0, 0.0, ps; jac = true, sparse = true)

# Compute sparse jacobian.
J = deepcopy(prob.f.jac_prototype)
prob.f.jac(J, prob.u0, prob.p, 1.0)
@test J == prob.f.jac(prob.u0, prob.p, 1.0) # Don't pass anymore.
@test J ≈ prob.f.jac(prob.u0, prob.p, 1.0) # Passes (difference between jacobians is very small).

# Compute normal Jacobian.
probd = ODEProblem(rn, u0, 0.0, ps; jac = true, sparse = false)
Jd = zeros(length(u0), length(u0))
probd.f.jac(Jd, prob.u0, prob.p, 1.0)
@test Jd == probd.f.jac(prob.u0, prob.p, 1.0) # This works for dense jacobians.
@test Matrix(J) == Jd # Now fails.
@test Matrix(J) ≈ Jd # works.

[TorkelE]

Some potentially useful outputs:

julia> J .- prob.f.jac(prob.u0, prob.p, 1.0)
13×13 SparseArrays.SparseMatrixCSC{Float64, Int64} with 1 stored entry:
 -1.11022e-16   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅            ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅            ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅            ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅            ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅            ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅            ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅            ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅            ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅            ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅            ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅            ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅            ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 

julia> Matrix(J) .- Jd
13×13 Matrix{Float64}:
 -1.11022e-16  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
  0.0          0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> prob.f.jac_prototype
13×13 SparseArrays.SparseMatrixCSC{Float64, Int64} with 39 stored entries:
 0.0       3.0e-323  6.0e-323   ⋅         ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         ⋅         ⋅             ⋅ 
 5.0e-324  4.0e-323  4.0e-323   ⋅         ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         ⋅         ⋅             ⋅ 
 1.0e-323  4.4e-323  7.4e-323   ⋅         ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         ⋅         ⋅             ⋅ 
 1.5e-323  5.0e-323   ⋅        9.4e-323   ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         ⋅         ⋅             ⋅ 
 1.5e-323   ⋅        8.0e-323   ⋅        1.0e-322    ⋅          ⋅         1.43e-322   ⋅          ⋅         ⋅         ⋅             ⋅ 
  ⋅         ⋅         ⋅         ⋅         ⋅         1.1e-322   1.3e-322    ⋅          ⋅          ⋅         ⋅         ⋅             ⋅ 
  ⋅         ⋅         ⋅         ⋅         ⋅         1.1e-322   1.33e-322   ⋅          ⋅          ⋅         ⋅         ⋅             ⋅ 
  ⋅         ⋅         ⋅         ⋅        1.04e-322  1.24e-322  1.4e-322   1.5e-322    ⋅          ⋅         ⋅         ⋅             ⋅ 
 1.5e-323  5.4e-323  8.0e-323   ⋅         ⋅          ⋅          ⋅          ⋅         1.43e-322   ⋅         ⋅         ⋅             ⋅ 
  ⋅         ⋅         ⋅         ⋅        1.1e-322    ⋅          ⋅         1.5e-322    ⋅         1.7e-322   ⋅         ⋅             ⋅ 
  ⋅         ⋅         ⋅         ⋅         ⋅          ⋅          ⋅          ⋅         1.63e-322   ⋅        1.4e-322   ⋅             ⋅ 
 3.0e-323  5.4e-323   ⋅         ⋅         ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         ⋅        1.0e-323       ⋅ 
  ⋅         ⋅        7.4e-323   ⋅         ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         ⋅        6.44477e-310  0.0

julia> ModelingToolkit.jacobian_sparsity(prob.f.sys)
13×13 SparseArrays.SparseMatrixCSC{Bool, Int64} with 38 stored entries:
 1  1  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 1  1  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 1  1  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 1  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 1  ⋅  1  ⋅  1  ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  1  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  1  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  1  1  1  1  ⋅  ⋅  ⋅  ⋅  ⋅
 1  1  1  ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  1  ⋅  ⋅  1  ⋅  1  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅  1  ⋅  ⋅
 1  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅
 ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  1

[ChrisRackauckas]

Show prob.f.jac(prob.u0, prob.p, 1.0)

[ChrisRackauckas]

Note that [4,4] is the hard zero, but it's entry 1 that's different so this is not due to #3535 related issues. This shows that is properly setting [4,4] as hard zero in the generated jac

[TorkelE]

As pointed out by Chris, setting cse = false in the problems solves this.

[TorkelE]

Some further debug info

probd = ODEProblem(rn, u0, 0.0, ps; jac = true, sparse = false, cse=true)
j1cse = probd.f.jac(Jd, prob.u0, prob.p, 1.0)
j2cse = probd.f.jac(prob.u0, prob.p, 1.0)

probd = ODEProblem(rn, u0, 0.0, ps; jac = true, sparse = false, cse=false)
j1ncse = probd.f.jac(Jd, prob.u0, prob.p, 1.0)
j2ncse = probd.f.jac(prob.u0, prob.p, 1.0)

j1cse == j1ncse # true
j2cse == j2ncse # false

prob = ODEProblem(rn, u0, 0.0, ps; jac = true, sparse = false, cse=true)
j1cse = probd.f.jac(J, prob.u0, prob.p, 1.0)
j2cse = probd.f.jac(prob.u0, prob.p, 1.0)

prob = ODEProblem(rn, u0, 0.0, ps; jac = true, sparse = false, cse=false)
j1ncse = probd.f.jac(J, prob.u0, prob.p, 1.0)
j2ncse = probd.f.jac(prob.u0, prob.p, 1.0)

j1cse == j1ncse # true
j2cse == j2ncse # true

[ChrisRackauckas]

@AayushSabharwal in-place and out of place Jacobians aren't CSEing the same, and in-place looks like it doesn't change its numerics with CSE true/false, can you check the function expression and see if it's actually CSEing?

[ChrisRackauckas]

I have a feeling that we might be putting out the same equation in both cases, but then LLVM's CSE pass might do something differently with sparse Jacobian allocations around and that might be a thing here, but we should check our expressions are exactly the same first

[AayushSabharwal]

I'm guessing this is due to

ModelingToolkit.jl/src/systems/diffeqs/abstractodesystem.jl

Lines 142 to 148 in 8e69dc9

	function assert_jac_length_header(sys)
	W = W_sparsity(sys)
	identity, expr -> Func([expr.args...], [], LiteralExpr(quote
	@assert $(findnz)($(expr.args[1]))[1:2] == $(findnz)($W)[1:2]
	$(expr.body)
	end))
	end

. It turns the entire function body into an Expr before CSE gets its hands on it, so CSE doesn't actually do anything.