format

Author: vyudu

src/network_analysis.jl

@@ -199,7 +199,8 @@ end
     Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap. 
 """
 function laplacianmat(rn::ReactionSystem, pmap = Dict(); sparse = false) 
-    D = incidencemat(rn; sparse); K = fluxmat(rn, pmap; sparse)
+    D = incidencemat(rn; sparse)
+    K = fluxmat(rn, pmap; sparse)
     D*K
 end
 
@@ -213,7 +214,9 @@ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false)
     rates = reactionrates(rn)
 
     !isempty(pmap) && (rates = substitutevals(rn, pmap, parameters(rn), rates))
-    rcmap = reactioncomplexmap(rn); nc = length(rcmap); nr = length(rates)
+    rcmap = reactioncomplexmap(rn)
+    nc = length(rcmap)
+    nr = length(rates)
     mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates)
     if sparse
         return fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates)
@@ -239,7 +242,8 @@ function fluxmat(::Type{SparseMatrixCSC{T, Int}}, rcmap, rates) where T
 end
 
 function fluxmat(::Type{Matrix{T}}, rcmap, rates) where T
-    nr = length(rates); nc = length(rcmap)
+    nr = length(rates)
+    nc = length(rcmap)
     K = zeros(T, nr, nc)
     for (i, (complex, rxs)) in enumerate(rcmap)
         for (rx, dir) in rxs
@@ -260,7 +264,7 @@ function fluxmat(rn::ReactionSystem, pmap::Tuple; sparse = false)
 end
 
 # Helper to substitute values into a (vector of) symbolic expressions. The syms are the symbols to substitute and the symexprs are the expressions to substitute into.
-function substitutevals(rn::ReactionSystem, map::Dict, syms, symexprs) 
+function substitutevals(rn::ReactionSystem, map::Dict, syms, symexprs)
     length(map) != length(syms) && error("Incorrect number of parameter-value pairs were specified.")
     map = symmap_to_varmap(rn, map)
     map = Dict(ModelingToolkit.value(k) => v for (k, v) in map)
@@ -272,11 +276,13 @@ end
 
     Return the vector whose entries correspond to the "mass action products" of each complex. For example, given the complex A + B, the corresponding entry of the vector would be ``A*B``, and for the complex 2X + Y, the corresponding entry would be ``X^2*Y``. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the [`laplacianmat`](@ref). This utility returns ``Φ(x)``.
     Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap.
-    If the `combinatoric_ratelaws` option is set, will include prefactors for that (see [introduction to Catalyst's rate laws](@ref introduction_to_catalyst_ratelaws).
+    If the `combinatoric_ratelaws` option is set, will include prefactors for that (see [introduction to Catalyst's rate laws](@ref introduction_to_catalyst_ratelaws). Will default to the default for the system.
 """
-function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric_ratelaws = true)
-    r = numreactions(rn); rxs = reactions(rn)
-    sm = speciesmap(rn); specs = species(rn)
+function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric_ratelaws = rn.combinatoric_ratelaws)
+    r = numreactions(rn)
+    rxs = reactions(rn)
+    sm = speciesmap(rn)
+    specs = species(rn)
 
     if !all(r -> ismassaction(r, rn), rxs)
         error("The supplied ReactionSystem has reactions that are not ismassaction. The mass action vector is only defined for pure mass action networks.")
@@ -298,12 +304,12 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric
     Φ 
 end
 
-function massactionvector(rn::ReactionSystem, scmap::Tuple; combinatoric_ratelaws = true) 
+function massactionvector(rn::ReactionSystem, scmap::Tuple; combinatoric_ratelaws = rn.combinatoric_ratelaws) 
     sdict = Dict(scmap)
     massactionvector(rn, sdict; combinatoric_ratelaws)
 end
 
-function massactionvector(rn::ReactionSystem, scmap::Vector; combinatoric_ratelaws = true) 
+function massactionvector(rn::ReactionSystem, scmap::Vector; combinatoric_ratelaws = rn.combinatoric_ratelaws) 
     sdict = Dict(scmap)
     massactionvector(rn, sdict; combinatoric_ratelaws)
 end

test/network_analysis/network_properties.jl

@@ -381,7 +381,8 @@ let
 
     K = fluxmat(MAPK)
     # Construct matrix from incidence matrix
-    mat = zeros(Num, 30, 26); D = incidencemat(MAPK)
+    mat = zeros(Num, 30, 26)
+    D = incidencemat(MAPK)
     rates = reactionrates(MAPK)
     for (i, col) in enumerate(eachcol(D))
         sub = findfirst(==(-1), col)
@@ -403,7 +404,8 @@ let
     A_k = Catalyst.laplacianmat(MAPK; sparse = true)
     @test Catalyst.issparse(A_k)
 
-    S = netstoichmat(MAPK); Y = complexstoichmat(MAPK)
+    S = netstoichmat(MAPK)
+    Y = complexstoichmat(MAPK)
     @test isequal(S*K, Y*A_k)
 
     eqs = Catalyst.assemble_oderhs(MAPK, specs)
@@ -440,7 +442,8 @@ let
         (k2, k3), 2Y + 2Z <--> 3X
     end
 
-    Φ = Catalyst.massactionvector(rn); specs = species(rn)
+    Φ = Catalyst.massactionvector(rn)
+    specs = species(rn)
     crvec = [rn.X^2/2 * rn.Y * rn.Z^3/6,
              1.,
              rn.Y^2/2 * rn.Z^2/2,
@@ -455,7 +458,10 @@ let
 
     # Test that the ODEs generated are the same.
     eqs = Catalyst.assemble_oderhs(rn, specs)
-    S = netstoichmat(rn); Y = complexstoichmat(rn); K = fluxmat(rn); A_k = laplacianmat(rn)
+    S = netstoichmat(rn)
+    Y = complexstoichmat(rn)
+    K = fluxmat(rn)
+    A_k = laplacianmat(rn)
     @test all(iszero, simplify(eqs - S*K*Φ))
     @test all(iszero, simplify(eqs - Y*A_k*Φ))
 
@@ -467,7 +473,8 @@ let
     k = rand(rng, numparams(rn))
     ratevec = collect(zip(parameters(rn), k))
     ratemap = Dict(ratevec)
-    K = fluxmat(rn, ratemap); A_k = laplacianmat(rn, ratemap)
+    K = fluxmat(rn, ratemap)
+    A_k = laplacianmat(rn, ratemap)
 
     numeqs = similar(eqs)
     for i in 1:length(eqs)