Merge pull request #109 from giordano/mg/unitful

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,Number}()
     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,Number}
     compartments::Dict{String,Compartment}
     species::Dict{String,Species}
     reactions::Dict{String,Reaction}

src/utils.jl

@@ -200,3 +200,64 @@ 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.  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)
+)
+
+# Get a `Unitful` quantity out of a `Unit_t`.
+get_unit(u::VPtr) =
+    UNIT_KIND_STRINGS_TO_UNIT[unsafe_string(
+        ccall(
+            sbml(:UnitKind_toString),
+            Cstring,
+            (Cint,),
+            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/data/units.xml

+<?xml version='1.0' encoding='UTF-8'?>
+<sbml xmlns:fbc="http://www.sbml.org/sbml/level3/version1/fbc/version2" xmlns="http://www.sbml.org/sbml/level3/version1/core" level="3" version="1" sboTerm="SBO:0000624" fbc:required="false">
+  <model fbc:strict="true" id="e_coli_core">
+    <listOfUnitDefinitions>
+      <unitDefinition id="non_existent">
+        <listOfUnits>
+          <unit kind="blahblah" scale="-3" multiplier="3.14" exponent="1"/>
+        </listOfUnits>
+      </unitDefinition>
+      <unitDefinition id="no_dimensions">
+        <listOfUnits>
+          <unit kind="dimensionless" scale="3" multiplier="0.001" exponent="1"/>
+          <unit kind="item" scale="1" multiplier="2" exponent="1"/>
+        </listOfUnits>
+      </unitDefinition>
+    </listOfUnitDefinitions>
+  </model>
+</sbml>

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,17 @@ 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"
+
+    m = readSBML(joinpath(@__DIR__, "data", "units.xml"))
+    @test m.units["non_existent"] == 0.00314
+    @test m.units["no_dimensions"] == 20.0
+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")