[WIP]- `subnetworks`, `linkagedeficiency`

[yewalenikhil65]

Creating here separate issue to keep track of stuff in addition to #394,
I think it's useful to know which reactions form subnetworks in the reaction network based on linkage class.

"""
   Given a `ReactionSystem` , returns a Vector of  Vector of Reactions that form   
subnetworks based on linkage class
"""
function subnetworks( rs::ReactionSystem;sparse::Bool=false,
               rmap::Vector{Vector{Pair{Int64, Int64}}} = collect(values(reactioncomplexmap(rs))),
               r::Vector{Reaction} = reactions(rs),
               ig::SimpleDiGraph{Int64} = incidencematgraph(reactioncomplexes(rs;sparse=sparse)[2]),
               lc::Vector{Vector{Int64}} = linkageclasses(ig) )

   subrs = Vector{Vector{Reaction}}()
   for i in 1:length(lc)
      rxind = unique(vcat([map(first, rmap[lc[i]][j]) for j in 1:length(lc[i])]...))
      push!(subrs , r[rxind])
   end
   subrs
end

rn = @reaction_network begin
    k₁, 2A --> B
    k₂, A --> C
    k₃, C --> D
    k₄, B + D --> E
end k₁ k₂ k₃ k₄

julia> subnetworks(rn)
3-element Vector{Vector{Reaction}}:
 [Reaction{Any, Int64}: k₁, 2A --> B]
 [Reaction{Any, Int64}: k₂, A --> C, Reaction{Any, Int64}: k₃, C --> D]
 [Reaction{Any, Int64}: k₄, B + D --> E]

Implying there are three subnetworks that do not share reaction-complexes.. This is also evident from

complexgraph(rn)

[yewalenikhil65]

This should form a basis for finding deficiency of each linkage class, (PR #394 involves deficiency of whole network).

These studies could create a good API for future work in steady-state analysis of certain reaction network.

[yewalenikhil65]

"""
   linkagedeficiency(network ;sparse ,subrn,ig ,lc )

Given a ReactionSystem, returns deficiency of each linkage class associated
with reaction network defined as,
deficiency = (number of reaction complexes in linkage class - 1 -
                  dim. of stochiometric subspace of linkage class)

Note: The kwargs as follows
- subrn,  subnetworks corresponding to each linkage class
- ig, incidencegraph from incidence matrix of the network
- lc , linkage classes from incidencegraph of the network
- sparse, a Bool for constructing a dense or sparse netstoichmat ,incidence matrix, etc.
"""

function linkagedeficiency(rs::ReactionSystem,sparse::Bool = false,
            subrn::Vector{Vector{Reaction}} = subnetworks(rs,sparse=sparse),
            ig::SimpleDiGraph{Int64} = incidencematgraph(reactioncomplexes(rs;sparse=sparse)[2]),
            lc::Vector{Vector{Int64}} = linkageclasses(ig))

   t  =ModelingToolkit.get_iv(rs)
   δ = zeros(Int,length(subrn))
   for i in 1:length(subrn)
      specs = unique(vcat([[subrn[i][j].substrates;subrn[i][j].products] for j in 1:length(subrn[i])]...))
      @named rxs = ReactionSystem(subrn[i], t, specs, Any[])
      ns_sub = netstoichmat(rxs;sparse=sparse)
      δ[i] = length(lc[i]) - 1 - AbstractAlgebra.rank(AbstractAlgebra.matrix(AbstractAlgebra.ZZ, ns_sub))
   end
   δ
end

[yewalenikhil65]

done with #396