utils.jl 10.2 KB
Newer Older
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
1
"""
2
    function stoichiometry_matrix(m::SBML.Model)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
3
4

Extract the vector of species (aka metabolite) identifiers, vector of reaction
5
6
7
identifiers, and a sparse stoichiometry matrix (of type `SparseMatrixCSC` from
`SparseArrays` package) from an existing `SBML.Model`. Returns a 3-tuple with
these values.
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
8
"""
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
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
44
    end
45
    return rows, cols, SparseArrays.sparse(SI, RI, SV, length(rows), length(cols))
46
47
end

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

51
52
53
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
54
"""
55
flux_bounds(m::SBML.Model)::NTuple{2,Vector{Tuple{Float64,String}}} =
56
    (broadcast(x -> x.lb, values(m.reactions)), broadcast(x -> x.ub, values(m.reactions)))
57

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
58
"""
59
    flux_objective(m::SBML.Model)::Vector{Float64}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
60
61
62

Extract the vector of objective coefficients of each reaction.
"""
63
flux_objective(m::SBML.Model)::Vector{Float64} = broadcast(x -> x.oc, values(m.reactions))
64
65
66
67

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

68
69
70
71
72
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.
73
74
75

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

76
77
78
79
80
81
82
83
84
# 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))
```
"""
85
initial_amounts(m::SBML.Model; convert_concentrations = false) = (
86
87
88
89
90
91
92
    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
93
    else
94
95
96
        nothing
    end for (k, s) in m.species
)
97
98
99
100

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

101
102
Return initial concentrations of the species in the model. Refer to work-alike
[`initial_amounts`](@ref) for details.
103
"""
104
initial_concentrations(m::SBML.Model; convert_amounts = false) = (
105
106
107
108
109
110
111
    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
112
    else
113
114
115
        nothing
    end for (k, s) in m.species
)
116
117
118


"""
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
119
120
121
122
123
124
125
126
127
    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)'",
            ),
        ),
    )
128
129
130
131
132

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
133
134
135
136
137
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`.

138
139
Handling of units in the conversion process is ignored in this version.
"""
140
141
142
143
144
145
146
147
148
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)'",
        ),
    ),
)
149
150
151
152
153
    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
154
        isnothing(sz) && (sz = handle_empty_compartment_size(x.id))
155
156
157
158
159
160
161
        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
162
163

"""
164
    get_error_messages(doc::VPtr, error::Exception, report_severities)
165
166
167

Show the error messages reported by SBML in the `doc` document and throw the
`error` if they are more than 1.
168
169
170

`report_severities` switches the reporting of certain error types defined by
libsbml; you can choose from `["Fatal", "Error", "Warning", "Informational"]`.
171
"""
172
function get_error_messages(doc::VPtr, error::Exception, report_severities)
173
    n_errs = ccall(sbml(:SBMLDocument_getNumErrors), Cuint, (VPtr,), doc)
174
    do_throw = false
175
    for i = 1:n_errs
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
176
        err = ccall(sbml(:SBMLDocument_getError), VPtr, (VPtr, Cuint), doc, i - 1)
177
        msg = string(strip(get_string(err, :XMLError_getMessage)))
178
179
180
        sev = string(strip(get_string(err, :XMLError_getSeverityAsString)))
        # keywords from `libsbml/src/sbml/xml/XMLError.cpp` xmlSeverityStringTable:
        if sev == "Fatal"
181
            sev in report_severities && @error "SBML reported fatal error: $(msg)"
182
183
            do_throw = true
        elseif sev == "Error"
184
            sev in report_severities && @error "SBML reported error: $(msg)"
185
186
            do_throw = true
        elseif sev == "Warning"
187
188
189
            sev in report_severities && @warn "SBML reported warning: $(msg)"
        elseif sev == "Informational"
            sev in report_severities && @info "SBML reported: $(msg)"
190
        end
191
    end
192
    do_throw && throw(error)
193
194
195
196
    nothing
end

"""
197
198
199
200
201
202
    check_errors(
        success::Integer,
        doc::VPtr,
        error::Exception,
        report_severities = ["Fatal", "Error"],
    )
203
204
205
206
207
208
209

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.
"""
210
211
212
213
214
215
check_errors(
    success::Integer,
    doc::VPtr,
    error::Exception,
    report_severities = ["Fatal", "Error"],
) = Bool(success) || get_error_messages(doc, error, report_severities)
216
217

# NOTE: this mapping is valid for Level 3/Version 2, it *may* not be valid for
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
# 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)
258
259
260
)

# Get a `Unitful` quantity out of a `Unit_t`.
Mosè Giordano's avatar
Mosè Giordano committed
261
262
263
function get_unit(u::VPtr)
    multiplier = ccall(sbml(:Unit_getMultiplier), Cdouble, (VPtr,), u)
    unit = UNIT_KIND_STRINGS_TO_UNIT[unsafe_string(
264
265
266
267
        ccall(
            sbml(:UnitKind_toString),
            Cstring,
            (Cint,),
268
269
            ccall(sbml(:Unit_getKind), Cint, (VPtr,), u),
        ),
Mosè Giordano's avatar
Mosè Giordano committed
270
271
272
273
274
    )]
    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
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
275
    return (multiplier * unit * exp10(scale))^exponent
Mosè Giordano's avatar
Mosè Giordano committed
276
end
277
278

# Get `Unitful` quantity out of a `UnitDefinition_t`.
279
280
281
282
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)
)