From 9a43432a93b38643a79f52221c3d02f3c1bfc4ed Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Tue, 3 Sep 2024 17:43:02 -0400
Subject: [PATCH 01/15] equations fix

---
 src/dsl.jl | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index 803b999d1b..5ef781cec2 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -297,7 +297,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
 
     # Get macro options.
     if length(unique(arg.args[1] for arg in option_lines)) < length(option_lines)
-        error("Some options where given multiple times.")
+        error("Some options were given multiple times.")
     end
     options = Dict(map(arg -> Symbol(String(arg.args[1])[2:end]) => arg,
         option_lines))
@@ -376,7 +376,8 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         push!(rxexprs.args, get_rxexprs(reaction))
     end
     for equation in equations
-        push!(rxexprs.args, equation)
+        equation = expand_equation_RHS!(equation)
+        push!(rxexprs.args, (equation))
     end
 
     quote
@@ -915,6 +916,13 @@ function recursive_expand_functions!(expr::ExprValues)
     expr
 end
 
+# Recursively expand the right-hand-side of an equation written using user-defined functions and special function calls like "hill(...)" with the actual corresponding expression. 
+function expand_equation_RHS!(eq::Expr) 
+    rhs = recursive_expand_functions!(eq.args[3])
+    eq.args[3] = rhs
+    eq
+end
+
 # Returns the length of a expression tuple, or 1 if it is not an expression tuple (probably a Symbol/Numerical).
 function tup_leng(ex::ExprValues)
     (typeof(ex) == Expr && ex.head == :tuple) && (return length(ex.args))

From a46df6dbd750914e5d8b75051ae5f4425135ebf5 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Tue, 3 Sep 2024 17:48:50 -0400
Subject: [PATCH 02/15] format

---
 src/dsl.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index 5ef781cec2..946c5b298a 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -917,7 +917,7 @@ function recursive_expand_functions!(expr::ExprValues)
 end
 
 # Recursively expand the right-hand-side of an equation written using user-defined functions and special function calls like "hill(...)" with the actual corresponding expression. 
-function expand_equation_RHS!(eq::Expr) 
+function expand_equation_RHS!(eq::Expr)
     rhs = recursive_expand_functions!(eq.args[3])
     eq.args[3] = rhs
     eq

From af712c30ec626c8763704f3a82223871264267c2 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Tue, 3 Sep 2024 22:49:04 -0400
Subject: [PATCH 03/15] parameter inferring inn equations

---
 Project.toml            |  3 ++-
 src/dsl.jl              | 15 +++++++++++-
 test/dsl/dsl_options.jl | 51 +++++++++++++++++++++++++++++++++++++++++
 3 files changed, 67 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 034f475666..8b7dd43d51 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "Catalyst"
 uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
-version = "14.4"
+version = "14.4.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -16,6 +16,7 @@ Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
diff --git a/src/dsl.jl b/src/dsl.jl
index 946c5b298a..a4722393d9 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -347,6 +347,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         variables)))
     species_extracted, parameters_extracted = extract_species_and_parameters!(reactions,
         declared_syms)
+    parameters_extracted_eq = extract_parameters_eq(equations, declared_syms)
+    parameters_extracted = vcat(parameters_extracted, parameters_extracted_eq)
+
     species = vcat(species_declared, species_extracted)
     parameters = vcat(parameters_declared, parameters_extracted)
 
@@ -377,7 +380,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     end
     for equation in equations
         equation = expand_equation_RHS!(equation)
-        push!(rxexprs.args, (equation))
+        push!(rxexprs.args, equation)
     end
 
     quote
@@ -526,6 +529,16 @@ function extract_species_and_parameters!(reactions, excluded_syms)
     collect(species), collect(parameters)
 end
 
+# Find undeclared parameters in the RHS of the equations. 
+function extract_parameters_eq(equations, excluded_syms)
+    parameters = OrderedSet{Union{Symbol, Expr}}()
+    for eq in equations
+        add_syms_from_expr!(parameters, eq.args[3], excluded_syms)
+    end
+
+    collect(parameters)
+end
+
 # Function called by extract_species_and_parameters!, recursively loops through an
 # expression and find symbols (adding them to the push_symbols vector).
 function add_syms_from_expr!(push_symbols::AbstractSet, rateex::ExprValues, excluded_syms)
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 2aecc8fffb..1ebe66f16a 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -13,6 +13,7 @@ seed = rand(rng, 1:100)
 
 # Sets the default `t` to use.
 t = default_t()
+D = default_time_deriv()
 
 ### Tests `@parameters`, `@species`, and `@variables` Options ###
 
@@ -152,6 +153,56 @@ let
     @test issetequal(parameters(rn5), @parameters k B X)
 end
 
+# Test whether user-defined functions are properly expanded in equations. 
+# TODO
+let
+    f(A, t) = 2*A*t
+
+    # Test user-defined function
+    rn = @reaction_network begin
+        @equations D(A) ~ f(A, t)
+    end
+    @test length(equations(rn)) == 1
+    @test equations(rn)[1] isa Equation
+    @species A(t)
+    @test isequal(equations(rn)[1], D(A) ~ 2*A*t)
+
+
+    # Test Catalyst function
+    rn2 = @reaction_network begin
+        @equations D(A) ~ hill(A, v, K, n)
+    end
+    @test length(equations(rn2)) == 1
+    @test equations(rn2)[1] isa Equation
+    @parameters K v n
+    @test isequal(Catalyst.expand_registered_functions(equations(rn2)[1]), D(A) ~ v*(A^n) / (A^n + K^n))
+
+
+    rn3 = @reaction_network begin
+        @species Iapp(t)
+        @equations begin
+            D(A) ~ Iapp
+            Iapp ~ f(A,t)
+        end
+    end
+    @test length(equations(rn3)) == 2
+    @test equations(rn3)[1] isa Equation
+    @test equations(rn3)[2] isa Equation
+    @variables Iapp(t)
+    @test isequal(equations(rn3)[1], D(A) ~ Iapp)
+    @test isequal(equations(rn3)[2], Iapp ~ 2*A*t)
+
+
+    g(A, K, n) = A^n + K^n
+    rn4 = @reaction_network begin
+        @equations D(A) ~ hill(A, v, K, n)*g(A, K, n)
+    end
+    @test length(equations(rn4)) == 1
+    @test equations(rn4)[1] isa Equation
+    @parameters v n
+    @test isequal(Catalyst.expand_registered_functions(equations(rn4)[1]), D(A) ~ v*(A^n))
+end
+
 # Test inferring with stoichiometry symbols and interpolation.
 let
     @parameters k g h gg X y [isconstantspecies = true]

From d07d7bbd3b032c44fd8b8e0deacc8345abed92be Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Tue, 3 Sep 2024 22:49:59 -0400
Subject: [PATCH 04/15] up

---
 Project.toml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 8b7dd43d51..d0a7c2730c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,7 +16,6 @@ Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
-OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"

From e20d90da62caee851dcf65fa5cd187340da10f47 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 4 Sep 2024 07:51:38 -0400
Subject: [PATCH 05/15] undo param-inferring

---
 src/dsl.jl | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index a4722393d9..c0f65e84e5 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -315,12 +315,12 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     parameters_declared = extract_syms(options, :parameters)
     variables_declared = extract_syms(options, :variables)
 
-    # Reads more options.
+    # Reads equations. 
     vars_extracted, add_default_diff, equations = read_equations_options(
         options, variables_declared)
     variables = vcat(variables_declared, vars_extracted)
 
-    # handle independent variables
+    # Handle independent variables
     if haskey(options, :ivs)
         ivs = Tuple(extract_syms(options, :ivs))
         ivexpr = copy(options[:ivs])
@@ -339,17 +339,15 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         combinatoric_ratelaws = true
     end
 
-    # Reads more options.
+    # Reads observables.
     observed_vars, observed_eqs, obs_syms = read_observed_options(
         options, [species_declared; variables], all_ivs)
 
+    # Collect species and parameters, including ones inferred from the reactions. 
     declared_syms = Set(Iterators.flatten((parameters_declared, species_declared,
         variables)))
     species_extracted, parameters_extracted = extract_species_and_parameters!(reactions,
         declared_syms)
-    parameters_extracted_eq = extract_parameters_eq(equations, declared_syms)
-    parameters_extracted = vcat(parameters_extracted, parameters_extracted_eq)
-
     species = vcat(species_declared, species_extracted)
     parameters = vcat(parameters_declared, parameters_extracted)
 
@@ -383,6 +381,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         push!(rxexprs.args, equation)
     end
 
+    # Output code corresponding to the reaction system. 
     quote
         $ivexpr
         $ps

From 05f1a7248a5b9f21c91fe42e134c9f48123e58b4 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 4 Sep 2024 08:52:11 -0400
Subject: [PATCH 06/15] up

---
 src/dsl.jl              | 20 +++++++++++++-------
 test/dsl/dsl_options.jl |  9 +++++++--
 2 files changed, 20 insertions(+), 9 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index c0f65e84e5..fa87a4f2d0 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -346,8 +346,14 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     # Collect species and parameters, including ones inferred from the reactions. 
     declared_syms = Set(Iterators.flatten((parameters_declared, species_declared,
         variables)))
