DSL issues with equations

[isaacsas]

1. Using an external function in a reaction works fine but has issues in equations

# works
f(A,t) = 2*A*t
rn = @reaction_network begin
    f(A,t), A --> 0
end 

# does not work
rn = @reaction_network begin
    @equations D(A) ~ f(A,t)
end

with the latter giving the error

ERROR: UndefVarError: `f` not defined
Stacktrace:
 [1] top-level scope
   @ ~/.julia/dev/Catalyst/src/dsl.jl:391

2. This doesn't work

rn = @reaction_network begin
    @equations begin
        D(A) ~ Iapp
        Iapp ~ f(A,t)
    end
end

with the error

ERROR: UndefVarError: `Iapp` not defined
Stacktrace:
 [1] top-level scope
   @ ~/.julia/dev/Catalyst/src/dsl.jl:391

3. This also doesn't work

rn = @reaction_network begin
    @variables Iapp(t)
    @equations begin
        D(A) ~ Iapp
        Iapp ~ f(A,t)
    end
end

with a similar error to case 1.

[isaacsas]

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

[isaacsas]

Also the (Catalyst.Catalyst).remake_ReactionSystem_internal seems quite weird.

[TorkelE]

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)

[isaacsas]

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.

[TorkelE]

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.

[isaacsas]

I removed it there in #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.

[vyudu]

Is case 2 above expected to work?

[isaacsas]

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.

[TorkelE]

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 #815

[TorkelE]

(the function thing is still a problem though)