ModelingToolkit.jl icon indicating copy to clipboard operation
ModelingToolkit.jl copied to clipboard

Deprication of `all_dimensionless` and how to adapt unit check in Catalyst?

Open TorkelE opened this issue 1 year ago • 3 comments

When I try to run

using Catalyst
rn = @reaction_network begin
    (p,d), 0 <--> X
end

on MTK#master I get an error:

ERROR: MethodError: no method matching all_dimensionless(::Vector{SymbolicUtils.BasicSymbolic{Real}})
Stacktrace:
 [1] ReactionSystem(eqs::Vector{…}, rxs::Vector{…}, iv::SymbolicUtils.BasicSymbolic{…}, sivs::Vector{…}, states::Vector{…}, spcs::Vector{…}, ps::Vector{…}, var_to_name::Dict{…}, observed::Vector{…}, name::Symbol, systems::Vector{…}, defaults::Dict{…}, connection_type::Nothing, nps::Catalyst.NetworkProperties{…}, cls::Bool, cevs::Vector{…}, devs::Vector{…}, complete::Bool; checks::Bool)
   @ Catalyst C:\Users\Torkel\.julia\packages\Catalyst\quz1b\src\reactionsystem.jl:532
 [2] ReactionSystem
   @ Catalyst C:\Users\Torkel\.julia\packages\Catalyst\quz1b\src\reactionsystem.jl:516 [inlined]
 [3] ReactionSystem(eqs::Vector{…}, iv::SymbolicUtils.BasicSymbolic{…}, states::Vector{…}, ps::Vector{…}; observed::Vector{…}, systems::Vector{…}, name::Symbol, default_u0::Dict{…}, default_p::Dict{…}, defaults::Dict{…}, connection_type::Nothing, checks::Bool, networkproperties::Nothing, combinatoric_ratelaws::Bool, balanced_bc_check::Bool, spatial_ivs::Nothing, continuous_events::Nothing, discrete_events::Nothing)
   @ Catalyst C:\Users\Torkel\.julia\packages\Catalyst\quz1b\src\reactionsystem.jl:648
 [4] make_ReactionSystem_internal(rxs_and_eqs::Vector{…}, iv::Num, sts_in::Vector{…}, ps_in::Vector{…}; spatial_ivs::Nothing, kwargs::@Kwargs{…})
   @ Catalyst C:\Users\Torkel\.julia\packages\Catalyst\quz1b\src\reactionsystem.jl:721
 [5] top-level scope
   @ C:\Users\Torkel\.julia\packages\Catalyst\quz1b\src\reaction_network.jl:414
Some type information was truncated. Use `show(err)` to see complete types.

Investigating, it seems like the MTK all_dimensionless function (which Catalyst depended on) has been removed on master(?), causing the error.

This makes sense I guess, given that MTK is redoing how it handles units. However, we need to replace it in Catalyst. What is the way to do this in MTK v9?

Here are the specific lines of code where the error is encountered:

if isempty(sivs) && (checks == true || (checks & MT.CheckUnits) > 0)
    nonrx_eqs = Equation[eq for eq in eqs if eq isa Equation]
    MT.all_dimensionless([states; ps; iv]) || check_units(nonrx_eqs)
end

Generally, all_dimensionless does not seem to be used anymore, and MTK seems to have gone from

 if checks == true || (checks & CheckUnits) > 0
    all_dimensionless([dvs; ps; iv]) || check_units(deqs)
end

to

if checks == true || (checks & CheckUnits) > 0
    u = __get_unit_type(dvs, ps, iv)
    check_units(u, deqs)
end

How should I adapt this to Catalyst?

TorkelE avatar Feb 13 '24 03:02 TorkelE

@TorkelE have you tried just using the new check that MTK uses?

isaacsas avatar Feb 13 '24 19:02 isaacsas

So, if I change

if isempty(sivs) && (checks == true || (checks & MT.CheckUnits) > 0)
    nonrx_eqs = Equation[eq for eq in eqs if eq isa Equation]
    MT.all_dimensionless([states; ps; iv]) || check_units(nonrx_eqs)
end

to

if isempty(sivs) && (checks == true || (checks & MT.CheckUnits) > 0)
    nonrx_eqs = Equation[eq for eq in eqs if eq isa Equation]
    u = MT.__get_unit_type(iv, sivs, unknowns, spcs, ps)
    check_units(u, nonrx_eqs)
end

and then change

# no units for species, time or parameters then assume validated
(specunits in (MT.unitless, nothing)) && (timeunits in (MT.unitless, nothing)) && MT.all_dimensionless(get_ps(rs)) && return true

to

# no units for species, time or parameters then assume validated
(specunits in (MT.unitless, nothing)) && (timeunits in (MT.unitless, nothing)) &&  return true

I can create ReactionSystems again. However, in the second step I am just skipping the all_dimensionless check, and there are some decent difference going on in the first case as well. I feel like it would be useful to know what is actually going on here.

TorkelE avatar Feb 13 '24 19:02 TorkelE

I have now replace

ModelingToolkit.all_dimensionless(vars)

with

all(u == 1.0 for u in ModelingToolkit.get_unit(vars))

If this seems equivalent (which I think it is), I am happy to close this.

TorkelE avatar Feb 15 '24 15:02 TorkelE