SciML / Catalyst.jl

Chemical reaction network and systems biology interface for scientific machine learning (SciML). High performance, GPU-parallelized, and O(1) solvers in open source software.
https://docs.sciml.ai/Catalyst/stable/
Other
463 stars 78 forks source link

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

Open isaacsas opened 1 month ago

isaacsas commented 1 month ago
    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))))
TorkelE commented 1 month ago

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.

isaacsas commented 1 month ago

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 commented 2 weeks ago

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

vyudu commented 1 week ago

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

isaacsas commented 1 week ago

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

vyudu commented 1 week ago

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 commented 1 week ago

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.

isaacsas commented 1 week ago

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)