Use `Unitful` for units, instead of a custom struct

Project.toml

@@ -9,10 +9,12 @@ Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 SBML_jll = "bb12108a-f4ef-5f88-8ef3-0b33ff7017f1"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 IfElse = "0.1"
 Symbolics = "0.1.21, 1"
+Unitful = "1"
 julia = "1.5"
 
 [extras]

src/SBML.jl

@@ -4,6 +4,7 @@ using SBML_jll, Libdl
 using SparseArrays
 using Symbolics
 using IfElse
+using Unitful
 
 include("types.jl")
 include("structs.jl")

src/readsbml.jl

@@ -177,30 +177,12 @@ function extractModel(mdl::VPtr)::SBML.Model
     end
 
     # parse out the unit definitions
-    units = Dict{String,Vector{SBML.UnitPart}}()
+    units = Dict{String,Quantity}()
     for i = 1:ccall(sbml(:Model_getNumUnitDefinitions), Cuint, (VPtr,), mdl)
         ud = ccall(sbml(:Model_getUnitDefinition), VPtr, (VPtr, Cuint), mdl, i - 1)
         id = get_string(ud, :UnitDefinition_getId)
-        units[id] = [
-            begin
-                u = ccall(sbml(:UnitDefinition_getUnit), VPtr, (VPtr, Cuint), ud, j - 1)
-                SBML.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 = 1:ccall(sbml(:UnitDefinition_getNumUnits), Cuint, (VPtr,), ud)
-        ]
+        units[id] = get_units(ud)
     end
-
     # parse out compartment names
     compartments = Dict{String,Compartment}()
     for i = 1:ccall(sbml(:Model_getNumCompartments), Cuint, (VPtr,), mdl)

src/structs.jl

 
-"""
-Part of a measurement unit definition that corresponds to the SBML definition
-of `Unit`. For example, the unit "per square megahour", Mh^(-2), is written as:
-
-    SBML.UnitPart("second",  # base SI unit, this says we are measuring time
-             -2,        # exponent, says "per square"
-             6,         # log-10 scale of the unit, says "mega"
-             1/3600)    # second-to-hour multiplier
-
-Compound units (such as "volt-amperes" and "dozens of yards per ounce") are
-built from multiple `UnitPart`s; see the definition of field `units` in
-[`SBML.Model`](@ref).
-"""
-struct UnitPart
-    kind::String
-    exponent::Int
-    scale::Int
-    multiplier::Float64
-    UnitPart(k, e, s, m) = new(k, e, s, m)
-end
-
-
 """
 Abstract type for all kinds of gene product associations
 """