-    species_extracted, parameters_extracted = extract_species_and_parameters!(reactions,
-        declared_syms)
+    species_extracted, parameters_extracted = extract_species_and_parameters!(reactions, declared_syms)
+
+    # Checks that all the symbols in the equation are declared. 
+    eq_syms = extract_syms_eq(equations, union(declared_syms, species_extracted, parameters_extracted))
+    if !isempty(eq_syms)
+        error("Symbols in equations must be explicitly defined using @variable, @species, or @parameter within the @reaction_network input.")
+    end
+
     species = vcat(species_declared, species_extracted)
     parameters = vcat(parameters_declared, parameters_extracted)
 
@@ -528,14 +534,14 @@ function extract_species_and_parameters!(reactions, excluded_syms)
     collect(species), collect(parameters)
 end
 
-# Find undeclared parameters in the RHS of the equations. 
-function extract_parameters_eq(equations, excluded_syms)
-    parameters = OrderedSet{Union{Symbol, Expr}}()
+# Find undeclared symbols in the RHS of the equations. 
+function extract_syms_eq(equations, excluded_syms)
+    syms = OrderedSet{Union{Symbol, Expr}}()
     for eq in equations
-        add_syms_from_expr!(parameters, eq.args[3], excluded_syms)
+        add_syms_from_expr!(syms, eq.args[3], excluded_syms)
     end
 
-    collect(parameters)
+    collect(syms)
 end
 
 # Function called by extract_species_and_parameters!, recursively loops through an
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 1ebe66f16a..b01e053291 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -154,7 +154,6 @@ let
 end
 
 # Test whether user-defined functions are properly expanded in equations. 
-# TODO
 let
     f(A, t) = 2*A*t
 
@@ -167,14 +166,19 @@ let
     @species A(t)
     @test isequal(equations(rn)[1], D(A) ~ 2*A*t)
 
+    # Tests that error is declared. 
+    @test_throws LoadError @eval @reaction_network begin
+        @equations D(A) ~ hill(A, v, K, n)
+    end
 
     # Test Catalyst function
     rn2 = @reaction_network begin
+        @parameters v K n
         @equations D(A) ~ hill(A, v, K, n)
     end
     @test length(equations(rn2)) == 1
     @test equations(rn2)[1] isa Equation
-    @parameters K v n
+    @parameters v K n
     @test isequal(Catalyst.expand_registered_functions(equations(rn2)[1]), D(A) ~ v*(A^n) / (A^n + K^n))
 
 
@@ -195,6 +199,7 @@ let
 
     g(A, K, n) = A^n + K^n
     rn4 = @reaction_network begin
+        @parameters v K n
         @equations D(A) ~ hill(A, v, K, n)*g(A, K, n)
     end
     @test length(equations(rn4)) == 1

From 0ab37c3814f2638130de88d2e83dcf5c649cc7ef Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 4 Sep 2024 09:02:53 -0400
Subject: [PATCH 07/15] revert, splitting off params-related changes

---
 src/dsl.jl              | 16 ----------------
 test/dsl/dsl_options.jl |  4 ----
 2 files changed, 20 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index fa87a4f2d0..0f6514fbb7 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -348,12 +348,6 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         variables)))
     species_extracted, parameters_extracted = extract_species_and_parameters!(reactions, declared_syms)
 
-    # Checks that all the symbols in the equation are declared. 
-    eq_syms = extract_syms_eq(equations, union(declared_syms, species_extracted, parameters_extracted))
-    if !isempty(eq_syms)
-        error("Symbols in equations must be explicitly defined using @variable, @species, or @parameter within the @reaction_network input.")
-    end
-
     species = vcat(species_declared, species_extracted)
     parameters = vcat(parameters_declared, parameters_extracted)
 
@@ -534,16 +528,6 @@ function extract_species_and_parameters!(reactions, excluded_syms)
     collect(species), collect(parameters)
 end
 
-# Find undeclared symbols in the RHS of the equations. 
-function extract_syms_eq(equations, excluded_syms)
-    syms = OrderedSet{Union{Symbol, Expr}}()
-    for eq in equations
-        add_syms_from_expr!(syms, eq.args[3], excluded_syms)
-    end
-
-    collect(syms)
-end
-
 # Function called by extract_species_and_parameters!, recursively loops through an
 # expression and find symbols (adding them to the push_symbols vector).
 function add_syms_from_expr!(push_symbols::AbstractSet, rateex::ExprValues, excluded_syms)
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index b01e053291..35778c97b9 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -166,10 +166,6 @@ let
     @species A(t)
     @test isequal(equations(rn)[1], D(A) ~ 2*A*t)
 
