Add writeSBML

Project.toml

@@ -12,17 +12,20 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
+ConstructionBase = "1.3"
 DocStringExtensions = "0.8, 0.9"
 IfElse = "0.1"
+SBML_jll = "5.19.5"
 Symbolics = "3, 4"
 Unitful = "1"
 julia = "1.6"
 
 [extras]
+ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Downloads", "SHA", "Symbolics", "Test"]
+test = ["ConstructionBase", "Downloads", "SHA", "Symbolics", "Test"]

src/SBML.jl

@@ -18,6 +18,7 @@ include("converters.jl")
 include("interpret.jl")
 include("math.jl")
 include("readsbml.jl")
+include("writesbml.jl")
 include("unitful.jl")
 include("utils.jl")
 
@@ -29,6 +30,10 @@ A shortcut that loads a function symbol from `SBML_jll`.
 sbml(sym::Symbol)::VPtr = dlsym(SBML_jll.libsbml_handle, sym)
 
 export readSBML, readSBMLFromString, stoichiometry_matrix, flux_bounds, flux_objective
+export writeSBML
 export set_level_and_version, libsbml_convert, convert_simplify_math
 
+# Read a file at precompile time, to improve time-to-first `readSBML`.
+writeSBML(readSBML(joinpath(@__DIR__, "..", "test", "data", "Dasgupta2020-written.xml")))
+
 end # module

src/math.jl

@@ -28,13 +28,73 @@ const relational_opers = Dict{Int32,String}(
     312 => "lt",
     313 => "neq",
 )
+# Inverse mapping, needed for creating `ASTNode_t` pointers from `MathApply` objects.
+const inv_relational_opers = Dict(val => key for (key, val) in relational_opers)
 
 function relational_oper(t::Int)
     haskey(relational_opers, t) ||
         throw(DomainError(t, "Unknown ASTNodeType value for relational operator"))
     relational_opers[t]
 end
 
