Closed isaacsas closed 6 days ago
Expanding the macro we see that function is becoming Catalyst.f
when used in an equation vs. when used in a reaction, i.e.
f(A,t) = 2*A*t
@macroexpand @reaction_network begin
f(A,t), A --> 0
end
gives
julia> @macroexpand @reaction_network begin
f(A,t), A --> 0
end
:(Catalyst.complete(begin
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:383 =#
var"#402###ps#349" = Catalyst.Num[]
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:384 =#
var"#403#t" = Catalyst.default_t()
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:385 =#
var"#404###vars#350" = Catalyst.Num[]
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:386 =#
begin
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:591 =#
var"#405###352" = begin
var"#406#A" = (identity)((Symbolics.wrap)((SymbolicUtils.setmetadata)(((Sym){(SymbolicUtils.FnType){Catalyst.NTuple{1, Catalyst.Any}, Real}}(:A))((Symbolics.value)(var"#403#t")), Symbolics.VariableSource, (:variables, :A))))
var"#406#A" = (Catalyst.ModelingToolkit).wrap(Catalyst.setmetadata((Catalyst.ModelingToolkit).value(var"#406#A"), (Catalyst.Catalyst).VariableSpecies, true))
begin
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:96 =#
Catalyst.all(!((Catalyst.Catalyst).isconstant) ∘ (Catalyst.ModelingToolkit).value, [var"#406#A"]) || Catalyst.throw(Catalyst.ArgumentError("isconstantspecies metadata can only be used with parameters."))
end
[var"#406#A"]
end
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:592 =#
var"#407###specs#351" = Catalyst.reduce(Catalyst.vcat, (Catalyst.Symbolics).scalarize(var"#405###352"))
end
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:387 =#
()
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:388 =#
var"#408###comps#353" = Catalyst.Num[]
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:389 =#
begin
end
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:391 =#
(Catalyst.Catalyst).remake_ReactionSystem_internal((Catalyst.Catalyst).make_ReactionSystem_internal((Catalyst.Catalyst).CatalystEqType[Catalyst.Reaction(f(var"#406#A", var"#403#t"), [var"#406#A"], nothing, [1], nothing, metadata = [:only_use_rate => false])], var"#403#t", Catalyst.setdiff(Catalyst.union(var"#407###specs#351", var"#404###vars#350", var"#408###comps#353"), []), var"#402###ps#349"; name = Catalyst.gensym(:ReactionSystem), spatial_ivs = nothing, observed = [], continuous_events = [], discrete_events = [], combinatoric_ratelaws = true); default_reaction_metadata = [])
end))
vs.
@macroexpand @reaction_network begin
@equations D(A) ~ f(A,t)
end
giving
@macroexpand @reaction_network begin
@equations D(A) ~ f(A,t)
end
:(Catalyst.complete(begin
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:383 =#
var"#409###ps#354" = Catalyst.Num[]
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:384 =#
var"#410#t" = Catalyst.default_t()
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:385 =#
begin
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:591 =#
var"#411###356" = begin
var"#412#A" = (identity)((Symbolics.wrap)((SymbolicUtils.setmetadata)(((Sym){(SymbolicUtils.FnType){Catalyst.NTuple{1, Catalyst.Any}, Real}}(:A))((Symbolics.value)(var"#410#t")), Symbolics.VariableSource, (:variables, :A))))
[var"#412#A"]
end
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:592 =#
var"#413###vars#355" = Catalyst.reduce(Catalyst.vcat, (Catalyst.Symbolics).scalarize(var"#411###356"))
end
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:386 =#
var"#414###specs#357" = Catalyst.Num[]
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:387 =#
()
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:388 =#
var"#415###comps#358" = Catalyst.Num[]
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:389 =#
begin
var"#416#D" = Catalyst.Differential(var"#410#t")
end
#= /Users/isaacsas/.julia/dev/Catalyst/src/dsl.jl:391 =#
(Catalyst.Catalyst).remake_ReactionSystem_internal((Catalyst.Catalyst).make_ReactionSystem_internal((Catalyst.Catalyst).CatalystEqType[var"#416#D"(var"#412#A") ~ Catalyst.f(var"#412#A", var"#410#t")], var"#410#t", Catalyst.setdiff(Catalyst.union(var"#414###specs#357", var"#413###vars#355", var"#415###comps#358"), []), var"#409###ps#354"; name = Catalyst.gensym(:ReactionSystem), spatial_ivs = nothing, observed = [], continuous_events = [], discrete_events = [], combinatoric_ratelaws = true); default_reaction_metadata = [])
end))
Also the (Catalyst.Catalyst).remake_ReactionSystem_internal
seems quite weird.
It is the same with the make_ReactionSystem_internal((Catalyst.Catalyst)
, which I think I just copied the former of. Maybe here:
Catalyst.remake_ReactionSystem_internal(
Catalyst.make_ReactionSystem_internal(
$rxexprs, $tiv, setdiff(union($spssym, $varssym, $compssym), $obs_syms),
$pssym; name = $name, spatial_ivs = $sivs, observed = $observed_eqs,
continuous_events = $continuous_events_expr,
discrete_events = $discrete_events_expr,
combinatoric_ratelaws = $combinatoric_ratelaws);
default_reaction_metadata = $default_reaction_metadata
)
The Catalyst.
is redundant? There are some more cases, like (Catalyst.Catalyst).isconstant
. Here we also use Catalyst.isconstant
in the code I think.
The equation bit is more serious, good catch. Should have been in the tests, will add that when I have the opportunity. (and try to fix as well)
The repeated Catalyst stuff is definitely weird and we should try to figure out all the places it is appearing. I'm actually surprised it doesn't cause issues.
In a brief test, removing the Catalys.
to:
remake_ReactionSystem_internal(
make_ReactionSystem_internal(
$rxexprs, $tiv, setdiff(union($spssym, $varssym, $compssym), $obs_syms),
$pssym; name = $name, spatial_ivs = $sivs, observed = $observed_eqs,
continuous_events = $continuous_events_expr,
discrete_events = $discrete_events_expr,
combinatoric_ratelaws = $combinatoric_ratelaws);
default_reaction_metadata = $default_reaction_metadata
)
you get expected b behaviour with Catalyst.remake_ReactionSystem_internal(Catalyst.make_ReactionSystem_internal ...
. Haven't investigated any wider implications though.
I removed it there in https://github.com/SciML/Catalyst.jl/pull/1024 and it hasn't seemed to cause any issues. The macro should be running within the Catalyst module's scope, so we shouldn't need to qualify internal Catalyst functions.
Is case 2 above expected to work?
I would think that an @variable
for Iapp
should be auto-created in case two -- @TorkelE is that correct? But maybe we require a tweaked form with it moved to the RHS:
julia> rn = @reaction_network begin
@equations begin
D(A) ~ Iapp
0 ~ -Iapp + f(A,t)
end
end
But note this still gives the same error about Iapp
not being defined.
Note that even this
julia> rn = @reaction_network begin
@variables Iapp(t)
@equations begin
D(A) ~ Iapp
0 ~ -Iapp + f(A,t)
end
end
ERROR: UndefVarError: `f` not defined
Stacktrace:
[1] top-level scope
@ ~/.julia/packages/Catalyst/Dgrwf/src/dsl.jl:392
errors, where I've explicitly declared the variable. (Though at least now it is back to the other error about not seeing f
.)
So I guess there are two issues within this issue; not auto-discovering new variables from algebraic equations, and not allowing the external function. The latter is higher priority to fix I'd say since we can always tell people they need to declare variables within algebraic equations for now.
The discussion is almost a year old know, but I think we decided to be (almost) maximally stringent with inferring variables, and only (some cases of) D(V) ~ ... Gave V being infered. Right now we do not do it for V ~ ... Although I would not mind at all adding that.
There is a section here that goes into all the details https://github.com/SciML/Catalyst.jl/pull/815
(the function thing is still a problem though)
does not work
rn = @reaction_network begin @equations D(A) ~ f(A,t) end
ERROR: UndefVarError:
f
not defined Stacktrace: [1] top-level scope @ ~/.julia/dev/Catalyst/src/dsl.jl:391with the error
with a similar error to case 1.