readsbml.jl 7.34 KB
Newer Older
1

2
const VPtr = Ptr{Cvoid}
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
3

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
4
5
6
7
8
"""
    function readSBML(fn::String)::Model

Read the SBML from a XML file in `fn` and return the contained `Model`.
"""
9
10
11
12
function readSBML(fn::String)::Model
    doc = ccall(sbml(:readSBML), VPtr, (Cstring,), fn)
    try
        n_errs = ccall(sbml(:SBMLDocument_getNumErrors), Cuint, (VPtr,), doc)
13
14
15
16
17
        for i = 0:n_errs-1
            err = ccall(sbml(:SBMLDocument_getError), VPtr, (VPtr, Cuint), doc, i)
            msg = unsafe_string(ccall(sbml(:XMLError_getMessage), Cstring, (VPtr,), err))
            @warn "SBML reported error: $msg"
        end
18
        if n_errs > 0
19
            throw(AssertionError("Opening SBML document has reported errors"))
20
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
21

22
        if 0 == ccall(sbml(:SBMLDocument_isSetModel), Cint, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
23
            throw(AssertionError("SBML document contains no model"))
24
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
25

26
        model = ccall(sbml(:SBMLDocument_getModel), VPtr, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
27

28
29
30
        return extractModel(model)
    finally
        ccall(sbml(:SBMLDocument_free), Nothing, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
31
    end
32
end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
33

34
function extractModel(mdl::VPtr)::Model
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
35
36
    mdl_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), mdl, "fbc")

St. Elmo's avatar
St. Elmo committed
37
38
39
40
41
42
43
44
    parameters = Dict{String,Float64}()
    for i = 1:ccall(sbml(:Model_getNumParameters), Cuint, (VPtr,), mdl)
        p = ccall(sbml(:Model_getParameter), VPtr, (VPtr, Cuint), mdl, i - 1)
        id = unsafe_string(ccall(sbml(:Parameter_getId), Cstring, (VPtr,), p))
        v = ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p)
        parameters[id] = v
    end

45
    units = Dict{String,Vector{UnitPart}}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
46
47
    for i = 1:ccall(sbml(:Model_getNumUnitDefinitions), Cuint, (VPtr,), mdl)
        ud = ccall(sbml(:Model_getUnitDefinition), VPtr, (VPtr, Cuint), mdl, i - 1)
48
49
50
        id = unsafe_string(ccall(sbml(:UnitDefinition_getId), Cstring, (VPtr,), ud))
        units[id] = [
            begin
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
51
                u = ccall(sbml(:UnitDefinition_getUnit), VPtr, (VPtr, Cuint), ud, j - 1)
52
53
54
55
56
57
58
59
60
61
62
63
64
                UnitPart(
                    unsafe_string(
                        ccall(
                            sbml(:UnitKind_toString),
                            Cstring,
                            (Cint,),
                            ccall(sbml(:Unit_getKind), Cint, (VPtr,), u),
                        ),
                    ),
                    ccall(sbml(:Unit_getExponent), Cint, (VPtr,), u),
                    ccall(sbml(:Unit_getScale), Cint, (VPtr,), u),
                    ccall(sbml(:Unit_getMultiplier), Cdouble, (VPtr,), u),
                )
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
65
            end for j = 1:ccall(sbml(:UnitDefinition_getNumUnits), Cuint, (VPtr,), ud)
66
        ]
67
68
    end

69
70
71
72
73
74
    compartments = [
        unsafe_string(
            ccall(
                sbml(:Compartment_getId),
                Cstring,
                (VPtr,),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
75
                ccall(sbml(:Model_getCompartment), VPtr, (VPtr, Cuint), mdl, i - 1),
76
            ),
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
77
        ) for i = 1:ccall(sbml(:Model_getNumCompartments), Cuint, (VPtr,), mdl)
78
79
80
    ]

    species = Dict{String,Species}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
81
82
    for i = 1:ccall(sbml(:Model_getNumSpecies), Cuint, (VPtr,), mdl)
        sp = ccall(sbml(:Model_getSpecies), VPtr, (VPtr, Cuint), mdl, i - 1)
83
84
85
86
        species[unsafe_string(ccall(sbml(:Species_getId), Cstring, (VPtr,), sp))] = Species(
            unsafe_string(ccall(sbml(:Species_getName), Cstring, (VPtr,), sp)),
            unsafe_string(ccall(sbml(:Species_getCompartment), Cstring, (VPtr,), sp)),
        )
