SciML / Catalyst.jl

Chemical reaction network and systems biology interface for scientific machine learning (SciML). High performance, GPU-parallelized, and O(1) solvers in open source software.
https://docs.sciml.ai/Catalyst/stable/
Other
455 stars 75 forks source link

non-SI units not working #1037

Open isaacsas opened 2 weeks ago

isaacsas commented 2 weeks ago

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.

isaacsas commented 2 weeks ago

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).

ctessum commented 1 week ago

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.

isaacsas commented 1 week ago

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/).

isaacsas commented 1 week ago

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.

isaacsas commented 1 week ago

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

ctessum commented 1 week ago

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 :)

isaacsas commented 1 week ago

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).

ctessum commented 1 week ago

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