Sparse `stoichmat` dispatches?

[yewalenikhil65]

Does it make sense to add Sparse stoichmat dispatches, now that we have many stoich functions ?
If yes, how do we decide for sparse or dense version as all these stoichmat functions take reaction system as argument ?
Maybe numreactions(rn) could be deciding factor ?(If yes, how much threshold ?)

[yewalenikhil65]

Or else, we can also pass a Bool as flag for sparsity in the function. default set to true ?

[isaacsas]

I guess we could make a three argument form for the various matrix functions. Maybe something like

function netstoichmat(::Type{SparseMatrixCSC}, rn::ReactionSystem; smap=speciesmap(rn))
...
end
function netstoichmat(::Type{Matrix}, rn::ReactionSystem; smap = speciesmap(rn))
# the existing function code
end
function netstoichmat(rn::ReactionSystem, smap=speciesmap(rn))
   netstoichmat(Matrix, rn; smap=smap)
end

[yewalenikhil65]

function netstoichmat(rn::ReactionSystem, smap=speciesmap(rn))
   netstoichmat(Matrix, rn; smap=smap)
end

This part I didn't understand properly.

But
I was thinking along following lines

function netstoichmat(rn::ReactionSystem, sparsity::Bool;smap=species(rn)) 

if sparsity == true
  #sparse code
else
  # existing function code
end


end

Also we should not use spzeros way in # sparse code part., but try constructing using sparse(Is, Js, Vs) way as it's more effective @btime wise

[yewalenikhil65]

hey @isaacsas , did a little benchmarking on large reaction mentioned here #378 (comment)
using spzeros vs sparse(I,J,K)

function substoichmat(rn::ReactionSystem; smap=speciesmap(rn))
    smat = spzeros(Int,numreactions(rn),numspecies(rn))
    for (k,rx) in enumerate(reactions(rn))
        stoich = rx.substoich
        for (i,sub) in enumerate(rx.substrates)
            smat[k,smap[sub]] = stoich[i]
        end
    end
    smat
end

julia> @benchmark substoichmat(rn)
BenchmarkTools.Trial: 5530 samples with 1 evaluation.
 Range (min … max):  759.505 μs …  14.391 ms  ┊ GC (min … max): 0.00% … 90.41%
 Time  (median):     799.192 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   894.657 μs ± 949.304 μs  ┊ GC (mean ± σ):  7.46% ±  6.61%

  █                                                              
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁  ▁ ▁
  1.21 ms          Histogram: frequency by time          770 μs <

 Memory estimate: 464.20 KiB, allocs estimate: 8247.

function substoichmat(::Type{SparseMatrixCSC}, rn::ReactionSystem; smap=speciesmap(rn))
    Is=Int64[];  Js=Int64[];  Vs=Int64[];
    for (k,rx) in enumerate(reactions(rn))
        stoich = rx.substoich
        for (i,sub) in enumerate(rx.substrates)
            push!(Is,k)
            push!(Js, smap[sub])
            push!(Vs, stoich[i])
        end
    end
    smat = sparse(Is,Js,Vs,numreactions(rn),numspecies(rn))
end

julia> @benchmark substoichmat(SparseMatrixCSC,rn)
BenchmarkTools.Trial: 6021 samples with 1 evaluation.
 Range (min … max):  706.090 μs …  13.417 ms  ┊ GC (min … max): 0.00% … 88.78%
 Time  (median):     738.844 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   822.599 μs ± 876.015 μs  ┊ GC (mean ± σ):  7.49% ±  6.60%

  █                                                              
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▁
  1.09 ms          Histogram: frequency by time          794 μs <

 Memory estimate: 477.27 KiB, allocs estimate: 7681.

[yewalenikhil65]

I will make a PR based on this now!

[yewalenikhil65]

hey @isaacsas ,
I am ready with this PR, but it needs SparseArrays as dependency, i tried to add in Catalyst.jl and Project.toml file.. But cant seem to figure this out how to add dependency

(@v1.6) pkg> activate .
  Activating environment at `~/Catalyst.jl-myCatalyst/Project.toml`
ERROR: Could not parse project: TOML Parser error:
/home/chetan/Catalyst.jl-myCatalyst/Project.toml:20:14 error: key already has a value
  SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[isaacsas]

Take a look at the Project.toml for OrdinaryDiffEq.jl to see how to add it as a dependency.

[yewalenikhil65]

Take a look at the Project.toml for OrdinaryDiffEq.jl to see how to add it as a dependency.

Thanks 👍 That helped me

[yewalenikhil65]

#387