Catalyst.jl

Catalyst.jl copied to clipboard

SI units works, i.e.

    @test_nowarn @reaction_network begin
        @ivs t [unit=u"s"]
        @species begin
            X1(t), [unit=u"mol/m^3"]
            Z1(t), [unit=u"mol/m^3"]
            X2(t), [unit=u"mol/m^3"]
            Z2(t), [unit=u"mol/m^3"]
            X3(t), [unit=u"mol/m^3"]
            Y3(t), [unit=u"mol/m^3"]
            Z3(t), [unit=u"mol/m^3"]
        end
        @parameters begin
            k1, [unit=u"(m^6)/(s*mol^2)"]
            v2, [unit=u"(m^6)/(s*mol^2)"]
            K2, [unit=u"mol/m^3"]
            v3, [unit=u"(m^3)/(s*mol)"]
            K3, [unit=u"mol/m^3"]
            n3
        end
        k1*X1, 2X1 --> Z1
        mm(X2, v2, K2), 3X2 --> Z2
        hill(X3, v3, K3, n3), X3 + Y3--> Z3
    end

passes, but non-SI units still have issues right now:

    @test_nowarn @reaction_network begin
        @ivs t [unit=u"s"]
        @species begin
            X1(t), [unit=u"μM"]
            Z1(t), [unit=u"μM"]
            X2(t), [unit=u"μM"]
            Z2(t), [unit=u"μM"]
            X3(t), [unit=u"μM"]
            Y3(t), [unit=u"μM"]
            Z3(t), [unit=u"μM"]
        end
        @parameters begin
            k1, [unit=u"1/(s*(μM)^2)"]
            v2, [unit=u"1/(s*(μM)^2)"]
            K2, [unit=u"μM"]
            v3, [unit=u"1/(s*μM)"]
            K3, [unit=u"μM"]
            n3
        end
        k1*X1, 2X1 --> Z1
        mm(X2, v2, K2), 3X2 --> Z2
        hill(X3, v3, K3, n3), X3 + Y3--> Z3
    end

gives

┌ Warning: Reaction rate laws are expected to have units of 0.0009999999999999998 m⁻³ s⁻¹ mol however, k1*X1(t), 2*X1 --> Z1 has units of 9.999999999999995e-7 m⁻³ s⁻¹ mol.
└ @ Catalyst ~/.julia/dev/Catalyst/src/reactionsystem.jl:1492
┌ Warning: Reaction rate laws are expected to have units of 0.0009999999999999998 m⁻³ s⁻¹ mol however, Catalyst.mm(X2(t), v2, K2), 3*X2 --> Z2 has units of 9.999999999999992e-10 m⁻³ s⁻¹ mol.
└ @ Catalyst ~/.julia/dev/Catalyst/src/reactionsystem.jl:1492
┌ Warning: Reaction rate laws are expected to have units of 0.0009999999999999998 m⁻³ s⁻¹ mol however, Catalyst.hill(X3(t), v3, K3, n3), X3 + Y3 --> Z3 has units of 9.999999999999995e-7 m⁻³ s⁻¹ mol.
└ @ Catalyst ~/.julia/dev/Catalyst/src/reactionsystem.jl:1492
Test Failed at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-R17H3W25T9.0/build/default-honeycrisp-R17H3W25T9-0/julialang/julia-release-1-dot-10/usr/share/julia/stdlib/v1.10/Test/src/Test.jl:903
  Expression: isempty(stderr_content)
   Evaluated: isempty("┌ Warning: Reaction rate laws are expected to have units of 0.0009999999999999998 m⁻³ s⁻¹ mol however, k1*X1(t), 2*X1 --> Z1 has units of 9.999999999999995e-7 m⁻³ s⁻¹ mol.\n└ @ Catalyst ~/.julia/dev/Catalyst/src/reactionsystem.jl:1492\n┌ Warning: Reaction rate laws are expected to have units of 0.0009999999999999998 m⁻³ s⁻¹ mol however, Catalyst.mm(X2(t), v2, K2), 3*X2 --> Z2 has units of 9.999999999999992e-10 m⁻³ s⁻¹ mol.\n└ @ Catalyst ~/.julia/dev/Catalyst/src/reactionsystem.jl:1492\n┌ Warning: Reaction rate laws are expected to have units of 0.0009999999999999998 m⁻³ s⁻¹ mol however, Catalyst.hill(X3(t), v3, K3, n3), X3 + Y3 --> Z3 has units of 9.999999999999995e-7 m⁻³ s⁻¹ mol.\n└ @ Catalyst ~/.julia/dev/Catalyst/src/reactionsystem.jl:1492\n")

ERROR: There was an error during testing

So it seems there are still issues with unit conversion to SI going on.

This really requires someone to dig into MTK and DynamicQuantites to understand how they are handling numeric coefficients that arise when converting and comparing units, and ensure that only numeric components arising from unit conversions are getting compared (i.e. numeric components that are part of expressions like X^2/2 should get dropped probably before conversion, including when tracing through registered functions).

MTK version 9.34 supports non-SI units as described here: https://github.com/SciML/ModelingToolkit.jl/pull/2621 . However Catalyst is currently pinned to MTK v9.32. So it may just be as simple as updating that dependency.

Catalyst isn't pinned to 9.32; are you sure this isn't something you did locally? I currently have ModelingToolkit 9.35.1 with Catalyst 14.4 locally.

Note that the line

ModelingToolkit = "9.32"

in the Project.toml means support is declared for any 9.X release with X >= 32 (see https://pkgdocs.julialang.org/v1/compatibility/).

I suspect the issue here is that ModelingToolkit converts non-SI units to SI internally it seems, and that can lead to expressions where numeric coefficients in a symbolic expression get merged with numeric coefficients generated from the unit conversions (like in converting micromolar to mol/m^3). There is then no effective way to unit check after this point. I think expressions need to have numeric coefficients removed prior to unit conversion to avoid this issue.

Or perhaps a way is needed for users to indicate how registered functions transform units to avoid having to try to trace through them.

Catalyst isn't pinned to 9.32; are you sure this isn't something you did locally? I currently have ModelingToolkit 9.35.1 with Catalyst 14.4 locally.

Note that the line

ModelingToolkit = "9.32"

in the Project.toml means support is declared for any 9.X release with X >= 32 (see https://pkgdocs.julialang.org/v1/compatibility/).

Oh, sorry, I guess I should have read the manual more carefully :)

No worries! I still get messed up on aspects of semver in Julia (and finding out where to even read about it can take a bit of effort).

This may be related to https://github.com/SciML/ModelingToolkit.jl/issues/3017