-    # Tests that error is declared. 
-    @test_throws LoadError @eval @reaction_network begin
-        @equations D(A) ~ hill(A, v, K, n)
-    end
 
     # Test Catalyst function
     rn2 = @reaction_network begin

From 44bb0c39cdf2b3747ed567b1b0c374606bfbea40 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 4 Sep 2024 09:22:34 -0400
Subject: [PATCH 08/15] format

---
 src/dsl.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index 0f6514fbb7..d932475559 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -346,7 +346,8 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     # Collect species and parameters, including ones inferred from the reactions. 
     declared_syms = Set(Iterators.flatten((parameters_declared, species_declared,
         variables)))
-    species_extracted, parameters_extracted = extract_species_and_parameters!(reactions, declared_syms)
+    species_extracted, parameters_extracted = extract_species_and_parameters!(
+        reactions, declared_syms)
 
     species = vcat(species_declared, species_extracted)
     parameters = vcat(parameters_declared, parameters_extracted)

From 05dde891e5f0647a9c2dc2ed97acb624f016b194 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 4 Sep 2024 11:52:56 -0400
Subject: [PATCH 09/15] up

---
 src/dsl.jl              | 24 ++++++++++++++----------
 test/dsl/dsl_options.jl | 17 ++++++++++++++---
 2 files changed, 28 insertions(+), 13 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index d932475559..684cdc1d92 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -378,7 +378,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         push!(rxexprs.args, get_rxexprs(reaction))
     end
     for equation in equations
-        equation = expand_equation_RHS!(equation)
+        equation = escape_equation_RHS!(equation)
         push!(rxexprs.args, equation)
     end
 
@@ -576,7 +576,7 @@ function get_rxexprs(rxstruct)
     subs_stoich_init = deepcopy(subs_init)
     prod_init = isempty(rxstruct.products) ? nothing : :([])
     prod_stoich_init = deepcopy(prod_init)
