Unverified Commit 2043139d authored by Miroslav Kratochvil's avatar Miroslav Kratochvil Committed by GitHub
Browse files

Merge pull request #147 from LCSB-BioCore/mk-sparse-stoi

simplify stoichiometry_matrix API & make it more efficient
parents 1bd9216f 3c2c3afe
"""
function stoichiometry_matrix(m::SBML.Model; zeros=spzeros)::Tuple{Vector{String},Vector{String},AbstractMatrix{Float64}}
function stoichiometry_matrix(m::SBML.Model)
Extract the vector of species (aka metabolite) identifiers, vector of reaction
identifiers, and the (dense) stoichiometry matrix from an existing `SBML.Model`.
Returns a tuple with these values.
The matrix is sparse by default (initially constructed by
`SparseArrays.spzeros`). You can fill in a custom empty matrix constructed to
argument `zeros`; e.g. running with `zeros=zeros` will produce a dense matrix.
identifiers, and a sparse stoichiometry matrix (of type `SparseMatrixCSC` from
`SparseArrays` package) from an existing `SBML.Model`. Returns a 3-tuple with
these values.
"""
function stoichiometry_matrix(
m::SBML.Model;
zeros = spzeros,
)::Tuple{Vector{String},Vector{String},AbstractMatrix{Float64}}
rows = [k for k in keys(m.species)]
cols = [k for k in keys(m.reactions)]
rowsd = Dict(k => i for (i, k) in enumerate(rows))
S = zeros(Float64, length(rows), length(cols))
for col = 1:length(cols)
s = m.reactions[cols[col]].reactants
S[getindex.(Ref(rowsd), keys(s)), col] .-= values(s)
s = m.reactions[cols[col]].products
S[getindex.(Ref(rowsd), keys(s)), col] .+= values(s)
function stoichiometry_matrix(m::SBML.Model)
rows = collect(keys(m.species))
cols = collect(keys(m.reactions))
row_idx = Dict(k => i for (i, k) in enumerate(rows))
col_idx = Dict(k => i for (i, k) in enumerate(cols))
nnz = 0
for (_, r) in m.reactions
for _ in r.reactants
nnz += 1
end
for _ in r.products
nnz += 1
end
end
SI = Int[]
RI = Int[]
SV = Float64[]
sizehint!(SI, nnz)
sizehint!(RI, nnz)
sizehint!(SV, nnz)
for (rid, r) in m.reactions
ridx = col_idx[rid]
for (sid, stoi) in r.reactants
push!(SI, row_idx[sid])
push!(RI, ridx)
push!(SV, -stoi)
end
for (sid, stoi) in r.products
push!(SI, row_idx[sid])
push!(RI, ridx)
push!(SV, stoi)
end
end
return rows, cols, S
return rows, cols, SparseArrays.sparse(SI, RI, SV, length(rows), length(cols))
end
"""
......@@ -253,7 +272,7 @@ function get_unit(u::VPtr)
exponent = ccall(sbml(:Unit_getExponent), Cint, (VPtr,), u)
# See page 44 of
# http://sbml.org/Special/specifications/sbml-level-3/version-2/core/release-2/sbml-level-3-version-2-release-2-core.pdf
return (multiplier * unit * exp10(scale)) ^ exponent
return (multiplier * unit * exp10(scale))^exponent
end
# Get `Unitful` quantity out of a `UnitDefinition_t`.
......
......@@ -28,9 +28,7 @@ end
mets, rxns, S = stoichiometry_matrix(mdl)
@test typeof(S) <: AbstractMatrix{Float64}
@test typeof(stoichiometry_matrix(mdl; zeros = spzeros)[3]) <: SparseMatrixCSC{Float64}
@test typeof(stoichiometry_matrix(mdl; zeros = zeros)[3]) == Matrix{Float64}
@test typeof(S) <: SparseMatrixCSC{Float64}
@test length(mets) == 77
@test length(rxns) == 77
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment