Catalyst.jl
Catalyst.jl copied to clipboard
DSL variables created from observable eqs are not equivalent to normal variables
rn1 = @reaction_network rn_observed begin
@variables X1(t)
@observables X2(t) ~ X1(t)
k, A --> 0
end
when expanded gives
:(Catalyst.complete(begin
var"#7616#t" = Catalyst.default_t()
begin
var"#7617###1945" = begin
var"#7618#k" = (ModelingToolkit.toparam)((Symbolics.wrap)((SymbolicUtils.setmetadata)((Sym){Real}(:k), Symbolics.VariableSource, (:parameters, :k))))
[var"#7618#k"]
end
var"#7619###ps#1944" = Catalyst.reduce(Catalyst.vcat, (Catalyst.Symbolics).scalarize(var"#7617###1945"))
end
begin
var"#7620###1947" = begin
var"#7621#X1" = (identity)((Symbolics.wrap)((SymbolicUtils.setmetadata)(((Sym){(SymbolicUtils.FnType){Catalyst.NTuple{1, Catalyst.Any}, Real}}(:X1))((Symbolics.value)(var"#7616#t")), Symbolics.VariableSource, (:variables, :X1))))
[var"#7621#X1"]
end
var"#7622###vars#1946" = Catalyst.reduce(Catalyst.vcat, (Catalyst.Symbolics).scalarize(var"#7620###1947"))
end
begin
var"#7623###1949" = begin
var"#7624#A" = (identity)((Symbolics.wrap)((SymbolicUtils.setmetadata)(((Sym){(SymbolicUtils.FnType){Catalyst.NTuple{1, Catalyst.Any}, Real}}(:A))((Symbolics.value)(var"#7616#t")), Symbolics.VariableSource, (:variables, :A))))
var"#7624#A" = (Catalyst.ModelingToolkit).wrap(Catalyst.setmetadata((Catalyst.ModelingToolkit).value(var"#7624#A"), (Catalyst.Catalyst).VariableSpecies, true))
begin
Catalyst.all(!((Catalyst.Catalyst).isconstant) ∘ (Catalyst.ModelingToolkit).value, [var"#7624#A"]) || Catalyst.throw(Catalyst.ArgumentError("isconstantspecies metadata can only be used with parameters."))
end
[var"#7624#A"]
end
var"#7625###specs#1948" = Catalyst.reduce(Catalyst.vcat, (Catalyst.Symbolics).scalarize(var"#7623###1949"))
end
begin
begin
var"#7626#X2" = (identity)((Symbolics.wrap)((SymbolicUtils.setmetadata)((Symbolics.CallWithMetadata)((Sym){(SymbolicUtils.FnType){Tuple, Real}}(:X2)), Symbolics.VariableSource, (:variables, :X2))))
[var"#7626#X2"]
end
var"#7626#X2" = var"#7626#X2"(Catalyst.sort(Catalyst.unique(Catalyst.reduce(Catalyst.vcat, [Catalyst.sorted_arguments((Catalyst.MT).unwrap(var"#7627#dep")) for var"#7627#dep" = Catalyst.filter(!((Catalyst.MT).isparameter), (Catalyst.Symbolics).get_variables(var"#7621#X1"(var"#7616#t")))])); by = ((var"#7636#iv",)->begin
Catalyst.findfirst(((Catalyst.MT).getname(var"#7636#iv") == var"#7635#ivs" for var"#7635#ivs" = [:t]))
end))...)
end
var"#7628###comps#1950" = Catalyst.Num[]
begin
end
var"#7629#sivs_vec" = nothing
var"#7630#rx_eq_vec" = Catalyst.CatalystEqType[Catalyst.Reaction(var"#7618#k", [var"#7624#A"], nothing, [1], nothing, metadata = [:only_use_rate => false])]
var"#7631#vars" = Catalyst.setdiff(Catalyst.union(var"#7625###specs#1948", var"#7622###vars#1946", var"#7628###comps#1950"), [var"#7626#X2"])
var"#7632#obseqs" = [var"#7626#X2" ~ var"#7621#X1"(var"#7616#t")]
var"#7633#cevents" = []
var"#7634#devents" = []
Catalyst.remake_ReactionSystem_internal(Catalyst.make_ReactionSystem_internal(var"#7630#rx_eq_vec", var"#7616#t", var"#7631#vars", var"#7619###ps#1944"; name = :rn_observed, spatial_ivs = var"#7629#sivs_vec", observed = var"#7632#obseqs", continuous_events = var"#7633#cevents", discrete_events = var"#7634#devents", combinatoric_ratelaws = true); default_reaction_metadata = [])
end))
Notice the difference in X1's definition vs. X2:
var"#7621#X1" = (identity)((Symbolics.wrap)((SymbolicUtils.setmetadata)(((Sym){(SymbolicUtils.FnType){Catalyst.NTuple{1, Catalyst.Any}, Real}}(:X1))((Symbolics.value)(var"#7616#t")), Symbolics.VariableSource, (:variables, :X1))))
vs
var"#7626#X2" = (identity)((Symbolics.wrap)((SymbolicUtils.setmetadata)((Symbolics.CallWithMetadata)((Sym){(SymbolicUtils.FnType){Tuple, Real}}(:X2)), Symbolics.VariableSource, (:variables, :X2))))
I wrote a comment on it here: https://github.com/SciML/Catalyst.jl/pull/1056#issuecomment-2364056963, but I really need someone who knows more about Symobilcs to explain, as stuff like these just isn't documented at all.
This seems to be the culprit: https://github.com/SciML/Catalyst.jl/blob/005a58fa70b8a5a85073a60c788cf55b45136ce7/src/dsl.jl#L788
We need to instead directly build the observable variables with the independent variables specified.
@TorkelE it would be good to fix this for V15 since it technically changes the semantics of these variables.
What's the expected behavior here? should isequal(rn1.X1, rn1.X2) return true?
If a symbolic is created via an observable the generated symbolic stood be the exact same as if it was created with @variables.
Hmm what would be the way to check that? I don't think there's a way that they could generate the same code when expanded if the observable IVs are inferred automatically
I think I got observables to generate the same way as variables but not sure what sort of tests to write to check they are actually the same.
Maybe something like this:
rn1 = @reaction_network rn_observed begin
@variables X1(t)
@observables X2 ~ X1
k, A --> 0
end
t = Catalyst.default_t()
@variables X2(t)
# this is false on the last release, but should be true:
isequal(rn1.X2, X2)
@TorkelE @vyudu this is fixed now right?
yes