@@ -182,7 +160,7 @@ objects.
 """
 struct Model
     parameters::Dict{String,Float64}
-    units::Dict{String,Vector{UnitPart}}
+    units::Dict{String,Quantity}
     compartments::Dict{String,Compartment}
     species::Dict{String,Species}
     reactions::Dict{String,Reaction}

src/utils.jl

@@ -200,3 +200,57 @@ check_errors(
     error::Exception,
     report_severities = ["Fatal", "Error"],
 ) = Bool(success) || get_error_messages(doc, error, report_severities)
+
+# NOTE: this mapping is valid for Level 3/Version 2, it *may* not be valid for
+# other versions
+const KIND_TO_UNIT = Dict(
+    0  => 1.0 * u"A", # UNIT_KIND_AMPERE
+    1  => ustrip(u"mol^-1", Unitful.Na), # UNIT_KIND_AVOGADRO
+    2  => 1.0 * u"Bq", # UNIT_KIND_BECQUEREL
+    3  => 1.0 * u"cd", # UNIT_KIND_CANDELA
+    4  => 1.0 * u"°C", # UNIT_KIND_CELSIUS
+    5  => 1.0 * u"C", # UNIT_KIND_COULOMB
+    6  => 1, # UNIT_KIND_DIMENSIONLESS
+    7  => 1.0 * u"F", # UNIT_KIND_FARAD
+    8  => 1.0 * u"g", # UNIT_KIND_GRAM
+    9  => 1.0 * u"Gy", # UNIT_KIND_GRAY
+    10 => 1.0 * u"H", # UNIT_KIND_HENRY
+    11 => 1.0 * u"Hz", # UNIT_KIND_HERTZ
+    12 => 1, # UNIT_KIND_ITEM
+    13 => 1.0 * u"J", # UNIT_KIND_JOULE
+    14 => 1.0 * u"kat", # UNIT_KIND_KATAL
+    15 => 1.0 * u"K", # UNIT_KIND_KELVIN
+    16 => 1.0 * u"kg", # UNIT_KIND_KILOGRAM
+    17 => 1.0 * u"L", # UNIT_KIND_LITER
+    18 => 1.0 * u"L", # UNIT_KIND_LITRE
+    19 => 1.0 * u"lm", # UNIT_KIND_LUMEN
+    20 => 1.0 * u"lx", # UNIT_KIND_LUX
+    21 => 1.0 * u"m", # UNIT_KIND_METER
+    22 => 1.0 * u"m", # UNIT_KIND_METRE
+    23 => 1.0 * u"mol", # UNIT_KIND_MOLE
+    24 => 1.0 * u"N", # UNIT_KIND_NEWTON
+    25 => 1.0 * u"Ω", # UNIT_KIND_OHM
+    26 => 1.0 * u"Pa", # UNIT_KIND_PASCAL
+    27 => 1.0 * u"rad", # UNIT_KIND_RADIAN
+    28 => 1.0 * u"s", # UNIT_KIND_SECOND
+    29 => 1.0 * u"S", # UNIT_KIND_SIEMENS
+    30 => 1.0 * u"Sv", # UNIT_KIND_SIEVERT
+    31 => 1.0 * u"sr", # UNIT_KIND_STERADIAN
+    32 => 1.0 * u"T", # UNIT_KIND_TESLA
+    33 => 1.0 * u"V", # UNIT_KIND_VOLT
+    34 => 1.0 * u"W", # UNIT_KIND_WATT
+    35 => 1.0 * u"W", # UNIT_KIND_WEBER
+    36 => 1, # UNIT_KIND_INVALID (let's treat is as a dimensionless quantity)
+)
+
+# Get a `Unitful` quantity out of a `Unit_t`.
+get_unit(u::VPtr) =
+    KIND_TO_UNIT[ccall(sbml(:Unit_getKind), Cint, (VPtr,), u)] ^
+    ccall(sbml(:Unit_getExponent), Cint, (VPtr,), u) *
+    exp10(ccall(sbml(:Unit_getScale), Cint, (VPtr,), u)) *
+    ccall(sbml(:Unit_getMultiplier), Cdouble, (VPtr,), u)
+
+# Get `Unitful` quantity out of a `UnitDefinition_t`.
+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))

test/ecoli_flux.jl

@@ -22,6 +22,8 @@ end
         sbmlfile * ".does.not.really.exist",
     )
 
+    @test mdl.units["mmol_per_gDW_per_hr"] == 2.7777e-7 * u"mol * g^-1 * s^-1"
+
     @test length(mdl.compartments) == 2
 
     mets, rxns, S = getS(mdl)

test/loadmodels.jl

@@ -93,6 +93,13 @@ end
     @test all(contains_time.(r.kinetic_math for (_, r) in m.reactions))
 end
 
+@testset "Units" begin
+    m = readSBML(joinpath(@__DIR__, "data", "sbml00852.xml"))
+    @test m.units["volume"] == 1 * u"L"
+    @test m.units["time"] == 1 * u"s"
+    @test m.units["substance"] == 1 * u"mol"
+end
+
 @testset "Initial amounts and concentrations" begin
     m = readSBML(joinpath(@__DIR__, "data", "sbml00852.xml"))
 

test/runtests.jl

@@ -3,6 +3,7 @@ using Test, SHA, SparseArrays, Downloads
 using SBML
 using SBML: Model, Reaction, Species
 using Symbolics
+using Unitful
 
 @testset "SBML test suite" begin
     include("version.jl")