87
88
    end

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
    objectives_fbc = Dict{String,Float64}()
    if mdl_fbc != C_NULL
        for i = 1:ccall(sbml(:FbcModelPlugin_getNumObjectives), Cuint, (VPtr,), mdl_fbc)
            o = ccall(
                sbml(:FbcModelPlugin_getObjective),
                VPtr,
                (VPtr, Cuint),
                mdl_fbc,
                i - 1,
            )
            # this part seems missing in C API docs...
            for j = 1:ccall(sbml(:Objective_getNumFluxObjectives), Cuint, (VPtr,), o)
                fo = ccall(sbml(:Objective_getFluxObjective), VPtr, (VPtr, Cuint), o, j - 1)
                objectives_fbc[unsafe_string(
                    ccall(sbml(:FluxObjective_getReaction), Cstring, (VPtr,), fo),
                )] = ccall(sbml(:FluxObjective_getCoefficient), Cdouble, (VPtr,), fo)
            end
        end
    end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
108

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
109
    reactions = Dict{String,Reaction}()
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
110
111
    for i = 1:ccall(sbml(:Model_getNumReactions), Cuint, (VPtr,), mdl)
        re = ccall(sbml(:Model_getReaction), VPtr, (VPtr, Cuint), mdl, i - 1)
112
113
114
        lb = (-Inf, "")
        ub = (Inf, "")
        oc = 0.0
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

        kl = ccall(sbml(:Reaction_getKineticLaw), VPtr, (VPtr,), re)
        if kl != C_NULL
            for j = 1:ccall(sbml(:KineticLaw_getNumParameters), Cuint, (VPtr,), kl)
                p = ccall(sbml(:KineticLaw_getParameter), VPtr, (VPtr, Cuint), kl, j - 1)
                id = unsafe_string(ccall(sbml(:Parameter_getId), Cstring, (VPtr,), p))
                pval = () -> ccall(sbml(:Parameter_getValue), Cdouble, (VPtr,), p)
                punit =
                    () ->
                        unsafe_string(ccall(sbml(:Parameter_getUnits), Cstring, (VPtr,), p))
                if id == "LOWER_BOUND"
                    lb = (pval(), punit())
                elseif id == "UPPER_BOUND"
                    ub = (pval(), punit())
                elseif id == "OBJECTIVE_COEFFICIENT"
                    oc = pval()
                end
132
            end
133
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
134

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
135
        # TRICKY: SBML spec is completely silent about what should happen if
St. Elmo's avatar
St. Elmo committed
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
        # someone specifies both the above and below formats of the flux bounds
        # for one reaction. Notably, these do not really specify much
        # interaction with units. In this case, we'll just set a special
        # "[fbc]" unit that has no specification in `units`, and hope the users
        # can make something out of it.
        re_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), re, "fbc")
        if re_fbc != C_NULL
            fbcb =
                ccall(sbml(:FbcReactionPlugin_getLowerFluxBound), Cstring, (VPtr,), re_fbc)
            if fbcb != C_NULL && haskey(parameters, unsafe_string(fbcb))
                lb = (parameters[unsafe_string(fbcb)], "[fbc]")
            end
            fbcb =
                ccall(sbml(:FbcReactionPlugin_getUpperFluxBound), Cstring, (VPtr,), re_fbc)
            if fbcb != C_NULL && haskey(parameters, unsafe_string(fbcb))
                ub = (parameters[unsafe_string(fbcb)], "[fbc]")
            end
        end

155
156
157
158
159
160
161
162
        stoi = Dict{String,Float64}()
        add_stoi =
            (sr, factor) ->
                stoi[unsafe_string(
                    ccall(sbml(:SpeciesReference_getSpecies), Cstring, (VPtr,), sr),
                )] =
                    ccall(sbml(:SpeciesReference_getStoichiometry), Cdouble, (VPtr,), sr) *
                    factor
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
163

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
164
165
        for j = 1:ccall(sbml(:Reaction_getNumReactants), Cuint, (VPtr,), re)
            sr = ccall(sbml(:Reaction_getReactant), VPtr, (VPtr, Cuint), re, j - 1)
166
167
            add_stoi(sr, -1)
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
168

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
169
170
        for j = 1:ccall(sbml(:Reaction_getNumProducts), Cuint, (VPtr,), re)
            sr = ccall(sbml(:Reaction_getProduct), VPtr, (VPtr, Cuint), re, j - 1)
171
172
            add_stoi(sr, 1)
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
173

Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
174
175
176
        reid = unsafe_string(ccall(sbml(:Reaction_getId), Cstring, (VPtr,), re))
        reactions[reid] =
            Reaction(stoi, lb, ub, haskey(objectives_fbc, reid) ? objectives_fbc[reid] : oc)
177
    end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
178

St. Elmo's avatar
St. Elmo committed
179
    return Model(parameters, units, compartments, species, reactions)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
180
end