utils.jl 10.1 KB
Newer Older
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
1
"""
2
    function stoichiometry_matrix(m::SBML.Model; zeros=spzeros)::Tuple{Vector{String},Vector{String},AbstractMatrix{Float64}}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
3
4

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

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

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
29
"""
30
    flux_bounds(m::SBML.Model)::NTuple{2, Vector{Tuple{Float64,String}}}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
31

32
33
34
Extract the vectors of lower and upper 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).
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
35
"""
36
flux_bounds(m::SBML.Model)::NTuple{2,Vector{Tuple{Float64,String}}} =
37
    (broadcast(x -> x.lb, values(m.reactions)), broadcast(x -> x.ub, values(m.reactions)))
38

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
39
"""
40
    flux_objective(m::SBML.Model)::Vector{Float64}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
41
42
43

Extract the vector of objective coefficients of each reaction.
"""
44
flux_objective(m::SBML.Model)::Vector{Float64} = broadcast(x -> x.oc, values(m.reactions))
45
46
47
48

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

49
50
51
52
53
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.
54
55
56

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

57
58
59
60
61
62
63
64
65
# 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))
```
"""
66
initial_amounts(m::SBML.Model; convert_concentrations = false) = (
67
68
69
70
71
72
73
    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
74
    else
75
76
77
        nothing
    end for (k, s) in m.species
)
78
79
80
81

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

82
83
Return initial concentrations of the species in the model. Refer to work-alike
[`initial_amounts`](@ref) for details.
84
"""
85
initial_concentrations(m::SBML.Model; convert_amounts = false) = (
86
87
88
89
90
91
92
    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
93
    else
94
95
96
        nothing
    end for (k, s) in m.species
)
97
98
99


"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
100
101
102
103
104
105
106
107
108
    function extensive_kinetic_math(
        m::SBML.Model,
        formula::SBML.Math;
        handle_empty_compartment_size = (id::String) -> throw(
            DomainError(
                "Non-substance-only-unit reference to species `\$id' in an unsized compartment `\$(m.species[id].compartment)'",
            ),
        ),
    )
109
110
111
112
113

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.

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
114
115
116
117
118
If the data is missing, you can supply a function that adds them. A common way
to handle errors is to assume that unsized compartments have volume 1.0 (of
whatever units), you can specify that behavior by supplying
`handle_empty_compartment_size = _ -> 1.0`.

119
120
Handling of units in the conversion process is ignored in this version.
"""
121
122
123
124
125
126
127
128
129
function extensive_kinetic_math(
    m::SBML.Model,
    formula::SBML.Math;
    handle_empty_compartment_size = (id::String) -> throw(
        DomainError(
            "Non-substance-only-unit reference to species `$id' in an unsized compartment `$(m.species[id].compartment)'",
        ),
    ),
)
130
131
132
133
134
    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
135
        isnothing(sz) && (sz = handle_empty_compartment_size(x.id))
136
137
138
139
140
141
142
        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
143
144

"""
145
    get_error_messages(doc::VPtr, error::Exception, report_severities)
146
147
148

Show the error messages reported by SBML in the `doc` document and throw the
`error` if they are more than 1.
149
150
151

`report_severities` switches the reporting of certain error types defined by
libsbml; you can choose from `["Fatal", "Error", "Warning", "Informational"]`.
152
"""
153
function get_error_messages(doc::VPtr, error::Exception, report_severities)
154
    n_errs = ccall(sbml(:SBMLDocument_getNumErrors), Cuint, (VPtr,), doc)
155
    do_throw = false
156
    for i = 1:n_errs
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
157
        err = ccall(sbml(:SBMLDocument_getError), VPtr, (VPtr, Cuint), doc, i - 1)
158
        msg = string(strip(get_string(err, :XMLError_getMessage)))
159
160
161
        sev = string(strip(get_string(err, :XMLError_getSeverityAsString)))
        # keywords from `libsbml/src/sbml/xml/XMLError.cpp` xmlSeverityStringTable:
        if sev == "Fatal"
162
            sev in report_severities && @error "SBML reported fatal error: $(msg)"
163
164
            do_throw = true
        elseif sev == "Error"
165
            sev in report_severities && @error "SBML reported error: $(msg)"
166
167
            do_throw = true
        elseif sev == "Warning"
168
169
170
            sev in report_severities && @warn "SBML reported warning: $(msg)"
        elseif sev == "Informational"
            sev in report_severities && @info "SBML reported: $(msg)"
171
        end
172
    end
173
    do_throw && throw(error)
174
175
176
177
    nothing
end

"""
178
179
180
181
182
183
    check_errors(
        success::Integer,
        doc::VPtr,
        error::Exception,
        report_severities = ["Fatal", "Error"],
    )
184
185
186
187
188
189
190

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.
"""
191
192
193
194
195
196
check_errors(
    success::Integer,
    doc::VPtr,
    error::Exception,
    report_severities = ["Fatal", "Error"],
) = Bool(success) || get_error_messages(doc, error, report_severities)
197
198

