readsbml.jl 4.56 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
13
function readSBML(fn::String)::Model
    doc = ccall(sbml(:readSBML), VPtr, (Cstring,), fn)
    try
        n_errs = ccall(sbml(:SBMLDocument_getNumErrors), Cuint, (VPtr,), doc)
        if n_errs > 0
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
14
            throw(SystemError("Opening SBML document has failed"))
15
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
16

17
        if 0 == ccall(sbml(:SBMLDocument_isSetModel), Cint, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
18
            throw(AssertionError("SBML document contains no model"))
19
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
20

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

23
24
25
        return extractModel(model)
    finally
        ccall(sbml(:SBMLDocument_free), Nothing, (VPtr,), doc)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
26
    end
27
end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
28

29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
function extractModel(mdl::VPtr)::Model
    units = Dict{String,Vector{UnitPart}}()
    for i = 0:ccall(sbml(:Model_getNumUnitDefinitions), Cuint, (VPtr,), mdl)-1
        ud = ccall(sbml(:Model_getUnitDefinition), VPtr, (VPtr, Cuint), mdl, i)
        id = unsafe_string(ccall(sbml(:UnitDefinition_getId), Cstring, (VPtr,), ud))
        units[id] = [
            begin
                u = ccall(sbml(:UnitDefinition_getUnit), VPtr, (VPtr, Cuint), ud, j)
                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),
                )
            end for j = 0:ccall(sbml(:UnitDefinition_getNumUnits), Cuint, (VPtr,), ud)-1
        ]
52
53
    end

54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
    compartments = [
        unsafe_string(
            ccall(
                sbml(:Compartment_getId),
                Cstring,
                (VPtr,),
                ccall(sbml(:Model_getCompartment), VPtr, (VPtr, Cuint), mdl, i),
            ),
        ) for i = 0:ccall(sbml(:Model_getNumCompartments), Cuint, (VPtr,), mdl)-1
    ]

    species = Dict{String,Species}()
    for i = 0:ccall(sbml(:Model_getNumSpecies), Cuint, (VPtr,), mdl)-1
        sp = ccall(sbml(:Model_getSpecies), VPtr, (VPtr, Cuint), mdl, i)
        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)),
        )
72
73
    end

74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
    reactions = Dict{String,Reaction}()
    for i = 0:ccall(sbml(:Model_getNumReactions), Cuint, (VPtr,), mdl)-1
        re = ccall(sbml(:Model_getReaction), VPtr, (VPtr, Cuint), mdl, i)
        kl = ccall(sbml(:Reaction_getKineticLaw), VPtr, (VPtr,), re)
        lb = (-Inf, "")
        ub = (Inf, "")
        oc = 0.0
        for j = 0:ccall(sbml(:KineticLaw_getNumParameters), Cuint, (VPtr,), kl)-1
            p = ccall(sbml(:KineticLaw_getParameter), VPtr, (VPtr, Cuint), kl, j)
            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
94
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
95

96
97
98
99
100
101
102
103
        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
104

105
106
107
108
        for j = 0:ccall(sbml(:Reaction_getNumReactants), Cuint, (VPtr,), re)-1
            sr = ccall(sbml(:Reaction_getReactant), VPtr, (VPtr, Cuint), re, j)
            add_stoi(sr, -1)
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
109

110
111
112
113
        for j = 0:ccall(sbml(:Reaction_getNumProducts), Cuint, (VPtr,), re)-1
            sr = ccall(sbml(:Reaction_getProduct), VPtr, (VPtr, Cuint), re, j)
            add_stoi(sr, 1)
        end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
114

115
116
117
        reactions[unsafe_string(ccall(sbml(:Reaction_getId), Cstring, (VPtr,), re))] =
            Reaction(stoi, lb, ub, oc)
    end
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
118

119
    return Model(units, compartments, species, reactions)
Miroslav Kratochvil's avatar
Miroslav Kratochvil committed
120
end