+# Mapping of AST node type value subset to mathematical operations. Depends on
+# `ASTNodeType.h` (also see below the case with AST_NAME_TIME)
+const math_opers = Dict{Int32,String}(43 => "+", 45 => "-", 42 => "*", 47 => "/", 94 => "^")
+# Inverse mapping, needed for creating `ASTNode_t` pointers from `MathApply` objects.
+const inv_math_opers = Dict(val => key for (key, val) in math_opers)
+
+# Mapping of AST node type value subset to logical operations. Depends on
+# `ASTNodeType.h` (also see below the case with AST_NAME_TIME)
+const logical_opers =
+    Dict{Int32,String}(304 => "and", 305 => "not", 306 => "or", 307 => "xor")
+# Inverse mapping, needed for creating `ASTNode_t` pointers from `MathApply` objects.
+const inv_logical_opers = Dict(val => key for (key, val) in logical_opers)
+
+# Mapping of AST node type value subset to mathematical functions. Depends on
+# `ASTNodeType.h` (also see below the case with AST_NAME_TIME)
+const math_funcs = Dict{Int32,String}(
+    269 => "abs",
+    270 => "arccos",
+    271 => "arccosh",
+    272 => "arccot",
+    273 => "arccoth",
+    274 => "arccsc",
+    275 => "arccsch",
+    276 => "arcsec",
+    277 => "arcsech",
+    278 => "arcsin",
+    279 => "arcsinh",
+    280 => "arctan",
+    281 => "arctanh",
+    282 => "ceiling",
+    283 => "cos",
+    284 => "cosh",
+    285 => "cot",
+    286 => "coth",
+    287 => "csc",
+    288 => "csch",
+    289 => "delay",
+    290 => "exp",
+    291 => "factorial",
+    292 => "floor",
+    293 => "ln",
+    294 => "log",
+    295 => "piecewise",
+    296 => "power",
+    297 => "root",
+    298 => "sec",
+    299 => "sech",
+    300 => "sin",
+    301 => "sinh",
+    302 => "tan",
+    303 => "tanh",
+)
+# Inverse mapping, needed for creating `ASTNode_t` pointers from `MathApply` objects.
+const inv_math_funcs = Dict(val => key for (key, val) in math_funcs)
+# Everybody wants to be a map!
+const all_inv_function_mappings =
+    merge(inv_relational_opers, inv_math_opers, inv_logical_opers, inv_math_funcs)
+
 """
 $(TYPEDSIGNATURES)
 
@@ -92,3 +152,95 @@ function parse_math(ast::VPtr)::Math
         return MathIdent("?unsupported?")
     end
 end
+
+## Inverse of `parse_math`: create `ASTNode_t` pointers from `Math` objects.
+
+function get_astnode_ptr(m::MathTime)::VPtr
+    astnode = ccall(sbml(:ASTNode_create), VPtr, ())
+    ccall(sbml(:ASTNode_setName), Cint, (VPtr, Cstring), astnode, m.id)
+    # Same comment as in `parse_math`: this constant must be kept in-sync with
+    # the value of `AST_NAME_TIME` in `libsbml/src/sbml/math/ASTNodeType.h`.
+    ccall(sbml(:ASTNode_setType), Cint, (VPtr, Cuint), astnode, 262)
+    ccall(sbml(:ASTNode_canonicalize), Cint, (VPtr,), astnode)
+    return astnode
+end
+
+function get_astnode_ptr(m::Union{MathIdent,MathConst})::VPtr
+    m.id in ("?invalid?", "?unsupported?") &&
+        error("Cannot get a pointer for `MathIdent` with ID \"$(m.id)\"")
+    astnode = ccall(sbml(:ASTNode_create), VPtr, ())
+    ccall(sbml(:ASTNode_setName), Cint, (VPtr, Cstring), astnode, m.id)
+    ccall(sbml(:ASTNode_canonicalize), Cint, (VPtr,), astnode)
+    return astnode
+end
+
+function get_astnode_ptr(m::MathVal{<:Integer})::VPtr
+    astnode = ccall(sbml(:ASTNode_create), VPtr, ())
+    ccall(sbml(:ASTNode_setInteger), Cint, (VPtr, Clong), astnode, m.val)
+    ccall(sbml(:ASTNode_canonicalize), Cint, (VPtr,), astnode)
+    return astnode
+end
+
+function get_astnode_ptr(m::MathVal{<:Rational})::VPtr
+    astnode = ccall(sbml(:ASTNode_create), VPtr, ())
+    # Note: this can be in principle a lossy reconstruction as `Rational`s in
+    # Julia are automatically simplified (e.g., 5//10 -> 1//2).
+    ccall(
+        sbml(:ASTNode_setRational),
+        Cint,
+        (VPtr, Clong, Clong),
+        astnode,
+        numerator(m.val),
+        denominator(m.val),
+    )
+    ccall(sbml(:ASTNode_canonicalize), Cint, (VPtr,), astnode)
+    return astnode
+end
+
+function get_astnode_ptr(m::MathVal{<:Real})::VPtr
+    astnode = ccall(sbml(:ASTNode_create), VPtr, ())
+    ccall(sbml(:ASTNode_setReal), Cint, (VPtr, Cdouble), astnode, m.val)
+    ccall(sbml(:ASTNode_canonicalize), Cint, (VPtr,), astnode)
+    return astnode
+end
+
+function get_astnode_ptr(m::MathApply)::VPtr
+    astnode = ccall(sbml(:ASTNode_create), VPtr, ())
+    # Set the name
+    ccall(sbml(:ASTNode_setName), Cint, (VPtr, Cstring), astnode, m.fn)
+    # Set the type
+    if m.fn in keys(all_inv_function_mappings)
+        ccall(
+            sbml(:ASTNode_setType),
+            Cint,
+            (VPtr, Cuint),
+            astnode,
+            all_inv_function_mappings[m.fn],
+        )
+    else
+        ccall(sbml(:ASTNode_setType), Cint, (VPtr, Cuint), astnode, 268) # 268 == AST_FUNCTION
+    end
+    # Add children
+    for child in m.args
+        child_ptr = get_astnode_ptr(child)
+        ccall(sbml(:ASTNode_addChild), Cint, (VPtr, VPtr), astnode, child_ptr)
+    end
+    ccall(sbml(:ASTNode_canonicalize), Cint, (VPtr,), astnode)
+    return astnode
+end
+
+function get_astnode_ptr(m::MathLambda)::VPtr
+    astnode = ccall(sbml(:ASTNode_create), VPtr, ())
+    # All arguments
+    for child in MathIdent.(m.args)
+        child_ptr = get_astnode_ptr(child)
+        ccall(sbml(:ASTNode_addChild), Cint, (VPtr, VPtr), astnode, child_ptr)
+    end
+    # Add the body
+    body = get_astnode_ptr(m.body)
+    ccall(sbml(:ASTNode_addChild), Cint, (VPtr, VPtr), astnode, body)
+    # Set the type
+    ccall(sbml(:ASTNode_setType), Cint, (VPtr, Cuint), astnode, 267) # 267 == AST_LAMBDA
+    # Done
+    return astnode
+end

src/readsbml.jl

@@ -35,6 +35,20 @@ end
 """
 $(TYPEDSIGNATURES)
 