-    reaction_func = :(Reaction($(recursive_expand_functions!(rxstruct.rate)), $subs_init,
+    reaction_func = :(Reaction($(recursive_escape_functions!(rxstruct.rate)), $subs_init,
         $prod_init, $subs_stoich_init, $prod_stoich_init,
         metadata = $(rxstruct.metadata)))
     for sub in rxstruct.substrates
@@ -908,10 +908,10 @@ end
 
 ### Generic Expression Manipulation ###
 
-# Recursively traverses an expression and replaces special function call like "hill(...)" with the actual corresponding expression.
-function recursive_expand_functions!(expr::ExprValues)
+# Recursively traverses an expression and escapes all the user-defined functions. Special function calls like "hill(...)" are not expanded. 
+function recursive_escape_functions!(expr::ExprValues)
     (typeof(expr) != Expr) && (return expr)
-    foreach(i -> expr.args[i] = recursive_expand_functions!(expr.args[i]),
+    foreach(i -> expr.args[i] = recursive_escape_functions!(expr.args[i]),
         1:length(expr.args))
     if expr.head == :call
         !isdefined(Catalyst, expr.args[1]) && (expr.args[1] = esc(expr.args[1]))
@@ -919,11 +919,15 @@ function recursive_expand_functions!(expr::ExprValues)
     expr
 end
 
-# Recursively expand the right-hand-side of an equation written using user-defined functions and special function calls like "hill(...)" with the actual corresponding expression. 
-function expand_equation_RHS!(eq::Expr)
-    rhs = recursive_expand_functions!(eq.args[3])
-    eq.args[3] = rhs
-    eq
+# Recursively escape functions in the right-hand-side of an equation written using user-defined functions. Special function calls like "hill(...)" are not expanded.
+function escape_equation_RHS!(eqexpr::Expr)
+    lhs = recursive_escape_functions!(eqexpr.args[2])
+    eqexpr.args[2] = lhs
+
+    rhs = recursive_escape_functions!(eqexpr.args[3])
+    eqexpr.args[3] = rhs
+
+    eqexpr
 end
 
 # Returns the length of a expression tuple, or 1 if it is not an expression tuple (probably a Symbol/Numerical).
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 35778c97b9..4a69e54974 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -167,15 +167,26 @@ let
     @test isequal(equations(rn)[1], D(A) ~ 2*A*t)
 
 
-    # Test Catalyst function
+    hill_unregistered(A, v, K, n) = v*(A^n) / (A^n + K^n)
     rn2 = @reaction_network begin
         @parameters v K n
-        @equations D(A) ~ hill(A, v, K, n)
+        @equations D(A) ~ hill_unregistered(A, v, K, n)
     end
     @test length(equations(rn2)) == 1
     @test equations(rn2)[1] isa Equation
     @parameters v K n
-    @test isequal(Catalyst.expand_registered_functions(equations(rn2)[1]), D(A) ~ v*(A^n) / (A^n + K^n))
+    @test isequal(equations(rn2)[1], D(A) ~ v*(A^n) / (A^n + K^n))
+
+    @register_symbolic hill2(A, v, K, n)
+    # Test Catalyst function
+    rn2r = @reaction_network begin
+        @parameters v K n
+        @equations D(A) ~ hill2(A, v, K, n)
+    end
+    @test length(equations(rn2r)) == 1
+    @test equations(rn2r)[1] isa Equation
+    @parameters v K n
+    @test isequal(equations(rn2r)[1], D(A) ~ hill2(A, v, K, n))
 
 
     rn3 = @reaction_network begin

From a7a33562558fa2cbdcefab469f903ed2ace8bdc9 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 4 Sep 2024 12:10:29 -0400
Subject: [PATCH 10/15] test up

---
 test/dsl/dsl_options.jl | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 4a69e54974..336241341a 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -177,6 +177,7 @@ let
     @parameters v K n
     @test isequal(equations(rn2)[1], D(A) ~ v*(A^n) / (A^n + K^n))
 
+    hill2(A, v, K, n) = v*(A^n) / (A^n + K^n)
     @register_symbolic hill2(A, v, K, n)
     # Test Catalyst function
     rn2r = @reaction_network begin
@@ -213,6 +214,16 @@ let
     @test equations(rn4)[1] isa Equation
     @parameters v n
     @test isequal(Catalyst.expand_registered_functions(equations(rn4)[1]), D(A) ~ v*(A^n))
+
+
+    rn5 = @reaction_network begin
+        @species A(t) B(t)
+        @equations D(A) + D(B) ~ f(A, t)
+    end
+    @test length(equations(rn5)) == 1
+    @test equations(rn5)[1] isa Equation
+    @species B(t)
+    @test isequal(equations(rn5)[1], D(A) + D(B) ~ 2*A*t)
 end
 
 # Test inferring with stoichiometry symbols and interpolation.

From bc4e62072a9e630f1a195fd01e416671b8c9b562 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 4 Sep 2024 13:25:08 -0400
Subject: [PATCH 11/15] up

---
 src/dsl.jl                                               | 6 +++---
 test/reactionsystem_core/coupled_equation_crn_systems.jl | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index 684cdc1d92..a3adeddf16 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -686,7 +686,7 @@ function read_equations_options(options, variables_declared)
     # When this is the case, the variable X and differential D are extracted (for automatic declaration).
     # Also performs simple error checks.
     vars_extracted = Vector{Symbol}()
-    add_default_diff = false
+    add_default_diff = true
     for eq in equations
         if (eq.head != :call) || (eq.args[1] != :~)
             error("Malformed equation: \"$eq\". Equation's left hand and right hand sides should be separated by a \"~\".")
@@ -921,8 +921,8 @@ end
 
 # Recursively escape functions in the right-hand-side of an equation written using user-defined functions. Special function calls like "hill(...)" are not expanded.
 function escape_equation_RHS!(eqexpr::Expr)
-    lhs = recursive_escape_functions!(eqexpr.args[2])
-    eqexpr.args[2] = lhs
+    # lhs = recursive_escape_functions!(eqexpr.args[2])
+    # eqexpr.args[2] = lhs
 
     rhs = recursive_escape_functions!(eqexpr.args[3])
     eqexpr.args[3] = rhs
diff --git a/test/reactionsystem_core/coupled_equation_crn_systems.jl b/test/reactionsystem_core/coupled_equation_crn_systems.jl
index db7e14b887..d6c8bd5dc6 100644
--- a/test/reactionsystem_core/coupled_equation_crn_systems.jl
+++ b/test/reactionsystem_core/coupled_equation_crn_systems.jl
@@ -1042,4 +1042,4 @@ let
     u0 = [S1 => 1.0, S2 => 2.0, V1 => 0.1]
     ps = [p1 => 2.0, p2 => 3.0]
     @test_throws Exception ODEProblem(rs, u0, (0.0, 1.0), ps; structural_simplify = true)
-end
\ No newline at end of file
+end

From 3f2a9a99e4d842374795f1137fe8093cbfabf777 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Fri, 6 Sep 2024 09:38:40 -0400
Subject: [PATCH 12/15] revert

---
 src/dsl.jl              |  2 +-
 test/dsl/dsl_options.jl | 10 ----------
 test/runtests.jl        | 12 ++++++------
 3 files changed, 7 insertions(+), 17 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index a3adeddf16..bbaf2ff4ba 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -686,7 +686,7 @@ function read_equations_options(options, variables_declared)
     # When this is the case, the variable X and differential D are extracted (for automatic declaration).
     # Also performs simple error checks.
     vars_extracted = Vector{Symbol}()
-    add_default_diff = true
+    add_default_diff = false 
     for eq in equations
         if (eq.head != :call) || (eq.args[1] != :~)
             error("Malformed equation: \"$eq\". Equation's left hand and right hand sides should be separated by a \"~\".")
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 336241341a..8718f3498d 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -214,16 +214,6 @@ let
     @test equations(rn4)[1] isa Equation
     @parameters v n
     @test isequal(Catalyst.expand_registered_functions(equations(rn4)[1]), D(A) ~ v*(A^n))
-
-
-    rn5 = @reaction_network begin
-        @species A(t) B(t)
-        @equations D(A) + D(B) ~ f(A, t)
-    end
-    @test length(equations(rn5)) == 1
-    @test equations(rn5)[1] isa Equation
-    @species B(t)
-    @test isequal(equations(rn5)[1], D(A) + D(B) ~ 2*A*t)
 end
 
 # Test inferring with stoichiometry symbols and interpolation.
diff --git a/test/runtests.jl b/test/runtests.jl
index 6ad68bf0b6..f786e60612 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -11,12 +11,12 @@ using SafeTestsets, Test
 @time begin
 
     # Tests the `ReactionSystem` structure and its properties.
-    @time @safetestset "Reaction Structure" begin include("reactionsystem_core/reaction.jl") end
-    @time @safetestset "ReactionSystem Structure" begin include("reactionsystem_core/reactionsystem.jl") end
-    @time @safetestset "Higher Order Reactions" begin include("reactionsystem_core/higher_order_reactions.jl") end
-    @time @safetestset "Symbolic Stoichiometry" begin include("reactionsystem_core/symbolic_stoichiometry.jl") end
-    @time @safetestset "Parameter Type Designation" begin include("reactionsystem_core/parameter_type_designation.jl") end
-    @time @safetestset "Custom CRN Functions" begin include("reactionsystem_core/custom_crn_functions.jl") end
+    # @time @safetestset "Reaction Structure" begin include("reactionsystem_core/reaction.jl") end
+    # @time @safetestset "ReactionSystem Structure" begin include("reactionsystem_core/reactionsystem.jl") end
+    # @time @safetestset "Higher Order Reactions" begin include("reactionsystem_core/higher_order_reactions.jl") end
+    # @time @safetestset "Symbolic Stoichiometry" begin include("reactionsystem_core/symbolic_stoichiometry.jl") end
+    # @time @safetestset "Parameter Type Designation" begin include("reactionsystem_core/parameter_type_designation.jl") end
+    # @time @safetestset "Custom CRN Functions" begin include("reactionsystem_core/custom_crn_functions.jl") end
     @time @safetestset "Coupled CRN/Equation Systems" begin include("reactionsystem_core/coupled_equation_crn_systems.jl") end
     @time @safetestset "Events" begin include("reactionsystem_core/events.jl") end
 

From 52e7753a74531ddf7f6f7dfae603ffdad7e1d852 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Fri, 6 Sep 2024 09:42:48 -0400
Subject: [PATCH 13/15] format

---
 src/dsl.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index bbaf2ff4ba..16bd9e4617 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -686,7 +686,7 @@ function read_equations_options(options, variables_declared)
     # When this is the case, the variable X and differential D are extracted (for automatic declaration).
     # Also performs simple error checks.
     vars_extracted = Vector{Symbol}()
-    add_default_diff = false 
+    add_default_diff = false
     for eq in equations
         if (eq.head != :call) || (eq.args[1] != :~)
             error("Malformed equation: \"$eq\". Equation's left hand and right hand sides should be separated by a \"~\".")

From 7790f87e1f885020809e3c462795b64bc91f35aa Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Fri, 6 Sep 2024 10:02:49 -0400
Subject: [PATCH 14/15] test update

---
 test/dsl/dsl_options.jl |  8 ++++++++
 test/runtests.jl        | 12 ++++++------
 2 files changed, 14 insertions(+), 6 deletions(-)

diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 8718f3498d..e471e07c70 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -204,6 +204,14 @@ let
     @test isequal(equations(rn3)[1], D(A) ~ Iapp)
     @test isequal(equations(rn3)[2], Iapp ~ 2*A*t)
 
+    @species Iapp(t) A(t)
+    eq = [D(A) ~ Iapp, Iapp ~ f(A, t)]
+    @named rn3_sym = ReactionSystem(eq, t)
+    rn3_sym = complete(rn3_sym)
+    @test isequivalent(rn3, rn3_sym)
+    
+    
+
 
     g(A, K, n) = A^n + K^n
     rn4 = @reaction_network begin
diff --git a/test/runtests.jl b/test/runtests.jl
index f786e60612..6ad68bf0b6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -11,12 +11,12 @@ using SafeTestsets, Test
 @time begin
 
     # Tests the `ReactionSystem` structure and its properties.
-    # @time @safetestset "Reaction Structure" begin include("reactionsystem_core/reaction.jl") end
-    # @time @safetestset "ReactionSystem Structure" begin include("reactionsystem_core/reactionsystem.jl") end
-    # @time @safetestset "Higher Order Reactions" begin include("reactionsystem_core/higher_order_reactions.jl") end
-    # @time @safetestset "Symbolic Stoichiometry" begin include("reactionsystem_core/symbolic_stoichiometry.jl") end
-    # @time @safetestset "Parameter Type Designation" begin include("reactionsystem_core/parameter_type_designation.jl") end
-    # @time @safetestset "Custom CRN Functions" begin include("reactionsystem_core/custom_crn_functions.jl") end
+    @time @safetestset "Reaction Structure" begin include("reactionsystem_core/reaction.jl") end
+    @time @safetestset "ReactionSystem Structure" begin include("reactionsystem_core/reactionsystem.jl") end
+    @time @safetestset "Higher Order Reactions" begin include("reactionsystem_core/higher_order_reactions.jl") end
+    @time @safetestset "Symbolic Stoichiometry" begin include("reactionsystem_core/symbolic_stoichiometry.jl") end
+    @time @safetestset "Parameter Type Designation" begin include("reactionsystem_core/parameter_type_designation.jl") end
+    @time @safetestset "Custom CRN Functions" begin include("reactionsystem_core/custom_crn_functions.jl") end
     @time @safetestset "Coupled CRN/Equation Systems" begin include("reactionsystem_core/coupled_equation_crn_systems.jl") end
     @time @safetestset "Events" begin include("reactionsystem_core/events.jl") end
 

From 672225ab7812e9d8deaf500e24c82032e467443a Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Mon, 9 Sep 2024 12:46:44 -0400
Subject: [PATCH 15/15] up

---
 src/dsl.jl              |   4 --
 test/dsl/dsl_options.jl | 144 ++++++++++++++++++++--------------------
 2 files changed, 73 insertions(+), 75 deletions(-)

diff --git a/src/dsl.jl b/src/dsl.jl
index 16bd9e4617..bb3fbd6df8 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -921,12 +921,8 @@ end
 
 # Recursively escape functions in the right-hand-side of an equation written using user-defined functions. Special function calls like "hill(...)" are not expanded.
 function escape_equation_RHS!(eqexpr::Expr)
-    # lhs = recursive_escape_functions!(eqexpr.args[2])
-    # eqexpr.args[2] = lhs
-
     rhs = recursive_escape_functions!(eqexpr.args[3])
     eqexpr.args[3] = rhs
-
     eqexpr
 end
 
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index e471e07c70..b92a9bf0e2 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -153,77 +153,6 @@ let
     @test issetequal(parameters(rn5), @parameters k B X)
 end
 
-# Test whether user-defined functions are properly expanded in equations. 
-let
-    f(A, t) = 2*A*t
-
-    # Test user-defined function
-    rn = @reaction_network begin
-        @equations D(A) ~ f(A, t)
-    end
-    @test length(equations(rn)) == 1
-    @test equations(rn)[1] isa Equation
-    @species A(t)
-    @test isequal(equations(rn)[1], D(A) ~ 2*A*t)
-
-
-    hill_unregistered(A, v, K, n) = v*(A^n) / (A^n + K^n)
-    rn2 = @reaction_network begin
-        @parameters v K n
-        @equations D(A) ~ hill_unregistered(A, v, K, n)
-    end
-    @test length(equations(rn2)) == 1
-    @test equations(rn2)[1] isa Equation
-    @parameters v K n
-    @test isequal(equations(rn2)[1], D(A) ~ v*(A^n) / (A^n + K^n))
-
-    hill2(A, v, K, n) = v*(A^n) / (A^n + K^n)
-    @register_symbolic hill2(A, v, K, n)
-    # Test Catalyst function
-    rn2r = @reaction_network begin
-        @parameters v K n
-        @equations D(A) ~ hill2(A, v, K, n)
-    end
-    @test length(equations(rn2r)) == 1
-    @test equations(rn2r)[1] isa Equation
-    @parameters v K n
-    @test isequal(equations(rn2r)[1], D(A) ~ hill2(A, v, K, n))
-
-
-    rn3 = @reaction_network begin
-        @species Iapp(t)
-        @equations begin
-            D(A) ~ Iapp
-            Iapp ~ f(A,t)
-        end
-    end
-    @test length(equations(rn3)) == 2
-    @test equations(rn3)[1] isa Equation
-    @test equations(rn3)[2] isa Equation
-    @variables Iapp(t)
-    @test isequal(equations(rn3)[1], D(A) ~ Iapp)
-    @test isequal(equations(rn3)[2], Iapp ~ 2*A*t)
-
-    @species Iapp(t) A(t)
-    eq = [D(A) ~ Iapp, Iapp ~ f(A, t)]
-    @named rn3_sym = ReactionSystem(eq, t)
-    rn3_sym = complete(rn3_sym)
-    @test isequivalent(rn3, rn3_sym)
-    
-    
-
-
-    g(A, K, n) = A^n + K^n
-    rn4 = @reaction_network begin
-        @parameters v K n
-        @equations D(A) ~ hill(A, v, K, n)*g(A, K, n)
-    end
-    @test length(equations(rn4)) == 1
-    @test equations(rn4)[1] isa Equation
-    @parameters v n
-    @test isequal(Catalyst.expand_registered_functions(equations(rn4)[1]), D(A) ~ v*(A^n))
-end
-
 # Test inferring with stoichiometry symbols and interpolation.
 let
     @parameters k g h gg X y [isconstantspecies = true]
@@ -1024,3 +953,76 @@ let
     @unpack k1, A = rn3
     @test isequal(rl, k1*A^2)
 end
+
+# Test whether user-defined functions are properly expanded in equations. 
+let
+    f(A, t) = 2*A*t
+
+    # Test user-defined function
+    rn = @reaction_network begin
+        @equations D(A) ~ f(A, t)
+    end
+    @test length(equations(rn)) == 1
+    @test equations(rn)[1] isa Equation
+    @species A(t)
+    @test isequal(equations(rn)[1], D(A) ~ 2*A*t)
+
+
+    # Test whether expansion happens properly for unregistered/registered functions. 
+    hill_unregistered(A, v, K, n) = v*(A^n) / (A^n + K^n)
+    rn2 = @reaction_network begin
+        @parameters v K n
+        @equations D(A) ~ hill_unregistered(A, v, K, n)
+    end
+    @test length(equations(rn2)) == 1
+    @test equations(rn2)[1] isa Equation
+    @parameters v K n
+    @test isequal(equations(rn2)[1], D(A) ~ v*(A^n) / (A^n + K^n))
+
+    hill2(A, v, K, n) = v*(A^n) / (A^n + K^n)
+    @register_symbolic hill2(A, v, K, n)
+    # Registered symbolic function should not expand. 
+    rn2r = @reaction_network begin
+        @parameters v K n
+        @equations D(A) ~ hill2(A, v, K, n)
+    end
+    @test length(equations(rn2r)) == 1
+    @test equations(rn2r)[1] isa Equation
+    @parameters v K n
+    @test isequal(equations(rn2r)[1], D(A) ~ hill2(A, v, K, n))
+
+
+    rn3 = @reaction_network begin
+        @species Iapp(t)
+        @equations begin
+            D(A) ~ Iapp
+            Iapp ~ f(A,t)
+        end
+    end
+    @test length(equations(rn3)) == 2
+    @test equations(rn3)[1] isa Equation
+    @test equations(rn3)[2] isa Equation
+    @variables Iapp(t)
+    @test isequal(equations(rn3)[1], D(A) ~ Iapp)
+    @test isequal(equations(rn3)[2], Iapp ~ 2*A*t)
+
+    # Test whether the DSL and symbolic ways of creating the network generate the same system
+    @species Iapp(t) A(t)
+    eq = [D(A) ~ Iapp, Iapp ~ f(A, t)]
+    @named rn3_sym = ReactionSystem(eq, t)
+    rn3_sym = complete(rn3_sym)
+    @test isequivalent(rn3, rn3_sym)
+    
+    
+    # Test more complicated expression involving both registered function and a user-defined function. 
+    g(A, K, n) = A^n + K^n
+    rn4 = @reaction_network begin
+        @parameters v K n
+        @equations D(A) ~ hill(A, v, K, n)*g(A, K, n)
+    end
+    @test length(equations(rn4)) == 1
+    @test equations(rn4)[1] isa Equation
+    @parameters v n
+    @test isequal(Catalyst.expand_registered_functions(equations(rn4)[1]), D(A) ~ v*(A^n))
+end
+