utils.jl 5.66 KB
Newer Older
1
2
const VPtr = Ptr{Cvoid}

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
3
"""
anand jain's avatar
anand jain committed
4
    function getS(m::SBML.Model; zeros=spzeros)::Tuple{Vector{String},Vector{String},AbstractMatrix{Float64}}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
5
6

Extract the vector of species (aka metabolite) identifiers, vector of reaction
anand jain's avatar
anand jain committed
7
identifiers, and the (dense) stoichiometry matrix from an existing `SBML.Model`.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
8
Returns a tuple with these values.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
9
10
11
12
13
14

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.
"""
function getS(
anand jain's avatar
anand jain committed
15
    m::SBML.Model;
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
16
17
18
    zeros = spzeros,
)::Tuple{Vector{String},Vector{String},AbstractMatrix{Float64}}
    rows = [k for k in keys(m.species)]
19
    cols = [k for k in keys(m.reactions)]
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
20
    rowsd = Dict(k => i for (i, k) in enumerate(rows))
21
    S = zeros(Float64, length(rows), length(cols))
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
22
23
24
    for col = 1:length(cols)
        stoi = m.reactions[cols[col]].stoichiometry
        S[getindex.(Ref(rowsd), keys(stoi)), col] .= values(stoi)
25
26
27
28
    end
    return rows, cols, S
end

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
29
"""
30
    getLBs(m::SBML.Model)::Vector{Tuple{Float64,String}}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
31
32
33
34
35

Extract a vector of lower bounds of reaction rates from the model. All bounds
are accompanied with the unit of the corresponding value (this behavior is
based on SBML specification).
"""
36
37
getLBs(m::SBML.Model)::Vector{Tuple{Float64,String}} =
    broadcast(x -> x.lb, values(m.reactions))
38

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
39
"""
40
    getUBs(m::SBML.Model)::Vector{Tuple{Float64,String}}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
41

42
Likewise to [`getLBs`](@ref), extract the upper bounds.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
43
"""
44
45
getUBs(m::SBML.Model)::Vector{Tuple{Float64,String}} =
    broadcast(x -> x.ub, values(m.reactions))
46

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
47
"""
48
    getOCs(m::SBML.Model)::Vector{Float64}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
49
50
51

Extract the vector of objective coefficients of each reaction.
"""
52
getOCs(m::SBML.Model)::Vector{Float64} = broadcast(x -> x.oc, values(m.reactions))
53
54
55
56

"""
    initial_amounts(m::SBML.Model; convert_concentrations = false)

57
58
59
60
61
Return initial amounts for each species as a generator of pairs
`species_name => initial_amount`; the amount is set to `nothing` if not
available. If `convert_concentrations` is true and there is information about
initial concentration available together with compartment size, the result is
computed from the species' initial concentration.
62
63
64

In the current version, units of the measurements are completely ignored.

65
66
67
68
69
70
71
72
73
# Example
```
# get the initial amounts as dictionary
Dict(initial_amounts(model, convert_concentrations = true))

# remove the empty entries
Dict(k => v for (k,v) in initial_amounts(model) if !isnothing(v))
```
"""
74
initial_amounts(m::SBML.Model; convert_concentrations = false) = (
75
76
77
78
79
80
81
    k => if !isnothing(s.initial_amount)
        s.initial_amount[1]
    elseif convert_concentrations &&
           !isnothing(s.initial_concentration) &&
           haskey(m.compartments, s.compartment) &&
           !isnothing(m.compartments[s.compartment].size)
        s.initial_concentration[1] * m.compartments[s.compartment].size
82
    else
83
84
85
        nothing
    end for (k, s) in m.species
)
86
87
88
89

"""
    initial_concentrations(m::SBML.Model; convert_amounts = false)

90
91
Return initial concentrations of the species in the model. Refer to work-alike
[`initial_amounts`](@ref) for details.
92
"""
93
initial_concentrations(m::SBML.Model; convert_amounts = false) = (
94
95
96
97
98
99
100
    k => if !isnothing(s.initial_concentration)
        s.initial_concentration[1]
    elseif convert_amounts &&
           !isnothing(s.initial_amount) &&
           haskey(m.compartments, s.compartment) &&
           !isnothing(m.compartments[s.compartment].size)
        s.initial_amount[1] / m.compartments[s.compartment].size
101
    else
102
103
104
        nothing
    end for (k, s) in m.species
)
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134


"""
    extensive_kinetic_math(m::SBML.Model, formula::SBML.Math)

Convert a SBML math `formula` to "extensive" kinetic laws, where the references
to species that are marked as not having only substance units are converted
from amounts to concentrations.

Handling of units in the conversion process is ignored in this version.
"""
function extensive_kinetic_math(m::SBML.Model, formula::SBML.Math)
    conv(x::SBML.MathIdent) = begin
        haskey(m.species, x.id) || return x
        sp = m.species[x.id]
        sp.only_substance_units && return x
        sz = m.compartments[sp.compartment].size
        isnothing(sz) && throw(
            DomainError(
                formula,
                "Non-substance-only-unit reference to species `$(x.id)' in an unsized compartment `$(sp.compartment)'.",
            ),
        )
        SBML.MathApply("/", [x, SBML.MathVal(sz)])
    end
    conv(x::SBML.MathApply) = SBML.MathApply(x.fn, conv.(x.args))
    conv(x::SBML.Math) = x

    conv(formula)
end
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165

"""
    get_error_messages(doc::Ptr{Cvoid}, error::Exception)

Show the error messages reported by SBML in the `doc` document and throw the
`error` if they are more than 1.
"""
function get_error_messages(doc::VPtr, error::Exception)
    n_errs = ccall(sbml(:SBMLDocument_getNumErrors), Cuint, (VPtr,), doc)
    for i = 1:n_errs
        err = ccall(sbml(:SBMLDocument_getError), VPtr, (VPtr, Cuint), doc, i-1)
        msg = string(strip(get_string(err, :XMLError_getMessage)))
        @error "SBML reported error: $(msg)"
    end
    if n_errs > 0
        throw(error)
    end
    nothing
end

"""
    check_errors(success::Integer, doc::Ptr{Cvoid}, error::Exception)

If success is a 0-valued `Integer` (a logical `false`), then call
[`get_error_messages`](@ref) to show the error messages reported by SBML in the
`doc` document and throw the `error` if they are more than 1.  `success` is
typically the value returned by an SBML C function operating on `doc` which
returns a boolean flag to signal a successful operation.
"""
check_errors(success::Integer, doc::VPtr, error::Exception) =
    Bool(success) || get_error_messages(doc, error)