Merge 3a92b99 into faf7189

src/Catalyst.jl

@@ -44,7 +44,8 @@ export species, params, reactions, speciesmap, paramsmap, numspecies, numreactio
 export make_empty_network, addspecies!, addparam!, addreaction!
 export dependants, dependents, substoichmat, prodstoichmat, netstoichmat
 export conservationlaws, conservedquantities
-
+export reaction_complexes, reaction_rates, complex_stoich_matrix, complex_incidence_matrix, complex_outgoing_matrix
+
 # for Latex printing of ReactionSystems
 include("latexify_recipes.jl")
 

src/networkapi.jl

@@ -183,6 +183,108 @@ function netstoichmat(rn; smap=speciesmap(rn))
 end
 
 
+######################## some reacion complexes matrices and reaction rates ###############################
+"""
+    reaction_complexes(network)
+Given a [`ReactionSystem`](@ref), return a vector of set of pairs ,
+where pair consists on (species indicies) => (stoichiometric coefficients)
+denoting the composition of reaction complexes
+Notes:
+- Empty pair denotes 𝛷 complex.
+"""
+function reaction_complexes(rn)
+    r = reactions(rn);
+    smap = speciesmap(rn);
+    complexes_mat = Set{Pair{Int64, Int64}}[]
+    for i in 1:numreactions(rn)
+        push!(complexes_mat, Set(Pair.([getindex(smap,
+                r[i].substrates[j]) for j in 1:length(r[i].substrates)],
+                    r[i].substoich)));
+
+        push!(complexes_mat, Set(Pair.([getindex(smap,
+                r[i].products[j]) for j in 1:length(r[i].products)],
+                    r[i].prodstoich)))
+    end
+    return unique!(complexes_mat) # deleting duplicate complexes
+end
+
+
+"""
+    reaction_rates(network)
+Given a [`ReactionSystem`](@ref), returns a vector of rate of reactions
+"""
+function reaction_rates(rn)
+    return [r.rate for r in reactions(rn)]
+end
+
+
+"""
+    complex_stoich_matrix(network)
+Given a [`ReactionSystem`](@ref), return a matrix with positive entries of size,
+num_of_species x num_of_complexes, where the non-zero positive entries in the k'th
+column denote stoichiometric coefficients of the species participating in
+compositon of the reaction complex
+"""
+function complex_stoich_matrix(rn; cmp_mat = reaction_complexes(rn))
+    Z = zeros(Int64, numspecies(rn), length(cmp_mat));
+    for i in 1:length(cmp_mat)
+        for j in 1:length(cmp_mat[i])
+            if collect(cmp_mat[i])[j].first != Int64[] # NOT A NULL COMPLEX
+                Z[collect(cmp_mat[i])[j].first, i] = collect(cmp_mat[i])[j].second
+            end
+        end
+    end
+    return Z
+end
+
+
+
+"""
+    complex_incidence_matrix(network)
+Given a [`ReactionSystem`](@ref), return a matrix of size,
+num_of_complexes x num_of_reactions, where ,
+Bᵢⱼ = -1, if the i'th complex is the substrate of the j'th reaction,
+       1, if the i'th complex is the product of the j'th reaction,
+       0, otherwise
+"""
+function complex_incidence_matrix(rn; cmp_mat = reaction_complexes(rn))
+    r = reactions(rn);
+    smap = speciesmap(rn);
+    B = zeros(Int64, length(cmp_mat), numreactions(rn));
+    for i in 1:numreactions(rn)
+        substrates_pair = Set(Pair.([getindex(smap,
+                    r[i].substrates[j]) for j in 1:length(r[i].substrates)],
+                    r[i].substoich))
+
+        products_pair = Set(Pair.([getindex(smap,
+                    r[i].products[j]) for j in 1:length(r[i].products)],
+                    r[i].prodstoich))
+
+        substrates_cmp_index = findall( x -> x == substrates_pair, cmp_mat)
+        products_cmp_index = findall(x -> x == products_pair, cmp_mat)
+        B[substrates_cmp_index,i] .= -1
+        B[products_cmp_index,i] .= 1
+    end
+    return B
+end
+
+
+"""
+    complex_incidence_matrix(network)
+Given a [`ReactionSystem`](@ref), return a matrix of size,
+num_of_complexes x num_of_reactions, where ,
+where,  (Bᵢⱼ being complex_incidence_matrix(network))
+Δᵢⱼ = 0,    if Bᵢⱼ =1
+      Bᵢⱼ,  otherwise
+"""
+function complex_outgoing_matrix(rn)
+    Δ = copy(complex_incidence_matrix(rn))
+    Δ[Δ .== 1] .= 0
+    return Δ
+end
+
+
+################################################################################################
 ######################## conservation laws ###############################
 
 """