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

DSL variables created from observable eqs are not equivalent to normal variables

Open isaacsas opened this issue 1 year ago • 2 comments

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

isaacsas avatar Sep 20 '24 16:09 isaacsas

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.

TorkelE avatar Sep 21 '24 14:09 TorkelE

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.

isaacsas avatar Sep 21 '24 20:09 isaacsas

@TorkelE it would be good to fix this for V15 since it technically changes the semantics of these variables.

isaacsas avatar Nov 01 '24 19:11 isaacsas

What's the expected behavior here? should isequal(rn1.X1, rn1.X2) return true?

vyudu avatar Nov 08 '24 17:11 vyudu

If a symbolic is created via an observable the generated symbolic stood be the exact same as if it was created with @variables.

isaacsas avatar Nov 08 '24 17:11 isaacsas

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

vyudu avatar Nov 08 '24 17:11 vyudu

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.

vyudu avatar Nov 08 '24 17:11 vyudu

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)

isaacsas avatar Nov 08 '24 20:11 isaacsas

@TorkelE @vyudu this is fixed now right?

isaacsas avatar Jan 13 '25 15:01 isaacsas

yes

TorkelE avatar Jan 13 '25 15:01 TorkelE