# NOTE: this mapping is valid for Level 3/Version 2, it *may* not be valid for
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
# other versions.  See
# https://github.com/sbmlteam/libsbml/blob/d4bc12abc4e72e451a0a0f2be4b0b6101ac94160/src/sbml/UnitKind.c#L46-L85
const UNIT_KIND_STRINGS_TO_UNIT = Dict(
    "ampere" => 1.0 * u"A", # UNIT_KIND_AMPERE
    "avogadro" => ustrip(u"mol^-1", Unitful.Na), # UNIT_KIND_AVOGADRO
    "becquerel" => 1.0 * u"Bq", # UNIT_KIND_BECQUEREL
    "candela" => 1.0 * u"cd", # UNIT_KIND_CANDELA
    "Celsius" => 1.0 * u"°C", # UNIT_KIND_CELSIUS
    "coulomb" => 1.0 * u"C", # UNIT_KIND_COULOMB
    "dimensionless" => 1, # UNIT_KIND_DIMENSIONLESS
    "farad" => 1.0 * u"F", # UNIT_KIND_FARAD
    "gram" => 1.0 * u"g", # UNIT_KIND_GRAM
    "gray" => 1.0 * u"Gy", # UNIT_KIND_GRAY
    "henry" => 1.0 * u"H", # UNIT_KIND_HENRY
    "hertz" => 1.0 * u"Hz", # UNIT_KIND_HERTZ
    "item" => 1, # UNIT_KIND_ITEM
    "joule" => 1.0 * u"J", # UNIT_KIND_JOULE
    "katal" => 1.0 * u"kat", # UNIT_KIND_KATAL
    "kelvin" => 1.0 * u"K", # UNIT_KIND_KELVIN
    "kilogram" => 1.0 * u"kg", # UNIT_KIND_KILOGRAM
    "liter" => 1.0 * u"L", # UNIT_KIND_LITER
    "litre" => 1.0 * u"L", # UNIT_KIND_LITRE
    "lumen" => 1.0 * u"lm", # UNIT_KIND_LUMEN
    "lux" => 1.0 * u"lx", # UNIT_KIND_LUX
    "meter" => 1.0 * u"m", # UNIT_KIND_METER
    "metre" => 1.0 * u"m", # UNIT_KIND_METRE
    "mole" => 1.0 * u"mol", # UNIT_KIND_MOLE
    "newton" => 1.0 * u"N", # UNIT_KIND_NEWTON
    "ohm" => 1.0 * u"Ω", # UNIT_KIND_OHM
    "pascal" => 1.0 * u"Pa", # UNIT_KIND_PASCAL
    "radian" => 1.0 * u"rad", # UNIT_KIND_RADIAN
    "second" => 1.0 * u"s", # UNIT_KIND_SECOND
    "siemens" => 1.0 * u"S", # UNIT_KIND_SIEMENS
    "sievert" => 1.0 * u"Sv", # UNIT_KIND_SIEVERT
    "steradian" => 1.0 * u"sr", # UNIT_KIND_STERADIAN
    "tesla" => 1.0 * u"T", # UNIT_KIND_TESLA
    "volt" => 1.0 * u"V", # UNIT_KIND_VOLT
    "watt" => 1.0 * u"W", # UNIT_KIND_WATT
    "weber" => 1.0 * u"W", # UNIT_KIND_WEBER
    "(Invalid UnitKind)" => 1, # UNIT_KIND_INVALID (let's treat is as a dimensionless quantity)
239
240
241
)

# Get a `Unitful` quantity out of a `Unit_t`.
Mosè Giordano's avatar
Mosè Giordano committed
242
243
244
function get_unit(u::VPtr)
    multiplier = ccall(sbml(:Unit_getMultiplier), Cdouble, (VPtr,), u)
    unit = UNIT_KIND_STRINGS_TO_UNIT[unsafe_string(
245
246
247
248
        ccall(
            sbml(:UnitKind_toString),
            Cstring,
            (Cint,),
249
250
            ccall(sbml(:Unit_getKind), Cint, (VPtr,), u),
        ),
Mosè Giordano's avatar
Mosè Giordano committed
251
252
253
254
255
256
257
    )]
    scale = ccall(sbml(:Unit_getScale), Cint, (VPtr,), u)
    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
end
258
259

# Get `Unitful` quantity out of a `UnitDefinition_t`.
260
261
262
263
get_units(ud::VPtr) = prod(
    get_unit(ccall(sbml(:UnitDefinition_getUnit), VPtr, (VPtr, Cuint), ud, j - 1)) for
    j = 1:ccall(sbml(:UnitDefinition_getNumUnits), Cuint, (VPtr,), ud)
)