+Like [`get_string`](@ref), but returns `nothing` instead of throwing an
+exception. Also returns values only if `fn_test` returns true.
+"""
+function get_optional_string(x::VPtr, fn_test, fn_sym)::Maybe{String}
+    if ccall(sbml(fn_test), Cint, (VPtr,), x) == 0
+        return nothing
+    else
+        get_optional_string(x, fn_sym)
+    end
+end
+
+"""
+$(TYPEDSIGNATURES)
+
 Helper for getting out boolean flags.
 """
 function get_optional_bool(x::VPtr, is_sym, get_sym)::Maybe{Bool}
@@ -329,6 +343,7 @@ function get_model(mdl::VPtr)::SBML.Model
                 :Species_getHasOnlySubstanceUnits,
             ),
             constant = get_optional_bool(sp, :Species_isSetConstant, :Species_getConstant),
+            metaid = get_optional_string(sp, :SBase_getMetaId),
             notes = get_notes(sp),
             annotation = get_annotation(sp),
         )
@@ -337,7 +352,7 @@ function get_model(mdl::VPtr)::SBML.Model
     # parse out the flux objectives (these are complementary to the objectives
     # that appear in the reactions, see comments lower)
     objectives = Dict{String,Objective}()
-    active_objective = nothing
+    active_objective = ""
     if mdl_fbc != C_NULL
         for i = 1:ccall(sbml(:FbcModelPlugin_getNumObjectives), Cuint, (VPtr,), mdl_fbc)
             flux_objectives = Dict{String,Float64}()
@@ -356,8 +371,7 @@ function get_model(mdl::VPtr)::SBML.Model
             end
             objectives[get_string(o, :Objective_getId)] = Objective(type, flux_objectives)
         end
-        active_objective =
-            get_optional_string(mdl_fbc, :FbcModelPlugin_getActiveObjectiveId)
+        active_objective = get_string(mdl_fbc, :FbcModelPlugin_getActiveObjectiveId)
     end
 
     # reactions!
@@ -386,8 +400,16 @@ function get_model(mdl::VPtr)::SBML.Model
 
         re_fbc = ccall(sbml(:SBase_getPlugin), VPtr, (VPtr, Cstring), re, "fbc")
         if re_fbc != C_NULL
-            lower_bound = get_optional_string(re_fbc, :FbcReactionPlugin_getLowerFluxBound)
-            upper_bound = get_optional_string(re_fbc, :FbcReactionPlugin_getUpperFluxBound)
+            lower_bound = get_optional_string(
+                re_fbc,
+                :FbcReactionPlugin_isSetLowerFluxBound,
+                :FbcReactionPlugin_getLowerFluxBound,
+            )
+            upper_bound = get_optional_string(
+                re_fbc,
+                :FbcReactionPlugin_isSetUpperFluxBound,
+                :FbcReactionPlugin_getUpperFluxBound,
+            )
         end
 
         # extract stoichiometry
@@ -440,6 +462,7 @@ function get_model(mdl::VPtr)::SBML.Model
             gene_product_association = association,
             kinetic_math = math,
             reversible,
+            metaid = get_optional_string(re, :SBase_getMetaId),
             notes = get_notes(re),
             annotation = get_annotation(re),
         )
@@ -461,8 +484,9 @@ function get_model(mdl::VPtr)::SBML.Model
 
             if id != nothing
                 gene_products[id] = GeneProduct(
+                    label = get_string(gp, :GeneProduct_getLabel),
                     name = get_optional_string(gp, :GeneProduct_getName),
-                    label = get_optional_string(gp, :GeneProduct_getLabel),
+                    metaid = get_optional_string(gp, :SBase_getMetaId),
                     notes = get_notes(gp),
                     annotation = get_annotation(gp),
                 )
@@ -519,18 +543,31 @@ function get_model(mdl::VPtr)::SBML.Model
             )
         end
 
-        trig_math_ptr = ccall(
-            sbml(:Trigger_getMath),
-            VPtr,
-            (VPtr,),
-            ccall(sbml(:Event_getTrigger), VPtr, (VPtr,), ev),
+        trigger_ptr = ccall(sbml(:Event_getTrigger), VPtr, (VPtr,), ev)
+        trig_math_ptr = ccall(sbml(:Trigger_getMath), VPtr, (VPtr,), trigger_ptr)
+        trigger = Trigger(;
+            persistent = ccall(sbml(:Trigger_getPersistent), Bool, (VPtr,), trigger_ptr),
+            initial_value = ccall(
+                sbml(:Trigger_getInitialValue),
+                Bool,
+                (VPtr,),
+                trigger_ptr,
+            ),
+            math = trig_math_ptr == C_NULL ? nothing : parse_math(trig_math_ptr),
         )
 
-        events[unsafe_string(ccall(sbml(:Event_getId), Cstring, (VPtr,), ev))] = SBML.Event(
-            get_optional_string(ev, :Event_getName),
-            trig_math_ptr == C_NULL ? nothing : parse_math(trig_math_ptr),
-            event_assignments,
-        )
+        events[unsafe_string(ccall(sbml(:Event_getId), Cstring, (VPtr,), ev))] =
+            SBML.Event(;
+                use_values_from_trigger_time = ccall(
+                    sbml(:Event_getUseValuesFromTriggerTime),
+                    Cint,
+                    (VPtr,),
+                    ev,
+                ),
+                name = get_optional_string(ev, :Event_getName),
+                trigger,
+                event_assignments,
+            )
     end
 
     # Rules
@@ -591,6 +628,7 @@ function get_model(mdl::VPtr)::SBML.Model
         events,
         name = get_optional_string(mdl, :Model_getName),
         id = get_optional_string(mdl, :Model_getId),
+        metaid = get_optional_string(mdl, :SBase_getMetaId),
         conversion_factor = get_optional_string(mdl, :Model_getConversionFactor),
         area_units = get_optional_string(mdl, :Model_getAreaUnits),
         extent_units = get_optional_string(mdl, :Model_getExtentUnits),
@@ -602,6 +640,3 @@ function get_model(mdl::VPtr)::SBML.Model
         annotation = get_annotation(mdl),
     )
 end
-
-# Precompilation statements
-precompile(readSBML, (String,))