Commit c1e55991 authored by Miroslav Kratochvil's avatar Miroslav Kratochvil
Browse files

Revive "Use `Unitful` for units"

This reverts^H^H^H^H^H revives commit dc0571e7.

Technically, we can't easily use unitful units for storing the SBML data
because of various semantic concerns, but we can make the unitful value
handling easily accessible by users. This commit exposes the unitful
functionality using a single `unitful` function, similarly as we did with
conversion to Symbolics.
parent 8179525d
......@@ -19,6 +19,7 @@ include("converters.jl")
include("math.jl")
include("readsbml.jl")
include("symbolics.jl")
include("unitful.jl")
include("utils.jl")
sbml(sym::Symbol)::VPtr = dlsym(SBML_jll.libsbml_handle, sym)
......
......@@ -205,12 +205,30 @@ function extractModel(mdl::VPtr)::SBML.Model
end
# parse out the unit definitions
units = Dict{String,Number}()
units = Dict{String,Vector{SBML.UnitPart}}()
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] = get_units(ud)
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)
]
end
# parse out compartment names
compartments = Dict{String,Compartment}()
for i = 1:ccall(sbml(:Model_getNumCompartments), Cuint, (VPtr,), mdl)
......
"""
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
"""
......@@ -186,7 +208,7 @@ $(TYPEDFIELDS)
"""
struct Model
parameters::Dict{String,Float64}
units::Dict{String,Number}
units::Dict{String,Vector{UnitPart}}
compartments::Dict{String,Compartment}
species::Dict{String,Species}
reactions::Dict{String,Reaction}
......
# 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 UNITFUL_KIND_STRING = 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)
)
"""
unitful(u::UnitPart)
Converts an [`UnitPart`](@ref) to a corresponding Unitful unit.
"""
unitful(u::UnitPart) =
(u.multiplier * UNITFUL_KIND_STRING[u.kind] * exp10(u.scale))^u.exponent
"""
unitful(units::Vector{UnitPart})
Converts a SBML unit (i.e., a vector of [`UnitPart`](@ref)s) to a corresponding
Unitful unit.
"""
unitful(u::Vector{UnitPart}) = prod(unitful.(u))
"""
unitful(m::Model, val::Tuple{Float64, String}; default_unit::Number=1.0)
Computes a properly unitful value from a value-unit pair stored in the model
`m`. If the unit is not found in the model, `default_unit` is used; scalar
value is used by default.
"""
unitful(m::Model, val::Tuple{Float64,String}; default_unit::Number = 1.0) =
get(unitful(m.units), val[2], default_unit) * val[1]
"""
unitful(m::Model, val::Tuple{Float64, String}; default_unit::String)
Overload of [`unitful`] that allows specification of the `default_unit` by
string ID.
"""
unitful(m::Model, val::Tuple{Float64,String}; default_unit::String) =
unitful(m, val, default_unit = unitful(m.units[default_unit]))
......@@ -301,73 +301,6 @@ check_errors(
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`.
function get_unit(u::VPtr)
multiplier = ccall(sbml(:Unit_getMultiplier), Cdouble, (VPtr,), u)
unit = UNIT_KIND_STRINGS_TO_UNIT[unsafe_string(
ccall(
sbml(:UnitKind_toString),
Cstring,
(Cint,),
ccall(sbml(:Unit_getKind), Cint, (VPtr,), u),
),
)]
scale = ccall(sbml(:Unit_getScale), Cint, (VPtr,), u)
exponent = ccall(sbml(:Unit_getExponent), Cint, (VPtr,), u)
# See page 44 of
# http://sbml.org/Special/specifications/sbml-level-3/version-2/core/release-2/sbml-level-3-version-2-release-2-core.pdf
return (multiplier * unit * exp10(scale))^exponent
end
# 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)
)
function Base.show(io::IO, ::MIME"text/plain", m::SBML.Model)
print(
io,
......
......@@ -22,7 +22,8 @@ end
sbmlfile * ".does.not.really.exist",
)
@test mdl.units["mmol_per_gDW_per_hr"] 3.6001008028224795 * u"mol * g^-1 * s^-1"
@test SBML.unitful(mdl.units["mmol_per_gDW_per_hr"])
3.6001008028224795 * u"mol * g^-1 * s^-1"
@test length(mdl.compartments) == 2
......
......@@ -131,13 +131,13 @@ 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"
@test SBML.unitful(m.units["volume"]) == 1 * u"L"
@test SBML.unitful(m.units["time"]) == 1 * u"s"
@test SBML.unitful(m.units["substance"]) == 1 * u"mol"
m = readSBML(joinpath(@__DIR__, "data", "custom.xml"))
@test m.units["non_existent"] == 0.00314
@test m.units["no_dimensions"] == 20.0
@test SBML.unitful(m.units["non_existent"]) == 0.00314
@test SBML.unitful(m.units["no_dimensions"]) == 20.0
end
@testset "Initial amounts and concentrations" begin
......
......@@ -3,7 +3,6 @@ using Test, SHA, SparseArrays, Downloads
using SBML
using SBML: Model, Reaction, Species
using Symbolics
using Unitful
@testset "SBML test suite" begin
include("version.jl")
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment